In this work, we propose a novel intrinsic approach to the approximation of ergodic SDEs on Riemannian manifolds, which include Riemannian Langevin dynamics. In opposition to the standard extrinsic approaches such as penalization methods and projection methods, our methodology does not use embeddings or coordinates and only relies on natural geometric operations: geodesics, parallel transport,... We give a criterion for high order of accuracy for the invariant measure, develop new intrinsic numerical methods designed solely for sampling the invariant measure, and derive high order conditions using a new algebraic operation on exotic Lie-Butcher series. In the spirit of the Leimkuhler-Matthews method, our approach prioritizes long time sampling efficiency over finite time accuracy, and outperforms the previous extrinsic and intrinsic approaches in terms of cost for a given accuracy, which we illustrate with several numerical experiments.