Adrien Laurent

Gilles Vilmart

Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs

Keywords: stochastic differential equations, invariant measure, ergodicity, exotic aromatic trees, order conditions. AMS subject classification (2010): 60H35, 37M25, 65L06, 41A58

We introduce a new algebraic framework based on a modification (called exotic) of aromatic Butcher-series for the systematic study of the accuracy of numerical integrators for the invariant measure of a class of ergodic stochastic differential equations (SDEs) with additive noise. The proposed analysis covers Runge-Kutta type schemes including the cases of partitioned methods and postprocessed methods. We also show that the introduced exotic aromatic B-series satisfy an isometric equivariance property.

Introduction

We consider a class of stochastic systems of differential equations of the form dXptq f pXptqqdt σdW ptq, (1.1) where Xptq R d is the solution with initial condition X 0 assumed deterministic for simplicity, the vector field f : R d Ñ R d is assumed smooth and globally Lipschitz, σ ¡ 0 is a constant, and W ptq is a standard d-dimensional Wiener process fulfilling the usual assumptions.

We say that problem (1.1) is ergodic if it has a unique invariant measure µ satisfying for all deterministic initial conditions X 0 and all smooth test functions φ,

lim T ÑV 1 T » T 0 φpXpsqqds » R d φpxqdµpxq, almost surely. (1.2)
Under appropriate smoothness and growth assumptions on the vector field f , the above ergodicity property is automatically satisfied, and in addition one has in general the following exponential convergence for all initial conditions X 0 and all appropriate test functions φ,

§ § § § ErφpXptqqs ¡ » R d φpxqdµpxq § § § § ¤ Ce ¡λt (1.3)
for all t ¡ 0 where C Cpφ, X 0 q and λ ¡ 0 are independent of time t ¡ 0. We refer to [START_REF] Bally | The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function[END_REF][START_REF] Mattingly | Convergence of numerical timeaveraging and stationary measures via Poisson 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] for details. A special case of the class of problems (1.1) is the well-studied overdamped Langevin equation, also called Brownian dynamics, 1 dXptq ¡∇V pXptqqdt σdW ptq (1.4) where V : R d Ñ R is a smooth potential, f ¡∇V and σ ¡ 0. Under appropriate growth conditions on the potential V (growing at least quadratically), (1.4) is ergodic and (1.2) and

(1.3) hold with the invariant measure dµpxq ρ V pxqdx where its density ρ V with respect to the Lebesgue measure is given by the Gibbs distribution

ρ V pxq Ze ¡ 2 σ 2 V pxq (1.5)
where Z is a normalization constant such that

³ R d ρ V pxqdx 1.
Computing integrals with respect to the invariant measure µ in high dimensions is in general very costly using a deterministic quadrature rule, and one can take advantage of the above ergodicity property to compute such integrals using numerical approximations of (1.1), the exact solution of (1.1) not being available in general. Let us mention that a natural way to sample without any bias from the invariant measure with the Gibbs density (1.5) is to apply Markov Chain Monte-Carlo methods, in particular Metropolis-Hastings algorithms, see for instance the survey [START_REF] Sanz-Serna | Markov chain Monte Carlo and numerical differential equations[END_REF]. However, as highlighted in [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF], "the Metropolis-Hastings algorithm is rather expensive due to the need of accept/reject steps and does not admit the use of powerful weak methods", which can be combined with the methodology of "rejecting exploding trajectories" [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF]. We thus focus in this paper on classes of weak methods and consider a one step numerical integrator for the approximation of (1.1) at time t n nh of the form X n 1 ΨpX n , h, ξ n q, (1.6) where h is a fixed time step size and ξ n are independent random vectors. We say that the numerical scheme has local weak order r if the weak error after one step satisfies

|ErφpX 1 q|X 0 xs ¡ ErφpXphqq|Xp0q xs| ¤ Ch r 1 , (1.7)
where C Cpφ, xq is independent of h assumed small enough and φ is a test function.

Note that under appropriate assumptions on the numerical scheme (to achieve in particular bounded numerical moments along time), one can in general deduce a global weak order r, |EpφpX n qq ¡ EpφpXpt n qqq| ¤ Ch r , as shown in [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2.2]). The numerical method is called ergodic if it has a unique invariant probability law µ h with finite moments of any order and lim

N ÑV 1 N 1 N ņ0 φpX n q » R d
φpxqdµ h pxq, almost surely, for all deterministic initial conditions X 0 x and all test functions φ. We say that the numerical method has order p with respect to the invariant measure of (1.1)

if § § § § » R d φpxqdµ h pxq ¡ » R d φpxqdµpxq § § § § ¤ Ch p ,
(1.8)

1 Up to a time transformation, one could for simplicity fix the value of σ (e.g. σ 1 or ?

2), however we choose not to fix it to better distinguish the various order conditions.

where C is independent of h assumed small enough. Under appropriate assumptions on the numerical scheme (see for instance [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF]), one also obtains the following exponential estimate similar to (1.3) (possibly with a different constant λ ¡ 0),

§ § § § ErφpX n qs ¡ » R d φpxqdµpxq § § § § ¤ Ke ¡λtn Ch p , (1.9)
where K Kpφ, xq, C Cpφq are independent of n and h assumed small enough. A simple way to achieve high order p for the invariant measure is to consider a numerical scheme with high standard weak order r, and it is known for large classes of SDEs that p ¥ r, see in particular [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF] in the context of locally Lipschitz vector fields with multiplicative noise. Note analogously that the strong order q of convergence, which corresponds to the numerical approximation of individual trajectories of (1.1), is in general lower than or equal to the weak order r of convergence. There are interestingly many schemes in the literature for which p ¡ r and a high order p for the invariant measure is obtained, while the standard weak order of accuracy remains low, typically of order r 1, i.e. the scheme is consistent in the weak convergence sense. This is the case in particular for the Langevin equation [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] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF][START_REF] Abdulle | Long time accuracy of Lie-Trotter splitting methods for Langevin dynamics[END_REF]. In [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], a methodology for the analysis and design of high order integrators for the invariant measure is introduced and serves as a crucial ingredient in this paper. The approach combines the usual Talay-Tubaro methodology [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] and recent developments of the theory of backward error analysis and modified differential equations in the stochastic context [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], a major tool in the area of deterministic geometric numerical integration [START_REF] Hairer | Geometric numerical integration[END_REF]. In [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF] for finite dimensions and in [START_REF] Bréhier | High order integrator for sampling the invariant distribution of a class of parabolic stochastic PDEs with additive space-time noise[END_REF] in the context of parabolic stochastic partial differential equations, this approach is combined with the idea of processing from Butcher [START_REF] Butcher | The effective order of Runge-Kutta methods[END_REF], to design efficient postprocessed integrators with high order for the invariant measure at a negligible overcost compared to standard low order schemes. The postprocessor methodology is extended in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] for a class of explicit stabilized schemes of order two for the invariant measure and with optimally large stability domains. The aim of this paper is to provide a unified algebraic framework based on aromatic trees and B-series, with a set of trees independent of the dimension d of the problem, for the systematic study of the order conditions for the invariant measure of a class of numerical integrators that includes Runge-Kutta type schemes for problems of the form (1.1). We show that the new framework permits to recover some schemes and simplify the calculations in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] and for postprocessed integrators in [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF][START_REF] Bréhier | High order integrator for sampling the invariant distribution of a class of parabolic stochastic PDEs with additive space-time noise[END_REF][START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]. Analogously to [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] (we study here the additive noise case), we consider in this paper Runge-Kutta methods of the form2

Y i X n h s °j1 a ij f pY j q l °k1 d pkq i σ c hξ pkq n , i 1, . . . , s, X n 1 X n h s °i1 b i f pY i q σ c hξ p1q n , (1.10) where a ij , b i , d pkq i
are the coefficients defining the Runge-Kutta scheme, and ξ pkq n N p0, I d q are independent Gaussian random vectors. We highlight once again that we focus in this paper on the high order p of accuracy for the invariant measure, while the order r of accuracy in the weak sense can remain low (typically r 1). The analysis in this paper applies to the class of methods (1.10) for any number l of random vectors in the internal stages. However we shall often consider l 1 random vector per internal stage3 , which is sufficient to achieve order p 2 or 3 for the invariant measure. In particular, we shall consider the θ-method as an illustrative example in Sections 4 and 5 and recover known results on its accuracy. It is defined for θ fixed as

X n 1 X n hp1 ¡ θqf pX n q hθf pX n 1 q σ c hξ n , (1.11)
For θ 0, we get the explicit Euler-Maruyama method while the scheme is implicit for θ $ 0. It can be put in Runge-Kutta form (1.10) for l 1 with the following coefficients.

c A d b 0 0 0 0 1 1 ¡ θ θ 1 1 ¡ θ θ
The usage of trees and B-series4 is known as a powerful standard tool for the numerical analysis of differential equations. B-series where introduced by Hairer and Wanner in [START_REF] Hairer | On the Butcher group and general multi-value methods[END_REF] based on the work of Butcher [START_REF] Butcher | An algebraic theory of integration methods[END_REF], and are now exposed in many articles and books [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF], see also the presentation in [START_REF] Chartier | Algebraic structures of B-series[END_REF][START_REF] Sanz-Serna | Formal series and numerical integrators: some history and some new techniques[END_REF]. In the last decades, several works extended trees and 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] 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, and [START_REF] Burrage | Order conditions of stochastic Runge-Kutta methods by B-series[END_REF][START_REF] Rößler | Stochastic Taylor expansions for the expectation of functionals of diffusion processes[END_REF][START_REF] Rößler | Runge-Kutta methods for Itô stochastic differential equations with scalar noise[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] Debrabant | B-series analysis of stochastic Runge-Kutta methods that use an iterative scheme to compute their internal stage values[END_REF][START_REF] Rößler | Strong and weak approximation methods for stochastic differential equations-some recent developments[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. Tree series were also used to describe schemes preserving quadratic invariants [START_REF] Anmarkrud | Order conditions for stochastic Runge-Kutta methods preserving quadratic invariants of Stratonovich SDEs[END_REF]. The recent work [START_REF] Alamo | A technique for studying strong and weak local errors of splitting stochastic integrators[END_REF] also studies algebraically strong and weak errors, but instead of using tree series, it uses word series because these are well suited in the context of splitting stochastic integrators.

In this paper, we focus on the long time accuracy of numerical integrators and derive in a systematic manner the order conditions for sampling the invariant measure of an ergodic system of the form (1.1). Additionally, in Section 5.5, we allow the inclusion of nonreversible perturbation as in [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]. The proposed algebraic framework relies on aromatic B-series, a generalisation of B-series introduced in [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF][START_REF] Mclachlan | B-series methods are exactly the affine equivariant methods[END_REF] (see also the presentation in [START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF]) to characterize all the schemes that are affine equivariant, i.e. that behave transparently with respect to an affine change of coordinates. These aromatic B-series rely on aromatic trees, which were first introduced in [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF] to represent the divergence of B-series in the context of deterministic value or first integral preserving ordinary differential equations.

This paper is organized as follows. In Section 2, we described the general setting and assumptions needed in our analysis. In Section 3, we introduce a new generalization of B-series, called "exotic aromatic B-series" by considering an additional new type of edge called "liana" compared to standard aromatic B-series. We also show that these new exotic aromatic B-series satisfy an isometric equivariance property. In Section 4, we explain how this new algebraic framework applies for the long time accuracy analysis of stochastic integrators for ergodic problems. In Section 5, we derive order conditions for integrators expandable as aromatic B-series methods, with special emphasis on Runge-Kutta type integrators and post-processed integrators. In particular, we show that the orders 2 and 3 for the invariant measure of Brownian dynamics (1.4) yield respectively 2 and 6 order conditions (see Table 1), compared to the 3 and 10 more restrictive conditions for the standard weak order of convergence (see Table 2).

Preliminaries

We first state the following smoothness and growth Assumptions 2.1 and 2.2 on the vector field of (1.1) which automatically yield that (1.1) satifies the ergodicity properties (1.2) and (1.3) (see [START_REF] Hasminskii | Stochastic stability of differential equations[END_REF] in the more general context of SDEs with multiplicative noise).

Assumption 2.1. The vector field f : R d Ñ R d is globally Lipchitz and C V , and there exist

C 1 , C 2 ¡ 0 such that for all x R d , x T f pxq ¤ ¡C 1 x T x C 2 .
The following stronger assumption yields the important special case of Brownian dynamics (1.4).

Assumption 2.2. The vector field f is a globally Lipschitz gradient, i.e. there exists a C V potential V : R d Ñ R such that f pxq ¡∇V pxq is globally Lipschitz and there exist

C 1 ¡ 0 and C 2 such that for all x R d , V pxq ¥ C 1 x T x ¡ C 2 .
We recall that under Assumption 2.2, the density of the unique invariant measure is

given by ρ V Z exp ¡ 2V σ 2 ¨where Z is such that ³ R d ρ V pxqdx 1 but Z is not numerically known in general. We note that ∇ρ V ρ V 2 σ 2 f , equivalently ∇plog ρ V q 2 σ 2 f .
We denote C V P pR d , Rq the vector space of C V functions such that all partial derivatives φ up to all orders have a polynomial growth of the form |φpxq| ¤ Cp1 |x| s q for some constants s and C independent of x (but depending on the order of differentiation).

For

where the generator L is defined as

Lφ f ¤ ∇φ σ 2 2 ∆φ (2.2)
where ∆φ °d i1 f 2 φ fx 2 i denotes the Laplace operator. We recall that, under Assumption 2.1 or 2.2, the density of the invariant measure satisfies

L ¦ ρ V 0 where L ¦ φ ¡divpfφq σ 2 2 ∆φ is the L 2 -adjoint of L.
We make the following natural assumptions on the numerical integrator (1.6). In the following theorem, we recall the characterisation of order p for the invariant measure of an ergodic integrator in terms of adjoints of the linear differential operators of Assumption 2.4. It was shown in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] in the more general context of ergodic SDEs with multiplicative noise based on backward error analysis results in [START_REF] Debussche | Weak backward error analysis for SDEs[END_REF] on the torus and generalizations [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] then the scheme has order p for the invariant measure and (1.8)-(1.9) hold. Remark 2.7. Assuming in addition that the scheme has weak order p ¡ 1 of accuracy, i.e. assuming in addition the stronger assumption A j L j 1 {pj 1q!, j 2, . . . , p ¡ 2, then Theorem 2.6 is an immediate consequence of the Talay-Tubaro expansion of the error [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF][START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] (see also [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF]Chap. 2.2,2.3]) given by

» R d φpxqdµ h pxq ¡ » R d φpxqdµpxq λ p h p Oph p 1 q with λ p » V 0 » R d ¢ A p ¡ 1 pp 1q! L p 1 upy, tqρ V pyqdydt
where upx, tq is the solution of (2.1). Indeed, considering the L 2 -adjoint of the operator A p ¡ 1 pp 1q! L p 1 and using L ¦ ρ V 0 and A ¦ p ρ V 0 then yield λ p 0 and the scheme has (at least) order p for the invariant measure.

The following extension of Theorem 2.6 permits to combine an integrator (1.6) with a postprocessor to achieve high order for the invariant measure at a negligible overcost compared to a standard scheme. Theorem 2.8. [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF] Assume the hypotheses of Theorem 2.6 and consider a postprocessor X n G n pX n q that admits the following weak Taylor expansion for all φ C V P pR d , Rq,

ErφpG n pxqqs φpxq p¡1 i1 α i h i L i φpxq h p A p φpxq ¤ ¤ ¤ , (2.5)
for some constants α i and a linear differential operator A p . Assume further that

pA p rL, A p sq ¦ ρ V 0 (2.6)
where rL, A p s LA p ¡ A p L is the Lie bracket. Then X n yields an approximation of order p 1 for the invariant measure and it satisfies (1.8) and (1.9) with p replaced by p 1 and X n replaced by X n .

Theorem 2.8 is stated and proved in [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF] in the special case α i 0, i 1, . . . , p ¡ 1 in (2.5). However, the proof for non zero α i 's is nearly identical and thus is omitted. Notice that the order conditions (2.4) and (2.6) are respectively equivalent to the identities

» R d pA j φqρ V dx 0, j 2, . . . , p ¡ 1, » R d pA p φ rL, A p sφqρ V dx 0, for all test function φ C V P pR d , Rq.
In the following section, we introduce the suitable algebraic framework based on exotic aromatic trees and B-series for the systematic study of these order conditions of accuracy for the invariant measure.

Exotic aromatic trees and forests

We first recall the known framework of aromatic B-series before introducing a modification well suited for invariant measure order conditions and called exotic aromatic B-series. We rely on the aromatic trees and forests introduced in [START_REF] Chartier | Preserving first integrals and volume forms of additively split systems[END_REF] and rely on the presentation in [START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF].

Aromatic trees and forests

We first consider directed graphs γ pV, Eq with V a finite set of nodes and E V ¢ V the set of directed edges. If pv, wq E, we say that the edge is going from v to w, and v is called a predecessor of w. Two directed graphs pV 1 , E 1 q and pV 2 , E 2 q are equivalent if there exists a bijection ϕ :

V 1 Ñ V 2 with pϕ ¢ ϕqpE 1 q E 2 .
For brevity of notation, to avoid drawing arrows on the forests, an edge linking two nodes goes from the top node to the bottom one. If there is an eventual cycle, the arrows on it are going in the clockwise direction. For example, .

We call aromatic forests the equivalence classes of directed graphs where each node has at most one outgoing edge. The connected components making an aromatic forest are called aromatic trees. According to the above definition, there are two types of trees:

• aromas are aromatic trees5 with exactly one cycle: , , , , . . .

• rooted trees do not have a cycle ; they have a unique node that has no outgoing edge and that is called the root, graphically represented at the bottom: , , , , . . . Thus, an aromatic forest is a collection of aromas and rooted trees. We call AT t , , . . . u the set of aromatic forests containing exactly one rooted tree, and we name its elements the aromatic rooted forests.

Definition 3.1 (Elementary differentials).

Let γ pV, Eq AT , and let f : R d Ñ R d be a smooth function. We denote πpvq tw V, pw, vq Eu the set of all predecessors of the node v V and r the root of γ. We also call V 0 V tru tv 1 , . . . , v m u the other nodes of γ. Finally we introduce the notation I πpvq pi q 1 , . . . , i qs q where the q k are the predecessors of v, and

f I πpvq f f s f fx iq 1 . . . fx iq s .
Then F pγq is defined as

F pγqpfq d iv 1 ,...,iv m 1 £ ¹ vV 0 f I πpvq f iv f I πprq f. Example. Let γ k l j i m and r γ j i m l k
in AT where we added indices to apply the formula of Definition 3.1. Note that there is no index for the root. Then the associated differ-

entials are respectively F pγqpfq °d i,j,k,l,m1 f m f m f i f i f j f k f j,k f l f l f divpf q ¤ f I f P pf I f, f q and F pr γqpf q °d i,j,k,l,m1 f l f m f m,k f l f k f i f i f j f j f °d m1 f I m ppf m f q I pfqq ¤ f I f I f .

Exotic aromatic trees and forests

We now introduce a new kind of edge, called a liana, for the aromatic forests. The corresponding generalization is called exotic aromatic forests. Let pV, Eq be an aromatic forest and L be a finite list of pairs of elements of V (possibly with duplicates), then γ pV, E, Lq is an exotic aromatic forest. The elements of L are called lianas and correspond to nonoriented edges between any two nodes of the forest. We graphically represent them with a dashed edge linking the two given nodes. As we authorize duplicates, there can be several lianas between two given nodes. Also lianas can link a node to itself. For a node v, Γpvq denotes the list of the lianas (also with possible duplicates) linked to v. The predecessors of v only take in account the edges of E. An exotic aromatic tree of an exotic aromatic forest γ pV, E, Lq is a connected component of the associated aromatic forest pV, Eq. We call

EAT the set of exotic aromatic forests with exactly one rooted tree, and name its elements exotic aromatic rooted forests.

Example. The lianas can link different trees of an aromatic forest and thus yield an exotic aromatic forest. For instance, linking the aroma and the rooted tree gives .

The definition of elementary differentials is extended as follows.

Definition 3.2. Let γ pV, E, Lq EAT , and let f : R d Ñ R d be a smooth function. We name r the root of γ and V 0 V tru tv 1 , . . . , v m u the other nodes of γ. We denote l 1 ,. . . ,l s the elements of L and for v V , J Γpvq the multiindex pj lx 1 , . . . , j lx t q where Γpvq tl x 1 , . . . , l xt u. Then F pγq is defined as

F pγqpfq d iv 1 ,...,iv m 1 d jl 1 ,...,j ls 1 £ ¹ vV 0 f I πpvq f J Γpvq f iv f I πprq f J Γprq f.
Examples. The differential that corresponds to the rooted tree with a single node and a single liana is F p qpfq ∆f . We can also represent as exotic aromatic forest more complicated derivatives. For instance, let

γ i j k l , then F pγqpfq d i,j,k1
divpf i f q ¤ f I ppf kl f q I pf P pf ijj f, f kl f qqq.

Grafted exotic aromatic trees

For the study of the order for the invariant measure of numerical integrators, we introduce an extension of exotic aromatic forests. The root now symbolizes a test function φ, and it has leafs (nodes without predecessors) that represent a random standard normal vector ξ.

Note that these new trees can be seen as bi-coloured trees in the context of P-series (see [START_REF] Hairer | Geometric numerical integration[END_REF]Chap. 3]), where the nodes represented with crosses cannot have predecessors.

Definition 3.3. A grafted node is a new type of node graphically represented by a cross.

Let V be a set of nodes whose subset of grafted nodes is V g , let E be a set of edges such that each node in V g has exactly one outgoing edge and no ingoing edge, and let L be a set of lianas that link nodes in V V g , then γ pV, E, Lq is a grafted exotic aromatic forest. We define as before the grafted exotic aromatic trees and grafted exotic aromatic rooted forests, that we denote EAT g .

If γ pV, E, Lq is a grafted exotic aromatic rooted forest, φ : R d Ñ R a smooth function, and ξ a random vector of R d whose components are independent and follow a standard normal law, the associated elementary differential of γ is, with the same notation as Definition 3.2 and V 0 V pV g truq,

F pγqpf, φ, ξq d iv 1 ,...,iv m 1 d jl 1 ,...,j ls 1 £ ¹ vV 0 f I πpvq f J Γpvq f iv ¤ ¥ ¹ vVg ξ iv f I πprq f J Γprq φ.
Example. The differential associated to the forest is F p qpf, φ, ξq φ I pf P pξ, ξqq.

If γ is such that V g is empty, we recover the exotic aromatic forests of Definition 3.2, where φ is replaced by f . For the rest of the paper (except Section 3.5), we update the definition of the elementary differential of an exotic aromatic forests so that the root is associated to the function φ. This definition can be straightforwardly extended on nonrooted exotic aromatic forests. For brevity of notation, we also write F pγqpφq instead of F pγqpf, φ, ξq. We note that φ Ñ F pγqpφq is a linear differential operator (dependent on f and ξ).

Grafted exotic aromatic B-series

In this section, we adapt the formalism of aromatic B-series of [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF] to grafted exotic aromatic forests, in order to use it as a numerical tool for weak Taylor expansions in the next sections.

We define the order |γ| of a tree γ EAT g . We denote N pγq the number of nodes, N l pγq the number of lianas, N c pγq the number of grafted nodes and N v pγq N pγq ¡ N c pγq ¡ 1 the number of nodes that are non grafted and different from the root, then We extend the definition of F on VectpEAT g q by writing

|γ| N v pγq N l pγq N c pγq
F ¢ γEATg h |γ| apγqγ pφq Bpaqpφq.
The variable h is formal and thus can be chosen to be equal to 1. If the series is indexed only on (exotic) aromatic rooted forests, then it is called an (exotic) aromatic B-series. In Section 3.5, we shall focus on exotic aromatic B-series.

Remark 3.5. The coefficients apγq of standard B-series are sometimes renormalized as apγq ρpγq where ρ is a function determined by the symmetries of the associated forest. If ρ is appropriately chosen, it greatly simplifies the composition laws of (aromatic) B-series (see [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Chartier | Algebraic structures of B-series[END_REF][START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF]). Finding the best definition of ρ for this exotic extension of B-series is out of the scope of this paper.

Isometric equivariance of exotic aromatic rooted forests

In this subsection, we show that the exotic aromatic B-series satisfy an isometric equivariance property in the spirit of [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF][START_REF] Mclachlan | B-series methods are exactly the affine equivariant methods[END_REF]. We consider exotic aromatic rooted forests γ where the differential associated to the root is f . As the function f is no longer fixed, we denote the associated differential F pγqpfq. Also we adapt the definition of exotic aromatic B-series to this change. First we add a new tree: the empty tree ∅. The function F is then extended on EAT g t∅u by F p∅qpfq Id R d . Then, for a function a : EAT t∅u Ñ R, the associated exotic aromatic B-series is Bpaqpf q γEAT t∅u apγqF pγqpfq.

We study (exotic) aromatic B-series Bpaq with ap∅q 1. We call these (exotic) aromatic B-series methods. Let G be a subgroup of GL d pRq R d , let the action of an element pA, bq G on R d be x Þ Ñ Ax b, and let the action on a vector field f : R d Ñ R d be ppA, bq ¦ f qpxq : Af pA ¡1 px ¡ bqq.

We simplify the notation by writing A ¦f : pA, 0q ¦f. We recall the definition of equivariance from [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF]. The property of equivariance means the method is unchanged by an affine coordinate transformation. Let Φ be a differential operator and let G be a subgroup of

GL d pRq R d , then Φ is called G-equivariant if dpA, bq G, df C V pR d , R d q, ΦppA, bq ¦ f q pA, bq ¥ Φpf q ¥ pA, bq ¡1 .
In particular, Φ is said to be affine equivariant if

G GL d pRqR d and isometric equivariant if G O d pRq R d .
Theorem 3.6. Consider an exotic aromatic B-series method Bpaq (with ap∅q 1), then

Bpaq is isometric equivariant.

Remark 3.7. It is proved in [START_REF] Munthe-Kaas | Aromatic Butcher series[END_REF] that standard B-series methods are exactly the affine equivariant methods. Analogously, it would be interesting to characterize the isometric equivariant maps.

For the sake of brevity, we omit the proof of Theorem 3.6. The proof can be made in the spirit of the result [39, Prop. 2.1] for affine equivariant B-series.

Analysis of invariant measure order conditions using exotic aromatic forests

In this section, we show how the framework of Section 3 applies for the study of order conditions for the invariant measure of numerical integrators.

Weak Taylor expansion using exotic aromatic forests

Let us begin this section with the example of the θ-method (1.11). We apply the usual methodology to expand in Taylor series ErφpX 1 q|X 0 xs as h Ñ 0. We refer to [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] for other examples of analogous calculations performed without exotic aromatic forests. Under X 0 x, we have

X 1 x c hσξ hf h c hθσf I ξ h 2 θf I f h 2 θσ 2 2 f P pξ, ξq ¤ ¤ ¤
Then we deduce ErφpX 1 q|X 0 xs φpxq hLφpxq h 2 A 1 φpxq ¤ ¤ ¤ , where

A 1 φ Erθφ I f I f 1 2 φ P pf, f q θσ 2 2 φ I f P pξ, ξq θσ 2 φ P pf I ξ,
f i φf jk f i Erξ j ξ k s i,j f i φf jj f i F p qpφq,
where Erξ j ξ k s 0 for j $ k because of the independence of the ξ i 's. We see that taking the expectation of the differential associated to a grafted tree amounts to linking the grafted nodes with lianas in all possible manners. For instance, for the following example with four grafted nodes:

ErF p qpφqs Erφ p4q pξ, ξ, ξ, ξqs i,j,k,l f i,j,k,l φErξ i ξ j ξ k ξ l s i f i,i,i,i φErξ 4 i s 3 i,j i$j f i,i,j,j φErξ 2 i sErξ 2 j s 3 i,j f i,i,j,j φ 3F p qpφq.
Let us now comment this computation. The interesting fact is that Erξ 4 i s 3 corresponds exactly to the number of ways to gather the indices i, j, k and l in pairs. This observation makes an exotic aromatic tree naturally appear. However, here we took only four grafted nodes and the differential form was symmetric in the arguments ξ. We need to study the expectation of general exotic aromatic forest elementary differentials. This is the aim of the following theorem.

Theorem 4.1. Let γ EAT g be a grafted exotic aromatic rooted forest with an even number of grafted nodes 2n, let φ : R d Ñ R be a smooth function, and let V ¢ tc 1 , . . . , c 2n u be the set of grafted nodes of γ. We call P 2 p2nq the set of partitions by pair of t1, . . . , 2nu, i.e. the set of surjections p : t1, . . . , 2nu Ñ t1, . . . , nu such that the preimage of each singleton has exactly two elements and the minima of those preimages follow an ascending order (minpp ¡1 ptiuqq minpp ¡1 ptjuqq for i j). Finally we define ϕ γ : P 2 p2nq Ñ EAT the application that maps the partition p of γ to the aromatic forest where the grafted nodes are linked by lianas according to p. Then, the expectation of F pγqpφq is given by ErF pγqpφqs pP 2 p2nq

F pϕ γ ppqqpφq.

This theorem states that taking the expectation of the differential associated to a forest amounts to sum the forests obtained by linking the grafted nodes together pairwise using lianas in all possible manners and take the associated differential.

Example. Let us take γ t , , , . . . u the tree with only a root, 2n grafted nodes and no liana, then ErF pγqpφqs p2nq!

2 n n! ∆ n φ
The integer p2nq! Theorem 4.1 follows from the following lemma, which is an extension of the Isserlis theorem [START_REF] Isserlis | On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables[END_REF] to the case of multilinear mappings. The Isserlis theorem states that if χ is a 2n-dimensional Gaussian random vector with mean zero and arbitrary covariance, then

E 2n ¹ i1 χ i ' pP 2 p2nq ¹ i j ppiqppjq Erχ i χ j s. For n 2, it gives Erχ 1 χ 2 χ 3 χ 4 s Erχ 1 χ 2 sErχ 3 χ 4 s Erχ 1 χ 3 sErχ 2 χ 4 s Erχ 1 χ 4 sErχ 2 χ 3 s. Lemma 4.2. Let B : R d ¢ ¤ ¤ ¤ ¢ R d R 2nd
Ñ R be a 2n-multilinear form, and let ξ be a Gaussian vector N p0, I d q, then ErBpξ, . . . , ξqs

pP 2 p2nq d i1 ,...,in1
Bpe ip q, (4.3) with e ip pe i pp1q , . . . , e i pp2nq q, and we recall that e 1 , . . . , e d denotes the canonical basis of R d . Proof. For the particular case of an elementary multilinear form B σ : px 1 , . . . , x 2n q Þ Ñ ± 2n j1 px j q σpjq where σ : t1, . . . , 2nu Ñ t1, . . . , du is a given mapping and px j q σpjq denotes the σpjq's component of x j R d , the identity (4.3) reduces to the Isserlis theorem. As any multilinear form can be decomposed as a linear combination of such elementary multilinear forms, the result (4. B γ,φ pe ip q pP 2 p2nq

F pϕ γ ppqqpφq, because B γ,φ pe ip q is the differential F pγq where we differentiate the non-grafted nodes linked to grafted nodes in the directions given by e ip , and F pϕ γ ppqqpφq is obtained by summing B γ,φ pe ip q over all the indices. ¥

Integration by parts of the exotic aromatic forests

The goal of this section is to integrate by parts ³ R d F pγqpφqρ V dx, for γ an exotic aromatic rooted forest, in order to write it in the form

³ R d φ I r f ρ V dx for a certain sum of elementary differentials r
f . The idea is to transform a high order differential operator A : φ Ñ F pγqpφq into a differential operator φ Ñ φ I r f of order 1 such that A ¦ ρ V ¡ divp r f ρ V q. The tree formalism previously defined makes this task systematic and very convenient. This also serves as a crucial ingredient in the next section. Let us first begin with an example.

» R d F p qpφqρ V dx i,j » R d f 2 φ fx i fx j f i f j ρ V dx ¡ i,j » R d fφ fx j ff i fx i f j ρ V dx » R d fφ fx j f i ff j fx i ρ V dx » R d fφ fx j f i f j fρ V fx i dx & , ( 4.4)
where we integrated by parts; note that the boundary term vanishes using the growth assumptions on φ.

Notation. We denote g logpρ V q, then ∇ρ V p∇gqρ V .

We have

» R d F p qpφqρ V dx ¡ » R d divpf qφ I f ρ V dx ¡ » R d φ I f I f ρ V dx ¡ » R d g I f φ I f ρ V dx. By writing r f ¡pdivpfqf f I f g I f f q, we deduce » R d F p qpφqρ V dx » R d φ I r f ρ V dx.
We see that even for a simple forest, the integration by parts requires some calculations and a new term appears: the function g logpρ V q, and its derivatives. We use below the exotic aromatic forests to make this task easier.

Definition 4.3 (Aromatic root and elementary differential). An aromatic root is a new type of node represented by a square that has no outgoing edge. An exotic aromatic tree that has an aromatic root is considered as an aroma. The definition of the sets EAT is extended to include these new aromas.

If γ pV, E, Lq EAT is an exotic aromatic rooted forest whose set of aromatic roots is V 1 V , if V 0 V pV 1 truq and V tr, v 1 , . . . , v m u where r is the root of γ, if l 1 ,. . . ,l s are the elements of L, then F pγq is defined as

F pγqpφq d iv 1 ,...,iv m 1 d jl 1 ,...,j ls 1 £ ¹ vV 1 f I πpvq f J Γpvq g £ ¹ vV 0 f I πpvq f J Γpvq f iv f I πprq f J Γprq φ.
One can also extend the definition of EAT g in the same way so that it includes aromatic roots.

Examples. F p qpφq g I f φ I f and F p qpφq °i,j pf j gq 2 f i gf i φ.

Then the equality (4.4) can be rewritten using forests:

» R d F p qpφqρ V dx ¡ » R d F p qpφqρ V dx ¡ » R d F p qpφqρ V dx ¡ » R d F p qpφqρ V dx.
We notice that integrating by parts F p qpφq amounts to unplugging an edge from the root and to replug it either to all the other nodes of the forests or to an aromatic root and then to sum over all possibilities. This intuition is made rigorous in the following theorem.

Theorem 4.4. Let γ EAT , we choose a direction to integrate by parts, i.e. an edge or a liana e connected to the root r. Then

» R d F pγqpφqρ V dx ¡ ŗ γU pγ,eq » R d F pr γqpφqρ V dx,
where U pγ, eq is the set of exotic aromatic rooted forests obtained by unplugging the chosen edge/liana e and by linking it either to a node different from r or to a new aromatic root. Remark 4.5. Theorem 4.4 can be extended in the following way: if n is a node of γ and e an edge ingoing to n or a liana connected to n, then the same result holds if we replace U pγ, eq by U pγ, n, eq, the set of exotic aromatic rooted forests obtained by unplugging the chosen edge/liana e and by linking it either to a node different from n or to a new aromatic root.

Proof of Theorem 4.4. We call r the root and V 0 tv 1 , . . . , v m u the other nodes of γ. We suppose for simplicity that γ does not have any aromatic root, otherwise the proof can be adapted straightforwardly. We denote l 1 ,. . . ,l s the elements of L and v k 1 ,. . . ,v kp the elements of πprq. We choose to integrate by parts in the direction of the edge

x iv k 1 . » R d F pγqpφqρ V dx d iv 1 ,...,iv m 1 d jl 1 ,...,j ls 1 » R d £ ¹ vV 0 f I πpvq f J Γpvq f iv f iv k 1 ...iv kp f J Γprq φρ V dx ¡ d iv 1 ,...,iv m 1 d jl 1 ,...,j ls 1 u V 0 » R d ¤ ¥ ¹ vV 0 tuu f I πpvq f J Γpvq f iv f iv k 1 f I πpuq f J Γpuq f iu f iv k 2 ...iv kp f J Γprq φρ V dx » R d £ ¹ vV 0 f I πpvq f J Γpvq f iv f iv k 2 ...iv kp f J Γprq φf iv k 1 gρ V dx &
Each term of the sum on u V 0 is the differential associated to the forest γ u . This forest is obtained by unplugging the root of its edge linking it to v k 1 , and sticking it to u. The last term of the computation is the differential of the forest obtained by linking the unplugged edge to an aromatic root. These terms are exactly what we expected, thus the theorem is proved for the case of edges. For the case of integrating in the direction of a liana, the proof is nearly identical. We just need to develop J Γprq instead of I πprq . ¥ Definition 4.6. Let γ 1 and γ 2 be two exotic aromatic B-series, we define the equivalence relation and write γ 1 γ 2 if we can transform γ 1 into γ 2 by integrating by parts the associated differentials according to the procedure presented in the previous theorem.

Examples. The integration by parts (4.4) can now be simply rewritten as

¡ ¡ ¡ .

One can also iterate the process of integration by parts to fully simplify the trees. For example, we have

¡ ¡ ¡ ¡ ¡

and summing up yields

¡ ¡ 2 .
Finally here is a last example that will be used in Section 5.1. We apply the procedure of integration by parts to the forest .

¡ ¡ ¡ ¡ 2 ¡ Using analytic formulas, ³ R d F p qpφqρ V dx ³ R d φ I r f ρ V dx, where r f i ¡f i p∆gq ¡ ∆gf i g ¡ 2 d j1 f j gf i,j g ¡ d j1
pf j gq 2 f i g.

Order conditions using exotic aromatic forests

In this section, we adapt Theorem 2.6 in the context of exotic aromatic forests. In the spirit of traditional B-series, we give, under Assumption 2.2, the general simplification for orders up to three of a general numerical method expandable in exotic aromatic B-series. With these, one can improve a method order as presented in Section 5.1, or derive conditions on the method to achieve high order for the invariant measure as we do in Sections 5.2, 5.3, 5.4 and 5.5. It is worth noting that with only integration by parts, we can formally derive numerical methods (see Section 5.1 for an example), but under Assumption 2.2, the methods can be simplified. In the first and second cases, one can replace the nodes A, B, C with aromatic roots and the result remains. For the third case, the node B also be replaced.

Proof. For the first pair of patterns, the associated differentials have the respective forms

°d i1 f j f k f i f j f i f l and °d j1 f j f k f j f i f i f l .
As f is a gradient, f I is a symmetric matrix and f i f j f j f i . The two differentials are then equal.

The second point is proved in the same way. For the third and fourth points, we just use that ∇g 2 σ 2 f . ¥

Example. We have F p q 2 σ 2 F p q and F p q 2 σ 2 F p q. The equivalence relation of Definition 4.6 is extended to include the simplification rules of Proposition 4.7. [START_REF] Hairer | Geometric numerical integration[END_REF]Chap. 3,Exercise 3]. Indeed, for a given exotic aromatic forest γ pV, E, Lq, we take a bijective numbering n V : V tru Ñ t1, 2, . . . , |V | ¡ 1u of the nodes and another one for the lianas n L : L Ñ t|V | , . . . , |γ|u, then we define for p N. We denote A i F pγ i q. If γ i r γ i and F pr γ i q 0 for all 1 ¤ i p, then the method is of order (at least) p for the invariant measure. In particular, under Assumption 2.2, we obtain for order 1

Remark 4.8. If f is a general vector field not assumed to be a gradient, we can prove that the elementary differentials of exotic aromatic forests are independent. The proof is an extension of the result in

φpxq ¹ vπprq x n V pvq ¹ lΓprq x n L plq , f i pxq ¹ vπpn ¡1 v piqq x n V pvq ¹ lΓpn ¡1 v piqq x n L plq , i 1, . . . , |V | ¡
r γ 0 ¢ a ¡ © ¡ 2 σ 2 a ¡ ©
, and in addition we have for order 2,

r γ 1 ¢ ap q ¡ 2 σ 2 ap q 2 σ 2 ap q ¡ 4 σ 4 ap q ¢ ap q ¡ ap q ap q ¡ 2 σ 2 ap q ¢ ap q ¡ 2 σ 2 ap q 4 σ 4 ap q .
For order 3, the expression of r γ 2 can be found in the appendix.

Remark 4.10. In contrast to Section 4.2, we choose in Theorem 4.9 not to reduce the tree for γ 1 and the three trees , and for γ 2 in forests with exactly one edge linked to the root. The reason is these trees do not simplify well, and the coefficient multiplying them vanishes for most methods (for all consistent Runge-Kutta methods for example). If one wants to compute the f i of Section 5.1, one should integrate by parts these trees first.

Construction of high order integrators

Improvement of a method order via a modified equation

In [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF], a recursive method to obtain integrators of any order for the convergence to the invariant measure is presented. Let us suppose we have an integrator of order exactly p ¥ 1 for the invariant measure, then, using Theorem 2.6, for all j p, A ¦ j ρ V 0 and A ¦ p ρ V $ 0.

By integrating by parts, we can write

³ R d A p φρ V dx as ³ R d φ I f p ρ V dx.
We then consider the same numerical integrator but for the modified equation where we replaced f by f ¡ h p f p .

Applying Theorem 2.6 to the new context, we see that this integrator is at least of order p 1 for the original equation.

In this section, we give tools to simplify the computation of those modified integrators, in particular to calculate simply the operators A j , and to find the function f p .

Example. For the θ-method (1.11), we have A 1 F pγq where γ is given by (4.2). Applying integration by parts as described in Section 4.2, we obtain

γ ¢ θ ¡ 1 2 ¡ 1 2 ¡ 1 2 σ 2 2 p1 ¡ θq σ 2 p1 ¡ θq σ 2 2 σ 2 2 ¡ σ 4 8 ¡ σ 4 8 ¡ σ 4 4 ¡ σ 4 8 .
Then f 1 is given by

pf 1 q i ¢ θ ¡ 1 2 f I i f ¡ 1 2 g I f f i ¡ 1 2 divpf qf i σ 2 2 p1 ¡ θq∆f i σ 2 p1 ¡ θq d j1 f j f i f j g σ 2 2 ∆gf i σ 2 2 d j1 pf j gq 2 f i ¡ σ 4 8 f i ∆g ¡ σ 4 8 ∆gf i g ¡ σ 4 4 d j1 f j gf i,j g ¡ σ 4 8 d j1
pf j gq 2 f i g.

Remark 5.1. Note that the above computation does not suppose Assumption 2.2 to be satisfied. But as the function g is not explicitly known in the general case in practice, we cannot implement easily the corresponding scheme. Under Assumption 2.2, we can directly apply Theorem 4.9 and find a simple expression of f 1 :

f 1 ¢ 1 2 ¡ θ pf I f σ 2 2 ∆f q.
Application. Let us calculate the modified integrator f 2 , where f follows Assumption 2.2. First we rewrite the differential operator A 2 for the modified equation where we replaced f by f ¡ hf 1 . We call it A p1q

2 . We find A p1q

2 F pγq where

γ θp3θ ¡ 1q θp3θ ¡ 1q σ 2 2 θp4θ ¡ 1q σ 2 2 θ 2 θ 2 σ 4 4 θ 2 σ 2 θp4θ ¡ 1q σ 2 2 θp2θ 1q σ 2 2 θp2θ 1q σ 4 4 4θ ¡ 1 2 p4θ ¡ 1q σ 2 4 θ 2 σ 2 2 p4θ ¡ 1q σ 2 4 p4θ ¡ 1q σ 4 8 θ σ 4 2 θσ 2 1 6 θ σ 4 2 σ 2 4 σ 4 8 σ 6 48 .
Using Theorem 4.9, we find

γ ¢ ¡2θ 2 2θ ¡ 1 2 ¢ ¡θ 2 θ ¡ 1 4 σ 2 ¢ ¡ 3θ 2 2 3θ 2 ¡ 1 3 σ 2 ¢ ¡θ 2 θ ¡ 1 6 ¢ ¡ θ 2 4 θ 4 ¡ 1 24 σ 4 ¢ ¡θ 2 θ ¡ 1 6 σ 2 .
Thus we define

f 2 ¢ ¡2θ 2 2θ ¡ 1 2 f I f I f ¢ ¡θ 2 θ ¡ 1 4 σ 2 f I ∆f ¢ ¡ 3θ 2 2 3θ 2 ¡ 1 3 σ 2 i f P pe i , f I pe i qq ¢ ¡θ 2 θ ¡ 1 6 f P pf, f q ¢ ¡ θ 2 4 θ 4 ¡ 1 24 σ 4 ∆ 2 f ¢ ¡θ 2 θ ¡ 1 6 σ 2 p∆fq I pfq,
and, if the θ-scheme applied to dX pf ¡hf 1 ¡h 2 f 2 qdt σdW is ergodic, then it has order 3 the invariant measure.

For f 1 , we recover the formula of [3, Prop. 5.1]. The computation of f 2 was first done for θ 0 in [3, Prop. 5.2], which reveals a typographical error. This method can give numerical integrators of any order, but it comes with a high computing price if the partial derivatives of f are difficult to compute. In the following sections, we present order conditions for certain classes of numerical schemes, in order to obtain high order methods avoiding derivatives and unnecessary evaluations of f .

Order conditions for stochastic Runge-Kutta schemes

We consider stochastic Runge-Kutta schemes (1.10) for the overdamped Langevin equation (1.1). We set c i s °j1 a ij . In this section, we also assume Assumption 2.2 to simplify the computation. Using the proposed framework, our goal is to find algebraic conditions on the coefficients A pa ij q, b pb i q and d pd i q to achieve a given order condition for the invariant measure.

First, we suppose °bi 1 in order for A 0 L in Assumption 2.4 to be satisfied. Then A 1 φ F pγ 1 qpφq where

γ 1 ¸bi c i σ 2 2 ¸bi d 2 i σ 2 ¸bi d i 1 2 σ 2 2 σ 4 8 .
Theorem 4.9 yields

γ 1 ¢ ¸bi c i 1 2 ¡ 2 ¸bi d i σ 2 2 ¢ ¸bi d 2 i 1 2 ¡ 2 ¸bi d i . Thus if we suppose ¸bi 1, ¸bi c i 1 2 ¡ 2 ¸bi d i 0, ¸bi d 2 i 1 2 ¡ 2 ¸bi d i 0, then A 1 φ 0.
We have a Runge-Kutta scheme of order 2.

By continuing this methodology, we obtain the order conditions of order 3, and our analysis allows us to obtain the conditions for any order. The following theorem states the order conditions for Runge-Kutta methods.

Theorem 5.2. Assume Assumption 2.2 and consider an ergodic Runge-Kutta method (1.10) with °bi 1. Using the same notation as in Theorem 4.9, if A, b and d are chosen such that F pr γ i q 0 for all 1 ¤ i p, then the method has at least order p for the invariant measure. In particular, Table 1 gives sufficient conditions to have consistency and order 2 or 3 for the invariant measure for Runge-Kutta schemes.

Order Tree τ F pτqpφq

Order condition

1 φ I f °bi 1 2 φ I f I f °bi c i ¡ 2 °bi d i ¡ 1 2 φ I ∆f °bi d 2 i ¡ 2 °bi d i ¡ 1 2 3 φ I f I f I f °bi a ij c j ¡ 2 °bi a ij d j °bi c i ¡ p °bi d i q 2 0 φ I f I ∆f °bi a ij d 2 j ¡ 2 °bi a ij d j °bi c i ¡ p °bi d i q 2 0 1 2 °bi c 2 i ¡ 2 °bi d i c i φ I f P pf, f q ¡2 °bi d i 2 °bi d 2 i °bi c i ¡ 1 3 °bi d i a ij d j ¡ °bi c i d i ¡ °bi d i °bi d 2 i °φI f P pf I pe i q, e i q °bi c i ¡ °bi a ij d j ¡ 1 2 p °bi d i q 2 ¡ 1 6 1 2 °bi c i d 2 i ¡ °bi d 3 i ¡ 2 °bi d i φ I pp∆fq I pfqq 5 2 °bi d 2 i ¡ °bi c i d i 1 2 °bi c i ¡ 1 3 1 8 °bi d 4 i ¡ 1 2 °bi d 3 i φ I ∆ 2 f ¡ 1 2 °bi d i 3 4 °bi d 2 i ¡ 1 12
Table 1: Runge-Kutta order conditions for the invariant measure (See Theorem 5.2). The sums are over all involved indices. Remark 5.3. To check if a scheme has weak order p, one can develop L i φ, i ¤ p, in exotic aromatic B-series (a simple method for this computation is proposed in Section 5.3) and prove that the forest coefficients of 1 i! L i φ and A i¡1 φ are equal, yielding the order p estimate (1.7). In Table 2, we collect the corresponding order conditions up to order p ¤ 3. We recover exactly the same order conditions as first derived in [START_REF] Debrabant | Runge-Kutta methods for third order weak approximation of SDEs with multidimensional additive noise[END_REF]Thm. 4] using different types of trees and B-series. We recall that the conditions for weak order 3 have no solution for a method of the form (1.10) with only l 1 noise. Indeed, fixing d p2q 0 in Table 2, we obtain the incompatible order conditions °bi d p1q i ¨2 1 2 (third line of Table 2) and °bi d p1q i ¨2 1 3 (last line of Table 2). Taking l 2 noises in the method (1.10) is sufficient for reaching weak order 3 for general f satisfying Assumption 2.2.

Remark 5.4. Notice that Assumption 2.2 permits us to identify the differentials F p q F p q and the corresponding order conditions, using Proposition 4.7. Thus, under Assumption 2.2, we can replace the two conditions °bi a ij d j 1

6 and

¡ °bi r d i © 2 1
12 by one

new condition °bi a ij d p1q j 1 2 ¡ °bi d p2q i © 2 5
24 . Remark 5.5. As explained in the introduction, the study of weak order conditions using rooted trees is already well documented in the literature, but the framework of exotic aromatic B-series has the advantage to involve rooted forests that do not depend on the dimension d, which permits us to compute integration by parts and hence derive the order conditions for the invariant measure. Since weak convergence implies convergence with at least the same order for the invariant measure, the weak order conditions (Table 2) imply the order conditions for the invariant measure (Table 1). In particular, comparing Table 1 and Table 2, we observe that there is a lower number of order conditions for the convergence to the invariant measure compared to the standard weak convergence.

Order conditions for postprocessed integrators

In this section, we extend our analysis to the case of integrators combined with postprocessors [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF]. As stated in Theorem 2.8, it permits us to increase the order for the invariant measure of a given method while maintaining a low number of function evaluations per time step. We show that exotic aromatic B-series simplify this approach, but one issue remains: the computation of the Lie bracket rL, A p sφ. This is done by the following theorem for the composition of exotic aromatic forests and based on the Leibniz rule. Theorem 5.6. Let γ 1 and γ 2 be two exotic aromatic rooted forests, and let φ : R d Ñ R be a smooth function. For ϕ : πpr 2 q Ñ V 1 and ψ : Γpr 2 q Ñ V 1 , we build γ ϕ,ψ by plugging all the edges connected to r 2 to the nodes of γ 1 according to ϕ, and all the lianas (counting multiplicity) according to ψ. Then the composition of forests is given by

F pγ 2 qpFpγ 1 qpφqq φ:πpr 2 qÑV 1 ψ:Γpr 2 qÑV 1 F pγ ϕ,ψ qpφq.
Various composition rules for B-series and aromatic B-series have been studied in the literature (see [START_REF] Hairer | Geometric numerical integration[END_REF][START_REF] Chartier | Algebraic structures of B-series[END_REF][START_REF] Bogfjellmo | Algebraic structure of aromatic B-series[END_REF] and the references therein). The main difference from these previous works is that we compose only the roots of exotic aromatic rooted forests, because this corresponds to composing linear differential operators.

Proof. Using Definition 3.2, we have

F pγ 1 qpφq iv p1q 1 ,...,i v p1q m 1 jl p1q 1 ,...,j l p1q s 1 ¤ ¥ ¹ vV 0 1 f I πpvq f J Γpvq f iv f I πpr 1 q f J Γpr 1 q φ. F pτqpφq Order condition 1 φ I f °bi 1 2 φ I f I f °bi c i 1 2 φ I ∆f °bi d p1q i ¨2 °bi d p2q i ¨2 1 2 °φP pe i , f I pe i qq °bi d p1q i 1 2 3 φ I f I f I f °bi a ij c j 1 6 φ I f I ∆f °bi a ij d p1q j ¨2 °bi a ij d p2q j ¨2 1 6 φ I f P pf, f q °bi c 2 i 1 3
°φI f P pf I pe i q, e i q °bi d p1q °φQ pe i , e j , f P pe i , e j qq °bi d p1q i ¨2 1 3

Table 2: Runge-Kutta standard weak order conditions for l 2 noises. The sums are over all involved indices. We recover the same conditions as in [START_REF] Debrabant | Runge-Kutta methods for third order weak approximation of SDEs with multidimensional additive noise[END_REF].

. Then we replace φ by F pγ 1 qpφq and use the Leibniz rule to distribute the partial derivatives. qpφq.

F pγ 2 qpFpγ 1 qpφqq iv p1q 1 ,...,i v p1q m 1 i v p2q 1 ,...,i v p2q m 2 jl p1q 1 ,...,j l p1q s 1 j l p2q 1 ,...,j l p2q s 2 ¤ ¥ ¹ vV 0 2 f I πpvq f J Γpvq f iv ¤ f I πpr 2 q f J Γpr 2 q ! ¤ ¥ ¹ vV 0 1 f I πpvq f J Γpvq f iv f I πpr 1 q f J Γpr 1 q φ ( ) φ:πpr 2 qÑV 1 ψ:Γpr 2 qÑV 1 iv p1q 1 ,...,i v p1q m 1 i v p2q 1 ,...,i v p2q m 2 jl p1q 1 ,...,j l p1q s 1 j l p2q 1 ,...,j l p2q s 2 ¤ ¥ ¹ vV 0 2 f I πpvq f J Γpvq f iv ¤ ¤ ¥ ¹ vV 0 1 f I πpvqϕ ¡1 ptvuq f J Γpvqψ ¡1 ptvuq f iv f I πpr 1 qϕ
Combining the two previous equalities, we deduce

L 2 φ F p σ 2 σ 2 2 σ 2 σ 4 4
qpφq.

Using Theorems 2.8 and 5.6, we obtain general conditions on postprocessors to increase by 1 the order of a given method. if we denote γ the exotic aromatic B-series such that F pγq pA p rL, A p sq and if γ 0, then X n is of order p 1 for the invariant measure. In particular, if the order 2 conditions in Table 3 are verified, then the method has order 2 for the invariant measure. If we suppose the order 2 for the invariant measure of the numerical method and ap q 2 σ 2 ap q 0 (in order to have A 0 αL), and if conditions in Table 4 are verified, then the method is of order 3 for the invariant measure.

Tree γ

Order condition ap q ¡ 2 σ 2 ap q 2 σ 2 ap q ¡ 4 σ 4 ap q ¡ 2ap q 4 σ 2 ap q 0 ap q ¡ ap q ap q ¡ 2 σ 2 ap q ¡ σ 2 2 ap q ap q 0 ap q ¡ 2 σ 2 ap q 4 σ 4 ap q 0

Table 3: General order 2 conditions with postprocessor (See Theorem 5.7).

Theorem 5.8. Consider an ergodic Runge-Kutta method of order p ¥ 1 for the invariant measure of the form (1.10) and the following associated postprocessor

Y i X n h s °j1 a ij f pY j q d i σ c h ξ n , i 1, . . . , s, X n X n h s °i1 b i f pY i q d 0 σ c h ξ n .
Assume Assumption 2.2. If γ is the exotic aromatic B-series such that F pγq pA p rL, A p sq and if A, b, d, A, b, d, d 0 are chosen such that γ 0 then the postprocessed method X n has at least order p 1 for the invariant measure. In particular, if the conditions of order 2 in Table 5 are verified, then the postprocessed integrator has order 2. If the Runge-Kutta method has order 2 for the invariant measure (see Table 1), if °bi d 0 2 and if the conditions of order 3 in Table 5 are verified, then the method has order 3.

Remark 5.9. The condition of in Table 3 is not modified by the postprocessor as it does not depend on any apγq. Thus if the scheme is a consistent Runge-Kutta method, this condition is automatically satisfied (see also Remark 4.10). If the postprocessor satisfies the equation ap q ¡ 2 σ 2 ap q 4 σ 4 ap q 0, (5.1)

then the conditions of , and of Table 4 are not modified by the postprocessor, and thus are automatically satisfied for consistent Runge-Kutta methods. Equation (5.1) is verified for the class of Runge-Kutta methods for postprocessors presented in Theorem 5.8 that satisfy A 0 αL. [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF], is of order 2 for the invariant measure of (1.4) (if it is ergodic).

Example. Under Assumption 2.2, the following Runge-Kutta method, introduced in

X n 1 X n hf pX n 1 ¡1 c 2 2 σ c hξ n q σ c hξ n , X n X n h c 2 2 f pX n q c 4 c 2¡1 2 σ c h ξ n . Tree τ Order condition ap q ¡ 2 σ 2 ap q ¡ 2 σ 2 ap q 2 σ 2 ap q 2 σ 2 ap q ¡ 4 σ 2 ap
q ¡ 4 σ 4 ap q 12 σ 4 ap q ¡ 24 σ 6 ap q ¡ 2ap q 4 σ 2 ap q 2ap q ¡ 8 σ 2 ap q 16 σ 4 ap q 0 ap q ¡ ap q ¡ ap q ap q ap q ¡ 2ap q ¡ 2 σ 2 ap q 6 σ 2 ap q ¡ 12 σ 4 ap q ¡ σ 2 2 ap q ¡ ap q 2ap q σ 2 ap q ¡ 4ap q 8 σ 2 ap q 0 ap q ¡ 2 σ 2 ap q 4 σ 4 ap q 2 σ 2 ap q ¡ 4 σ 4 ap q 4 σ 4 ap q ¡ 8 σ 6 ap q ¡ ap q ¡ ap q 2 σ 2 ap q 0 ap q ¡ ap q ¡ ap q ¡ ap q 2ap q 2 σ 2 ap q ap q ¡ 2ap q ¡ 4 σ 2 ap q 8 σ 2 ap q ¡ 16 σ 4 ap q ¡ σ 2 ap q σ 2 ap q 2ap q ¡ 4ap q 8 σ 2 ap q 0 ap q ¡ 2 σ 2 ap q ¡ ap q 2 σ 2 ap q ap q 4 σ 2 ap q ¡ 4 σ 2 ap q 4 σ 2 ap q ¡ 8 σ 4 ap q ¡ σ 2 2 ap q ¡ ap q 2ap q ¡ 2ap q 4 σ 2 ap q 0 ap q ¡ ap q ap q ap q ¡ ap q ap q ¡ 2 σ 2 ap q ¡ σ 2 2 ap q σ 2 2 ap q ¡ σ 2 2 ap q ap q 0 ap q ¡ 2 σ 2 ap q 4 σ 4 ap q ¡ 8 σ 6 ap q 0 ap q ¡ ap q 2ap q ¡ 2 σ 2 ap q 2 σ 2 ap q ¡ 6 σ 2 ap q 12 σ 4 ap q ¡ σ 2 ap q 2ap q ¡ 4 σ 2 ap q 0 ap q ¡ 2 σ 2 ap q ¡ 2 σ 2 ap q 4 σ 4 ap q 4 σ 2 ap q ¡ 12 σ 4 ap q 24 σ 6 ap q ¡ 4ap q 8 σ 2 ap q ¡ 16 σ 4 ap q 0 Table 4: General order 3 conditions with postprocessor (See Theorem 5.7).

Order Tree τ Order condition

2 °bi c i ¡ 2 °bi d i ¡ °bi 2d 0 2 ¡ 1 2 °bi d 2 i ¡ 2 °bi d i ¡ °bi d 0 2 ¡ 1 2 °bi a ij c j ¡ 2 °bi a ij d j °bi c i 3 ¡ p °bi d i q 2 ¡ 2 °bi c i 4d 0 °bi d i ¡ d 0 4 0 °bi a ij d 2 j ¡ 2 °bi a ij d j °bi c i ¡ p °bi d i q 2 ¡ °bi c i ¡ °bi d i 2 4d 0 °bi d i ¡ d 0 4 0 1 2 °bi c 2 i ¡ 2 °bi d i c i ¡ 2 °bi d i 2 °bi d 2 i °bi c i ¡ °bi c i 2d 0 °bi d i ¡ d 0 4 2 ¡ 1 3 °bi d i a ij d j ¡ °bi c i d i ¡ °bi d i °bi d 2 i °bi c i ¡ °bi a ij d j ¡ 1 2 p °bi d i q 2 ¡ °bi c i 2d 0 °bi d i ¡ d 0 4 2 ¡ 1 6 1 2 °bi c i d 2 i ¡ °bi d 3 i ¡ 2 °bi d i 5 2 °bi d 2 i ¡ °bi c i d i 1 2 °bi c i ¡ 1 2 °bi c i ¡ 1 2 °bi d i 2 2d 0 °bi d i ¡ d 0 4 2 ¡ 1 3 1 8 °bi d 4 i ¡ 1 2 °bi d 3 i ¡ 1 2 °bi d i 3 4 °bi d 2 i ¡ 1 4 °bi d i 2 1 2 d 0 °bi d i ¡ d 0 4 8 ¡ 1 12
Table 5: Order conditions for Runge-Kutta method with Runge-Kutta postprocessor (See Theorem 5.8). The sums are over all involved indices.

Order conditions for partitioned methods

In (1.1), we assume f f 1 f 2 and we consider partitioned integrators that apply different numerical treatments to each f i . We explain in this section how to extend the exotic aromatic B-series formalism to compute order conditions for such partitioned integrators. The advantage is to treat differently each part of f according to their properties. For example, if f 1 is stiff and f 2 is non-stiff, one would like to apply an implicit method to f 1 and an explicit method for f 2 (IMEX methods).

Here we follow the formalism of [23, Sect. III.2] for bicoloured B-series, called P-series. We introduce white nodes ; they represent the function f 2 . Black nodes now correspond to f 1 but the root still corresponds to φ. We call these new forests exotic aromatic P-forests. There are two slight changes in the computation rules compared to the non-partitioned case:

• Simplification rule: if f 1 ∇V 1 and f 2 ∇V 2 are both gradients, then

A B 2 σ 2 ¢ A B A B .
Furthermore, the node B can be replaced by an aromatic root or a white node.

• The operator L is now written as

L F p σ 2 2 q.
In addition to the partitioning of the one can also add a postprocessor. The results of Section 4.3, 5.2 and 5.3 are straightforwardly adapted to the P-forests. Theorem 5.10. Consider a Runge-Kutta method of order p of the form

Y i X n h s °j1 a ij f 1 pY j q x a ij f 2 pY j q d i σ c hξ n , i 1, . . . , s, X n 1 X n h s °i1 b i f 1 pY i q p b i f 2 pY i q σ c hξ n ,
together with the following Runge-Kutta postprocessor

Y i X n h s °j1 a ij f 1 pY j q x a ij f 2 pY j q d i σ c h ξ n , i 1, . . . , s, X n X n h s °i1 b i f 1 pY i q p b i f 2 pY i q d 0 σ c h ξ n .
Under the notation and assumptions of Theorem 5.8, if we suppose f 1 and f 2 are gradients, if we choose f 1 , f 2 and the coefficients of the method such that γ 0 then the method has at least order p 1 for the invariant measure. In particular, the conditions for consistency and order 2 are in Table 6.

Order Tree τ F pτqpφq

Order condition

1 φ I f 1 °bi 1 φ I f 2 °p b i 1 2 φ I f I 1 f 1 °bi c i ¡ 2 °bi d i ¡ 2 °bi 2d 0 2 ¡ 1 2 φ I f I 1 f 2 °bi p c i ¡ 2 °bi d i ¡ °bi ¡ °p b i 2d 0 2 ¡ 1 2 φ I f I 2 f 1 °p b i c i ¡ 2 °p b i d i ¡ °bi ¡ °p b i 2d 0 2 ¡ 1 2 φ I f I 2 f 2 °p b i p c i ¡ 2 °p b i d i ¡ 2 °p b i 2d 0 2 ¡ 1 2 φ I ∆f 1 °bi d 2 i ¡ 2 °bi d i ¡ °bi d 0 2 ¡ 1 2 φ I ∆f 2 °p b i d 2 i ¡ 2 °p b i d i ¡ °p b i d 0 2 ¡ 1 2
Table 6: Order conditions for partitioned Runge-Kutta method with postprocessor (See Theorem 5.10). The sums are over all involved indices.

Example. Using the previously introduced formalism, we see that, if f 1 and f 2 are gradients and f f 1 f 2 satisfies Assumption 2.2, the following method, adapted from [10, Lemma 2.9], is of order 2 for the invariant measure of (1.4) (if it is ergodic).

X n 1 X n h 2 f 1 pX n 1 1 2 σ c hξ n q h 2 f 1 pX n 1 3 2 σ c hξ n q 2 pX n 1 2 σ c hξ n q σ c hξ n , X n X n 1 2 σ c h ξ n .
It can be put in Runge-Kutta form with the coefficients below: 2 for the postprocessor.

If we add a family of independent noises pχ n q n independent of pξ n q n , then by extending Theorem 5.10, we can show that, under the same hypothesis as the previous example, the following IMEX method has order 2 for the invariant measure of (1.4):

X n 1 X n hf 1 pX n 1 1 2 σ c hχ n q hf 2 pX n 1 2 σ c hξ n q σ c hξ n , X n X n 1 2 σ c h ξ n .

Non-reversible perturbation

An interesting modification of (1.4) is to introduce a non gradient perturbation that preserves the invariant measure. It permits for some classes of problems to improve the rate of convergence to equilibrium [START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF], and it can also reduce the variance [START_REF] Duncan | Variance reduction using nonreversible Langevin samplers[END_REF]. As in Section 5.4, we consider the equation (1.1) where f f 1 f 2 and we use bicoloured forests. We suppose f 1 ¡∇V is a gradient, and f 2 is a perturbation of f 1 that satisfies div ¡ f 2 e ¡ 2 σ 2 V © 0.

(5.

2)

The perturbation f 2 does not modify the invariant measure. Indeed equation (5.2) implies that the adjoint of Bφ φ I pf 2 q satisfies B ¦ ρ V 0, and thus the invariant measure is preserved. A simple example of such non gradient perturbation is f 2 J∇V , with J a fixed antisymmetric matrix. We can now apply all the results of Section 4 that do not use Assumption 2.2. We have the following useful properties.

• We still have the simplification rule (see Proposition 4.7):

A B 2 σ 2 A B
. Furthermore, the node B can be replaced by an aromatic root or a white node.

• The generator reads L F p σ 2 2 q.

• We have F p q ¡Fp q, and these differentials vanish if f 2 J∇V .

The two first properties allow us to simplify lianas in the forests as we did in Section 5.4. Then we are left with forests with white nodes such as . This is where the last property comes in handy, as we can integrate by part this tree and obtain ¡ . We deduce the following theorem.

Theorem 5.11. Consider an ergodic Runge-Kutta method and a postprocessor as in Theorem 5.10, and suppose f 1 ¡∇V and f 2 satisfies (5.2). Then under the notation and assumptions of Theorem 5.8, if the coefficients of the method are chosen such that γ 0 then the method has at least order p 1 for the invariant measure. In particular, the conditions for order 1 and order 2 are in Table 7.

Order Tree τ F pτqpφq Order condition We note that if f 2 J∇V , we have σ 2

1 φ I f 1 °bi 1 2 φ I f I 1 f 1 °bi c i ¡ 2 °bi d i ¡ 2 °bi 2d 0 2 ¡ 1 2 φ I f I 1 f 2 °bi p c i °bi ¡ °p b i 0 φ I f I 2 f 1 °p b i c i ¡ 2 °p b i d i °p b i ¡ °bi ¡ °p b i 2d 0 2 0 φ I f I 2 f 2 °p b i p c i ¡ 1 2 ¡ °p b i © 2 0 φ I ∆f 1 °bi d 2 i ¡ 2 °bi d i ¡ °bi d 0 2 ¡ 1 2 φ I ∆f 2 °p b i d 2 i ¡ 2 °p b i d i °p b i ¡ °p b i d 0 2 0

2

¡ ¡ . In this case, the order condition for can be omitted and the two conditions of and are respectively replaced by

2 ¸p b i d i ¡ ¸p b i ¸bi p c i ¡ ¸p b i d 2 i ¸bi ¡ d 0 2 0, ¸p b i c i ¸p b i d 2
i ¡ ¸bi d 0 2 0. Remark 5.12. In order for the method to satisfy A 0 L, the condition °p b i 1 should be added in Table 7, but it is not necessary to achieve order 1 for the invariant measure.

Example. If f 1 satisfies Assumption 2.2 and f 2 satisfies (5.2), the following consistent postprocessed scheme has order 2 for the invariant measure (if it is ergodic):

X n 1 X n hf 1 pX n σ 2 c
hξ n q 5 4 hf 2 pX n σ

2 c hξ n q ¡ 1 4 hf 2 pX n ¡ 2hf 2 pX n σ 2 c hξ n q ¡ σ 2 c hξ n q σ c hξ n , X n X n σ 2 c h ξ n .
(5.3)

If f 2 J∇V , it needs two evaluations of ∇V per timestep similarly to a standard Runge-Kutta weak order 2 method. For f 2 0, note that the scheme (5.3) coincides with the one proposed in [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF], formulated in a different manner (See [START_REF] Vilmart | Postprocessed integrators for the high order integration of ergodic SDEs[END_REF]).

q ¡ 24 σ 6 ap q ¢ ap q ¡ ap q ¡ ap q ap q ap q ¡ 2ap q ¡ 2 σ 2 ap q 6 σ 2 ap q ¡ 12 σ 4 ap q ¢ ap q ¡ 2 σ 2 ap q 4 σ 4 ap q 2 σ 2 ap q ¡ 4 σ 4 ap q 4 σ 4 ap q ¡ 8 σ 6 ap q ¢ ap q ¡ ap q ¡ ap q ¡ ap q 2ap q 2 σ 2 ap q ap q ¡ 2ap q ¡ 4 σ 2 ap q 8 σ 2 ap q ¡ 16 σ 4 ap q ¢ ap q ¡ 2 σ 2 ap q ¡ ap q 2 σ 2 ap q ap q 4 σ 2 ap q ¡ 4 σ 2 ap q 4 σ 2 ap q ¡ 8 σ 4 ap q ¢ ap q ¡ ap q ap q ap q ¡ ap q ap q ¡ 2 σ 2 ap q ¢ ap q ¡ 2 σ 2 ap q 4 σ 4 ap q ¡ 8 σ 6 ap q ¢ ap q ¡ ap q 2ap q ¡ 2 σ 2 ap q 2 σ 2 ap q ¡ 6 σ 2 ap q 12 σ 4 ap q ¢ ap q ¡ 2 σ 2 ap q ¡ 2 σ 2 ap q 4 σ 4 ap q 4 σ 2 ap q ¡ 12 σ 4 ap q 24 σ 6 ap q .

Assumption 2 . 3 .

 23 The numerical scheme (1.6) has bounded moments of any order along time, i.e. for all integer k ¥ 0, sup n¥0 Er|X n | 2k s V.

Proposition 4 . 7 (

 47 Simplification rules). Under Assumption 2.2, the two forest patterns gathered in each of the following pairs represent the same differential:

i a ij d p1q j °bi d p2q i a ij d p2q j 1 6 φ3 1 6i © 2 1 12 1 3

 6111 °φP pe i , f I f I pe i qq °bi a ij d p1q j °φP pf I pe i q, f I pe i qq ¡ °bi d p2q °φP pe i , f P pe i , f qq °bi c i d p1q i °φP pe i , p∆fq I pe i qq °bi d p1q

Theorem 5 . 7 .

 57 Using notation and assumptions of Theorem 2.8, if the numerical scheme and the postprocessor can be developed in exotic aromatic B-series of the respective forms 4.5 andErφpG n pxqqs F p qpφq γEAT 1¤|γ|¤p h |γ| apγqF pγqpφq . . . ,

and s 0 and d 0 1

 1

Assumption 2.4. The

  numerical scheme(1.6) has a weak Taylor expansion of the formErφpX 1 q|X 0 xs φpxq hA 0 φpxq h 2 A 1 φpxq ¤ ¤ ¤ ¡ X n |X n s| ¤ Cp1 |X n |qh, |X n 1 ¡ X n | ¤ M n p1 |X n |q c

	(2.3)
	for all φ C V P pR d , Rq, where A i , i 0, 1, 2, . . . are linear differential operators with coeffi-
	cients depending smoothly on the drift f , and its derivatives (and depending on the choice
	of the integrator). In addition, we assume that A 0 coincides with the generator L given in
	(2.2), i.e. it is consistent and has (at least) local order one in the weak sense, A 0 L.
	Remark 2.5. A convenient sufficient condition to satisfy Assumption 2.3 is given in [37,
	Lemma 2.2.2]: if X 0 is deterministic or has bounded moments of all order and the Markov chain pX n q n satisfies
	|ErX n 1 h,
	for C a constant independent of h and M n a random variable whose moments are all bounded
	uniformly with respect to h small enough, then the numerical scheme satisfies Assumption
	2.3. Note that consistent Runge-Kutta type schemes such as (1.10) satisfy Assumption 2.3
	(using the global Lipschitz assumption of f ) and Assumption 2.4.

  in the space R d .

	Theorem 2.6. [3] Assume Assumption 2.1 or 2.2. Consider the one step integrator (1.6)
	and assume that it is ergodic when applied to (1.1). Assume further Assumptions 2.3 and
	2.4. If	A ¦ j ρ V 0, j 2, . . . , p ¡ 1,	(2.4)

  Let a : EAT g Ñ R a map, and let f : R d Ñ R d and φ : R d Ñ R be two smooth functions, then the grafted exotic aromatic B-series Bpaqpφq is a formal series indexed over EAT g defined by

	2	.
	Definition 3.4. Bpaqpφq

γEATg h |γ| apγqF pγqpφq.

  The operator A 1 is now convenient to write with exotic aromatic trees. Applying Theorem 4.1 to (4.1), we deduce A 1 F pγq with

	γ θ 1 2	θσ 2 2	θσ 2	σ 2 2	σ 4 8	.	(4.2)
	2 n n! is exactly the number of ways to gather the grafted nodes by pairs. An
	other example is
	ErF p		qpφqs 3F p q,
	where the coefficient 3 accounts for the number of choices for linking the grafted nodes
	pairwise.

Application. Using Theorem 4.1, we immediately obtain that Runge-Kutta methods (1.10) can be developed in exotic aromatic forests. As a special Runge-Kutta method, we get back to the θ-method

(1.11)

.

  3) is proved by linearity with respect to B.

	pP 2 p2nq

¥ Proof of Theorem 4.1. We consider F pγqpφq as a 2n-multilinear form B γ,φ evaluated in pξ, . . . , ξq (see Definition 3.3). Lemma 4.2 gives ErB γ,φ pξ, . . . , ξqs d i1 ,...,in1

  1, and f i 0, i |V | , . . . , |γ|. With this choice, F pr γqpf, φqp0q $ 0 if and only if r γ γ, thus giving the independence of elementary differentials. Note that this result does not hold if f is assumed gradient due to simplification rules (see Proposition 4.7). now adapt Theorem 2.6 to the context of exotic aromatic forests.

	We

Theorem 4.9. We assume Assumptions 2.1 and 2.3. We consider an ergodic numerical scheme that can be developed in exotic aromatic B-series (and thus has a development of the form (2.3))

ErφpX 1 q|X 0 xs F p qpφq

γEAT 1¤|γ|¤p h |γ| apγqF pγqpφq ¤ ¤ ¤ , (4.5)

  ¡1 ptr 1 uq f J Γpr 1 qψ ¡1 ptr 1 uq φ We recall L F p σ 2

	φ:πpr 2 qÑV 1 ψ:Γpr 2 qÑV 1
	2	q, then we can compute L 2 φ. Using Theorem 5.6, we
	obtain
	F p qpLφq F p		σ 2 2 qpφq
	and
	F p qpLφq F p	2	σ 2 2

F pγ ϕ,ψ qpφq. ¥ Example.

Table 7 :

 7 Order conditions for partitioned Runge-Kutta method with postprocessor for the perturbed equation (See Theorem 5.11). The sums are over all involved indices.

Ecole Normale Supérieure de Rennes and Université de Rennes 1, Campus de Ker Lann, avenue Robert Schumann, F-35170 Bruz, France, Adrien.Laurent@unige.ch

Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland, Gilles.Vilmart@unige.ch

Note that the internal stages Yi depend on n, but this dependence is omitted for brevity of the notation.

In this case, we denote d d p1q and ξn ξ p1q n .

originally named Butcher-series

Such graphs with one cycle are not strictly speaking "trees", they are however called aromatic trees in the literature as an analogy with carbon chemistry.

Acknowledgements

. The authors would like to thank Hans Munthe-Kaas and Olivier Verdier for helpful discussions about an earlier version of this work. This work was partially supported by the Swiss National Science Foundation, grants No. 200021_162404, 200020_178752 and 200020_144313/1.

Appendix

Expression of r γ 2 in Theorem 4.9.