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.