  We derive a new methodology for the construction of high order integrators
for sampling the invariant measure of ergodic stochastic differential equations
with dynamics constrained on a manifold. We obtain the order conditions for
sampling the invariant measure for a class of Runge-Kutta methods applied to
the constrained overdamped Langevin equation. The analysis is valid for
arbitrarily high order and relies on an extension of the exotic aromatic
Butcher-series formalism. To illustrate the methodology, a method of order two
is introduced, and numerical experiments on the sphere, the torus and the
special linear group confirm the theoretical findings.
