Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds

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.


Introduction
We consider systems of stochastic differential equations (SDEs) in R d subject to a smooth scalar constraint and a Stratonovich noise of the form dXptq " Π M pXptqqf pXptqqdt`Π M pXptqqΣpXptqq˝dW ptq, Xp0q " X 0 P M, (1.1) where Π M : R d Ñ R dˆd is the orthogonal projection on the tangent bundle of the manifold M " tx P R d , ζpxq " 0u of codimension q, ζ : R d Ñ R q is a given constraint, f : R d Ñ R d is a smooth drift, Σ : R d Ñ R dˆd is a smooth diffusion coefficient and W is a standard ddimensional Brownian motion in R d on a probability space equipped with a filtration and fulfilling the usual assumptions. For simplicity of the analysis, we assume that M is a compact smooth manifold of codimension q " 1. The smoothness and compactness of M guaranty in particular the existence and uniqueness of a solution to (1.1) with bounded moments for all times t ą 0. 2 In addition, thanks to the projection operator Π M , the solution Xptq lies on M for all t ą 0. In the additive noise case where Σpxq " σI d with σ ą 0, equation (1.1) can also be rewritten equivalently with a Lagrange multiplier (see [44,Sect. 3.2.4.1] or [45,Sect. 3.3]) as dXptq " f pXptqqdt`σdW ptq`gpXptqqdλ t , ζpXptqq " 0, Xp0q " X 0 P M, where g " ∇ζ and λ is an adapted stochastic process determined by the equation ζpXq " 0.
A major motivation of model (1.1) appears in computational problems in molecular dynamics with the constrained overdamped Langevin equation (obtained in the particular case where Σpxq " σI d is a constant homothety), dXptq " Π M pXptqqf pXptqqdt`σΠ M pXptqq˝dW ptq, Xp0q " X 0 P M, (1.2) with σ ą 0, f "´∇V and V : R d Ñ R is a smooth potential. The overdamped Langevin equation is widely used to model the motion of a set of particles subject to a potential V in a high friction regime. The possible constraints can be induced for example by strong covalent bonds between atoms, or fixed angles in molecules. Sampling from the constrained overdamped Langevin equation allows to compute the so-called free energy, which is a key quantity in thermodynamic (see, for instance, [18,44,45] and references therein). Equations of the form (1.1) appear naturally when studying conservative SDEs, that is, SDEs possessing an invariant H conserved almost surely by all realisations of (1.1). The solution of conservative SDEs are subject to the constraint ζpXq " 0 with ζpxq " Hpxq´HpX 0 q. Drawing samples on a manifold also has many applications in statistics (see [10,23] and references therein). Under regularity conditions on the generator of the SDE and on M, it was shown in [18,25] that the solution Xptq of the SDE (1.1) is ergodic, that is, there exists a unique invariant measure dµ 8 on M that has a density ρ 8 with respect to dσ M , the canonical measure on M induced by the Euclidean metric of R d , such that for all test functions φ P C 8 pR d , Rq, In the case of the overdamped Langevin equation (1.2) on M, the process is naturally ergodic and the invariant measure is given by dµ 8 " ρ 8 dσ M " 1 Z exp`´2 σ 2 V˘dσ M with Z " ş M exp`´2 σ 2 V˘dσ M . Approximating the quantity ş M φpxqdµ 8 pxq is a computational challenge when the dimension d is high, which is the case in the context of molecular dynamics where the dimension is proportional to the number of particles, because a standard quadrature formula becomes prohibitively expensive. We emphasize that µ 8 is singular with respect to the Lebesgue measure on R d . In addition, the integrator samples should remain on the manifold M. Hence, the order conditions for sampling the invariant measure in the Euclidean context of R d do not generalize straightforwardly to the manifold case. The main goal of this article is to build and analyse high order one-step integrators for approximating ş M φpxqdµ 8 pxq that lie on the manifold M and that have the form X n`1 " ΦpX n , h, ξ n q, (1.4) where the ξ n are standard independent random vectors and h is the numerical step.
There are different ways to approximate the solution of the SDE problem (1.1). A strong approximation focuses on approaching the realisation of a single trajectory of (1.1) for a given realisation of the Wiener process W . A weak approximation approaches the average of functionals of the solution. We focus here on the approximation for the invariant measure, that is approaching averages of functionals of the solution in the stationary state. This convergence is the numerical equivalent of (1.3). The integrator (1.4) is said to have order p for the invariant measure if for all φ P C 8 pR d , Rq, there exists a positive constant Cpφq independent of the initial condition X 0 such that epφ, hq ď Cpφqh p where epφ, hq "ˇˇˇˇlim (1.5) We recall that a scheme of weak order r immediately has order p ě r for the invariant measure. For the underdamped and overdamped Langevin dynamics in R d , the articles [9,40,3,4,41] proposed multiple schemes of high order for the invariant measure with low weak order (typically r " 1). We mention in particular the work [3] that introduced a methodology for the analysis and design of high order integrators for the invariant measure. This methodology, that relies on Talay-Tubaro expansions [63], backward error analysis and modified differential equations for SDEs [67,1,22,36,37], is generalised in the context of manifolds in the present paper.
A widely used and simple numerical scheme for sampling the invariant measure distribution on manifolds is the Euler scheme (see [17,42,44,45] for instance). Two variants exist for the overdamped Langevin equation (1.2), both of order one in the weak sense, or for sampling the invariant measure: the Euler integrator with explicit projection direction X n`1 " X n`h f pX n q`σ ? hξ n`λ gpX n q, ζpX n`1 q " 0, (1.6) and alternatively the Euler integrator with implicit projection direction To the best of our knowledge, no high order numerical integrators for sampling the invariant measure of the overdamped Langevin equation with constraints (1.2) have been proposed in the literature. In [46], an order two discretization based on the RATTLE integrator (see [58,5,31]) is applied to the underdamped Langevin equation, rather than to the overdamped Langevin dynamic (1.2). The previously described discretizations can be combined with Metropolis-Hastings rejection procedures [50,32]. We quote in particular the Markov-Chain Monte-Carlo (MCMC) methods [27,10,45] and the Hybrid Monte-Carlo methods [65,46], where the need for a reverse projection check is shown to be a key step. We also mention the integrators in [66,47] that are based on an Euler discretization and present new approaches for projecting on the manifold. The alternative approach of using Metropolis-Hastings rejection procedure allows to fully remove the bias on the invariant measure. Analogously to the Euclidean case, this procedure does not make high order discretizations obsolete because, in particular, the rejection rate depends on the quality of the discretization and the dimension of the problem in general, and in the case of stiff problems or problems in high dimension, it suffers from timestep restrictions. Note also that in the specific case where M is a Lie group, high order integrators can be naturally obtained using splitting methods, that are, however, typically limited to weak order two of accuracy (see [7] for further details in the context of ODEs). This article proposes new tools for constructing integrators of any high order for sampling the invariant measure of constrained SDEs of the form (1.1) and relies on the formalism of trees and Butcher-series. Originally introduced by Hairer and Wanner in [30], and based on the work of Butcher [13], B-series have proved to be a powerful standard tool for the numerical analysis of deterministic differential equations, as presented, for instance, in the textbooks [29,14]. In the last decades, several works extended B-series to the stochastic context. We mention in particular Burrage and Burrage [11,12] and Komori, Mitsui and Sugiura [35] who first introduced stochastic trees and B-series for studying the order conditions of strong convergence of SDEs, Rößler [53,54,55,56,57] and Debrabant and Kvaernø [20,19,21] for the design and analysis of high order weak and strong integrators on a finite time interval, [6] for creating schemes preserving quadratic invariants, and [38], where tree series were applied to a class of stochastic differential algebraic equations (SDAEs) for the computation of strong order conditions. Finally we mention the recent work [39], that introduced the exotic aromatic B-series for the computation of order conditions for sampling the invariant measure of ergodic SDEs in R d , and that we extend in this paper to the context of SDEs on manifolds.
This article is organized as follows. Section 2 is devoted to the analysis of the accuracy of integrators for sampling the invariant measure on a manifold M. In Section 3, we apply this methodology on a class of Runge-Kutta methods for solving the constrained overdamped Langevin equation (1.2), to derive arbitrary high order conditions for the invariant measure, with special emphasis on order two conditions, and to introduce a new order two scheme that uses only a few evaluations of f per step. The detailed calculations of the order conditions for the invariant measure are done in Section 4 with the help of an extension of the exotic aromatic B-series formalism [39]. We compare in Section 5 the new order two scheme with the Euler scheme (1.7) in numerical experiments on a sphere, a torus and the special linear group SLpmq to confirm its order of convergence for sampling the invariant measure.

High order ergodic approximation on a manifold
In this section, we present a new criterion for building integrators of any order for the invariant measure by extending the R d results in [22,3] to the context of manifolds. We first settle down a few notations and assumptions, before we recall the standard weak expansions of the exact and numerical solution using the backward Kolmogorov equation. For ζ : R d Ñ R a smooth map, we denote g " ∇ζ its gradient, and Gpxq " g T pxqgpxq " |gpxq| 2 the Gram function related to the manifold M " tx P R d , ζpxq " 0u, where we denote by |x| " px T xq 1{2 the Euclidean norm in R d . We assume in the rest of the article that M is a compact and smooth manifold of codimension one embedded in R d . We suppose in addition that the Gram function G is strictly positive on M, Gpxq ě α ą 0 for all x P M. With these notations, the projection Π M on the tangent bundle is given by Π M pxq " I´Gpxq´1gpxqgpxq T . We denote L the generator of the SDE (1.1). It is given, for φ P C 8 pR d , Rq, by where pe i q i"1,...,d is the canonical basis of R d and, for all vectors a 1 , . . . , a m P R d , we use the following notation for differentials in R d , For the overdamped Langevin equation (1.2), the generator (2.1) reduces to where ∇ M ψ :" Π M ∇ψ and div M pHq :" divpHq´G´1pg, H 1 pgqq. The adjoint L˚of the generator (2.1) in L 2 pdσ M q for the SDE (1.1), i.e. , the operator that satisfies for all test is given by Remark 2.1. As L is a self-adjoint operator in L 2 pdµ 8 q, but not in L 2 pdσ M q in general, it could be more natural to perform the analysis in the space L 2 pdµ 8 q. However, as we allow the substages of our numerical integrators to explore the open neighbourhood of M in R d , we shall work in this paper with differential operators that cannot be rewritten in general with intrinsic derivatives on the manifold M. In addition, performing directly the integration by parts calculations in L 2 pdµ 8 q with such operators is not straightforward, and this motivated the choice of L 2 pdσ M q for the analysis. A similar choice was done in [39] in the context of R d .
We follow the framework of [25]. In particular, we rely on the construction of the local orthogonal coordinates. In a neighbourhood N M of the manifold M, there exists an atlas of local orthogonal coordinate systems py, zq P pV Ă R d´1 qˆp´ε, εq for ε ą 0, with respect to local charts ψ : U Ă N M Ñ pV Ă R d´1 qˆp´ε, εq, such that if ψpxq " py, zq, then z " ζpxq. We make the following regularity assumption on the generator L.

Assumption 2.2. On an open neighbourhood
where r x P M is such that ψpr xq " py, 0q and, for k " 1, . . . , d, p r Σ ik py, zqq i P R d´1 is defined as the restriction of the vector pΠ M pr xqΣ ik pxqq i P R d to the tangent space T r x M of M, rewritten in the local orthogonal coordinate system.
This assumption is a variant in the manifold case of the uniform ellipticity property of the generator L in the Euclidean context of R d . In addition, Assumption 2.2 is automatically satisfied for the constrained overdamped Langevin equation (1.2) and yields that the function upx, tq " ErφpXptqq|Xp0q " xs satisfies the backward Kolmogorov equation (see [25]): We refer to [36,37] for similar results in the context of R d . The backward Kolmogorov equation (2.3) allows us to write the following expansion of upx, hq " ErφpXphqq|X 0 " xs for h small enough, 4) where N M is an open neighbourhood of M in R d and the remainder satisfies the esti-mateˇˇR h N pφ, xqˇˇď C N pφq where the constant C N pφq is independent of h and x. We now assume the existence and uniqueness of an invariant measure, as well as an additional regularity property on L, in the spirit of [22, in the context of R d . 2) (see [25,Sect. 2.3] for further details). Assumption 2.3 yields the ergodicity of the process Xptq solution of (1.1) with the unique invariant measure dµ 8 " ρ 8 dσ M on M. To proceed further, we shall assume that the integrator (1.4) is ergodic, that is, there exists a measure dµ h that has a density with respect to dσ M such that (2.5) We refer to [60,61,48,62] in the Euclidian case, and to [25] in the manifold case, and references therein, for further details on the ergodicity of numerical integrators. In addition, we suppose that ErφpX 1 q|X 0 " xs, the numerical analog of upx, hq, can be developed in powers of h as was done, for instance, in [63,3] in the context of R d .  Under Assumptions 2.2 and 2.4, by comparing the expansions (2.6) and (2.4), the integrator has at least weak order p if A j´1 " L j {j! for j " 1, . . . , p. However, as observed already in R d , high order for the invariant measure can be achieved in spite of a low weak order. This is the purpose of Theorem 2.5 where we present a new sufficient condition for a scheme to have order r for the invariant measure. This result, that relies on the powerful tool of backward error analysis for SDEs, is similar to [3,Thm. 3.3] in the context of smooth compact manifolds.

Theorem 2.5. Under Assumptions 2.2, 2.3 and 2.4, if the numerical scheme is consistent (that is, A 0 " L) and ergodic, and if it satisfies in
then it has order r for the invariant measure and the numerical error (1.5) satisfies, for h Ñ 0, where ρ r P C 8 pN M , Rq is the unique solution of the Poisson problem L˚ρ r "´Ar ρ 8 The proof of Theorem 2.5 is detailed in Appendix A for the sake of completeness. The idea is to write an expansion of the error in the spirit of [63], and to generalize the analysis in [22,3] on T d and in [25] to the context of smooth compact manifolds.
Theorem 2.5 states the result for times t Ñ 8. A bound of the error at finite time t n " nh is typically given by the following exponential estimate (see [25,22] where the constant µ ą 0 is in practice the spectral gap of a certain operator that depends on the numerical integrator. Reducing the error term Ke´µ tn is out of the scope of this paper, though the recent works [43,24,2] proposed numerical methods in R d that improve the rate of convergence at infinity, while sometimes also reducing the variance. Remark 2.6. One can consider possible generalisations of Theorem 2.5 in the case where M is not compact, or if M is a manifold of any dimension. We refer to [3] for the non-compact extension of Theorem 2.5 in the context of R d .

High order integrators for constrained Langevin dynamics
In this section, we propose a new class of Runge-Kutta methods for sampling the invariant measure of equation (1.2), and present the methodology for deriving the conditions of any order for the invariant measure using Theorem 2.5. In particular, we compute exactly the consistency and order two conditions for the invariant measure as they are the most relevant for the applications.

Runge-Kutta methods for constrained overdamped Langevin
When discretizing naively equation (1.2), one cannot ensure in general that the integrator stays on M. It is natural to discretize instead the equivalent formulation with Lagrange multipliers The class of numerical schemes we obtain is in the spirit of deterministic Runge-Kutta methods for differential algebraic problems such as the methods SHAKE and RATTLE (see [58,5,31]), introduced in the context of constrained Hamiltonian dynamics, or the SPARK class of methods for general DAEs (see [34]). Since evaluating f is in practical applications the most expensive part of the algorithm compared to evaluating g, we propose high order integrators that are implicit in g and explicit in f in the spirit of implicit-explicit (IMEX) integrators (see, e.g., [31]), so that there are only a few evaluations of f per step. We thus consider the following class of Runge-Kutta integrators where A " pa ij q, p A " pp a ij q P R sˆs and δ i " ř s j"1 p a ij P t0, 1u are the given Runge-Kutta coefficients, and the ξ n " N p0, I d q are independent standard Gaussian random vectors in R d (an alternative with discrete bounded random variables is discussed in Remark 3.2). We fix δ s " 1 so that X n`1 P M and we ask that if δ i " 0, then p a ij " 0 for j " 1, . . . , s (internal stages without projection, Y i R M a.s.). Ideally, one aims for IMEX integrators with a low number of evaluations of f , we hence assume in addition that p A is a lower triangular matrix and A is a strictly lower triangular matrix (in the spirit of DIRK methods). We represent the numerical integrators with their associated Butcher tableau [39] in the Euclidean case R d ) to study partitioned problems where f " f 1`f2 and, for example, to create IMEX schemes. In order to improve the order of the method without increasing its cost, one could also apply a postprocessor (in the spirit of [64] in R d ) or use multiple independent noises in (3.1) instead of only one random variable ξ n " N p0, I d q.
This last extension can increase the number of conditions but may also increase the set of solutions. We refer in particular to [19,39] in the context of R d , where it is shown for a class of stochastic Runge-Kutta method that the order conditions for weak order 3 cannot be satisfied in general, unless we use at least two independent noises. In addition, if we rewrite the internal stages of (3.1) as where g : R d Ñ R dˆq and λ i P R q , then the same class of methods is also fit for solving (1.2) with a multidimensional constraint ζ : R d Ñ R q . Note that the coefficients of the method do not depend on the dimension of the space d or the codimension q of the manifold. This will be studied in future work.

Remark 3.2. If ξ n is a Gaussian random variable, its realisations can be arbitrarily large, and the existence and uniqueness of the solution of the system (3.1) does not hold in general.
A standard remedy to ensure that the projection on M always exists for h ď h 0 small enough is to replace the standard Gaussian random vectors ξ in (3.1) by bounded discrete random vectors p ξ that have the same first moments in the spirit of [51,Chap. 2]. This way, the order of the method is preserved both in the weak sense and for the invariant measure, and the method is well-posed for all h small enough. For weak/ergodic order two, one can consider, for instance, the random vectors p ξ with independent components p ξ i that satisfy The following lemma guarantees the well-posedness of a method of the form (3.1) with bounded random variables p ξ n . The result is still true when A and p A are general matrices, but we consider only the lower triangular case for the sake of brevity. This result is in the spirit of [31, Chap. VII] for deterministic DAEs.  .1) where the ξ n are replaced by bounded random variables p ξ n , there exists h 0 ą 0 such that for all h ď h 0 , for any initial condition X n P M, there exists a unique solution X n`1 of (3.1) in a neighbourhood of X n . Furthermore, the internal stages satisfy Proof. We proceed by induction on i. We assume that for j ă i, the Y j are already defined and The result is straightforward if δ i " 0. We thus assume that δ i " 1 and prove the existence of a unique solution to the equations of the internal stage i: Using ζpX n q " 0, we rewrite equation (3.4) as Multiplying both sides of equation (3.3) by ş 1 0 g T pX n`τ pY i´Xn qqdτ´ř i j"1 p a ij gpY j q¯, and substituting λ i in (3.3) with its value from (3.6), we deduce that F pY i , hq " 0, where the function F : R dˆR Ñ R d is given by F py, tq " As F pX n , 0q " 0 and the partial differential B y F pX n , 0q " GpX n qI d is invertible, the implicit function theorem yields the existence and uniqueness of Y i in a ball of center X n for h ď h 0 small enough. As p ξ n is bounded and M is compact, there exists a deterministic h 0 that works for every initial condition X n P M. Now that Y i is well-posed, we deduce from the identity F pY i , hq " 0 that Y i " X n`O p ? hq and we derive from (3.6) that λ i is well-posed for h small enough and satisfies λ i " Op ? hq. Finally we observe that pY i , λ i q is indeed a solution to (3.3)-(3.4).

Remark 3.4.
In practice, one can solve numerically each internal stage of the set of equations (3.1) with a fixed point iterations or a Newton method starting from Y i " X n and λ i " 0.
As M is compact, if the ξ n are replaced by bounded random variables, these two methods converge for h ď h 0 where h 0 is small enough and independent of the initial condition. It is crucial to initialize the Y i in a neighbourhood of X n as (3.1) has multiple solutions in general. For example, the Euler scheme (1.7) always has two solutions if M is a sphere (the two intersections of M and a straight line going through the center of M).
Before looking at the consistency and the order conditions of the class of methods (3.1), we introduce a concise notation for multiplying vectors component-wise. We present below the detailed calculation of the consistency conditions of the class of methods (3.1) for the constrained overdamped Langevin equation (1.2). Similar proofs can be found in [44,Prop. 3.24] for the Euler schemes (1.6)-(1.7), and in [3,39] for Runge-Kutta methods in R d .

Proposition 3.6. For a Runge-Kutta method of the form
In

7)
then the method is consistent, that is, A 0 " L.
Proof. If we apply one step of a method of the form (3.1) with the initial condition X 0 " x, then the internal stages Y i satisfy the following expansion where the remainder satisfiesˇˇR hˇď Ch 3{2 , and where we used that λ i can be developed in powers of ? h as λ i " ? hλ 1{2,i`h λ 1,i`. . . in the spirit of [44,Lemma 3.25]. If δ i " 1, ζpY i q can also be expanded as where we omitted the dependency in x of G, g, g 1 and the λ k{2,j 's for brevity. We have ζpY i q " ζpxq " 0 (since x P M), thus by identifying each term of the expansion with zero, we get For φ a test function, the operator A 0 φ satisfies By replacing Y s with its expansion in powers of h 1{2 , and by identifying the first terms, we deduce that where we used that δ s " 1 and that all the terms containing an odd number of ξ vanish since odd moments of ξ are zero. Distributing the expectation on each term and using c s " b T 1 yields the desired expression of A 0 φ. We deduce the consistency conditions b T 1 " d s " 1 and p b T d " p b T pδ˛dq in order to get A 0 " L.

Order conditions for the invariant measure on manifolds
We now derive the methodology for getting the conditions of arbitrary high order for sampling the invariant measure of the constrained overdamped Langevin equation (1.2). In particular, the following theorem presents the Runge-Kutta conditions for order two for the invariant measure on M. Note that the number of conditions does not depend on the dimension of the space d.
In the particular case where we set δ " 1, the order two conditions reduce to the following: For simplicity, we used in Theorem 3.8 the notation˛of Definition 3.5. For instance, the The order conditions of Theorem 3.8 can be obtained from straightforward calculations with the following methodology. We compute the operator A 1 with the same method used for A 0 in Proposition 3.6. It is a differential operator of order four with the following first terms where B is a differential operator of order three. We present the complete expansion of A 1 in Section 4 by using a B-series approach. If we assume that p b T d " b T d, then we can integrate by parts to transform ş M A 1 φdµ 8 into an integral of the form ş M A 0 1 φdµ 8 where A 0 1 φ is a differential operator of order one in φ (in the spirit of [3,39]). On a manifold, the integration by parts is a corollary of the Green theorem (see, for instance, [59, Chap. II]). As we shall see below, it reveals a crucial tool for deriving order conditions for the invariant measure. To perform the calculations in a systematic manner, a formalization of the integration by parts process with trees and B-series is presented in Section 4. p2k`1qG´p k`2q pg, g 1 gqpg, Hqψ´G´k divpHqψ`2kG´p k`1q pg, g 1 Hqψ For example, let us integrate by part the terms of order four w.r.t. φ of the operator A 1 φ in equation (3.8). Applying identity (3.9) with ψ " σ 4 8 ∆φ 1 pe i q, H " e i and k " 0, and then summing on i " 1, . . . , d yields ż We apply again (3.9) with ψ " σ 4 8 φ p3q pg, g, e i q, H " e i and k " 1, and then sum on i " 1, . . . , d to get ż M " σ 4 8 G´1∆φ 2 pg, gq´σ 4 8 G´2φ p4q pg, g, g, gq G´2φ p3q pg, g, g 1 gq G´2pg, f qφ p3q pg, g, gq´σ 2 4 G´1φ p3q pg, g, f q ı dµ 8 .
Substracting (3.11) from (3.10) allows to express ş M A 1 φdµ 8 with derivatives of φ of order strictly less than 4. We iterate this method to obtain is an operator of order one in φ, and then find sufficient conditions such that A 0 1 " 0. This implies that A1ρ 8 " 0, and Theorem 2.5 then gives the order two for the invariant measure. The computation of A 0 1 is further detailed in Section 4. Although constructing methods of high weak order is not the main focus of this paper, considering the explicit formula for A 1 and comparing with L 2 {2 (see Section 4 for their detailed expansion in B-series), one immediately obtains the following theorem for weak order two of accuracy. Theorem 3.10 (Runge-Kutta conditions for weak order two). We consider a Runge-Kutta method of the form (3.1) and assume that it satisfies (3.7). If the following conditions are satisfied, then the integrator has weak order two: Remark 3.11. For δ " 1, the weak order two conditions of Theorem 3.10 have no solution, which is in contrast with the invariant measure case presented in Theorem 3.8. Indeed, the condition p b T App1´δq˛dq " 1 4 cannot be fulfilled if we fix δ " 1.

Illustrative examples of high order Runge-Kutta methods on manifolds
In hd 3 ξ n`λ3 rp a 31 gpY 1 q`p a 32 gpY 2 q`p a 33 gpY 3 qs , where λ 1 , λ 2 , λ 3 , λ 4 are such that ζpY 1 q " ζpY 2 q " ζpY 3 q " ζpX n`1 q " 0, and with the values of c i , d i , p a ij given in Appendix C. To implement one step of this scheme, we apply a few iterations of the Newton method to find the projections on M. We emphasize that if the stepsize h is not small enough, the fixed point problems of finding λ i such that ζpY i q " 0 may not be well defined, leading to diverging Newton iterations. Following Remark 3.2, we replace the standard Gaussian random vectors ξ in (3.1) by independent bounded discrete random vectors p ξ that satisfy (3.2). This way, the order two for the invariant measure is preserved and the method is well-posed for h small enough.
With the same methodology we used to obtain the order conditions of Theorem 3.8 and Theorem 3.10, and with the expressions of A 1 φ and A 0 1 φ (see Section 4 for further details), we also get classes of Runge-Kutta integrators and their order conditions for the following specific subproblems.
Euclidean case R d . Fixing g " 0 in the expressions of A 1 φ and A 0 1 φ yields the order two conditions in the weak sense and for the invariant measure in R d as given in [39, Tables 1-2].
Deterministic case. Fixing σ " 0 in the expression of A 1 φ yields the order conditions for approximating the solution of ODEs of the form 9 x " Π M pxqf pxq, where f is a gradient. Note that this equation can be rewritten as the following differential algebraic equation (DAE) of index two (see [31,Chap. VII]): We obtain a class of deterministic Runge-Kutta methods for solving DAEs of the form (3.13) by setting σ " 0 in (3.1). A Runge-Kutta method of this form is consistent if b T 1 " 1, and has order two if p b T c " p b T pδ˛cq " b T c " b T pδ˛cq " 1{2. For instance, an order two method for solving ODEs of the form (3.13) is X n`1 " X n`h f pX n q`f pX n`1 q 2`λ gpX n q`gpX n`1 q 2 , ζpX n`1 q " 0.
Spherical case. In the simple case where M is the unit sphere in R d (that is, when the constraint is of the form ζpxq " p|x| 2´1 q{2 and gpxq " x), the consistency conditions (3.7) reduce to b T 1 " d s " 1. The weak order two conditions of Theorem 3.10 reduce to the following conditions: On the other hand, the order two conditions for the invariant measure of Theorem 3.8 on the sphere are the following: For example, the following integrator has order two for the invariant measure if M is a sphere: Brownian motions on manifolds. Runge-Kutta methods of the form (3.1) can also be used for simulating a Brownian motion on a manifold (see [33,Chap. III]) by solving numerically We recall that in the context of R d , the Euler-Maruyama integrator is exact for approximating a Brownian motion in law. However, in the context of manifolds, there are no exact Runge-Kutta integrators for simulating a Brownian motion on M in general. In particular, the Euler scheme (1.7) only has weak order one for solving (3.14) in general. Fixing f " 0 in (3.1) yields a class of Runge-Kutta methods for solving (3.14). The consistency conditions are d s " 1 and p b T d " p b T pδ˛dq. The conditions for order two for the invariant measure (respectively for weak order two) of such a Runge-Kutta method are obtained by deleting the order conditions in Theorem 3.8 (respectively in Theorem 3.10) that involve A, b or c. In the specific case where M is a sphere, the consistency conditions (3.7) become d s " 1 and the weak order two conditions of Theorem 3.10 reduce to the two following conditions For example, a weak order two method for simulating a Brownian motion on a sphere is X n`1 " X n`? hξ n`λ 3X n`? hξ n`Xn`1 4 , ζpX n`1 q " 0.
In addition, there are no additional order two conditions for the invariant measure, that is, any consistent integrator, such as the Euler scheme (1.7), has at least order two for the invariant measure on the sphere.

Exotic aromatic B-series for computing order conditions
As described in the introduction, B-series were introduced to tackle the calculations of order conditions of ODEs by representing Taylor expansions with trees. In [16], an extension of the original B-series, called aromatic B-series, was used to study volume-preserving integrators. It allowed in particular to represent the divergence of a B-series. B-series and aromatic Bseries were also studied later in [52,49,26] for their geometric properties, and in [15,8] for their algebraic structure of Hopf algebras. In [39], a new formalism of B-series, called exotic aromatic B-series, was introduced for computing order conditions for sampling the invariant measure of SDEs in R d . It added a new kind of edge, called liana, to the aromatic trees in order to represent new terms such as the Laplacian of an aromatic B-series. In this section, we extend the formalism of exotic aromatic B-series by allowing the representation of scalar products, and show that the operators L j and A j can be represented conveniently in the form of B-series. We also rewrite the integration by parts formula (3.9) as a straightforward process on graphs, and apply it to compute A 0 1 .
We consider graphs γ " pV, E, Lq where V is the set of nodes, E the set of edges and L the set of lianas. We split the set of edges into E " E 0 Y E S where E 0 are the standard oriented edges as defined in [39], and where E S is a new set of non-oriented edges represented as double horizontal straight lines. If pv, wq " pw, vq P E S , we consider this edge as an outgoing edge for both v and w, but v and w are not predecessors of each other. If pv, wq P E S , we denote Spvq " w and Spvq " v otherwise. We consider graphs where each node has exactly one outgoing edge, except exactly one node, called the root r, that has none. If we consider only the graph pV, Eq, where we erase the lianas, it can be decomposed in two kinds of connected components: one that contains the root, that we name the rooted tree, and the other components that we name aromas. We decompose the set of nodes in V " V f Y V g Y tru where V f are the nodes representing the function f and are represented with black disks (respectively V g represent the function g and are drawn with white disks). We write N f pγq the number of elements of V f (respectively N g pγq the number of elements of V g ) and N l pγq the number of lianas. The order of a directed graph γ " pV, E, Lq is defined as For instance, the graph γ " pV, E, Lq with satisfies |γ| " 7 and is represented as We say that two directed graphs pV 1 , E 1 , L 1 q and pV 2 , E 2 , L 2 q are equivalent if there exists a bijection ϕ : We call exotic aromatic forests the equivalence classes of these directed graphs γ " pV, E, Lq, and we denote EAT the set of exotic aromatic forests. In addition, we need a different set of rooted forests where the root is in V f or V g . We call them exotic aromatic vector fields and gather them together in the set EAV. The elementary differential associated to an exotic aromatic forest is given by the following definition. Definition 4.1. Let γ " pV, E, Lq P EAT , and let f , g : R d Ñ R d and φ : R d Ñ R be smooth functions. We denote l 1 , . . . , l s the elements of L, v 1 , . . . , v m the elements of V tru and δ i,j the Kronecker symbol (δ i,j " 1 if i " j, δ i,j " 0 else). We use the notation for v P V , I πpvq " pi q 1 , . . . , i qs q where πpvq " tq 1 , . . . , q s u are the predecessors of v, and J Γpvq " pj lx 1 , . . . , j lx t q where Γpvq " tl x 1 , . . . , l xt u are the lianas linked to v. Then F pγq is defined as For example, the differential associated to the exotic aromatic forest γ given by (4.1) is We extend the definition of F on SpanpEAT q by linearity and write, for the sake of simplicity, F pγqpφq instead of F pγqpf, g, φq. An exotic aromatic B-series is a formal series indexed over EAT of the form Bpaqpφq " ÿ γPEAT h |γ| apγqF pγqpφq.

Remark 4.2.
As we assumed that the functions f and g are gradients, multiple exotic aromatic forests can represent the same differential. We do not detail here the method to identify two such forests as it is similar to [39,Prop. 4.7] in the context of R d .
The following result states that the operators L j {j! and A j can be written with exotic aromatic forests. We omit the proof for the sake of brevity as it is similar to [39,Thm. 4.1].

Proposition 4.3. Take a Runge-Kutta method of the form (3.1), then the expansions (2.4) and (2.6) can be formally written with exotic aromatic B-series, that is, there exists two maps e and a over EAT such that
ErφpXphqq|Xp0q " xs " Bpeqpφqpxq, ErφpX 1 q|X 0 " xs " Bpaqpφqpxq, and where the operators are given by If epγq " apγq for all γ P EAT with 1 ď |γ| ď p, then the integrator has at least weak order p.
For example, the operator L in (2.2) can be rewritten with exotic aromatic forests as We present in Table 2 (see Appendix D) the decomposition in exotic aromatic forests of the operators L 2 φ{2 " LpLφq{2 and A 1 φ under the consistency condition (3.7).

Remark 4.4.
If we replace the functions g and φ by f and fix σ " G " 1, the newly obtained exotic aromatic B-series satisfy an isometric equivariance property, that is, they stay unchanged when applying an isometric coordinate transformation. It was proved in [52] that, under a condition of locality, aromatic B-series are exactly the affine equivariant methods, that is, the maps that stay unchanged when applying an affine coordinate transformation. Analogously, it would be interesting to make a link between the isometric equivariant maps and the exotic aromatic B-series.
In the spirit of the Butcher product on trees [29,Chap. III], we introduce a few notations for writing with ease different operations on forests.
Notation. Let γ be an exotic aromatic forest/vector field, τ be an exotic aromatic vector field and v a node of γ, then we define the following operators on forests.
1. For example, let γ " , τ " and v " r the root of γ, then we get The integration by parts (3.9) can be rewritten conveniently with exotic aromatic forests.
We write γ " r γ if it is possible to go from γ P EAT to r γ P SpanpEAT q with the processes of integration by parts (4.2) or (4.3). We extend this relation by linearity on SpanpEAT q and make it symmetric so that " becomes an equivalence relation on SpanpEAT q. For example, the integrations by parts (3.10) and (3.11) can be rewritten with exotic aromatic B-series by using (4.3) with γ " and γ " . It yields For the sake of completeness, we present in Appendix B the integrations by parts for the order 3 terms of A 1 φ. The computations are similar for the terms of order two in φ.

Remark 4.6.
In the Euclidean case R d , that is, for a forest γ P EAT and a vector field τ P EAV with N g pγq " N g pτ q " 0 and g " 0, Lemma 4.5 reduces to the two following equations: We recover the process of integration by parts described in [39,Thm. 4.4] in the context of exotic aromatic B-series in R d .
We can now revisit the statement of Theorem 2.5 in terms of B-series.
Theorem 4.7. Take a consistent ergodic Runge-Kutta method of the form (3.1). We denote A i " F pγ i q with γ i P EAT . If γ i " γ 0 i and F pγ 0 i q " 0 for 1 ď i ă r, then the method has at least order r for the invariant measure.
By applying repeatedly the process of integrations by parts described in Lemma 4.5, one can simplify the operator A i " F pγ i q into an operator of the form A 0 i " F pγ 0 i q such that γ i " γ 0 i . The complete decomposition of A 0 1 into exotic aromatic forests is detailed in Table 3 (see Appendix D). According to Theorem 4.7, choosing the coefficients of the Runge-Kutta method such that γ 0 1 " 0 yields the order two conditions for the invariant measure, as stated in Theorem 3.8.

Remark 4.8.
We call EAT 0 the subset of exotic aromatic forests whose root has only one predecessor (that is, the forests associated to an order one operator) or that have a rooted tree of the form , , , . . . Then, if γ P EAT , there exists γ 0 P EAT 0 such that γ " γ 0 . For instance, for a consistent method of the form (3.1), the operator A 0 1 " F pγ 0 1 q has the form a 0 pγqγ, so that γ 0 1 P EAT 0 , and A 0 1 is a differential operator of order one if the condition b T d " p b T d holds.

Numerical experiments
In this section, we perform numerical experiments to confirm the theoretical findings, first on a sphere and a torus in R 3 , and then on the special linear group.

Invariant measure approximation on a sphere and a torus
To check the numerical order two of the Runge-Kutta integrator (3.12) presented in Section 3.3, we first compare it with the Euler scheme (1.7) on the unit sphere in R 3 , where the constraint is given by ζpxq " px 2 1`x 2 2`x 2 3´1 q{2. We choose the potential V pxq " 25p1´x 2 1´x 2 2 q, with σ " ? 2, φpxq " x 2 3 , f "´∇V , g " ∇ζ, M " 10 7 independent trajectories to have a small Monte-Carlo error and a final time T " 20. Observe that for the smaller final time T " 10 (not included in the figures for conciseness), the convergence curves reveal nearly identical to the case T " 20 considered in Figure 1, which suggests that the numerical solutions are already very close to equilibrium at these final times. Following Remark 3.2 and Lemma 3.3, we use discrete bounded random variables satisfying (3.2) in the implementation of the integrators. For both integrators, we compute the Monte-Carlo is the m-th realisation of the integrator at time t n " nh, and N is an integer satisfying N h " T . We compare this approximation with a reference value of ş M φdµ 8 computed via a standard quadrature formula, and we plot the error for the invariant measure (1.5) versus different timestep h. We also plot an estimate of the Monte-Carlo error by using the standard error of the mean estimator We observe in all convergence plots that the Monte-Carlo error prevails for small values of the timestep h. On Figure 1, we observe as expected order one for the Euler scheme (1.7) and order two for the Runge-Kutta scheme (3.12). We then apply the Euler scheme (1.7) and the Runge-Kutta integrator (3.12) on a torus defined by the constraint ζpxq " px 2 1`x 2 2 q with R " 3 and r " 1. The potential is V pxq " 25px 3´r q 2 and we choose σ " ? 2, φpxq " x 2 3 , f "´∇V , g " ∇ζ, a final time T " 20 and M " 10 7 independent trajectories. On Figure 2, we plot the error for the invariant measure versus the timestep h, by using a reference value for ş M φdµ 8 obtained with a standard quadrature formula. As expected, we observe order two for the proposed integrator. These curves confirm the theoretical findings presented in Section 3. In particular, the scheme (3.12) has order two of accuracy for the invariant measure on manifolds, according to Theorem 3.8. Note that if we had chosen a very short final time T , we would have observed the weak order one instead of the order two for the invariant measure as we would not have reached equilibrium.

Invariant measure approximation on the special linear group
Sampling on a manifold M is especially useful to compute integrals of the form ş M φpxqdµ 8 when M is a manifold of high dimension. The class of methods (3.1) is convenient as the number of order conditions does not increase with the dimension of the space increasing. We apply Method (3.12) on a Lie group (in the spirit of [65,66]) to see how it performs in high dimension. We choose the special linear group SLpmq " tM P R mˆm , detpM q " 1u, seen as a submanifold of R m 2 of codimension 1. As explained in Remark 2.6, our analysis still applies to SLpmq if we choose a potential V with appropriate growth assumptions, even if it is not a compact manifold. We compare the Euler scheme (1.7) and the Runge-Kutta integrator (3.12) on M " SLpmq for different m (that is, with the constraint ζpxq " detpxq´1), where we use does not converge for approximately 0.005% of the trajectories for m " 4. We choose to discard these trajectories, which induces a negligible bias in the expectation. This does not occur for the Runge-Kutta integrator (3.12). We recall that for a small enough timestep h, the Newton method would always converge (see also Remark 3.4). The reference solution for Jpmq " ş SLpmq φpxqdµ 8 pxq is computed with the Runge-Kutta method (3.12) with h ref " 2´1 4 T . With the factor 25 in the potential (5.1), the solution of (1.2) stays close to I m 2 , and Jpmq is close to φpI m 2 q " m. This choice of factor permits to explore a reasonably small area of SLpmq with moderate manifold curvature. We observe numerically that replacing the factor 25 by 1 in (5.1) induces a severe timestep restriction (results not included for conciseness). The computation of Jpmq could also be done via the parametrization given by the Iwasawa decomposition for SLpmq (see, for instance, [28,Chap. 1]) and the use of standard quadrature methods, but these methods have prohibitive costs in high dimension. We put together the numerical results in Table 1 Table 1: Numerical approximation of the integral Jpmq " ş SLpmq φpxqdµ 8 for 2 ď m ď 4 with the estimator s J " M´1 ř M k"1 φpX pkq N q where pX n q is given by the Euler scheme (1.7) for s J Euler and by the Runge-Kutta integrator (3.12) for s J 2 , with their respective errors. The average is taken over M " 10 6 trajectories, with the potential (5.1), φpxq " Trpxq, a final time T " 10 and a timestep h " 2´1 2 T .

Appendices A Proof of Theorem 2.5
In the spirit of backward error analysis for differential equations (see [29,67,1,22]), we build a modified generator L h such that U px, hq " ErφpX 1 q|X 0 " xs formally satisfies U px, hq " ÿ where the B l are the Bernoulli numbers (see [22,67,3] for similar expansions in T d or R d ).
Using Assumption 2.3, we build recursively a sequence of functions pρ n q such that L˚ρ n "´n ÿ l"1 Ll ρ n´l and ρ 0 " ρ 8 , and R is a differential operator of order two. Using Lemma 4.5 multiple times, we get the following integrations by parts of the order three terms of Bφ.

C Coefficients of the order two Runge-Kutta method
The coefficients of the Runge-Kutta method (3.12) used in Section 5 are