Adrien Laurent

Gilles Vilmart

C Rq

Order Conditions for Sampling the Invariant Measure of Ergodic Stochastic Differential Equations on Manifolds

Keywords: constrained stochastic differential equations, manifolds, invariant measure, ergodicity, exotic aromatic B-series, order conditions AMS subject classification (2020): 60H35, 70H45, 37M25, 65L06

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 pXptqqfpXptqqdt Π M pXptqqΣpXptqq ¥ dW ptq, Xp0q X 0 M, (1.1) where Π M : R d Ñ R d¢d is the orthogonal projection on the tangent bundle of the manifold M tx 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 guarantee 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 [START_REF] Lelièvre | Free energy computations[END_REF]Sect. 3.2.4.1] or [START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF]Sect. 3.3]) as dXptq f pXptqqdt σdW ptq gpXptqqdλ t , ζpXptqq 0, Xp0q X 0 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 pXptqqfpXptqqdt σΠ M pXptqq ¥ dW ptq, Xp0q X 0 M,

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, [START_REF] Ciccotti | Projection of diffusions on submanifolds: application to mean force computation[END_REF][START_REF] Lelièvre | Free energy computations[END_REF][START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF] 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 [START_REF] Brubaker | A family of MCMC methods on implicitly defined manifolds[END_REF][START_REF] Diaconis | Advances in modern statistical theory and applications: a Festschrift in honor of Morris L. Eaton[END_REF] and references therein). Under regularity conditions on the generator of the SDE and on M, it was shown in [START_REF] Ciccotti | Projection of diffusions on submanifolds: application to mean force computation[END_REF][START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF] that the solution Xptq of the SDE (1.1) is ergodic, that is, there exists a unique invariant measure dµ V on M that has a density ρ V with respect to dσ M , the canonical measure on M induced by the Euclidean metric of R d , such that for all test functions φ C V pR d , Rq,

lim T ÑV 1 T » T 0 φpXptqqdt » M φpxqdµ V pxq almost surely. (1.3)
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µ

V ρ V dσ M 1 Z exp ¡ 2 σ 2 V ¨dσ M with Z ³ M exp ¡ 2 σ 2 V ¨dσ M .
Approximating the quantity ³ M φpxqdµ V 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 µ V 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µ V 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 φ C V 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

N ÑV 1 N 1 N ņ0 φpX n q ¡ » M φdµ V § § § § . (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 [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF][START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF][START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for Langevin dynamics[END_REF][START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF] proposed multiple schemes of high order for the invariant measure with low weak order (typically r 1). We mention in particular the work [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] that introduced a methodology for the analysis and design of high order integrators for the invariant measure. This methodology, which relies on Talay-Tubaro expansions [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], backward error analysis and modified differential equations for SDEs [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Abdulle | High weak order methods for stochastic differential equations based on modified equations[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF][START_REF] Kopec | Weak backward error analysis for overdamped Langevin processes[END_REF], 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 [START_REF] Ciccotti | Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics[END_REF][START_REF] Lelièvre | Analyse de certains schémas de discrétisation pour des équations différentielles stochastiques contraintes[END_REF][START_REF] Lelièvre | Free energy computations[END_REF][START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF] 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 hf pX n q σ c hξ n λgpX n q, ζpX n 1 q 0, (1.6) and alternatively the Euler integrator with implicit projection direction

X n 1 X n hf pX n q σ c hξ n λgpX n 1 q, ζpX n 1 q 0.

(1.7)

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 [START_REF] Lelièvre | Hybrid Monte Carlo methods for sampling probability measures on submanifolds[END_REF], an order two discretization based on the RATTLE integrator (see [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF][START_REF] Andersen | Rattle: A "velocity" version of the Shake algorithm for molecular dynamics calculations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF]) 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 [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. We quote in particular the Markov-Chain Monte-Carlo (MCMC) methods [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Brubaker | A family of MCMC methods on implicitly defined manifolds[END_REF][START_REF] Lelièvre | Langevin dynamics with constraints and computation of free energy differences[END_REF] and the Hybrid Monte-Carlo methods [START_REF] Zappa | Monte Carlo on manifolds: sampling densities and integrating functions[END_REF][START_REF] Lelièvre | Hybrid Monte Carlo methods for sampling probability measures on submanifolds[END_REF], where the need for a reverse projection check is shown to be a key step. We also mention the integrators in [START_REF] Zhang | Ergodic SDEs on submanifolds and related numerical sampling schemes[END_REF][START_REF] Lelièvre | Multiple projection MCMC algorithms on submanifolds[END_REF] 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. Analogous 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 [START_REF] Blanes | On the necessity of negative coefficients for operator splitting schemes of order higher than two[END_REF] 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 [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF], and based on the work of Butcher [START_REF] Butcher | An algebraic theory of integration methods[END_REF], 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 [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF]. In the last decades, several works extended B-series to the stochastic context. We mention in particular Burrage and Burrage [START_REF] Burrage | High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations[END_REF][START_REF] Burrage | Order conditions of stochastic Runge-Kutta methods by B-series[END_REF] and Komori, Mitsui and Sugiura [START_REF] Komori | Rooted tree analysis of the order conditions of ROW-type scheme for stochastic differential equations[END_REF] who first introduced stochastic trees and B-series for studying the order conditions of strong convergence of SDEs, Rößler [START_REF] Rößler | Stochastic Taylor expansions for the expectation of functionals of diffusion processes[END_REF][START_REF] Rößler | Rooted tree analysis for order conditions of stochastic Runge-Kutta methods for the weak approximation of stochastic differential equations[END_REF][START_REF] Rößler | Runge-Kutta methods for Itô stochastic differential equations with scalar noise[END_REF][START_REF] Rößler | Stochastic Taylor expansions for functionals of diffusion processes[END_REF][START_REF] Rößler | Strong and weak approximation methods for stochastic differential equations-some recent developments[END_REF] and Debrabant and Kvaernø [START_REF] Debrabant | B-series analysis of stochastic Runge-Kutta methods that use an iterative scheme to compute their internal stage values[END_REF][START_REF] Debrabant | Runge-Kutta methods for third order weak approximation of SDEs with multidimensional additive noise[END_REF][START_REF] Debrabant | Composition of stochastic B-series with applications to implicit Taylor methods[END_REF] for the design and analysis of high order weak and strong integrators on a finite time interval, [START_REF] Anmarkrud | Order conditions for stochastic Runge-Kutta methods preserving quadratic invariants of Stratonovich SDEs[END_REF] for creating schemes preserving quadratic invariants, and [START_REF] Küpper | A Runge-Kutta method for index 1 stochastic differentialalgebraic equations with scalar noise[END_REF], 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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF], 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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]. 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 [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] 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 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 M. With these notations, the projection Π M on the tangent bundle is given by Π M pxq I ¡ Gpxq ¡1 gpxqgpxq T . We denote L the generator of the SDE (1.1). It is given, for φ C V pR d , Rq, by

Lφ φ I pΠ M f q 1 2 d i1 φ I ppΠ M Σe i q I pΠ M Σe i qq 1 2 d i1 φ P pΠ M Σe i , Π M Σe i q, (2.1)
where pe i q i1,...,d is the canonical basis of R d and, for all vectors a 1 , . . . , a m R d , we use the following notation for differentials in R d ,

φ pmq pa 1 , . . . , a m q d i1 ,...,im1 f i 1 ,...,im φ a 1 i 1 . . . a m im d i1 ,...,im1 f m φ fx i 1 . . . fx im a 1 i 1 . . . a m im .
For the overdamped Langevin equation (1.2), the generator (2.1) reduces to

Lφ φ I f ¡ G ¡1 pg, f qφ I g ¡ σ 2 2 G ¡1 divpgqφ I g σ 2 2 G ¡2 pg, g I gqφ I g σ 2 2 ∆φ (2.2) ¡ σ 2 2 G ¡1 φ P pg, gq σ 2 2 exp ¡ 2 σ 2 V is given by L ¦ φ ¡ div M pφfq 1 2 d i1
div M pdiv M pφΣe i qΣe i q.

Remark 2.1. As L is a self-adjoint operator in L 2 pdµ V 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µ V 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µ V 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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF] in the context of R d .

We follow the framework of [START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF]. 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 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. x M is such that ψpr xq py, 0q and, for k 1, . . . , d, p r Σ ik py, zqq i R d¡1 is defined as the restriction of the vector pΠ M pr xqΣ ik pxqq i R d to the tangent space T r

Assumption 2.2. On an open neighbourhood

N M of M in R d ,
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 [START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF]): fu ft px, tq Lupx, tq, upx, 0q φpxq, x N M , t ¡ 0.

(2.3)

We refer to [START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF][START_REF] Kopec | Weak backward error analysis for overdamped Langevin processes[END_REF]

N j1 h j j! L j φpxq h N 1 R h N pφ, xq, x N M , (2.4)
where N M is an open neighbourhood of M in R d and the remainder satisfies the estimate § § 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 [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF] in the context of R d .

Assumption 2.3. There exists an open neighbourhood N

M of M in R d and a unique positive function ρ V C V pN M , Rq satisfying ³ M ρ V dσ M 1 and L ¦ ρ V 0 on N M . Moreover, for all φ C V pN M , Rq such that ³ M φdσ M 0, there exists a unique solution ρ C V pN M , Rq to the Poisson problem L ¦ ρ φ that satisfies ³ M ρdσ M 0.
The existence and uniqueness of the invariant measure are in particular satisfied for the constrained overdamped Langevin equation (1.2) (see [START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF]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µ V ρ V 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 lim

N ÑV 1 N 1 N ņ0 φpX n q » M φdµ h almost surely.
(2.5)

We refer to [START_REF] Talay | Second-order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF][START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Talay | Stochastic Hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF] in the Euclidean case, and to [START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF] 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 [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] in the context of R d .

Assumption 2.4. For all φ C V pR d , Rq, the numerical integrator (1.4) has a weak Taylor expansion of the form

ErφpX 1 q|X 0 xs φpxq N j1 h j A j¡1 φpxq h N 1 R h N pφ, xq, x N M , (2.6)
for all h assumed small enough, and where N M is an open neighbourhood of M in R d and the remainder satisfies § § R h N pφ, xq § § ¤ C N pφq where the constant C N pφq is independent of h and x. The A j 's, j 0, 1, 2, . . . are linear differential operators with coefficients depending smoothly on f , g and their (high order) derivatives (and depending on the choice of the integrator).

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 [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF]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 L 2 pdσ M q A ¦ j ρ V 0, j 1, . . . , r ¡ 1, then it has order r for the invariant measure and the numerical error (1.5) satisfies, for h Ñ 0,

epφ, hq h r » M φpxqρ r pxqdσ M pxq Oph r 1 q h r » V 0 » M upx, tqA ¦ r ρ V pxqdσ M pxqdt Oph r 1 q, where ρ r C V pN M , Rq is the unique solution of the Poisson problem L ¦ ρ r ¡A ¦ r ρ V in N M that satisfies ³ M ρ r dσ M 0, with N M an open neighbourhood of M in R d .
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 [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], and to generalize the analysis in [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] on T d and in [START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF] to the context of smooth compact manifolds. Theorem 2.5 states the result for times t Ñ V. A bound of the error at finite time t n nh is typically given by the following exponential estimate (see [START_REF] Faou | Conservative stochastic differential equations: mathematical and numerical analysis[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF]

) § § § § ErφpX n qs ¡ » M φpxqdµ V pxq § § § § ¤ Ke ¡µtn Ch p ,
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 [START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF][START_REF] Duncan | Variance reduction using nonreversible Langevin samplers[END_REF][START_REF] Abdulle | Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations[END_REF] 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 [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] 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 dX f pXqdt σdW gpXqdλ t , ζpXq 0, Xp0q X 0 M.

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 [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF][START_REF] Andersen | Rattle: A "velocity" version of the Shake algorithm for molecular dynamics calculations[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF]), introduced in the context of constrained Hamiltonian dynamics, or the SPARK class of methods for general DAEs (see [START_REF] Jay | Structure preservation for constrained dynamics with super partitioned additive Runge-Kutta methods[END_REF]). 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., [START_REF] Hairer | Solving ordinary differential equations[END_REF]), so that there are only a few evaluations of f per step. We thus consider the following class of Runge-Kutta integrators

Y i X n h s j1 a ij f pY j q σ c hd i ξ n λ i s j1 p a ij gpY j q, i 1, . . . , s, ζpY i q 0 if δ i 1, i 1, . . . , s, (3.1) X n 1 Y s ,
where A pa ij q, p A pp a ij q R s¢s and δ i °s j1 p a ij 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 M and we ask that if δ i 0, then p a ij 0 for j 1, . . . , s (internal stages without projection, Y i 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, where b pa s,i q i , p b pp a s,i q i , c A1 and 1 p1, . . . , 1q T .

c A δ p A d b T p b T
For instance, the Euler schemes can be written as Runge-Kutta methods of the form (3.1)

with s 2 and the following Butcher tableaux.

Euler (1.6) :

0 0 0 0 0 0 0 1 1 0 1 1 0 1 1 0 1 0 Euler (1.7) : 0 0 0 0 0 0 0 1 1 0 1 0 1 1 1 0 0 1
Note that the class of methods (3.1) satisfies automatically Assumption 2.4.

Remark 3.1. The class of Runge-Kutta methods (3.1) can be straightforwardly generalized (as done in [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF] in the Euclidean case R d ) to study partitioned problems where f f 1 f 2 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 [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF] 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 [START_REF] Debrabant | Runge-Kutta methods for third order weak approximation of SDEs with multidimensional additive noise[END_REF][START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF] 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

Y i X n h s j1 a ij f pY j q σ c hd i ξ n ¢ s j1 p a ij gpY j q λ i ,
where g : R d Ñ R d¢q and λ i 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 [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]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

Pp p ξ i 0q 2 3 and Pp p ξ i ¨c3q 1 6 , i 1, . . . , d. (3.2)
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 [START_REF] Hairer | Solving ordinary differential equations[END_REF]Chap. VII] for deterministic DAEs. Lemma 3.3. For Runge-Kutta methods of the form (3.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 M, there exists a unique solution X n 1 of (3.1) in a neighbourhood of X n . Furthermore, the internal stages satisfy Y i X n Op c hq and λ i Op c hq for i 1, . . . , s.

Proof. We proceed by induction on i. We assume that for j i, the Y j are already defined and satisfy Y j X n Op c hq. 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:

Y i X n h i¡1 j1 a ij f pY j q σ c hd i p ξ n λ i i j1 p a ij gpY j q, (3.3
)

ζpY i q 0. (3.4)
Using ζpX n q 0, we rewrite equation (3.4) as

ζpY i q ¡ ζpX n q » 1 0 g T pX n τ pY i ¡ X n qqdτpY i ¡ X n q 0. (3.5)
Inserting (3.3) in (3.5) yields » 1 0 g T pX n τ pY i ¡ X n qqdτ h i¡1 j1 a ij f pY j q σ c hd i p ξ n λ i i j1
p a ij gpY j q % 0.

(3.6)

Multiplying both sides of equation (3.3) by

³ 1 0 g T pX n τ pY i ¡ X n qqdτ ¡ °i j1 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 » 1 0 g T pX n τ py ¡ X n qqdτ ¡ t i¡1 j1 a ij f pY j q σ c td i p ξ n ©¡ i¡1 j1 p a ij gpY j q p a ii gpyq © ¡ i¡1 j1 p a ij gpY j q p a ii gpyq ©¡ y ¡ X n ¡ t i¡1 j1 a ij f pY j q ¡ σ c td i p ξ n © & .
As F pX n , 0q 0 and the partial differential f 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 M. Now that Y i is well-posed, we deduce from the identity F pY i , hq 0 that Y i X n Op c hq and we derive from (3.6) that λ i is well-posed for h small enough and satisfies λ i Op c 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 [START_REF] Lelièvre | Free energy computations[END_REF]Prop. 3.24] for the Euler schemes (1.6)-(1.7), and in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF][START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF] for Runge-Kutta methods in R d . Proposition 3.6. For a Runge-Kutta method of the form (3.1), the operator

A 0 in (2.6) is given for φ C V pR d , Rq by A 0 φ b T 1φ I f ¡ b T 1G ¡1 pg, f qφ I g ¡ σ 2 2 d 2 s G ¡1 divpgqφ I g σ 2 2 d 2 s ∆φ ¡ σ 2 2 d 2 s G ¡1 φ P pg, gq σ 2 d s ¢ p b T d ¡ p b T pδ dq 1 2 d s G ¡2 pg, g I gqφ I g σ 2 d s ¡ p b T pδ dq ¡ p b T d © G ¡1 φ I g I g.
In particular, if

b T 1 d s 1 and p b T d p b T pδ dq, (3.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

Y i x σ c hd i ξ hc i f pxq R h , if δ i 0, Y i x c h σd i ξ λ 1{2,i pxqgpxq $ h c i f pxq λ 1,i pxqgpxq σλ 1{2,i pxq s j1 p a ij d j g I pxqξ λ 1{2,i pxq s j1 p a ij λ 1{2,j pxqδ j g I pxqgpxq % R h , if δ i 1,
where the remainder satisfies § § R h § § ¤ Ch 3{2 , and where we used that λ i can be developed in powers of c h as λ i c 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

ζpY i q ζpxq c h σd i pg, ξq λ 1{2,i G $ h c i pg, f q λ 1,i G λ 1{2,i s j1 p a ij λ 1{2,j δ j pg, g I gq σλ 1{2,i s j1 p a ij d j pg, g I ξq 1 2 σ 2 d 2 i pξ, g I ξq σλ 1{2,i d i pg, g I ξq 1 2 λ 2 1{2,i pg, g I gq % . . .
where we omitted the dependency in x of G, g, g I and the λ k{2,j 's for brevity. We have ζpY i q ζpxq 0 (since x M), thus by identifying each term of the expansion with zero, we get

λ 1{2,i ¡σδ i d i G ¡1 pg, ξq, λ 1,i ¡δ i c i G ¡1 pg, f q σ 2 δ i £ s j1 p a ij d i d j d 2 i G ¡2 pg, ξqpg, g I ξq ¡ σ 2 δ i £ s j1 p a ij δ j d i d j 1 2 d 2 i G ¡3 pg, ξq 2 pg, g I gq ¡ σ 2 2 δ i d 2 i G ¡1 pξ, g I ξq.
For φ a test function, the operator A 0 φ satisfies ErφpX 1 qs ErφpY s qs φpxq hA 0 φpxq h 2 A 1 φpxq . . . By replacing Y s with its expansion in powers of h 1{2 , and by identifying the first terms, we deduce that

A 0 φ E c s φ I f ¡ c s G ¡1 pg, f qφ I g ¡ σ 2 2 d 2 s G ¡1 pξ, g I ξqφ I g σ 2 d s p p b T d d s qG ¡2 pg, ξqpg, g I ξqφ I g ¡ σ 2 d s ¢ p b T pδ dq 1 2 d s G ¡3 pg, ξq 2 pg, g I gqφ I g σ 2 2 d 2 s φ P pξ, ξq ¡ σ 2 d 2 s G ¡1 pg, ξqφ P pg, ξq σ 2 2 d 2 s G ¡2 pg, ξq 2 φ P pg, gq σ 2 d s p b T pδ dqG ¡2 pg, ξq 2 φ I g I g ¡ σ 2 d s p b T dG ¡1 pg, ξqφ I g I ξ % ,
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 yield 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.

Remark 3.7. The analysis presented in Section 3.1 is conducted for the overdamped Langevin dynamics (1.2). It would be interesting to consider extensions with multiplicative noise or a non-gradient vector field f . The calculations would likely become more involved and we may get more order conditions (see, for instance, [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF]Thm. 3.3] and [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]Remark 5.1 and Sect. 5.5] in the context of R d , where many additional terms arise, in particular for the integration by parts calculations). This will be studied in future work.

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.

Theorem 3.8 (Runge-Kutta conditions for order two for the invariant measure). We consider a Runge-Kutta method of the form (3.1) and assume the consistency condition (3.7). If the method is ergodic and if the following conditions are satisfied, then the integrator has order two for the invariant measure:

p b T d b T d, b T c b T pδ cq b T d 2 b T pδ d 2 q 2 p b T d ¡ 1 2 , p b T c p b T pδ cq p b T d 2 p b T pδ d 2 q p b T d 3 p b T pδ d 3 q 2 p b T d ¡ 1 2 , p b T pc dq p b T pδ c dq, b T pd p App1 ¡ δq dqqq 0, p b T Appδ ¡ 1q dqq p b T pδ Appδ ¡ 1q dqqq p p b T dq 2 ¡ 2 p b T d 1 2 , p b T pd p Acq p b T pd p Ad 2 q p b T pd p Apδ d 2 qq 2 p b T pd p Adq p p b T dq 2 ¡ 2 p b T d 1 2 , p b T pd 2 p Adq p b T pd p Adq 1 2 p p b T dq 2 , p b T pc p Appδ ¡ 1q dq p b T pd p Appδ ¡ 3 ¤ 1q dq p b T pd p Apδ cq 2p p b T dq 2 ¡ 4 p b T d 1, p b T pd 2 p Apδ dqq p b T pd p Apδ dqq 2 p b T pd p Adq 3 2 p p b T dq 2 ¡ 2 p b T d 1 2 , p b T pd p p App1 ¡ δq dqq 2 q 0, p b T pd p p Adq 2 q 3 p b T pd p Apd p App1 ¡ δq dqqq p4 ¡ 2 p b T dq p b T pd p Adq 3p p b T dq 2 ¡ 4 p b T d 1.
In the particular case where we set δ 1, the order two conditions reduce to the following:

p p b T dq 2 ¡ 2 p b T d 1 2 0, p b T d b T d, b T c b T d 2 p b T c p b T d 2 p b T d 3 2 p b T d ¡ 1 2 , p b T pd p Acq p b T pd p Ad 2 q 2 p b T pd p Adq, p b T pd 2 p Adq p b T pd p Adq p b T d ¡ 1 4 , p b T pd p p Adq 2 q p4 ¡ 2 p b T dq p b T pd p Adq 2 p b T d ¡ 1 2 .
For simplicity, we used in Theorem 3.8 the notation of Definition 3.5. For instance, the

condition p b T pd 2 p Adq p b T pd p Adq 1 2 p p b T dq 2 rewrites into s i,j1 p b i d 2 i p a ij d j d i,j1 p b i d i p a ij d j 1 2 ¡ d i1 p b i d i © 2
.

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

A 1 φ σ 4 8 ∆ 2 φ ¡ σ 4 4 G ¡1 ∆φ P pg, gq σ 4 8 G ¡2 φ p4q pg, g, g, gq Bφ, (3.8)
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µ V into an integral of the form ³ M A 0 1 φdµ V where A 0
1 φ is a differential operator of order one in φ (in the spirit of [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF][START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]). 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.

» M p∇ M ψ, Hqdσ M ¡ » M ψ div M pΠ M Hqdσ M ,
where ∇ M ψ : Π M ∇ψ and div M pHq : divpHq ¡ G ¡1 pg, H I gq. In addition, with the invariant measure dµ V ρ V dσ M and k ¥ 0, we obtain

» M G ¡k ψ I H ¡ G ¡pk 1q pg, Hqψ I g % dµ V » M G ¡pk 1q pg, H I gqψ (3.9) ¡ p2k 1qG ¡pk 2q pg, g I gqpg, Hqψ ¡ G ¡k divpHqψ 2kG ¡pk 1q pg, g I Hqψ G ¡pk 1q divpgqpg, Hqψ 2 σ 2 G ¡pk 1q pg, f qpg, Hqψ ¡ 2 σ 2 G ¡k pf, Hqψ % dµ V .
For example, let us integrate by parts the terms of order four w.r.t. φ of the operator A 1 φ in equation (3.8). Applying identity (3.9) with ψ σ 4 8 ∆φ I pe i q, H e i and k 0, and then summing on i 1, . . . , d yields

» M σ 4 8 ∆ 2 φ ¡ σ 4 8 G ¡1 ∆φ P pg, gq % dµ V » M ¡ σ 4 8 G ¡2 pg, g I gq∆φ I g (3.10) σ 4 8 G ¡1 divpgq∆φ I g σ 2 4 G ¡1 pg, f q∆φ I g ¡ σ 2 4 ∆φ I f % dµ V .
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 ∆φ P pg, gq ¡ σ 4 8 G ¡2 φ p4q pg, g, g, gq % dµ V (3.11) » M ¡ σ 4 4 G ¡1 i φ p3q pg, f i g, eiq σ 4 2 G ¡2 φ p3q pg, g, g I gq ¡ 3σ 4 8 G ¡3 pg, g I gqφ p3q pg, g, gq σ 4 8 G ¡2 divpgqφ p3q pg, g, gq σ 2 4 G ¡2 pg, f qφ p3q pg, g, gq ¡ σ 2 4 G ¡1 φ p3q pg, g, f q % dµ V .
Subtracting (3.11) from (3.10) allows to express

³ M A 1 φdµ
V with derivatives of φ of order strictly less than 4. We iterate this method to obtain

³ M A 1 φdµ V ³ M A 0 1 φdµ V where A 0 1
is an operator of order one in φ, and then find sufficient conditions such that A 0 1 0. This implies that A ¦ 1 ρ V 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:

b T d b T c b T pδ cq b T d 2 b T pδ d 2 q 1 2 , p b T d p b T c p b T pδ cq p b T d 2 p b T pδ d 2 q p b T d 3 p b T pδ d 3 q 1 2 , p b T pc dq p b T pδ c dq, p b T pd p Adq 1 8 , b T pd p App1 ¡ δq dqq 0, p b T App1 ¡ δq dq p b T pδ App1 ¡ δq dqq 1 4 , p b T pd p Acq p b T pd p Ad 2 q p b T pd p Apδ d 2 qq 0, p b T pd 2 p Adq 1 4 , p b T pc p App1 ¡ δq dqq ¡ p b T pd p Apδ dqq ¡ p b T pd p Apδ cqq

Illustrative examples of high order Runge-Kutta methods on manifolds

In this section, we present several examples of high order Runge-Kutta methods of the form (3.1). The purpose of these examples is to illustrate our analysis, and deriving new integrators with small error constant, favourable stability properties, small variance and fast convergence to equilibrium is a challenging open question which is not addressed in the present paper. First, we introduce a method that has order two for sampling the invariant measure of the constrained Langevin dynamics (1.2). Since there are many solutions to the order conditions, we obtain this integrator by solving numerically an optimization problem: we minimize the absolute values of the coefficients of the method under the constraints given by the order conditions of Theorem 3.8. This method is explicit in f and uses only three evaluations of f per step. It is defined by the following Butcher tableau

0 0 0 0 0 1 1 0 0 0 d 1 c 2 c 2 0 0 0 1 p a 21 p a 22 0 0 d 2 c 3 0 c 3 0 0 1 p a 31 p
Y 1 X n σ c hd 1 ξ n λ 1 gpY 1 q, Y 2 X n hc 2 f pY 1 q σ c hd 2 ξ n λ 2 rp a 21 gpY 1 q p a 22 gpY 2 qs , Y 3 X n hc 3 f pY 2 q σ c hd 3 ξ n λ 3 rp a 31 gpY 1 q p a 32 gpY 2 q p a 33 gpY 3 qs , X n 1 X n h 3 j1 p a 4j f pY j q σ c hξ n λ 4 3 j1 p a 4j gpY j q, (3.12)
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 12]. 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 pxqfpxq, where f is a gradient. Note that this equation can be rewritten as the following differential algebraic equation (DAE) of index two (see [31,

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 ¡ 1q{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:

b T d b T c b T pδ cq b T d 2 b T pδ d 2 q p b T d p b T c 1 2 , p b T pd p Adq 1 8 , p b T App1 ¡ δq dq 1 4 , p b T pd p Acq 0.
On the other hand, the order two conditions for the invariant measure of Theorem 3.8 on the sphere are the following:

p b T d b T d, b T c b T pδ cq b T d 2 b T pδ d 2 q p b T c 2 p b T d ¡ 1 2 , p b T pd p Acq 2 p b T pd p Adq p p b T dq 2 ¡ 2 p b T d 1 2 , p b T Appδ ¡ 1q dq p p b T dq 2 ¡ 2 p b T d 1 2 .
For example, the following integrator has order two for the invariant measure if M is a sphere:

Y 1 X n h ¢ 3 2 ¡ c 2 f pY 2 q σ c h ¢ 1 ¡ c 2 2 ξ n λ 1 p2Y 1 ¡ Y 2 q, ζpY 1 q 0, Y 2 X n hf pY 1 q σ c hξ n λ 2 Y 1 , ζpY 2 q 0, X n 1 Y 2 .
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

dX Π M pXq ¥ dW, Xp0q X 0 M. (3.14)
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

p b T d 1 2 , p b T pd p Adq 1 8 .
For example, a weak order two method for simulating a Brownian motion on a sphere is

X n 1 X n c hξ n λ 3X n c hξ n X n 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 [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF], 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 [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF][START_REF] Mclachlan | B-series methods are exactly the affine equivariant methods[END_REF][START_REF] Fløystad | The universal pre-Lie-Rinehart algebras of aromatic trees[END_REF] for their geometric properties, and in [START_REF] Chartier | Algebraic structures of B-series[END_REF][START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF] for their algebraic structure of Hopf algebras. In [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF], 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 E S where E 0 are the standard oriented edges as defined in [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF], and where E S is a new set of non-oriented edges represented as double horizontal straight lines. If pv, wq pw, vq 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 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 V g 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

|γ| N f pγq N l pγq N g pγq 2 ¡ |E S | . For instance, the graph γ pV, E, Lq with V f tv 2 , v 5 , v 6 u, V g tv 1 , v 3 , v 4 , v 7 u, E S tpv 6 , v 7 qu, (4.1)
E 0 tpv 1 , rq, pv 2 , v 1 q, pv 3 , rq, pv 4 , v 4 q, pv 5 , v 4 qu, L tpv 2 , v 2 q, pv 3 , v 5 q, pv 5 , v 6 qu, 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 ϕ :

V 1 Ñ V 2 such that ϕpV 1 f q V 2 f , ϕpV 1 g q V 2 g , pϕ ¢ ϕqpE 1 q E 2 , pϕ ¢ ϕqpE 1 S q E 2 S , pϕ ¢ ϕqpL 1 q L 2 .
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 with an exotic aromatic forest is given by the following definition. Definition 4.1. Let γ pV, E, Lq 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 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

F pγqpf, g, φq σ 2p|γ|¡N f pγqq G ¡Ngpγq{2 d iv 1 ,...,iv m 1 d jl 1 ,...,j ls 1 ¤ ¥ ¹ vV f δ iv,i Spvq f I πpvq f J Γpvq f iv ¤ ¤ ¥ ¹ vVg δ iv,i Spvq f I πpvq f J Γpvq g iv f I πprq f J Γprq φ.
For example, the differential associated with the exotic aromatic forest γ given by (4.1) is

F pγqpf, g, φq σ 8 G ¡2 d iv 1 ,...,iv 7 1 d jl 1 ,...,j l 3 1 f j l 1 j l 1 f iv 2 f j l 2 j l 3 f iv 5 δ iv 6 ,iv 7 f j l 3 f iv 6 ¤ f iv 2 g iv 1 f j l 2 g iv 3 f iv 4 iv 5 g iv 4 δ iv 7 ,iv 6 g iv 7 f iv 1 iv 3 φ.
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 γEAT 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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]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 [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]Thm. 4 If epγq apγq for all γ 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

Lφ φ I f ¡ G ¡1 pg, f qφ I g ¡ σ 2 2 G ¡1 divpgqφ I g σ 2 2 G ¡2 pg, g I gqφ I g σ 2 2 ∆φ ¡ σ 2 2 G ¡1 φ P pg, gq F ¡ ¡ ¡ 1 2 1 2 1 2 ¡ 1 2 © pφq.
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 [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF] 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.

τ γ : sum of all exotic aromatic forests/vector fields obtained by linking the root of τ to a node of γ with a new edge in E 0

2.

τ τ (resp. , where operation 1 is always applied first.

For example, let γ , τ and v r the root of γ, then we get

τ γ 2 , τ τ , τ , γ γv 2 , v v γ τ τ .
The integration by parts (3.9) can be rewritten conveniently with exotic aromatic forests.

Lemma 4.5. Let γ EAT and τ EAV, then the process of integration by parts rewrites into

» M F ¡ τ γ ¡ τ τ γ © pφqdµ V » M F ¡ τ τ γ ¡ pN g pγq N g pτq 1q τ τ γ ¡ τ γ (4.2) pN g pγq N g pτqq τ τ γ τ τ γ 2 τ τ γ ¡ 2 τ τ γ © pφqdµ V , » M F ¡ γ γv ¡ v v γ © pφqdµ V » M F ¡ ¡ pN g pγq 1q v v γ (4.3) N g pγq v v γ v v γ 2 v v γ ¡ 2 v v γ © pφqdµ V .
We write γ r γ if it is possible to go from γ EAT to r γ 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

1 8 ¡ 1 8 ¡ 1 8 1 8 1 4 ¡ 1 4 , 1 8 ¡ 1 8 ¡ 1 4 1 2 ¡ 3 8 1 8 1 4 ¡ 1 4 .
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 γ EAT and a vector field τ EAV with N g pγq N g pτq 0 and g 0, Lemma 4.5 reduces to the two following equations:

τ γ ¡ τ γ ¡ 2 τ τ γ , γ γv ¡2 v v γ .
We recover the process of integration by parts described in [START_REF] Laurent | Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs[END_REF]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 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 with an order one operator) or that have a rooted tree of the form , , , . . . Then, if γ EAT , there exists γ 0 EAT 0 such that γ γ 0 .

For instance, for a consistent of the form (3.1), the operator A 0 1 F pγ 0 1 q has the form

γ 0 1 pb T d ¡ p b T dq |γ|2 |πprq|1
a 0 pγqγ, so that γ 0 1 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 ¡ 1q{2. We choose the potential V pxq 25p1 ¡ x 2 1 ¡ x 2 2 q, with σ c 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 estimator s J 1 M °M m1 φpX pmq N q ErφpX N qs, where X pmq n 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µ V 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

°M m1 pφpX pmq N q¡ s Jq 2 ¨1{2 { M pM ¡ 1q.
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).

Figure 1: A trajectory of the order two method (left) and the convergence curve for the sphere for the invariant measure (right) with the potential V pxq 25p1 ¡ x 2 1 ¡ x 2 2 q, φpxq x 2 3 , a final time T 20 and M 10 7 trajectories.

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 x 2 3 R 2 ¡r 2 q 2 ¡4R 2 px 2 1 x 2 2 q with R 3 and r 1. The potential is V pxq 25px 3 ¡ rq 2 and we choose σ c 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µ V 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µ V
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 [START_REF] Zappa | Monte Carlo on manifolds: sampling densities and integrating functions[END_REF][START_REF] Zhang | Ergodic SDEs on submanifolds and related numerical sampling schemes[END_REF]) to see how it performs in high dimension. We choose the special linear group SLpmq tM 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µ V pxq is computed with the Runge-Kutta method (3.12) with h ref 2 ¡14 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, [START_REF] Goldfeld | Automorphic forms and L-functions for the group GLpn[END_REF]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 and observe that the Runge-Kutta method (3.12) performs significantly better than the Euler scheme (1.7 J M ¡1 °M k1 φ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 ¡12 T . where ³ M ρ n dσ M 0 for n ¥ 1. We denote ρ h r °r n0 h n ρ n and dµ h r ρ h r dσ M and adapt on the manifold M the following result from [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF]Thm. 2.1] in the context of R d . Lemma A.1. Under Assumptions 2.2, 2.3 and 2.4, for all φ C V pR d , Rq, for every positive integer r, there exists a constant C r pφq independent of h such that, for all h small enough,

§ § § § » M φdµ h ¡ » M φdµ h r § § § § ¤ C r pφqh r 1 . (A.4)
We omit the proof of Lemma A.1 as it is exactly the same as in [22, Thm. 2.1] by replacing dx by dσ M and T d by M. We are able to prove Theorem 2.5.

Proof of Theorem 2.5. As A ¦ j ρ V 0 for j 1, . . . , r ¡ 1, we deduce recursively from (A.2) and (A.3) that ρ j 0 for j 1, . . . , r ¡ 1, which yields ρ h r ρ V h r ρ r . Using the definition of the error for the invariant measure (1.5) and the ergodicity of the integrator (2.5), equation

(A.4) becomes § § § § § epφ, hq ¡ h r » M φpxqρ r pxqdσ M pxq § § § § § ¤ Ch r 1 .
We are left to prove that

» M φpxqρ r pxqdσ M pxq » V 0 » M upx, tqA ¦ r ρ V pxqdσ M pxqdt.
By the backward Kolmogorov equation and ergodicity, u satisfies lim

T ÑV upx, T q φpxq » V 0 Lupx, tqdt » M φpyqdµ V pyq. Using L ¦ ρ r ¡L ¦ r ρ V ¡A ¦ r ρ V , we deduce » M φpxqρ r pxqdσ M pxq ¡ » V 0 » M Lupx, tqρ r pxqdσ M pxqdt » M φpyqdµ V pyq » M ρ r pxqdσ M pxq ¡ » V 0 » M upx, tqL ¦ ρ r pxqdσ M pxqdt » V 0 » M upx, tqA ¦ r ρ V pxqdσ M pxqdt,
where we used that ³ M ρ r pxqdσ M pxq 0. This concludes the proof of Theorem 2.5.

B Integration by parts using the tree formalism

We provide here the detailed calculations of the integrations by parts of the order three terms, that are needed for the proof of Theorem 3.8. After applying the operations (3.10) and (3.11),

³ M A 1 φdµ V is transformed into ³ M Bφdµ V
where B is a differential operator of order three given by

Bφ F ¡ 1 4 ¡ 1 4 1 2 ¡ 1 4 1 4 1 8 ¡ 3 8 ¡ 1 4 ¡ 1 8 1 8 © pφq Rφ,
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φ.

1 4 ¡ 1 4 ¡ 1 4 1 4 ¡ 1 4 1 4 1 2 ¡ 1 2 ¡ 1 4 1 4 1 4 1 4 ¡ 1 4 ¡ 1 4 3 4 ¡ 1 2 ¡ 1 4 ¡ 1 2 1 2 ¡ 1 4 1 4 1 4 1 4 1 4 ¡ 3 4 ¡ 1 4 ¡ 1 4 3 4 ¡ 1 4 ¡ 1 2 1 2 ¡ 1 8 1 8 1 8 1 8 ¡ 1 8 ¡ 3 8 3 8 ¡ 1 8 ¡ 1 4 1 4 1 8 ¡ 1 8 ¡ 1 8 ¡ 1 4 ¡ 1

C Coefficients of the order two Runge-Kutta method

The coefficients of the Runge-Kutta method (3.12) Terms of order 3 w.r.t. φ Terms of order 2 w.r.t. φ φ P pf, f q Terms of order 1 w.r.t. φ

σ 2 ∆φ I f 1 2 1 2 σ 2 G ¡1 φ p3q pg, g, f q ¡ 1 2 ¡ 1 2 σ 4 G ¡2 φ p3q pg, g, g I gq 1 1 σ 4 G ¡1 °φp3q pg, g I e i , e i q ¡ 1 2 ¡ 1 2 σ 2 G ¡2 pg, f qφ p3q pg, g, gq 1 2 1 2 σ 4 G ¡2 divpgqφ p3q pg, g, gq
1 2 1 2 σ 2 °φP pf I e i , e i q 1 2 b T d σ 2 G ¡1 φ P pg, g I f q ¡1 ¡1 σ 2 G ¡1 φ P pg, f I gq ¡1 ¡b T d ¡ p b T d σ 4 G ¡2 φ P pg, g I g I gq 3 2 ¡2 p b T pd p Adq ¡ p p b T dq 2 2 p b T d 1
¡ p b T pd p Adq ¡ 1 2 p p b T dq 2 p b T d σ 4 G ¡2 φ P pg, g P pg, gqq 1 2 2 p b T d 2 ¡ 2 p b T pδ d 2 q 1 2 σ 4 G ¡1 φ P pg, ∆gq ¡ 1 2 ¡ 1 2 σ 4 G ¡1 °φP pg I e i , g I e i q ¡ 1 4 p b T pd p Adq 1 2 p p b T dq 2 ¡ p b T d σ 4 G ¡1 °φP pg P pg, e i q, e i q 0 p b T pδ d 2 q ¡ p b T d 2 G ¡1 pg, f qφ P pg, f q ¡1 ¡1 σ 2 G ¡1 divpgqφ P pg, f q ¡ 1 2 ¡ 1 2 σ 2 G ¡2 pg, g I gqφ P pg, f q 1 2 1 2 σ 2 G ¡1 pg, f q °φP pe i , g I e i q ¡ 1 2 ¡ p b T d σ 4 G ¡1 divpgq °φP pe i , g I e i q ¡ 1 4 ¡ 1 2 p b T d σ 4 G ¡2 pg, g I gq °φP pe i , g I e i q 1 4 1 2 p b T d σ 2 G ¡2 pg, f qφ P pg, g I gq 2 2 p b T d 1 σ 4 G ¡2 divpgqφ P pg, g I gq 1 p b T d 1 2 σ 4 G ¡3 pg, g I gqφ P pg, g I gq ¡ 5 2 2 p b T pd p Adq p p b T dq 2 ¡ 3 p b T d ¡ 3 2 σ 2 G ¡2 pg, g I f qφ P pg, gq 1 1 σ 2 G ¡2 pg, f I gqφ P pg, gq 1 2 b T d σ 4 G ¡2 °fij g i g j φ P pg, gq 1 2 1 2 σ 4 G ¡2 °fj g i f i g j φ P pg, gq 1 4 1 4 σ 4 G ¡3 pg, g I g I gqφ P pg, gq ¡ 7 4 p b T pd p Adq 1 2 p p b T dq 2 ¡ p b T d ¡ 3 2 σ 4 G ¡3 pg, g P pg, gqqφ P pg, gq ¡ 1 2 p b T pδ d 2 q ¡ p b T d 2 ¡ 1 2 G ¡2 pg, f q 2 φ P pg, gq 1 2 1 2 σ 2 G ¡2 pg, f q divpgqφ P pg, gq 1 2 1 2 σ 2 G ¡3 pg, f qpg, g I gqφ P pg, gq ¡2 ¡ p b T d ¡ 3 2
1 2 b T c 1 4 1 2 b T d 2 ¡ 1 2 p b T pd p Acq ¡ p b T d ¡ 1 2 p b T Appδ ¡ 1q dq ¡ p b T db T d 0 b T pd p Appδ ¡ 1q dqq ¡ 1 4 ¡ 1 2 b T pδ d 2 q 0 p b T pδ c dq ¡ p b T pc dq ¡ 1 4 1 2 p b T pd p Ad 2 q ¡ 1 2 p b T d ¡ p b T pd p p Adq 2 q ¡ 3 p b T pd p Apd p App1 ¡ δq dqqq 1 2 ¡2 p b T d p b T pd p Adq p p b T dq 2 p b T d 1 4 p b T pd p Ad 2 q ¡ 3 2 p b T pd p Apδ d 2 qq p b T dp p b T d 2 ¡ p b T pδ d 2 qq 1 2 p b T d 0 ¡2 p b T pd 2 p Adq p p b T d 1q p b T d 2 ¡ p b T d p b T pδ d 2 q 0 p b T pd 2 p Adq ¡ 1 2 p b T d 2 0 1 2 p b T pδ d 3 q ¡ 1 2 p b T d 3 ¡ 1 2 ¡ p b T c
¡ 1 4 ¡ 1 2 p b T c 1 4 1 2 p b T c ¡ 1 2 ¡b T pδ cq ¡ 1 4 ¡ 1 2 b T pδ d 2 q 1 4 b T pd p App1 ¡ δq dqqq 1 2 b T pδ d 2 q p b T pc p App1 ¡ δq dqq p p b T dq 2 p b T d 1 p b T pd p App1 ¡ δq dqq ¡ p b T pd p Apδ cqq 1 2 p b T pd p App1 ¡ δq dqq 1 2 p p b T dq 2 1 2 p b T d 1 2 1 2 p b T pd 2 p App1 ¡ δq dqq ¡ 1 2 p b T pd p Apδ d 2 qq 3 p b T pd p Apd p App1 ¡ δq dqqq p2 p b T d ¡ 1 2 q p b T pd p Adq ¡3 p b T pd p p App1 ¡ δq dqq 2 q 1 2 p b T pd p Apδ dqq p b T pd p p Adq 2 q ¡1 1 2 p b T pd 2 p Appδ ¡ 1q dqq 1 2 p b T pd p Apδ d 2 qq ¡ 3 2 p p b T dq 2 ¡ 3 2 p b T d ¡ 1 4 ¡ 1 2 p b T d 2 ¡ 1 8 ¡ 1 4 p b T d 2 1 8 1 4 p b T d 2 1 4 1 2 p b T pδ d 2 q 1 8 1 4 p b T pδ d 2 q ¡ 1 8 p b T pd 2 p Adq ¡ p p b T d 1 2 q p b T d 2 p p b T d ¡ 1 4 q p b T pδ d 2 q 1 2 ¡ p b T pd p Acq p b T d 1 2 p b T pδ App1 ¡ δq dqq p b T db T d 1 4 ¡ 1 2 p b T pd p Ad 2 q 1 2 p b T d 0 ¡ p b T pd 2 p Adq 1 2 p b T pδ d 2 q
¡ δq dqqq ¡ 1 2 2 p b T pd 2 p Adq 2 p b T d p b T pd p Adq ¡ p b T pδ d 2 q ¡ p p b T dq 2 ¡ p b T d ¡ 1 4 1 2 p b T pd p App3δ ¡ 2 ¤ 1q d 2 qq p b T dp p b T pδ d 2 q ¡ p b T d 2 q ¡ 1 2 p b T d 1 2 p b T pδ cq 1 2 1 2 p b T pδ d 2 q 1 2 p b T pδ cq p b T pc p Appδ ¡ 1q dqq p b T pd p Appδ ¡ 1q dqq ¡ 3 2 p b T pd p Apδ cqq ¡ 1 2 p b T pδ d 2 q ¡ 1 2 p b T pδ cq ¡ p p b T dq 2 ¡ p b T d 1 8 1 4 p b T pδ d 2 q 1 2 p b T pd p Appδ ¡ 1q dqq 1 2 p b T pd p Apδ d 2 qq ¡ 3 4 1 2 p b T pd 2 p Appδ ¡ 1q dqq ¡ 1 2 p b T pδ d 2 q ¡ 1 2 p p b T dq 2 ¡ 1 2 p b T d 3 p b T pd p p App1 ¡ δq dqq 2 q ¡ 1 2 p b T pd 2 p Appδ 1q dqq ¡ p b T pd p p Adq 2 q ¡ 3 p b T pd p Apd p App1 ¡ δq dqqq p 1 2 ¡ 2 p b T dq p b T pd p Adq ¡ 1 2 p b T pd p Apδ dqq ¡ 1 2 p b T pd p Apδ d 2 qq 9 8 9 2 p b T pδ d 3 q ¡ 15 4 p b T pδ d 2 q 3 2 p p b T dq 2 3 2 p b T d ¡ 1 2 ¡b T c ¡ 1 2 ¡ 1 2 ¡ 1 2 ¡ 1 2 ¡ 1 4 ¡ 1 2 b T d 2 ¡ 1 2 ¡b T d 3 2 ¡ p b T pd p Acq p b T d 1 3 2 b T pd p App1 ¡ δq dqq p b T App1 ¡ δq dq p p b T d 2qb T d 1 2 p b T pc dq ¡ p b T pδ c dq 1 2 1 4 1 2 b T pδ d 2 q 3 4 ¡ 1 2 p b T pd p Ad 2 q 1 2 p b T d 1 2
p b T d 3 ¡ 1 2 p b T pδ d 3 q 1 4 1 2 ¡ p b T pd 2 p Adq 3 2 p b T d 2 ¡ p b T pδ d 2 q 1 2 1 4 ¡ p b T pd p Adq ¡ 1 2 p p b T dq 2 p b T d p b T pd p p Adq 2 q 3 p b T pd p Apd p App1 ¡ δq dqqq ¡ 9 4 p2 p b T d 3q p b T pd p Adq 1 2 p p b T dq 2 ¡ 4 p b T d ¡ 1 1 2 p b T pd p App3δ ¡ 2 ¤ 1q d 2 qq 2 p b T pd 2 p Adq ¡ 7 4 ¡p2 p b T d 3q p b T d 2 2p p b T d 1q p b T pδ d 2 q ¡ 1 2 p b T d ¡ 3 2 ¡ 1 8 ¡ 1 8 ¡ 1 8 ¡ 1 8 3 2 p b T c 1 3 4 1 2 p b T c 1 2 ¡ 9 4 p b T pd p Acq ¡ 1 2 p b T c ¡ p b T d ¡ 3 2 1 2 b T pδ cq 1 4 1 2 b T pδ d 2 q p b T pδ Appδ ¡ 1q dqq ¡ 5 4 b T pd p Appδ ¡ 1q dqq ¡ 1 2 b T pδ d 2 q ¡ p p b T d 1qb T d 3 4 1 2 p b T d 2 1 2 3 8 1 4 p b T d 2 1 4 ¡ 9 8 1 2 p b T pd p Ad 2 q ¡ 1 4 p b T d 2 ¡ 1 2 p b T d ¡ 3 4 1 2 p b T d 1 4 1 2 p b T d ¡ 1 2 p b T pd 2 p Adq ¡ 1 2 p b T pδ d 2 q ¡ 1 2 p b T d ¡ 1

Definition 3 . 5 .

 35 For y, y p1q , . . . , y pnq R d and m ¥ 0, we define the diamond product and the diamond power as the vectors in R d , y p1q . . . y pnq ¢ n ¹ k1 y pkq i i and y m py m i q i .

Lemma 3 . 9 (

 39 Integration by parts on M). If ψ : R d Ñ R and H : R d Ñ R d are smooth functions, then

a 32 p a 33 0 d 3 1 p

 1 a 41 p a 42 p a 43 0 1 p a 41 p a 42 p a 43 0 1 p a 41 p a 42 p a 43 0 p a 41 p a 42 p a 43 0 or by the associated set of equations

  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

Figure 2 :

 2 Figure 2: A trajectory of the order two method (left) and the convergence curve for the torus for the invariant measure (right) with the potential V pxq 25px 3 ¡ rq 2 , φpxq x 2 3 , a final time T 20 and M 10 7 trajectories.

σ 4 G 4 ¡ 3 4σ 2 G 2 ¡ 1 2σ 4 G 4 ¡ 1 4σ 4 G

 4432214414 ¡2 pg, g I gqφ p3q pg, g, gq ¡ 3 ¡1 pg, f q∆φ I pgq ¡ 1 ¡1 divpgq∆φ I pgq ¡ 1 ¡2 pg, g I gq∆φ I pgq

¡ p b T pd p Adq ¡ 1 2 p p b T dq 2 3 2 p b T d 15 8

 28

4 p 1 Table 2 (

 412 b T pc p Appδ ¡ 1q dqq p b T pd p Apδ cqq ¡3 p b T pd p Appδ ¡ 1q dqq ¡ p p b T dq 2 ¡ 3 p b T d ¡ Part 6/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1).

p 3 2p b T d ¡ 1 2 3 2 p b T pd 2 p 2 p 2 p b T pd p Apδ dqq 23 4 2 p b T d 3 ¡ 3 4 ¡ 1 2p b T pδ d 2 p b T pd p App2 ¤ 1 ¡ 8 p2 p b T d 3 2 q p b T d 2 ¡ p2 p b T d 3 4 q p b T pδ d 2 q 1 2 p b T d 5 4 ¡1 ¡ p b T pδ cq ¡ 1 2 ¡1 ¡ 1 2 p b T pδ cq ¡ 1 2 p b T pδ d 2 q ¡ 1 2 p 2 ¡ p b T pd p Apδ cqq 1 2 p b T pδ cq 1 2 p 2 ¡ 1 4 ¡ 1 4p b T pδ d 2 q ¡ 1 8 1 2 p 2 p b T pd p Apδ d 2 qq 7 4 1 2 p b T pd p App1 ¡ δq dqq 1 2 p 4 p 2 p b T pd 2 pδq dqq 1 2 p b T pd p Apδ d 2 qq ¡ 9 2 p b T pδ d 3 4 p 8 Table 2 ( 2 1 2 b 4 p 2 p 2 b 4 p 2 p 4 ¡ 1 1 2 p b T pd p App2 ¤ 1 ¡ 4 ¡2 p b T pd 2 p 2 p b T pd 2 p 2 p 4 1 2 p b T pδ d 3 q ¡ 1 2 p 2 ¡ 1 2p b T c p b T d ¡ 1 4 3 (

 3222224233412182242222222214182242242222348222422424121422224222113 b T pd p Appδ ¡ 1q dqq ¡ 1 2 p p b T dq 2 ¡ p b T pd p p App1 ¡ δq dqq 2 q ¡ 1 App3 ¤ 1 δq dqq ¡2 p b T pd p p Adq 2 q ¡ 6 p b T pd p Apd p App1 ¡ δq dqqq ¡ 1 b T pd p Apδ d 2 qq ¡ 1 ¡p4 p b T d 5 2 q p b T pd p Adq p b T pδ d 2 q p p b T dq 2 13 3δq d 2 qq ¡ p b T pd 2 pAdq 13 b T pc p App1 ¡ δq dqq p b T pd p App1 ¡ δq dqq7 b T pδ d 2 q p p b T dq 2 2 p b T d 3 b T pd 2 p App1 ¡ δq dqq ¡ 1 b T pδ d 2 q 1 2 p p b T dq 2 p b T d 3 b T pd p p Adq 2 q 3 p b T pd p Apd p App1 ¡ δq dqqq ¡3 p b T pd p p App1 ¡ δq dqq 2 q p2 p b T d 1 2 q p b T pd p Adq 1 App1 b T pδ d 2 q ¡ p p b T dq 2 ¡ 3 p b T d ¡ 15Part 7/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1).γ a 0 pγq b T d ¡ p b T d b T c ¡ 2b T d 1 T d 2 ¡ b T d 1 b T pd p Acq ¡ 2 p b T pd p Adq ¡ p p b T dq 2 2 p b T d ¡ 1 b T Appδ ¡ 1q dqq ¡ p b T db T d 2 p b T d ¡ 1 T pd p Appδ ¡ 1q dqq ¡ 1 2 b T pδ d 2 q b T d ¡ 1 b T pδ c dq ¡ p b T pc dq 2 p b T d 2 ¡ 2 p b T pδ d 2 q 1 b T pd p Ad 2 q ¡ p b T pd p Adq ¡ 1 2 p p b T dq 2 p b T d ¡ 1 p b T pd p p Adq 2 q ¡ 3 p b T pd p Apd p App1 ¡ δq dqqq p4 ¡ 2 p b T dq p b T pd p Adq 3p p b T dq 2 ¡ 4 p b T d 3δq d 2 qq p b T pd p Adq p p b T d ¡ 1qp p b T d 2 ¡ p b T pδ d 2 qq 1 2 p p b T dq 2 ¡ p b T d 1 Adq 2 p b T pd p Adq p p b T d ¡ 2q p b T d 2 p3 ¡ p b T dq p b T pδ d 2 q p p b T dq 2 ¡ 2 p b T d 1 Adq ¡ p b T pd p Adq 1 b T d 2 ¡ p b T pδ d 2 q ¡ 1 2 p p b T dq 2 p b T d ¡ 1 b T d 3 p b T d 2 ¡ p b T pδ d 2 q ¡ p b T d 2 p b T pδ d 2 q ¡ p b T c 2 p b T d ¡ 1Table Part 1/5): Coefficients in exotic aromatic B-series of the operator A 0 1 φ °a0 pγqFpγqpφq for consistent Runge-Kutta methods of the form (3.1).

γ a 0 pγq 1 2 p b T c ¡ p b T d 1 4 ¡b T pδ cq 2b T d ¡ 1 2 ¡ 1 2 b 4 b 4 p b T pd p App3 ¤ 1 ¡ 2 p 2 p b T pd p Apδ d 2 qq 1 2 p 2 p b T pd 2 p App1 ¡ δq dqq 1 2 p b T pd p Apδ dqq 1 2 p 2 ¡ 1 2p b T d 2 p b T d ¡ 1 4 ¡ 1 4p b T d 2 1 2 p b T d ¡ 1 8 1 4 p b T d 2 ¡ 1 2 p b T d 1 8 ¡2 p b T d 2 5 2 p 4 5 4 p 2 p b T d 1 8 p 5 2 1 2 p p b T dq 2 3 2 p b T d ¡ 3 8 ¡ 2 p 2 ¡ 1 2p b T pd p Ad 2 q p b T pd p Adq 1 2 p p b T dq 2 ¡ p b T d 1 4 ¡ p b T pd 2 p 2 p b T pδ d 2 q 1 2 p p b T dq 2 ¡ p b T d 1 4 3 (

 242244122222221412842824428512822112422213 T pδ d 2 q b T d ¡ 1 T pd p App1 ¡ δq dqq 1 2 b T pδ d 2 q ¡ b T d 1 δq dqq ¡ p b T pd p Apδ cqq p b T pc p App1 ¡ δq dqq 2p p b T dq 2 ¡ 4 p b T d 1 1 b T pd p App3 ¤ 1 ¡ δq dqq ¡ 1 b T pd 2 p App1 ¡ δq dqq p p b T dq 2 ¡ 2 p b T d 1 2 p b T pd p p Adq 2 q 3 p b T pd p Apd p App1 ¡ δq dqqq p2 p b T d ¡ 11 2 q p b T pd p Adq ¡3 p b T pd p p App1 ¡ δq dqq 2 q ¡ 1 b T pd p Apδ d 2 qq ¡ 4p p b T dq 2 6 p b T d ¡ 3 b T pδ d 2 q ¡ p b T d 1 b T pδ d 2 q ¡ p b T d 2 ¡ 1 ¡ p b T dq p b T d 2 p p b T d ¡ 13 4 q p b T pδ d 2 q p b T pd 2 p Adq ¡ p b T pd p Adq ¡ p b T pd p Acq 2 p b T pd p Adq p p b T dq 2 ¡ 2 p b T d 1 b T pδ App1 ¡ δq dqq p b T db T d ¡ 2 p b T d 1 Adq p b T pd pAdq 1 Table Part 2/5): Coefficients in exotic aromatic B-series of the operator A 0 1 φ °a0 pγqFpγqpφq for consistent Runge-Kutta methods of the form (3.1).

γ a 0 2 1 2 p b T pd p App3δ ¡ 2 ¤ 4 p 2 1 2 p b T pδ d 2 q 1 2 p 2 p 2 1 4 p 2 p b T d 1 8 1 2 p 2 p b T pd p Apδ d 2 qq 1 2 p b T pd 2 p 2 p 4 3App1 ¡ δq dqqq ¡ 1 2 p b T pd 2 p 2 p b T pd p Apδ d 2 qq 9 2 p b T pδ d 3 q ¡ 1 2 p b T pd p Apδ dqq ¡ 15 4 p b T pδ d 2 q 9 2 p p b T dq 2 ¡ 15 2 p b T d 15 8 ¡b T c 2 p b T d ¡ 1 2 ¡ 1 2 b 2 b 2 p 1 2 b T pδ d 2 q ¡ p b T d 1 4 ¡ 1 2p b T pd p Ad 2 q p b T pd p Adq 1 2 p p b T dq 2 ¡ p b T d 1 4 1 2 p b T d 3 ¡ 1 2 p

 022242222242822222422222428222214112422 pγq p b T pd p p Adq 2 q 3 p b T pd p Apd p App1 ¡ δq dqqq 2 p b T pd 2 p Adq p2 p bT d ¡ 6q p b T pd p Adq ¡ p b T pδ d 2 q ¡ 4p p b T dq 2 6 p b T d ¡ 3 1q d 2 qq ¡ p b T pd p Adq p1 ¡ p b T dqp p b T d 2 ¡ p b T pδ d 2 qq ¡ 1 2 p p b T dq 2 p b T d ¡ 1 b T pδ cq ¡ 2 p b T d 1 b T pδ cq ¡ 2 p b T d 1 b T pc p Appδ ¡ 1q dqq p b T pd p Appδ ¡ 3 ¤ 1q dqq p b T pd p Apδ cqq ¡ 1 2 p b T pδ d 2 q ¡ 1 2 p b T pδ cq ¡ 2p p b T dq 2 6 p b T d ¡ 3 b T pδ d 2 q ¡ 1 b T pd p Appδ ¡ 3 ¤ 1q dqq 1 Appδ ¡ 1q dqq ¡ 1 b T pδ d 2 q ¡ p p b T dq 2 3 p b T d ¡ 3 p b T pd p p App1 ¡ δq dqq 2 q p 13 2 ¡ 2 p b T dq p b T pd p Adq ¡ p b T pd p p Adq 2 q ¡ 3 p b T pd p Apd p Appδ 1q dqq ¡ 1 T d 2 p b T d ¡ 1 4 ¡b T d p b T d ¡ p b T pd p Acq 2 p b T pd p Adq p p b T dq 2 ¡ 2 p b T d 1 T pd p App1 ¡ δq dqq p b T App1 ¡ δq dq p p b T d 2qb T d ¡ 4 p b T d 1 b T pc dq ¡ p b T pδ c dq ¡ 2 p b T d 2 2 p b T pδ d 2 q b T pδ d 3 q ¡ p b T d 2 p b T pδ d 2 q Table 3 (Part 3/5): Coefficients in exotic aromatic B-series of the operator A 0 1 φ °a0 pγqFpγqpφq for consistent Runge-Kutta methods of the form (3.1).

  there exists a constant C ¡ 0 such that for all x N M and py, zq ψpxq, for all one-form field v :T M Ñ R on M of norm Σ ik py, zq r Σ jk py, zqv i pr xqv j pr xq ¥ C,

	one, we have
	d¡1 i,j1	d ķ1	r
	where r

  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,
	upx, hq φpxq

  App1 ¡ δq dqq 2 q 0, p b T pd p p Adq 2 q 3 p b T pd p

	p b T pd 2 p Apδ dqq p b T pd p Apδ dqq 1 8 ,	1 8 ,
	Apd p App1 ¡ δq dqqq 1 8 .

p b T pd p p 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.

  .1]. 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,

	Proposition 4.3. and where the operators are given by
	L j j! F	¢	γEAT ,|γ|j	epγqγ	, A j¡1 F	¢	γEAT ,|γ|j	apγqγ	.

  τ τ ): aroma obtained by linking the root of τ to a white node (resp. a black node) with a new edge in E S

	3.	τ : sum of all aromas obtained by linking the root of τ to a node of τ with a new edge
		in E 0
	4. γ γv : sum of all exotic aromatic forests/vector fields obtained by linking the node v to a
		node of γ with a new liana
		τ τ
	5. v v γ

: forest obtained by linking the root of τ to the node v of γ with a new edge in E 0

For simplicity, we combine multiple operations on a same forest as in v v γ and τ τ

Table 1 :

 1 ). Numerical approximation of the integral Jpmq ³

	m	dimpSLpmqq	Jpmq	s J Euler	error for s J Euler	s J 2	error for s J 2
	2 3 4	3 8 15	2.00967 3.01954 4.02930	2.01031 3.02068 4.03095	6.4 ¤ 10 ¡4 1.1 ¤ 10 ¡3 1.6 ¤ 10 ¡3	2.00962 3.01934 4.02907	4.4 ¤ 10 ¡5 2.0 ¤ 10 ¡4 2.3 ¤ 10 ¡4
					SLpmq φpxqdµ V for 2 ¤ m ¤ 4 with the
	estimator s

D Decomposition of the operators in exotic aromatic forests

  a 22 1 ¡ p a 21 , p a 43 1 ¡ p a 41 ¡ p a 42 . a 33 1 ¡ p a 31 ¡ p a 32 ,

	Forest γ	Differential F pγqpφq Exact epγq	Numerical approximation apγq
		Terms of order 4 w.r.t. φ
		σ 4 ∆ 2 φ	1 8	1 8
		σ 4 G ¡1 ∆φ P pg, gq	¡ 1 4	¡ 1 4
		σ 4 G ¡2 φ p4q pg, g, g, gq	1 8	1 8
				used in Section 5 are
		c 2 0.621729189582953540,		c 3 0.102032386582165330,
		d 1 ¡0.898931652839146019,		d 2 ¡1.66233102561284629,
		d 3 0.318924515019668897,	p a 21 0.584372887990673524,
		p a 31 0.887706593835748395, p a 41 0.0547449506054026516,	p a 32 ¡0.345018694936693742, p a 42 ¡0.0205123070437693053,

p p

Table 2 (

 2 Part 1/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1).

	γ	F pγqpφq	epγq	apγq
		σ 2 G ¡2 φ P pg I g, g I gq	1 4

Table 2 (

 2

Part 2/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1).

Table 2 (

 2 Part 3/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1).

	γ	epγq	apγq

Table 2 (

 2 Part 4/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1). T pd p p Adq 2 q 3 p b T pd p

	γ	epγq	apγq
			Apd p App1

p b

Table 2 (

 2 Part 5/7): Coefficients in exotic aromatic B-series of the operators L 2 φ{2 °epγqF pγqpφq and A 1 φ °apγqF pγqpφq for consistent Runge-Kutta methods of the form (3.1).

	γ	epγq	apγq
		1	1
		4	2

Université de Genève, Section de mathématiques, UNI DUFOUR,

rue du Général Dufour, CP 64, 1211 Geneva 4, Switzerland. Adrien.Laurent@unige.ch, Gilles.Vilmart@unige.ch.[START_REF] Abdulle | Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations[END_REF] The extension to manifolds of any codimension q ¥ 1 will be further discussed in Remarks 2.6 and 3.1.

Acknowledgements. This work was partially supported by the Swiss National Science Foundation, grants No. 200020_184614, No. 200021_162404 and No. 200020_178752. The computations were performed at the University of Geneva on the Baobab cluster using the Julia programming language.

Appendices A Proof of Theorem 2.5

In the spirit of backward error analysis for differential equations (see [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Abdulle | High weak order methods for stochastic differential equations based on modified equations[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF]), we build a modified generator L h such that U px, hq ErφpX 1 q|X 0 xs formally satisfies U px, hq j¥0 h j j! pL h q j φpxq. (A.1) Truncating this formal series yields an estimate of the form U px, hq φpxq

where

For this, we write formally L h L °n¥1 h n L n and compare the series expression in (2.6) and (A.1). By formally identifying the powers of h, we deduce the following rigorous definition of the

where the B l are the Bernoulli numbers (see [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] 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 γ a 0 pγq