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

Adrien Laurent, Gilles Vilmart

► To cite this version:

Adrien Laurent, Gilles Vilmart. Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs. Mathematics of Computation, 2020, 89 (321), pp.169-202. 10.1090/mcom/3455. hal-04347647

HAL Id: hal-04347647

https://hal.science/hal-04347647v1

Submitted on 15 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs

Adrien Laurent<fnref n="1" />,<fnref n="2" /> and Gilles Vilmart<fnref n="2" />

July 1, 2019

## Abstract

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.

**Keywords:** stochastic differential equations, invariant measure, ergodicity, exotic aromatic trees, order conditions.

AMS subject classification (2010): 60H35, 37M25, 65L06, 41A58

## 1 Introduction

We consider a class of stochastic systems of differential equations of the form

$$dX(t) = f(X(t))dt + \sigma dW(t), \quad (1.1)$$

where $X(t) \in \mathbb{R}^d$ is the solution with initial condition $X_0$ assumed deterministic for simplicity, the vector field $f: \mathbb{R}^d \to \mathbb{R}^d$ is assumed smooth and globally Lipschitz, $\sigma > 0$ is a constant, and $W(t)$ 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 $\mu$ satisfying for all deterministic initial conditions $X_0$ and all smooth test functions $\phi$,

$$\lim_{T \to \infty} \frac{1}{T} \int_0^T \phi(X(s))ds = \int_{\mathbb{R}^d} \phi(x)d\mu(x), \quad \text{almost surely.} \quad (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 $\phi$,

$$\left| \mathbb{E}[\phi(X(t))] - \int_{\mathbb{R}^d} \phi(x)d\mu(x) \right| \le Ce^{-\lambda t} \quad (1.3)$$

<footnote n="1">¹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</footnote>

<footnote n="2">²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</footnote>
for all $t > 0$ where $C = C(\phi, X_0)$ and $\lambda > 0$ are independent of time $t > 0$. We refer to [7, 34, 21, 28, 29] for details. A special case of the class of problems (1.1) is the well-studied overdamped Langevin equation, also called Brownian dynamics,<fnref n="1" />

$$dX(t) = -\nabla V(X(t))dt + \sigma dW(t) \qquad (1.4)$$

where $V: \mathbb{R}^d \to \mathbb{R}$ is a smooth potential, $f = -\nabla V$ and $\sigma > 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\mu(x) = \rho_\infty(x)dx$ where its density $\rho_\infty$ with respect to the Lebesgue measure is given by the Gibbs distribution

$$\rho_{\infty}(x) = Z e^{-\frac{2}{\sigma^2} V(x)} \qquad (1.5)$$

where $Z$ is a normalization constant such that $\int_{\mathbb{R}^d} \rho_\infty(x)dx = 1$. Computing integrals with respect to the invariant measure $\mu$ 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 [44]. However, as highlighted in [38], “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” [38]. 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} = \Psi(X_n, h, \xi_n), \qquad (1.6)$$

where $h$ is a fixed time step size and $\xi_n$ are independent random vectors. We say that the numerical scheme has local weak order $r$ if the weak error after one step satisfies

$$|\mathbb{E}[\phi(X_1)|X_0 = x] - \mathbb{E}[\phi(X(h))|X(0) = x]| \le Ch^{r+1}, \qquad (1.7)$$

where $C = C(\phi, x)$ is independent of $h$ assumed small enough and $\phi$ 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$, $|\mathbb{E}(\phi(X_n)) - \mathbb{E}(\phi(X(t_n)))| \le Ch^r$, as shown in [36] (see [37, Chap. 2.2]). The numerical method is called ergodic if it has a unique invariant probability law $\mu^h$ with finite moments of any order and

$$\lim_{N \to \infty} \frac{1}{N+1} \sum_{n=0}^{N} \phi(X_n) = \int_{\mathbb{R}^d} \phi(x) d\mu^h(x), \quad \text{almost surely,}$$

for all deterministic initial conditions $X_0 = x$ and all test functions $\phi$. We say that the numerical method has order $p$ with respect to the invariant measure of (1.1) if

$$\left| \int_{\mathbb{R}^d} \phi(x) d\mu^h(x) - \int_{\mathbb{R}^d} \phi(x) d\mu(x) \right| \le Ch^p, \qquad (1.8)$$

<footnote n="1">¹Up to a time transformation, one could for simplicity fix the value of $\sigma$ (e.g. $\sigma = 1$ or $\sqrt{2}$), however we choose not to fix it to better distinguish the various order conditions.</footnote>
where C is independent of h assumed small enough. Under appropriate assumptions on the numerical scheme (see for instance [21]), one also obtains the following exponential estimate similar to (1.3) (possibly with a different constant $\lambda > 0$),

$$ \left| \mathbb{E}[\phi(X_n)] - \int_{\mathbb{R}^d} \phi(x) d\mu(x) \right| \le K e^{-\lambda t_n} + Ch^p, \quad (1.9) $$

where $K = K(\phi, x)$, $C = C(\phi)$ 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 \ge r$, see in particular [33] 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 [9, 30, 31, 4]. In [3, 4], 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 [46] and recent developments of the theory of backward error analysis and modified differential equations in the stochastic context [48, 2, 21, 28, 29], a major tool in the area of deterministic geometric numerical integration [23]. In [47] for finite dimensions and in [10] in the context of parabolic stochastic partial differential equations, this approach is combined with the idea of processing from Butcher [13], 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 [1] 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 [3] and for postprocessed integrators in [47, 10, 1]. Analogously to [41] (we study here the additive noise case), we consider in this paper Runge-Kutta methods of the form<fnref n="2" />

$$ \begin{aligned} Y_i &= X_n + h \sum_{j=1}^{s} a_{ij} f(Y_j) + \sum_{k=1}^{l} d_i^{(k)} \sigma \sqrt{h} \xi_n^{(k)}, \quad i = 1, \dots, s, \\ X_{n+1} &= X_n + h \sum_{i=1}^{s} b_i f(Y_i) + \sigma \sqrt{h} \xi_n^{(1)}, \end{aligned} \quad (1.10) $$

where $a_{ij}, b_i, d_i^{(k)}$ are the coefficients defining the Runge-Kutta scheme, and $\xi_n^{(k)} \sim \mathcal{N}(0, I_d)$ 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

<footnote n="2">²Note that the internal stages $Y_i$ depend on $n$, but this dependence is omitted for brevity of the notation.</footnote>
stages. However we shall often consider $l = 1$ random vector per internal stage<fnref n="3" />, which is sufficient to achieve order $p = 2$ or 3 for the invariant measure. In particular, we shall consider the $\theta$-method as an illustrative example in Sections 4 and 5 and recover known results on its accuracy. It is defined for $\theta$ fixed as

$$X_{n+1} = X_n + h(1-\theta)f(X_n) + h\theta f(X_{n+1}) + \sigma\sqrt{h}\xi_n, \quad (1.11)$$

For $\theta = 0$, we get the explicit Euler-Maruyama method while the scheme is implicit for $\theta \neq 0$. It can be put in Runge-Kutta form (1.10) for $l = 1$ with the following coefficients.

$$\frac{c}{b} \begin{vmatrix} A & d \end{vmatrix} = \frac{\begin{vmatrix} 0 & 0 & 0 \\ 1 & 1-\theta & \theta \\ 1-\theta & \theta & 1 \end{vmatrix}}{\begin{vmatrix} 1 & 1-\theta & \theta \end{vmatrix}}$$

The usage of trees and B-series<fnref n="4" /> is known as a powerful standard tool for the numerical analysis of differential equations. B-series where introduced by Hairer and Wanner in [24] based on the work of Butcher [14], and are now exposed in many articles and books [23, 15], see also the presentation in [16, 45]. In the last decades, several works extended trees and B-series to the stochastic context, we mention in particular Burrage and Burrage [11] and Komori, Mitsui and Sugiura [27] who first introduced stochastic trees and B-series for studying the order conditions of strong convergence of SDEs, and [12, 40, 42, 41, 19, 43, 20] 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 [6]. The recent work [5] 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 non-reversible perturbation as in [32, 22]. The proposed algebraic framework relies on aromatic B-series, a generalisation of B-series introduced in [39, 35] (see also the presentation in [8]) 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 [17] 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).

<footnote n="3">³In this case, we denote $d = d^{(1)}$ and $\xi_n = \xi_n^{(1)}$.</footnote>

<footnote n="4">⁴originally named Butcher-series</footnote>
## 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) satisfies the ergodicity properties (1.2) and (1.3) (see [25] in the more general context of SDEs with multiplicative noise).

**Assumption 2.1.** The vector field $f: \mathbb{R}^d \to \mathbb{R}^d$ is globally Lipchitz and $C^\infty$, and there exist $C_1, C_2 > 0$ such that for all $x \in \mathbb{R}^d$,

$$x^T f(x) \leq -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^\infty$ potential $V: \mathbb{R}^d \to \mathbb{R}$ such that $f(x) = -\nabla V(x)$ is globally Lipschitz and there exist $C_1 > 0$ and $C_2$ such that for all $x \in \mathbb{R}^d$, $V(x) \geq C_1 x^T x - C_2$.

We recall that under Assumption 2.2, the density of the unique invariant measure is given by $\rho_\infty = Z \exp(-\frac{2V}{\sigma^2})$ where $Z$ is such that $\int_{\mathbb{R}^d} \rho_\infty(x)dx = 1$ but $Z$ is not numerically known in general. We note that $\nabla \rho_\infty = \rho_\infty \frac{2}{\sigma^2} f$, equivalently $\nabla(\log \rho_\infty) = \frac{2}{\sigma^2} f$.

We denote $C_p^\infty(\mathbb{R}^d, \mathbb{R})$ the vector space of $C^\infty$ functions such that all partial derivatives $\phi$ up to all orders have a polynomial growth of the form

$$|\phi(x)| \leq C(1 + |x|^s)$$

for some constants $s$ and $C$ independent of $x$ (but depending on the order of differentiation). For $\phi \in C_p^\infty(\mathbb{R}^d, \mathbb{R})$ we define $u(x, t) = E[\phi(X(t))|X(0) = x]$. A classical tool for the study of (1.1) is the backward Kolmogorov equation [37, Chap. 2], which states that $u(x, t)$ solves the following deterministic parabolic PDE in $\mathbb{R}^d$,

$$\frac{\partial u}{\partial t} = \mathcal{L}u, \quad u(x, 0) = \phi(x), \quad x \in \mathbb{R}^d, t > 0, \qquad (2.1)$$

where the generator $\mathcal{L}$ is defined as

$$\mathcal{L}\phi = f \cdot \nabla\phi + \frac{\sigma^2}{2}\Delta\phi \qquad (2.2)$$

where $\Delta\phi = \sum_{i=1}^d \frac{\partial^2\phi}{\partial x_i^2}$ denotes the Laplace operator. We recall that, under Assumption 2.1 or 2.2, the density of the invariant measure satisfies

$$\mathcal{L}^* \rho_\infty = 0$$

where $\mathcal{L}^*\phi = -\text{div}(f\phi) + \frac{\sigma^2}{2}\Delta\phi$ is the $L^2$-adjoint of $\mathcal{L}$.

We make the following natural assumptions on the numerical integrator (1.6).

**Assumption 2.3.** The numerical scheme (1.6) has bounded moments of any order along time, i.e. for all integer $k \geq 0$,

$$\sup_{n \geq 0} \mathbb{E}[|X_n|^{2k}] < \infty.$$
**Assumption 2.4.** The numerical scheme (1.6) has a weak Taylor expansion of the form

$$ \mathbb{E}[\phi(X_1)|X_0 = x] = \phi(x) + h\mathcal{A}_0\phi(x) + h^2\mathcal{A}_1\phi(x) + \dots \quad (2.3) $$

for all $\phi \in C_p^\infty(\mathbb{R}^d, \mathbb{R})$, where $\mathcal{A}_i$, $i = 0, 1, 2, \dots$ are linear differential operators with coefficients depending smoothly on the drift $f$, and its derivatives (and depending on the choice of the integrator). In addition, we assume that $\mathcal{A}_0$ coincides with the generator $\mathcal{L}$ given in (2.2), i.e. it is consistent and has (at least) local order one in the weak sense, $\mathcal{A}_0 = \mathcal{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 $(X_n)_n$ satisfies

$$ |\mathbb{E}[X_{n+1} - X_n|X_n]| \le C(1 + |X_n|)h, \quad |X_{n+1} - X_n| \le M_n(1 + |X_n|)\sqrt{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 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 [3] in the more general context of ergodic SDEs with multiplicative noise based on backward error analysis results in [21] on the torus and generalizations [28, 29] in the space $\mathbb{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

$$ \mathcal{A}_j^* \rho_\infty = 0, \quad j = 2, \dots, p-1, \quad (2.4) $$

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 $\mathcal{A}_j = \mathcal{L}^{j+1}/(j+1)!$, $j = 2, \dots, p-2$, then Theorem 2.6 is an immediate consequence of the Talay-Tubaro expansion of the error [36, 46] (see also [37, Chap. 2.2, 2.3]) given by

$$ \int_{\mathbb{R}^d} \phi(x) d\mu^h(x) - \int_{\mathbb{R}^d} \phi(x) d\mu(x) = \lambda_p h^p + \mathcal{O}(h^{p+1}) $$

with

$$ \lambda_p = \int_0^{+\infty} \int_{\mathbb{R}^d} \left( A_p - \frac{1}{(p+1)!} \mathcal{L}^{p+1} \right) u(y,t) \rho_\infty(y) dydt $$

where $u(x,t)$ is the solution of (2.1). Indeed, considering the $L^2$-adjoint of the operator $A_p - \frac{1}{(p+1)!} \mathcal{L}^{p+1}$ and using $\mathcal{L}^* \rho_\infty = 0$ and $\mathcal{A}_p^* \rho_\infty = 0$ then yield $\lambda_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.** [47] Assume the hypotheses of Theorem 2.6 and consider a postprocessor

$$ \overline{X}_n = G_n(X_n) $$

that admits the following weak Taylor expansion for all $\phi \in C_p^\infty(\mathbb{R}^d, \mathbb{R})$,

$$ \mathbb{E}[\phi(G_n(x))] = \phi(x) + \sum_{i=1}^{p-1} \alpha_i h^i \mathcal{L}^i \phi(x) + h^p \overline{\mathcal{A}}_p \phi(x) + \dots, \quad (2.5) $$

for some constants $\alpha_i$ and a linear differential operator $\overline{\mathcal{A}}_p$. Assume further that

$$ (\mathcal{A}_p + [\mathcal{L}, \overline{\mathcal{A}}_p])^* \rho_\infty = 0 \quad (2.6) $$

where $[\mathcal{L}, \overline{\mathcal{A}}_p] = \mathcal{L}\overline{\mathcal{A}}_p - \overline{\mathcal{A}}_p\mathcal{L}$ is the Lie bracket. Then $\overline{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 $\overline{X}_n$.

Theorem 2.8 is stated and proved in [47] in the special case $\alpha_i = 0$, $i = 1, \dots, p-1$ in (2.5). However, the proof for non zero $\alpha_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

$$ \int_{\mathbb{R}^d} (\mathcal{A}_j \phi) \rho_\infty dx = 0, \quad j = 2, \dots, p-1, $$

$$ \int_{\mathbb{R}^d} (\mathcal{A}_p \phi + [\mathcal{L}, \overline{\mathcal{A}}_p] \phi) \rho_\infty dx = 0, $$

for all test function $\phi \in C_p^\infty(\mathbb{R}^d, \mathbb{R})$. 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.

# 3 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 [17] and rely on the presentation in [8].

## 3.1 Aromatic trees and forests

We first consider directed graphs $\gamma = (V, E)$ with $V$ a finite set of nodes and $E \subset V \times V$ the set of directed edges. If $(v, w) \in E$, we say that the edge is going from $v$ to $w$, and $v$ is called a predecessor of $w$. Two directed graphs $(V_1, E_1)$ and $(V_2, E_2)$ are equivalent if there exists a bijection $\varphi: V_1 \to V_2$ with $(\varphi \times \varphi)(E_1) = 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,

<figure><img src="image_2.png" /></figure>

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 trees<fnref n="5" /> 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 $\mathcal{AT} = \{\bullet, \circ, \circ, \circ, \dots\}$ 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 $\gamma = (V, E) \in \mathcal{AT}$, and let $f: \mathbb{R}^d \to \mathbb{R}^d$ be a smooth function. We denote $\pi(v) = \{w \in V, (w, v) \in E\}$ the set of all predecessors of the node $v \in V$ and $r$ the root of $\gamma$. We also call $V^0 = V \setminus \{r\} = \{v_1, \dots, v_m\}$ the other nodes of $\gamma$. Finally we introduce the notation $I_{\pi(v)} = (i_{q_1}, \dots, i_{q_s})$ where the $q_k$ are the predecessors of $v$, and

$$ \partial_{I_{\pi(v)}} f = \frac{\partial^s f}{\partial x_{i_{q_1}} \dots \partial x_{i_{q_s}}}. $$

Then $F(\gamma)$ is defined as

$$ F(\gamma)(f) = \sum_{i_{v_1}, \dots, i_{v_m} = 1}^{d} \left( \prod_{v \in V^0} \partial_{I_{\pi(v)}} f_{i_v} \right) \partial_{I_{\pi(r)}} f. $$

**Example.** Let $\gamma = \circ \circ \circ \circ \circ$ and $\tilde{\gamma} = \circ \circ \circ \circ \circ$ in $\mathcal{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 differentials are respectively $F(\gamma)(f) = \sum_{i,j,k,l,m=1}^{d} \partial_m f_m f_i \partial_i f_j f_k \partial_{j,k} f_l \partial_l f = \text{div}(f) \cdot f' f''(f'f, f)$ and $F(\tilde{\gamma})(f) = \sum_{i,j,k,l,m=1}^{d} \partial_l f_m \partial_{m,k} f_l f_k f_i \partial_i f_j \partial_j f = \sum_{m=1}^{d} f'_m ((\partial_m f)'(f)) \cdot f' f' f$.

## 3.2 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 $(V, E)$ be an aromatic forest and $L$ be a finite list of pairs of elements of $V$ (possibly with duplicates), then $\gamma = (V, E, L)$ is an exotic aromatic forest. The elements of $L$ are called lianas and correspond to non-oriented 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$, $\Gamma(v)$ 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 $\gamma = (V, E, L)$ is a connected component of the associated aromatic forest $(V, E)$. We call $\mathcal{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 ◯.

<figure><img src="image_3.png" /></figure>

<footnote n="5">⁵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.</footnote>
The definition of elementary differentials is extended as follows.

**Definition 3.2.** Let $\gamma = (V, E, L) \in \mathcal{EAT}$, and let $f : \mathbb{R}^d \to \mathbb{R}^d$ be a smooth function. We name $r$ the root of $\gamma$ and $V^0 = V \setminus \{r\} = \{v_1, \dots, v_m\}$ the other nodes of $\gamma$. We denote $l_1, \dots, l_s$ the elements of $L$ and for $v \in V$, $J_{\Gamma(v)}$ the multiindex $(j_{l_1}, \dots, j_{l_s})$ where $\Gamma(v) = \{l_{x_1}, \dots, l_{x_t}\}$. Then $F(\gamma)$ is defined as

$$F(\gamma)(f) = \sum_{i_{v_1}, \dots, i_{v_m} = 1}^{d} \sum_{j_{l_1}, \dots, j_{l_s} = 1}^{d} \left( \prod_{v \in V^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{i_v} \right) \partial_{I_{\pi(r)}} \partial_{J_{\Gamma(r)}} f.$$

**Examples.** The differential that corresponds to the rooted tree with a single node and a single liana is $F(\gamma)(f) = \Delta f$. We can also represent as exotic aromatic forest more

complicated derivatives. For instance, let $\gamma = \gamma$, then

$$F(\gamma)(f) = \sum_{i,j,k=1}^{d} \mathrm{div}(\partial_i f) \cdot f'((\partial_{kl} f)'(f''(\partial_{ij} f, \partial_{kl} f))).$$

## 3.3 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 $\phi$, and it has leafs (nodes without predecessors) that represent a random standard normal vector $\xi$. Note that these new trees can be seen as bi-coloured trees in the context of P-series (see [23, 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 \setminus V_g$, then $\gamma = (V, E, L)$ is a grafted exotic aromatic forest. We define as before the grafted exotic aromatic trees and grafted exotic aromatic rooted forests, that we denote $\mathcal{EAT}_g$.

If $\gamma = (V, E, L)$ is a grafted exotic aromatic rooted forest, $\phi : \mathbb{R}^d \to \mathbb{R}$ a smooth function, and $\xi$ a random vector of $\mathbb{R}^d$ whose components are independent and follow a standard normal law, the associated elementary differential of $\gamma$ is, with the same notation as Definition 3.2 and $V^0 = V \setminus (V_g \cup \{r\})$,

$$F(\gamma)(f, \phi, \xi) = \sum_{i_{v_1}, \dots, i_{v_m} = 1}^{d} \sum_{j_{l_1}, \dots, j_{l_s} = 1}^{d} \left( \prod_{v \in V^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{i_v} \right) \left( \prod_{v \in V_g} \xi_{i_v} \right) \partial_{I_{\pi(r)}} \partial_{J_{\Gamma(r)}} \phi.$$

**Example.** The differential associated to the forest is $F(\gamma)(f, \phi, \xi) = \phi'(f''(\xi, \xi))$.

If $\gamma$ is such that $V_g$ is empty, we recover the exotic aromatic forests of Definition 3.2, where $\phi$ 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 $\phi$. This definition can be straightforwardly extended on non-rooted exotic aromatic forests. For brevity of notation, we also write $F(\gamma)(\phi)$ instead of $F(\gamma)(f, \phi, \xi)$. We note that $\phi \to F(\gamma)(\phi)$ is a linear differential operator (dependent on $f$ and $\xi$).
## 3.4 Grafted exotic aromatic B-series

In this section, we adapt the formalism of aromatic B-series of [39] 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 $|\gamma|$ of a tree $\gamma \in \mathcal{E}\mathcal{A}\mathcal{T}_g$. We denote $N(\gamma)$ the number of nodes, $N_l(\gamma)$ the number of lianas, $N_c(\gamma)$ the number of grafted nodes and $N_v(\gamma) = N(\gamma) - N_c(\gamma) - 1$ the number of nodes that are non grafted and different from the root, then

$$|\gamma| = N_v(\gamma) + N_l(\gamma) + \frac{N_c(\gamma)}{2}.$$

**Definition 3.4.** Let $a : \mathcal{E}\mathcal{A}\mathcal{T}_g \to \mathbb{R}$ a map, and let $f : \mathbb{R}^d \to \mathbb{R}^d$ and $\phi : \mathbb{R}^d \to \mathbb{R}$ be two smooth functions, then the grafted exotic aromatic B-series $B(a)(\phi)$ is a formal series indexed over $\mathcal{E}\mathcal{A}\mathcal{T}_g$ defined by

$$B(a)(\phi) = \sum_{\gamma \in \mathcal{E}\mathcal{A}\mathcal{T}_g} h^{|\gamma|}a(\gamma)F(\gamma)(\phi).$$

We extend the definition of $F$ on $\text{Vect}(\mathcal{E}\mathcal{A}\mathcal{T}_g)$ by writing

$$F\left(\sum_{\gamma \in \mathcal{E}\mathcal{A}\mathcal{T}_g} h^{|\gamma|}a(\gamma)\gamma\right)(\phi) = B(a)(\phi).$$

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 $a(\gamma)$ of standard B-series are sometimes renormalized as $\frac{a(\gamma)}{\rho(\gamma)}$ where $\rho$ is a function determined by the symmetries of the associated forest. If $\rho$ is appropriately chosen, it greatly simplifies the composition laws of (aromatic) B-series (see [23, 16, 8]). Finding the best definition of $\rho$ for this exotic extension of B-series is out of the scope of this paper.

## 3.5 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 [39, 35]. We consider exotic aromatic rooted forests $\gamma$ where the differential associated to the root is $f$. As the function $f$ is no longer fixed, we denote the associated differential $F(\gamma)(f)$. Also we adapt the definition of exotic aromatic B-series to this change. First we add a new tree: the empty tree $\emptyset$. The function $F$ is then extended on $\mathcal{E}\mathcal{A}\mathcal{T}_g \cup \{\emptyset\}$ by $F(\emptyset)(f) = \text{Id}_{\mathbb{R}^d}$. Then, for a function $a : \mathcal{E}\mathcal{A}\mathcal{T} \cup \{\emptyset\} \to \mathbb{R}$, the associated exotic aromatic B-series is

$$B(a)(f) = \sum_{\gamma \in \mathcal{E}\mathcal{A}\mathcal{T} \cup \{\emptyset\}} a(\gamma)F(\gamma)(f).$$

We study (exotic) aromatic B-series $B(a)$ with $a(\emptyset) = 1$. We call these (exotic) aromatic B-series methods. Let $G$ be a subgroup of $\text{GL}_d(\mathbb{R}) \times \mathbb{R}^d$, let the action of an element $(A, b) \in G$ on $\mathbb{R}^d$ be $x \mapsto Ax + b$, and let the action on a vector field $f : \mathbb{R}^d \to \mathbb{R}^d$ be

$$(A, b) * f)(x) := Af(A^{-1}(x - b)).$$
where $\mathbb{E}[\xi_j\xi_k] = 0$ for $j \neq k$ because of the independence of the $\xi_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:

$$ 
\begin{aligned}
\mathbb{E}[F(\begin{smallmatrix} \times & \times \\ \bullet & \bullet \end{smallmatrix}) (\phi)] &= \mathbb{E}[\phi^{(4)}(\xi, \xi, \xi, \xi)] = \sum_{i,j,k,l} \partial_{i,j,k,l} \phi \mathbb{E}[\xi_i \xi_j \xi_k \xi_l] \\
&= \sum_i \partial_{i,i,i} \phi \mathbb{E}[\xi_i^4] + 3 \sum_{\substack{i,j \\ i \neq j}} \partial_{i,i,j,j} \phi \mathbb{E}[\xi_i^2] \mathbb{E}[\xi_j^2] = 3 \sum_{i,j} \partial_{i,i,j,j} \phi \\
&= 3F(\begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix}) (\phi).
\end{aligned}
$$

Let us now comment this computation. The interesting fact is that $\mathbb{E}[\xi_i^4] = 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 $\xi$. 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 $\gamma \in \mathcal{E}\mathcal{A}\mathcal{T}_g$ be a grafted exotic aromatic rooted forest with an even number of grafted nodes $2n$, let $\phi: \mathbb{R}^d \to \mathbb{R}$ be a smooth function, and let $V^\times = \{c_1, \dots, c_{2n}\}$ be the set of grafted nodes of $\gamma$. We call $\mathcal{P}_2(2n)$ the set of partitions by pair of $\{1, \dots, 2n\}$, i.e. the set of surjections $p: \{1, \dots, 2n\} \to \{1, \dots, n\}$ such that the preimage of each singleton has exactly two elements and the minima of those preimages follow an ascending order $(\min(p^{-1}(\{i\})) < \min(p^{-1}(\{j\}))$ for $i < j$. Finally we define $\varphi_\gamma: \mathcal{P}_2(2n) \to \mathcal{E}\mathcal{A}\mathcal{T}$ the application that maps the partition $p$ of $\gamma$ to the aromatic forest where the grafted nodes are linked by lianas according to $p$. Then, the expectation of $F(\gamma)(\phi)$ is given by

$$ \mathbb{E}[F(\gamma)(\phi)] = \sum_{p \in \mathcal{P}_2(2n)} F(\varphi_\gamma(p))(\phi). $$

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 $\gamma \in \{\bullet, \bullet^*, \bullet^\times, \bullet^\times, \bullet^\times, \dots\}$ the tree with only a root, $2n$ grafted nodes and no liana, then

$$ \mathbb{E}[F(\gamma)(\phi)] = \frac{(2n)!}{2^n n!} \Delta^n \phi $$

The integer $\frac{(2n)!}{2^n n!}$ is exactly the number of ways to gather the grafted nodes by pairs. An other example is

$$ \mathbb{E}[F(\begin{smallmatrix} \times & \times \\ \bullet & \bullet \end{smallmatrix}) (\phi)] = 3F(\begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix}), $$

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). The operator $\mathcal{A}_1$ is now convenient to write with exotic aromatic trees. Applying Theorem 4.1 to (4.1), we deduce $\mathcal{A}_1 = F(\gamma)$ with

$$ \gamma = \theta^{\bullet} + \frac{1}{2} \bullet \bullet + \frac{\theta \sigma^2 \bullet \bullet}{2} + \theta \sigma^2 \bullet \bullet + \frac{\sigma^2 \bullet \bullet}{2} + \frac{\sigma^4 \bullet \bullet \bullet}{8} \quad (4.2) $$

Theorem 4.1 follows from the following lemma, which is an extension of the Isserlis theorem [26] to the case of multilinear mappings. The Isserlis theorem states that if $\chi$ is a 2n-dimensional Gaussian random vector with mean zero and arbitrary covariance, then

$$ \mathbb{E}\left[\prod_{i=1}^{2n} \chi_i\right] = \sum_{p \in \mathcal{P}_2(2n)} \prod_{\substack{i<j \\ p(i)=p(j)}} \mathbb{E}[\chi_i \chi_j]. $$

For $n=2$, it gives $\mathbb{E}[\chi_1\chi_2\chi_3\chi_4] = \mathbb{E}[\chi_1\chi_2]\mathbb{E}[\chi_3\chi_4] + \mathbb{E}[\chi_1\chi_3]\mathbb{E}[\chi_2\chi_4] + \mathbb{E}[\chi_1\chi_4]\mathbb{E}[\chi_2\chi_3]$.

**Lemma 4.2.** Let $B: \mathbb{R}^d \times \dots \times \mathbb{R}^d = \mathbb{R}^{2nd} \to \mathbb{R}$ be a 2n-multilinear form, and let $\xi$ be a Gaussian vector $\mathcal{N}(0, I_d)$, then

$$ \mathbb{E}[B(\xi, \dots, \xi)] = \sum_{p \in \mathcal{P}_2(2n)} \sum_{i_1, \dots, i_n = 1}^{d} B(e_{i_p}), \quad (4.3) $$

with $e_{i_p} = (e_{i_{p(1)}}, \dots, e_{i_{p(2n)}})$, and we recall that $e_1, \dots, e_d$ denotes the canonical basis of $\mathbb{R}^d$.

*Proof.* For the particular case of an elementary multilinear form $B_\sigma: (x_1, \dots, x_{2n}) \mapsto \prod_{j=1}^{2n} (x_j)_{\sigma(j)}$ where $\sigma: \{1, \dots, 2n\} \to \{1, \dots, d\}$ is a given mapping and $(x_j)_{\sigma(j)}$ denotes the $\sigma(j)$'s component of $x_j \in \mathbb{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.3) is proved by linearity with respect to $B$. $\square$

*Proof of Theorem 4.1.* We consider $F(\gamma)(\phi)$ as a 2n-multilinear form $B_{\gamma,\phi}$ evaluated in $(\xi, \dots, \xi)$ (see Definition 3.3). Lemma 4.2 gives

$$ \mathbb{E}[B_{\gamma,\phi}(\xi, \dots, \xi)] = \sum_{p \in \mathcal{P}_2(2n)} \sum_{i_1, \dots, i_n = 1}^{d} B_{\gamma,\phi}(e_{i_p}) = \sum_{p \in \mathcal{P}_2(2n)} F(\varphi_{\gamma}(p))(\phi), $$

because $B_{\gamma,\phi}(e_{i_p})$ is the differential $F(\gamma)$ where we differentiate the non-grafted nodes linked to grafted nodes in the directions given by $e_{i_p}$, and $F(\varphi_{\gamma}(p))(\phi)$ is obtained by summing $B_{\gamma,\phi}(e_{i_p})$ over all the indices. $\square$

## 4.2 Integration by parts of the exotic aromatic forests

The goal of this section is to integrate by parts $\int_{\mathbb{R}^d} F(\gamma)(\phi)\rho_\infty dx$, for $\gamma$ an exotic aromatic rooted forest, in order to write it in the form $\int_{\mathbb{R}^d} \phi' \tilde{f} \rho_\infty dx$ for a certain sum of elementary differentials $\tilde{f}$. The idea is to transform a high order differential operator $\mathcal{A}: \phi \to F(\gamma)(\phi)$ into a differential operator $\phi \to \phi' \tilde{f}$ of order 1 such that $\mathcal{A}^* \rho_\infty = -\operatorname{div}(\tilde{f} \rho_\infty)$. 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.

$$
\begin{align}
\int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx &= \sum_{i,j} \int_{\mathbb{R}^d} \frac{\partial^2 \phi}{\partial x_i \partial x_j} f_i f_j \rho_{\infty} dx \nonumber \\
&= -\sum_{i,j} \left[ \int_{\mathbb{R}^d} \frac{\partial \phi}{\partial x_j} \frac{\partial f_i}{\partial x_i} f_j \rho_{\infty} dx + \int_{\mathbb{R}^d} \frac{\partial \phi}{\partial x_j} f_i \frac{\partial f_j}{\partial x_i} \rho_{\infty} dx \right. \tag{4.4} \\
&\qquad \left. + \int_{\mathbb{R}^d} \frac{\partial \phi}{\partial x_j} f_i f_j \frac{\partial \rho_{\infty}}{\partial x_i} dx \right], \nonumber
\end{align}
$$

where we integrated by parts; note that the boundary term vanishes using the growth assumptions on $\phi$.

**Notation.** We denote $g = \log(\rho_\infty)$, then $\nabla \rho_\infty = (\nabla g) \rho_\infty$.

We have

$$
\int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx = - \int_{\mathbb{R}^d} \mathrm{div}(f) \phi' f \rho_{\infty} dx - \int_{\mathbb{R}^d} \phi' f' f \rho_{\infty} dx - \int_{\mathbb{R}^d} g' f \phi' f \rho_{\infty} dx.
$$

By writing $\tilde{f} = -(\mathrm{div}(f)f + f'f + g'ff)$, we deduce

$$
\int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx = \int_{\mathbb{R}^d} \phi' \tilde{f} \rho_{\infty} dx.
$$

We see that even for a simple forest, the integration by parts requires some calculations and a new term appears: the function $g = \log(\rho_\infty)$, 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 $\mathcal{EAT}$ is extended to include these new aromas.

If $\gamma = (V, E, L) \in \mathcal{EAT}$ is an exotic aromatic rooted forest whose set of aromatic roots is $V^1 \subset V$, if $V^0 = V \setminus (V^1 \cup \{r\})$ and $V = \{r, v_1, \dots, v_m\}$ where $r$ is the root of $\gamma$, if $l_1, \dots, l_s$ are the elements of $L$, then $F(\gamma)$ is defined as

$$
F(\gamma)(\phi) = \sum_{i v_1, \dots, i v_m = 1}^{d} \sum_{j l_1, \dots, j l_s = 1}^{d} \left( \prod_{v \in V^1} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} g \right) \left( \prod_{v \in V^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{i_v} \right) \partial_{I_{\pi(r)}} \partial_{J_{\Gamma(r)}} \phi.
$$

One can also extend the definition of $\mathcal{EAT}_g$ in the same way so that it includes aromatic roots.

**Examples.** $F(\bullet, \bullet)(\phi) = g'f\phi'f$ and $F(\circ, \circ)(\phi) = \sum_{i,j} (\partial_j g)^2 \partial_i g \partial_i \phi$.

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

$$
\int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx = - \int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx - \int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx - \int_{\mathbb{R}^d} F(\cdot, \cdot)(\phi) \rho_{\infty} dx.
$$

We notice that integrating by parts $F(\cdot, \cdot)(\phi)$ 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 $\gamma \in \mathcal{EAT}$, we choose a direction to integrate by parts, i.e. an edge or a liana e connected to the root r. Then

$$ \int_{\mathbb{R}^d} F(\gamma)(\phi) \rho_{\infty} dx = - \sum_{\tilde{\gamma} \in U(\gamma, e)} \int_{\mathbb{R}^d} F(\tilde{\gamma})(\phi) \rho_{\infty} dx, $$

where $U(\gamma, e)$ 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 $\gamma$ and e an edge ingoing to n or a liana connected to n, then the same result holds if we replace $U(\gamma, e)$ by $U(\gamma, n, e)$, 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 = \{v_1, \dots, v_m\}$ the other nodes of $\gamma$. We suppose for simplicity that $\gamma$ does not have any aromatic root, otherwise the proof can be adapted straightforwardly. We denote $l_1, \dots, l_s$ the elements of L and $v_{k_1}, \dots, v_{k_p}$ the elements of $\pi(r)$. We choose to integrate by parts in the direction of the edge $x_{i_{v_{k_1}}}$.

$$ \begin{aligned} \int_{\mathbb{R}^d} F(\gamma)(\phi) \rho_{\infty} dx &= \sum_{i_{v_1}, \dots, i_{v_m} = 1}^{d} \sum_{j_{l_1}, \dots, j_{l_s} = 1}^{d} \int_{\mathbb{R}^d} \left( \prod_{v \in V^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{iv} \right) \partial_{i_{v_{k_1}}} \dots \partial_{i_{v_{k_p}}} \partial_{J_{\Gamma(r)}} \phi \rho_{\infty} dx \\ &= - \sum_{i_{v_1}, \dots, i_{v_m} = 1}^{d} \sum_{j_{l_1}, \dots, j_{l_s} = 1}^{d} \left[ \sum_{u \in V^0} \int_{\mathbb{R}^d} \left( \prod_{v \in V^0 \setminus \{u\}} \partial_{I_{\pi(u)}} \partial_{J_{\Gamma(u)}} f_{iv} \right) \partial_{i_{v_{k_1}}} \partial_{I_{\pi(u)}} \partial_{J_{\Gamma(u)}} f_{iu} \partial_{i_{v_{k_2}}} \dots \partial_{i_{v_{k_p}}} \partial_{J_{\Gamma(r)}} \phi \rho_{\infty} dx \right. \\ & \qquad \left. + \int_{\mathbb{R}^d} \left( \prod_{v \in V^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{iv} \right) \partial_{i_{v_{k_2}}} \dots \partial_{i_{v_{k_p}}} \partial_{J_{\Gamma(r)}} \phi \partial_{i_{v_1}} g \rho_{\infty} dx \right] \end{aligned} $$

Each term of the sum on $u \in V^0$ is the differential associated to the forest $\gamma_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_{\Gamma(r)}$ instead of $I_{\pi(r)}$. $\square$

**Definition 4.6.** Let $\gamma_1$ and $\gamma_2$ be two exotic aromatic B-series, we define the equivalence relation $\sim$ and write $\gamma_1 \sim \gamma_2$ if we can transform $\gamma_1$ into $\gamma_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

$$ \bullet \sim -\bullet - \bullet \bullet - \bullet \circ \circ $$
**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 [23, Chap. 3, Exercise 3]. Indeed, for a given exotic aromatic forest $\gamma = (V, E, L)$, we take a bijective numbering $n_V : V \setminus \{r\} \to \{1, 2, \dots, |V| - 1\}$ of the nodes and another one for the lianas $n_L : L \to \{|V|, \dots, |\gamma|\}$, then we define

$$ \phi(x) = \prod_{v \in \pi(r)} x_{n_V(v)} \prod_{l \in \Gamma(r)} x_{n_L(l)}, \\ f_i(x) = \prod_{v \in \pi(n_v^{-1}(i))} x_{n_V(v)} \prod_{l \in \Gamma(n_v^{-1}(i))} x_{n_L(l)}, \quad i = 1, \dots, |V| - 1, $$

and $f_i = 0$, $i = |V|, \dots, |\gamma|$. With this choice, $F(\tilde{\gamma})(f, \phi)(0) \neq 0$ if and only if $\tilde{\gamma} = \gamma$, 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).

We now adapt Theorem 2.6 to the context of exotic aromatic forests.

**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))

$$ \mathbb{E}[\phi(X_1)|X_0 = x] = F(\bullet)(\phi) + \sum_{\substack{\gamma \in \mathcal{A}T \\ 1 \le |\gamma| \le p}} h^{|\gamma|} a(\gamma) F(\gamma)(\phi) + \dots, \quad (4.5) $$

for $p \in \mathbb{N}$. We denote $A_i = F(\gamma_i)$. If $\gamma_i \sim \tilde{\gamma}_i$ and $F(\tilde{\gamma}_i) = 0$ for all $1 \le 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

$$ \tilde{\gamma}_0 = \left( a\left(\bullet\right) - \frac{2}{\sigma^2} a\left(\bullet\right) \right)\bullet, $$

and in addition we have for order 2,

$$ \tilde{\gamma}_1 = \left( a(\bullet) - \frac{2}{\sigma^2} a(\bullet) + \frac{2}{\sigma^2} a(\bullet) - \frac{4}{\sigma^4} a(\bullet) \right) \bullet + \left( a(\bullet) - a(\bullet) + a(\bullet) - \frac{2}{\sigma^2} a(\bullet) \right) \bullet \\ + \left( a(\bullet) - \frac{2}{\sigma^2} a(\bullet) + \frac{4}{\sigma^4} a(\bullet) \right) \bullet. $$

For order 3, the expression of $\tilde{\gamma}_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 $\bullet\bullet$ for $\gamma_1$ and the three trees $\bullet\bullet\bullet$, $\bullet\bullet$ and $\bullet\bullet\bullet$ for $\gamma_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.
# 5 Construction of high order integrators

## 5.1 Improvement of a method order via a modified equation

In [3], 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 \ge 1$ for the invariant measure, then, using Theorem 2.6, for all $j < p$, $A_j^* \rho_\infty = 0$ and $A_p^* \rho_\infty \ne 0$. By integrating by parts, we can write $\int_{\mathbb{R}^d} A_p \phi \rho_\infty dx$ as $\int_{\mathbb{R}^d} \phi' f_p \rho_\infty 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 $\theta$-method (1.11), we have $A_1 = F(\gamma)$ where $\gamma$ is given by (4.2). Applying integration by parts as described in Section 4.2, we obtain

$$ \gamma \sim \left( \theta - \frac{1}{2} \right) \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} - \frac{1}{2} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} - \frac{1}{2} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \frac{\sigma^2}{2} (1 - \theta) \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \sigma^2 (1 - \theta) \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \frac{\sigma^2}{2} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} \\ + \frac{\sigma^2}{2} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} - \frac{\sigma^4}{8} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} - \frac{\sigma^4}{8} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} - \frac{\sigma^4}{4} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} - \frac{\sigma^4}{8} \begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix}. $$

Then $f_1$ is given by

$$ (f_1)_i = \left(\theta - \frac{1}{2}\right) f'_i f - \frac{1}{2} g' f f_i - \frac{1}{2} \operatorname{div}(f) f_i + \frac{\sigma^2}{2} (1 - \theta) \Delta f_i + \sigma^2 (1 - \theta) \sum_{j=1}^{d} \partial_j f_i \partial_j g \\ + \frac{\sigma^2}{2} \Delta g f_i + \frac{\sigma^2}{2} \sum_{j=1}^{d} (\partial_j g)^2 f_i - \frac{\sigma^4}{8} \partial_i \Delta g - \frac{\sigma^4}{8} \Delta g \partial_i g - \frac{\sigma^4}{4} \sum_{j=1}^{d} \partial_j g \partial_{i,j} g \\ - \frac{\sigma^4}{8} \sum_{j=1}^{d} (\partial_j g)^2 \partial_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 = \left(\frac{1}{2} - \theta\right) \left(f'f + \frac{\sigma^2}{2}\Delta f\right). $$

**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 - h f_1$. We call it $A_2^{(1)}$. We find $A_2^{(1)} = F(\gamma)$ where

$$ \gamma = \theta(3\theta - 1)\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta(3\theta - 1)\frac{\sigma^2}{2}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta(4\theta - 1)\frac{\sigma^2}{2}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta^2\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta^2\frac{\sigma^4}{4}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} \\ + \theta^2\sigma^2\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta(4\theta - 1)\frac{\sigma^2}{2}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta(2\theta + 1)\frac{\sigma^2}{2}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \theta(2\theta + 1)\frac{\sigma^4}{4}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} + \frac{4\theta - 1}{2}\begin{bmatrix} \bullet \\ \bullet \\ \bullet \\ \bullet \end{bmatrix} $$
$$ 
\begin{aligned}
&+ (4\theta - 1) \frac{\sigma^2}{4} \bullet \bullet \bullet + \theta^2 \frac{\sigma^2}{2} \bullet \bullet \bullet + (4\theta - 1) \frac{\sigma^2}{4} \bullet \bullet \bullet + (4\theta - 1) \frac{\sigma^4}{8} \bullet \bullet \bullet + \theta \frac{\sigma^4}{2} \bullet \bullet \\
&+ \theta \sigma^2 \bullet \bullet \bullet + \frac{1}{6} \bullet \bullet \bullet + \theta \frac{\sigma^4}{2} \bullet \bullet \bullet + \frac{\sigma^2}{4} \bullet \bullet \bullet + \frac{\sigma^4}{8} \bullet \bullet \bullet + \frac{\sigma^6}{48} \bullet \bullet \bullet.
\end{aligned}
$$

Using Theorem 4.9, we find

$$ 
\begin{aligned}
\gamma &\sim \left(-2\theta^2 + 2\theta - \frac{1}{2}\right) \bullet \bullet \bullet + \left(-\theta^2 + \theta - \frac{1}{4}\right) \sigma^2 \bullet \bullet \bullet + \left(-\frac{3\theta^2}{2} + \frac{3\theta}{2} - \frac{1}{3}\right) \sigma^2 \bullet \bullet \bullet \\
&+ \left(-\theta^2 + \theta - \frac{1}{6}\right) \bullet \bullet \bullet + \left(-\frac{\theta^2}{4} + \frac{\theta}{4} - \frac{1}{24}\right) \sigma^4 \bullet \bullet \bullet + \left(-\theta^2 + \theta - \frac{1}{6}\right) \sigma^2 \bullet \bullet \bullet.
\end{aligned}
$$

Thus we define

$$ 
\begin{aligned}
f_2 &= \left(-2\theta^2 + 2\theta - \frac{1}{2}\right) f'f'f + \left(-\theta^2 + \theta - \frac{1}{4}\right) \sigma^2 f'\Delta f \\
&\quad + \left(-\frac{3\theta^2}{2} + \frac{3\theta}{2} - \frac{1}{3}\right) \sigma^2 \sum_i f''(e_i, f'(e_i)) + \left(-\theta^2 + \theta - \frac{1}{6}\right) f''(f, f) \\
&\quad + \left(-\frac{\theta^2}{4} + \frac{\theta}{4} - \frac{1}{24}\right) \sigma^4 \Delta^2 f + \left(-\theta^2 + \theta - \frac{1}{6}\right) \sigma^2 (\Delta f)'(f),
\end{aligned}
$$

and, if the $\theta$-scheme applied to $dX = (f - hf_1 - h^2 f_2)dt + \sigma dW$ is ergodic, then it has order 3 for the invariant measure.

For $f_1$, we recover the formula of [3, Prop. 5.1]. The computation of $f_2$ was first done for $\theta = 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$.

## 5.2 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 = \sum_{j=1}^s 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 = (a_{ij})$, $b = (b_i)$ and $d = (d_i)$ to achieve a given order condition for the invariant measure.

First, we suppose $\sum b_i = 1$ in order for $A_0 = \mathcal{L}$ in Assumption 2.4 to be satisfied. Then $A_1\phi = F(\gamma_1)(\phi)$ where

$$ \gamma_1 = \sum b_i c_i \bullet \bullet \bullet + \frac{\sigma^2}{2} \sum b_i d_i^2 \bullet \bullet \bullet + \sigma^2 \sum b_i d_i \bullet \bullet \bullet + \frac{1}{2} \bullet \bullet \bullet + \frac{\sigma^2}{2} \bullet \bullet \bullet + \frac{\sigma^4}{8} \bullet \bullet \bullet. $$

Theorem 4.9 yields

$$ \gamma_1 \sim \left( \sum b_i c_i + \frac{1}{2} - 2 \sum b_i d_i \right) \bullet \bullet \bullet + \frac{\sigma^2}{2} \left( \sum b_i d_i^2 + \frac{1}{2} - 2 \sum b_i d_i \right) \bullet \bullet \bullet. $$
Thus if we suppose

$$ \sum b_i = 1, \quad \sum b_i c_i + \frac{1}{2} - 2 \sum b_i d_i = 0, \quad \sum b_i d_i^2 + \frac{1}{2} - 2 \sum b_i d_i = 0, $$

then $A_1\phi \sim 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 $\sum b_i = 1$. Using the same notation as in Theorem 4.9, if A, b and d are chosen such that $F(\tilde{\gamma}_i) = 0$ for all $1 \le 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.

<table><caption>Table 1: Runge-Kutta order conditions for the invariant measure (See Theorem 5.2). The sums are over all involved indices.</caption><thead><tr><th>Order</th><th>Tree τ</th><th>F(τ)(φ)</th><th>Order condition</th></tr></thead><tbody><tr><td>1</td><td>...</td><td>φ'f</td><td>∑bᵢ = 1</td></tr><tr><td rowspan="2">2</td><td>...</td><td>φ'f'f</td><td>∑bᵢcᵢ - 2∑bᵢdᵢ = -½</td></tr><tr><td>...</td><td>φ'Δf</td><td>∑bᵢdᵢ² - 2∑bᵢdᵢ = -½</td></tr><tr><td rowspan="6">3</td><td>...</td><td>φ'f'f'f</td><td>∑bᵢaᵢⱼcⱼ - 2∑bᵢaᵢⱼdⱼ + ∑bᵢcᵢ - (∑bᵢdᵢ)² = 0</td></tr><tr><td>...</td><td>φ'f'Δf</td><td>∑bᵢaᵢⱼdⱼ² - 2∑bᵢaᵢⱼdⱼ + ∑bᵢcᵢ - (∑bᵢdᵢ)² = 0</td></tr><tr><td>...</td><td>φ'f''(f, f)</td><td>½∑bᵢcᵢ² - 2∑bᵢdᵢcᵢ<br/>-2∑bᵢdᵢ + 2∑bᵢdᵢ² + ∑bᵢcᵢ = -½</td></tr><tr><td>...</td><td>∑φ'f''(f'(eᵢ), eᵢ)</td><td>∑bᵢdᵢaᵢⱼdⱼ - ∑bᵢcᵢdᵢ - ∑bᵢdᵢ + ∑bᵢdᵢ²<br/>+ ∑bᵢcᵢ - ∑bᵢaᵢⱼdⱼ - ½(∑bᵢdᵢ)² = -½</td></tr><tr><td>...</td><td>φ'((Δf)'(f))</td><td>½∑bᵢcᵢdᵢ² - ∑bᵢdᵢ³ - 2∑bᵢdᵢ<br/>+ ½∑bᵢdᵢ² - ∑bᵢcᵢdᵢ + ½∑bᵢcᵢ = -½</td></tr><tr><td>...</td><td>φ'Δ²f</td><td>½∑bᵢdᵢ⁴ - ½∑bᵢdᵢ³<br/>- ½∑bᵢdᵢ + ¾∑bᵢdᵢ² = -¾</td></tr></tbody></table>

**Remark 5.3.** To check if a scheme has weak order p, one can develop $L^i\phi$, $i \le 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 $\frac{1}{i!}L^i\phi$ and $A_{i-1}\phi$ are equal, yielding the order p estimate (1.7). In Table 2, we collect the corresponding order conditions up to order $p \le 3$. We
recover exactly the same order conditions as first derived in [18, 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^{(2)} = 0$ in Table 2, we obtain the incompatible order conditions $\sum b_i (d_i^{(1)})^2 = \frac{1}{2}$ (third line of Table 2) and $\sum b_i (d_i^{(1)})^2 = \frac{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(\cdot)$ and the corresponding order conditions, using Proposition 4.7. Thus, under Assumption 2.2, we can replace the two conditions $\sum b_i a_{ij} d_j = \frac{1}{6}$ and $(\sum b_i \tilde{d}_i)^2 = \frac{1}{12}$ by one new condition $\sum b_i a_{ij} d_j^{(1)} + \frac{1}{2} (\sum b_i d_i^{(2)})^2 = \frac{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.

## 5.3 Order conditions for postprocessed integrators

In this section, we extend our analysis to the case of integrators combined with postprocessors [47]. 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 $[\mathcal{L}, \bar{A}_p]\phi$. This is done by the following theorem for the composition of exotic aromatic forests and based on the Leibniz rule.

**Theorem 5.6.** Let $\gamma_1$ and $\gamma_2$ be two exotic aromatic rooted forests, and let $\phi : \mathbb{R}^d \to \mathbb{R}$ be a smooth function. For $\varphi : \pi(r_2) \to V_1$ and $\psi : \Gamma(r_2) \to V_1$, we build $\gamma_{\varphi,\psi}$ by plugging all the edges connected to $r_2$ to the nodes of $\gamma_1$ according to $\varphi$, and all the lianas (counting multiplicity) according to $\psi$. Then the composition of forests is given by

$$F(\gamma_2)(F(\gamma_1)(\phi)) = \sum_{\substack{\varphi:\pi(r_2) \to V_1 \\ \psi:\Gamma(r_2) \to V_1}} F(\gamma_{\varphi,\psi})(\phi).$$

Various composition rules for B-series and aromatic B-series have been studied in the literature (see [23, 16, 8] 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(\gamma_1)(\phi) = \sum_{v_1^{(1)}, \dots, v_m^{(1)}} \sum_{j_1^{(1)}, \dots, j_s^{(1)}} \left( \prod_{v \in V_1^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{iv} \right) \partial_{I_{\pi(r_1)}} \partial_{J_{\Gamma(r_1)}} \phi.$$
<table><caption>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 [18].</caption><thead><tr><th>Order</th><th>Tree τ</th><th>F(τ)(φ)</th><th>Order condition</th></tr></thead><tbody><tr><td>1</td><td>...</td><td>φ'f</td><td>Σb<sub>i</sub> = 1</td></tr><tr><td rowspan="3">2</td><td>...</td><td>φ'f'f</td><td>Σb<sub>i</sub>c<sub>i</sub> = 1/2</td></tr><tr><td>...</td><td>φ'Δf</td><td>Σb<sub>i</sub>(d<sub>i</sub><sup>(1)</sup>)<sup>2</sup> + Σb<sub>i</sub>(d<sub>i</sub><sup>(2)</sup>)<sup>2</sup> = 1/2</td></tr><tr><td>...</td><td>Σφ''(e<sub>i</sub>, f'(e<sub>i</sub>))</td><td>Σb<sub>i</sub>d<sub>i</sub><sup>(1)</sup> = 1/2</td></tr><tr><td rowspan="10">3</td><td>...</td><td>φ'f'f'f</td><td>Σb<sub>i</sub>a<sub>ij</sub>c<sub>j</sub> = 1/6</td></tr><tr><td>...</td><td>φ'f'Δf</td><td>Σb<sub>i</sub>a<sub>ij</sub>(d<sub>j</sub><sup>(1)</sup>)<sup>2</sup> + Σb<sub>i</sub>a<sub>ij</sub>(d<sub>j</sub><sup>(2)</sup>)<sup>2</sup> = 1/6</td></tr><tr><td>...</td><td>φ'f''(f, f)</td><td>Σb<sub>i</sub>c<sub>i</sub><sup>2</sup> = 1/3</td></tr><tr><td>...</td><td>Σφ'f''(f'(e<sub>i</sub>), e<sub>i</sub>)</td><td>Σb<sub>i</sub>d<sub>i</sub><sup>(1)</sup>a<sub>ij</sub>d<sub>j</sub><sup>(1)</sup> + Σb<sub>i</sub>d<sub>i</sub><sup>(2)</sup>a<sub>ij</sub>d<sub>j</sub><sup>(2)</sup> = 1/6</td></tr><tr><td>...</td><td>φ'((Δf)'(f))</td><td>Σb<sub>i</sub>c<sub>i</sub>(d<sub>i</sub><sup>(1)</sup>)<sup>2</sup> + Σb<sub>i</sub>c<sub>i</sub>(d<sub>i</sub><sup>(2)</sup>)<sup>2</sup> = 1/3</td></tr><tr><td>...</td><td>φ'Δ<sup>2</sup>f</td><td>Σb<sub>i</sub>(d<sub>i</sub><sup>(1)</sup>)<sup>2</sup> + (d<sub>i</sub><sup>(2)</sup>)<sup>2</sup> = 1/3</td></tr><tr><td>...</td><td>Σφ''(e<sub>i</sub>, f'f'(e<sub>i</sub>))</td><td>Σb<sub>i</sub>a<sub>ij</sub>d<sub>j</sub><sup>(1)</sup> = 1/6</td></tr><tr><td>...</td><td>Σφ''(f'(e<sub>i</sub>), f'(e<sub>i</sub>))</td><td>(Σb<sub>i</sub>d<sub>i</sub><sup>(2)</sup>)<sup>2</sup> = 1/12</td></tr><tr><td>...</td><td>Σφ''(e<sub>i</sub>, f''(e<sub>i</sub>, f))</td><td>Σb<sub>i</sub>c<sub>i</sub>d<sub>i</sub><sup>(1)</sup> = 1/3</td></tr><tr><td>...</td><td>Σφ''(e<sub>i</sub>, (Δf)'(e<sub>i</sub>))</td><td>Σb<sub>i</sub>(d<sub>i</sub><sup>(1)</sup>)<sup>3</sup> + Σb<sub>i</sub>d<sub>i</sub><sup>(1)</sup>(d<sub>i</sub><sup>(2)</sup>)<sup>2</sup> = 1/3</td></tr><tr><td>...</td><td>Σφ'''(e<sub>i</sub>, e<sub>j</sub>, f''(e<sub>i</sub>, e<sub>j</sub>))</td><td>Σb<sub>i</sub>(d<sub>i</sub><sup>(1)</sup>)<sup>2</sup> = 1/3</td></tr></tbody></table>
Then we replace $\phi$ by $F(\gamma_1)(\phi)$ and use the Leibniz rule to distribute the partial derivatives.

$$
\begin{align*}
F(\gamma_2)(F(\gamma_1)(\phi)) &= \sum_{\substack{i_1^{(1)}, \dots, i_{m_1}^{(1)} \\ i_1^{(2)}, \dots, i_{m_2}^{(2)}}} \sum_{\substack{j_1^{(1)}, \dots, j_{s_1}^{(1)} \\ j_1^{(2)}, \dots, j_{s_2}^{(2)}}} \left( \prod_{v \in V_2^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{iv} \right) \\
&\quad \cdot \partial_{I_{\pi(r_2)}} \partial_{J_{\Gamma(r_2)}} \left[ \left( \prod_{v \in V_1^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{iv} \right) \partial_{I_{\pi(r_1)}} \partial_{J_{\Gamma(r_1)}} \phi \right] \\
&= \sum_{\substack{\varphi:\pi(r_2) \to V_1 \\ \psi:\Gamma(r_2) \to V_1}} \sum_{\substack{i_1^{(1)}, \dots, i_{m_1}^{(1)} \\ i_1^{(2)}, \dots, i_{m_2}^{(2)}}} \sum_{\substack{j_1^{(1)}, \dots, j_{s_1}^{(1)} \\ j_1^{(2)}, \dots, j_{s_2}^{(2)}}} \left( \prod_{v \in V_2^0} \partial_{I_{\pi(v)}} \partial_{J_{\Gamma(v)}} f_{iv} \right) \\
&\quad \cdot \left( \prod_{v \in V_1^0} \partial_{I_{\pi(v)} \cup \varphi^{-1}(\{v\})} \partial_{J_{\Gamma(v)} \cup \psi^{-1}(\{v\})} f_{iv} \right) \partial_{I_{\pi(r_1)} \cup \varphi^{-1}(\{r_1\})} \partial_{J_{\Gamma(r_1)} \cup \psi^{-1}(\{r_1\})} \phi \\
&= \sum_{\substack{\varphi:\pi(r_2) \to V_1 \\ \psi:\Gamma(r_2) \to V_1}} F(\gamma_{\varphi, \psi})(\phi). \quad \square
\end{align*}
$$

**Example.** We recall $\mathcal{L} = F(\bullet + \frac{\sigma^2}{2}\bullet)$, then we can compute $\mathcal{L}^2\phi$. Using Theorem 5.6, we obtain

$$F(\bullet)(\mathcal{L}\phi) = F(\bullet + \bullet + \frac{\sigma^2}{2}\bullet)(\phi)$$

and

$$F(\bullet)(\mathcal{L}\phi) = F(\bullet + \bullet + 2\bullet + \frac{\sigma^2}{2}\bullet)(\phi).$$

Combining the two previous equalities, we deduce

$$\mathcal{L}^2\phi = F(\bullet + \bullet + \sigma^2\bullet + \frac{\sigma^2}{2}\bullet + \sigma^2\bullet + \frac{\sigma^4}{4}\bullet)(\phi).$$

Using Theorems 2.8 and 5.6, we obtain general conditions on postprocessors to increase by 1 the order of a given method.

**Theorem 5.7.** 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 and

$$\mathbb{E}[\phi(G_n(x))] = F(\bullet)(\phi) + \sum_{\substack{\gamma \in \mathcal{A} \\ 1 \le |\gamma| \le p}} h^{|\gamma|} \bar{a}(\gamma) F(\gamma)(\phi) + \dots,$$

if we denote $\gamma$ the exotic aromatic B-series such that $F(\gamma) = (\mathcal{A}_p + [\mathcal{L}, \overline{\mathcal{A}_p}])$ and if $\gamma \sim 0$, then $\overline{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 $\bar{a}(\bullet) + \frac{2}{\sigma^2}\bar{a}(\bullet) = 0$ (in order to have $\overline{\mathcal{A}}_0 = \alpha\mathcal{L}$), and if conditions in Table 4 are verified, then the method is of order 3 for the invariant measure.
<table><caption>Table 3: General order 2 conditions with postprocessor (See Theorem 5.7).</caption><thead><tr><th>Tree γ</th><th>Order condition</th></tr></thead><tbody><tr><td>⋮</td><td>a(•) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(•) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(•) - <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(•) - 2a(•) + <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(•) = 0</td></tr><tr><td>⋮</td><td>a(•) - a(•) + a(•) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(•) - <sup>σ<sup>2</sup></sup>&frasl;<sub>2</sub>a(•) + a(•) = 0</td></tr><tr><td>⋮</td><td>a(•) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(•) + <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(•) = 0</td></tr></tbody></table>

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

$$ 
\begin{aligned}
\bar{Y}_i &= X_n + h \sum_{j=1}^{s} \bar{a}_{ij} f(\bar{Y}_j) + \bar{d}_i \sigma \sqrt{h} \bar{\xi}_n, \quad i = 1, \dots, s, \\
\bar{X}_n &= X_n + h \sum_{i=1}^{s} \bar{b}_i f(\bar{Y}_i) + \bar{d}_0 \sigma \sqrt{h} \bar{\xi}_n.
\end{aligned}
$$

Assume Assumption 2.2. If $\gamma$ is the exotic aromatic B-series such that $F(\gamma) = (A_p + [\mathcal{L}, \overline{A_p}])$ and if $A, b, d, \overline{A}, \overline{b}, \overline{d}, \overline{d}_0$ are chosen such that $\gamma \sim 0$ then the postprocessed method $\overline{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 $\sum \bar{b}_i = \bar{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 $\cdot \cdot \cdot$ in Table 3 is not modified by the postprocessor as it does not depend on any $\bar{a}(\gamma)$. 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

$$ \bar{a}(\cdot \cdot \cdot) - \frac{2}{\sigma^2} \bar{a}(\cdot) + \frac{4}{\sigma^4} \bar{a}(\cdot \cdot \cdot) = 0, \qquad (5.1) $$

then the conditions of $\cdot \cdot \cdot, \cdot \cdot \cdot$ and $\cdot \cdot \cdot$ 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 $\overline{A}_0 = \alpha \mathcal{L}$.

**Example.** Under Assumption 2.2, the following Runge-Kutta method, introduced in [47], is of order 2 for the invariant measure of (1.4) (if it is ergodic).

$$ 
\begin{aligned}
X_{n+1} &= X_n + h f(X_{n+1} + \frac{-1+\sqrt{2}}{2}\sigma\sqrt{h}\xi_n) + \sigma\sqrt{h}\xi_n, \\
\overline{X}_n &= X_n + h \frac{\sqrt{2}}{2} f(\overline{X}_n) + \frac{\sqrt{4\sqrt{2}-1}}{2}\sigma\sqrt{h}\overline{\xi}_n.
\end{aligned}
$$

Indeed, its coefficients, placed in the following Butcher tableau, fulfil the conditions of order 2 of Theorem 5.8 (See Table 5).

$$ \frac{c}{b} \quad | \quad A \quad | \quad d \quad | \quad \bar{c} \quad | \quad \bar{A} \quad | \quad \bar{d} \quad = \frac{1}{1} \quad | \quad 1 \quad | \quad \frac{1+\sqrt{2}}{2} \quad | \quad \frac{\sqrt{2}}{2} \quad | \quad \frac{\sqrt{2}}{2} \quad | \quad \frac{\sqrt{4\sqrt{2}-1}}{2} \quad | \quad \frac{\sqrt{2}}{2} \quad | \quad \frac{\sqrt{4\sqrt{2}-1}}{2} $$
<table><caption>Table 4: General order 3 conditions with postprocessor (See Theorem 5.7).</caption><thead><tr><th>Tree τ</th><th>Order condition</th></tr></thead><tbody><tr><td>...</td><td>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) + <sup>12</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - <sup>24</sup>&frasl;<sub>σ<sup>6</sup></sub>a(●) - 2a(●) + <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + 2a(●) - <sup>8</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>16</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - a(●) - a(●) + a(●) + a(●) - 2a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>6</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>12</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - <sup>σ<sup>2</sup></sup>&frasl;<sub>2</sub>a(●) - a(●) + 2a(●) + σ<sup>2</sup>a(●) - 4a(●) + <sup>8</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) + <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - <sup>8</sup>&frasl;<sub>σ<sup>6</sup></sub>a(●) - a(●) - a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - a(●) - a(●) - a(●) + 2a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + a(●) - 2a(●) - <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>8</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - <sup>16</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - σ<sup>2</sup>a(●) + σ<sup>2</sup>a(●) + 2a(●) - 4a(●) + <sup>8</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + a(●) + <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - <sup>8</sup>&frasl;<sub>σ<sup>6</sup></sub>a(●) - <sup>σ<sup>2</sup></sup>&frasl;<sub>2</sub>a(●) - a(●) + 2a(●) - 2a(●) + <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - a(●) + a(●) + a(●) - a(●) + a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>σ<sup>2</sup></sup>&frasl;<sub>2</sub>a(●) + <sup>σ<sup>2</sup></sup>&frasl;<sub>2</sub>a(●) - <sup>σ<sup>2</sup></sup>&frasl;<sub>2</sub>a(●) + a(●) = 0</td></tr><tr><td>...</td><td>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - <sup>8</sup>&frasl;<sub>σ<sup>6</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - a(●) + 2a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>6</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>12</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) - σ<sup>2</sup>a(●) + 2a(●) - <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) = 0</td></tr><tr><td>...</td><td>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>2</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) + <sup>4</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) + <sup>4</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>12</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) + <sup>24</sup>&frasl;<sub>σ<sup>6</sup></sub>a(●) - 4a(●) + <sup>8</sup>&frasl;<sub>σ<sup>2</sup></sub>a(●) - <sup>16</sup>&frasl;<sub>σ<sup>4</sup></sub>a(●) = 0</td></tr></tbody></table>
<table><caption>Table 5: Order conditions for Runge-Kutta method with Runge-Kutta postprocessor (See Theorem 5.8). The sums are over all involved indices.</caption><thead><tr><th>Order</th><th>Tree τ</th><th>Order condition</th></tr></thead><tbody><tr><td rowspan="2">2</td><td>...</td><td>Σ b<sub>i</sub>c<sub>i</sub> - 2Σ b<sub>i</sub>d<sub>i</sub> - 2Σ b&#773;<sub>i</sub> + 2d&#773;<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>...</td><td>Σ b<sub>i</sub>d<sub>i</sub><sup>2</sup> - 2Σ b<sub>i</sub>d<sub>i</sub> - Σ b&#773;<sub>i</sub> + d&#773;<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td rowspan="6">3</td><td>...</td><td>Σ b<sub>i</sub>a<sub>ij</sub>c<sub>j</sub> - 2Σ b<sub>i</sub>a<sub>ij</sub>d<sub>j</sub> + Σ b<sub>i</sub>c<sub>i</sub> - (Σ b<sub>i</sub>d<sub>i</sub>)<sup>2</sup> - 2Σ b&#773;<sub>i</sub>c&#773;<sub>i</sub> + 4d&#773;<sub>0</sub>Σ b&#773;<sub>i</sub>d<sub>i</sub> - d&#773;<sub>0</sub><sup>4</sup> = 0</td></tr><tr><td>...</td><td>Σ b<sub>i</sub>a<sub>ij</sub>d<sub>j</sub><sup>2</sup> - 2Σ b<sub>i</sub>a<sub>ij</sub>d<sub>j</sub> + Σ b<sub>i</sub>c<sub>i</sub> - (Σ b<sub>i</sub>d<sub>i</sub>)<sup>2</sup> - Σ b&#773;<sub>i</sub>c&#773;<sub>i</sub> - Σ b&#773;<sub>i</sub>d<sub>i</sub><sup>2</sup> + 4d&#773;<sub>0</sub>Σ b&#773;<sub>i</sub>d<sub>i</sub> - d&#773;<sub>0</sub><sup>4</sup> = 0</td></tr><tr><td>...</td><td>1/2 Σ b<sub>i</sub>c<sub>i</sub><sup>2</sup> - 2Σ b<sub>i</sub>d<sub>i</sub>c<sub>i</sub> - 2Σ b<sub>i</sub>d<sub>i</sub> + 2Σ b<sub>i</sub>d<sub>i</sub><sup>2</sup> + Σ b<sub>i</sub>c<sub>i</sub> - Σ b<sub>i</sub>a<sub>ij</sub>d<sub>j</sub> - 1/2 (Σ b<sub>i</sub>d<sub>i</sub>)<sup>2</sup> - Σ b&#773;<sub>i</sub>c&#773;<sub>i</sub> + 2d&#773;<sub>0</sub>Σ b&#773;<sub>i</sub>d<sub>i</sub> - d&#773;<sub>0</sub><sup>4</sup> = -1/6</td></tr><tr><td>...</td><td>1/2 Σ b<sub>i</sub>c<sub>i</sub>d<sub>i</sub><sup>2</sup> - Σ b<sub>i</sub>d<sub>i</sub><sup>3</sup> - 2Σ b<sub>i</sub>d<sub>i</sub> + 5/2 Σ b<sub>i</sub>d<sub>i</sub><sup>2</sup> - Σ b<sub>i</sub>c<sub>i</sub>d<sub>i</sub> + 1/2 Σ b<sub>i</sub>c<sub>i</sub> - 1/2 Σ b&#773;<sub>i</sub>c&#773;<sub>i</sub> - 1/2 Σ b&#773;<sub>i</sub>d<sub>i</sub><sup>2</sup> + 2d&#773;<sub>0</sub>Σ b&#773;<sub>i</sub>d<sub>i</sub> - d&#773;<sub>0</sub><sup>4</sup> = -1/3</td></tr><tr><td>...</td><td>1/8 Σ b<sub>i</sub>d<sub>i</sub><sup>4</sup> - 1/2 Σ b<sub>i</sub>d<sub>i</sub><sup>3</sup> - 1/2 Σ b<sub>i</sub>d<sub>i</sub> + 3/4 Σ b<sub>i</sub>d<sub>i</sub><sup>2</sup> - 1/4 Σ b&#773;<sub>i</sub>d<sub>i</sub><sup>2</sup> + 1/2d&#773;<sub>0</sub>Σ b&#773;<sub>i</sub>d<sub>i</sub> - d&#773;<sub>0</sub><sup>4</sup> = -1/12</td></tr></tbody></table>

## 5.4 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 $\phi$. 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 = \nabla V_1$ and $f_2 = \nabla V_2$ are both gradients, then

$$A \Box \cdot B \sim \frac{2}{\sigma^2} \left( \begin{pmatrix} \bullet^A \\ \bullet^B \end{pmatrix} + \begin{pmatrix} \circ^A \\ \circ^B \end{pmatrix} \right).$$

Furthermore, the node $B$ can be replaced by an aromatic root or a white node.
• The operator $\mathcal{L}$ is now written as

$$ \mathcal{L} = F(\cdot + \cdot + \frac{\sigma^2}{2} \cdot). $$

In addition to the partitioning of the method, 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

$$ \begin{aligned} Y_i &= X_n + h \sum_{j=1}^{s} a_{ij} f_1(Y_j) + \widehat{a}_{ij} f_2(Y_j) + d_i \sigma \sqrt{h} \xi_n, \quad i = 1, \dots, s, \\ X_{n+1} &= X_n + h \sum_{i=1}^{s} b_i f_1(Y_i) + \widehat{b}_i f_2(Y_i) + \sigma \sqrt{h} \xi_n, \end{aligned} $$

together with the following Runge-Kutta postprocessor

$$ \begin{aligned} \bar{Y}_i &= X_n + h \sum_{j=1}^{s} \bar{a}_{ij} f_1(\bar{Y}_j) + \widehat{\bar{a}}_{ij} f_2(\bar{Y}_j) + \bar{d}_i \sigma \sqrt{h} \bar{\xi}_n, \quad i = 1, \dots, s, \\ \bar{X}_n &= X_n + h \sum_{i=1}^{s} \bar{b}_i f_1(\bar{Y}_i) + \widehat{\bar{b}}_i f_2(\bar{Y}_i) + \bar{d}_0 \sigma \sqrt{h} \bar{\xi}_n. \end{aligned} $$

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 $\gamma \sim 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.

<table><caption>Table 6: Order conditions for partitioned Runge-Kutta method with postprocessor (See Theorem 5.10). The sums are over all involved indices.</caption><thead><tr><th>Order</th><th>Tree τ</th><th>F(τ)(φ)</th><th>Order condition</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>●</td><td>φ'f<sub>1</sub></td><td>Σb<sub>i</sub> = 1</td></tr><tr><td>●</td><td>φ'f<sub>2</sub></td><td>Σb̂<sub>i</sub> = 1</td></tr><tr><td rowspan="6">2</td><td>●</td><td>φ'f'<sub>1</sub>f<sub>1</sub></td><td>Σb<sub>i</sub>c<sub>i</sub> - 2Σb<sub>i</sub>d<sub>i</sub> - 2Σb̂<sub>i</sub> + 2d̅<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>●</td><td>φ'f'<sub>1</sub>f<sub>2</sub></td><td>Σb<sub>i</sub>ĉ<sub>i</sub> - 2Σb<sub>i</sub>d<sub>i</sub> - Σb̂<sub>i</sub> - Σb̂<sub>i</sub> + 2d̅<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>●</td><td>φ'f'<sub>2</sub>f<sub>1</sub></td><td>Σb̂<sub>i</sub>c<sub>i</sub> - 2Σb̂<sub>i</sub>d<sub>i</sub> - Σb̂<sub>i</sub> - Σb̂<sub>i</sub> + 2d̅<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>●</td><td>φ'f'<sub>2</sub>f<sub>2</sub></td><td>Σb̂<sub>i</sub>ĉ<sub>i</sub> - 2Σb̂<sub>i</sub>d<sub>i</sub> - 2Σb̂<sub>i</sub> + 2d̅<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>●</td><td>φ'Δf<sub>1</sub></td><td>Σb<sub>i</sub>d<sub>i</sub><sup>2</sup> - 2Σb<sub>i</sub>d<sub>i</sub> - Σb̂<sub>i</sub> + d̅<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>●</td><td>φ'Δf<sub>2</sub></td><td>Σb̂<sub>i</sub>d<sub>i</sub><sup>2</sup> - 2Σb̂<sub>i</sub>d<sub>i</sub> - Σb̂<sub>i</sub> + d̅<sub>0</sub><sup>2</sup> = -1/2</td></tr></tbody></table>
**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).

$$ 
\begin{aligned}
X_{n+1} &= X_n + \frac{h}{2} f_1(X_{n+1} + \frac{1}{2}\sigma\sqrt{h\xi_n}) + \frac{h}{2} f_1(X_{n+1} + \frac{3}{2}\sigma\sqrt{h\xi_n}) \\
&\quad + h f_2(X_n + \frac{1}{2}\sigma\sqrt{h\xi_n}) + \sigma\sqrt{h\xi_n}, \\
\bar{X}_n &= X_n + \frac{1}{2}\sigma\sqrt{h\xi_n}.
\end{aligned}
$$

It can be put in Runge-Kutta form with the coefficients below:

$$ 
\begin{array}{c|c|c|c|c}
c & A & \hat{c} & \hat{A} & d \\
\hline
b & & \hat{b} & &
\end{array}
=
\begin{array}{c|c|c|c|c}
0 & 0 & 0 & 0 & 0 \\
1 & 0 & 1/2 & 1/2 & 1 \\
1 & 0 & 1/2 & 1/2 & 1 \\
\hline
0 & 1/2 & 1/2 & & 1
\end{array}
\begin{array}{c|c|c|c|c}
0 & 0 & 0 & 0 & 1/2 \\
1 & 0 & 0 & 0 & 1/2 \\
1 & 0 & 0 & 0 & 3/2 \\
\hline
0 & 1/2 & 1/2 & & 1
\end{array}
$$

and $s = 0$ and $\bar{d}_0 = \frac{1}{2}$ for the postprocessor.

If we add a family of independent noises $(\chi_n)_n$ independent of $(\xi_n)_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):

$$ 
\begin{aligned}
X_{n+1} &= X_n + h f_1(X_{n+1} + \frac{1}{2}\sigma\sqrt{h\chi_n}) + h f_2(X_n + \frac{1}{2}\sigma\sqrt{h\xi_n}) + \sigma\sqrt{h\xi_n}, \\
\bar{X}_n &= X_n + \frac{1}{2}\sigma\sqrt{h\xi_n}.
\end{aligned}
$$

## 5.5 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 [32], and it can also reduce the variance [22]. 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 = -\nabla V$ is a gradient, and $f_2$ is a perturbation of $f_1$ that satisfies

$$ \operatorname{div} \left( f_2 e^{-\frac{2}{\sigma^2} V} \right) = 0. \qquad (5.2) $$

The perturbation $f_2$ does not modify the invariant measure. Indeed equation (5.2) implies that the adjoint of $\mathcal{B}\phi = \phi'(f_2)$ satisfies $\mathcal{B}^*\rho_\infty = 0$, and thus the invariant measure is preserved. A simple example of such non gradient perturbation is $f_2 = J\nabla 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 \hat{\bullet} B \sim \frac{2}{\sigma^2} \bullet B^A$. Furthermore, the node $B$ can be replaced by an aromatic root or a white node.

* The generator reads $\mathcal{L} = F(\bullet + \bullet + \frac{\sigma^2}{2}\bullet)$.

* We have $F(\bullet) = -F(\bullet)$, and these differentials vanish if $f_2 = J\nabla 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 $\bullet \circ \circ$. This is where the last property comes in handy, as we can integrate by part this tree and obtain $\circ \circ \sim \bullet \circ$. 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 = -\nabla 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 $\gamma \sim 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.

<table><caption>Table 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.</caption><thead><tr><th>Order</th><th>Tree τ</th><th>F(τ)(φ)</th><th>Order condition</th></tr></thead><tbody><tr><td>1</td><td>●</td><td>φ'f<sub>1</sub></td><td>Σb<sub>i</sub> = 1</td></tr><tr><td rowspan="5">2</td><td>●●</td><td>φ'f'<sub>1</sub>f<sub>1</sub></td><td>Σb<sub>i</sub>c<sub>i</sub> - 2Σb<sub>i</sub>d<sub>i</sub> - 2Σb̄<sub>i</sub> + 2d̄<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td>●●</td><td>φ'f'<sub>1</sub>f<sub>2</sub></td><td>Σb<sub>i</sub>ĉ<sub>i</sub> + Σb̄<sub>i</sub> - Σb̄<sub>i</sub> = 0</td></tr><tr><td>●●</td><td>φ'f'<sub>2</sub>f<sub>1</sub></td><td>Σb̂<sub>i</sub>c<sub>i</sub> - 2Σb̂<sub>i</sub>d<sub>i</sub> + Σb̂<sub>i</sub> - Σb̄<sub>i</sub> - Σb̂<sub>i</sub> + 2d̄<sub>0</sub><sup>2</sup> = 0</td></tr><tr><td>●●</td><td>φ'f'<sub>2</sub>f<sub>2</sub></td><td>Σb̂<sub>i</sub>ĉ<sub>i</sub> - 1/2 (Σb̂<sub>i</sub>)<sup>2</sup> = 0</td></tr><tr><td>●</td><td>φ'Δf<sub>1</sub></td><td>Σb<sub>i</sub>d<sub>i</sub><sup>2</sup> - 2Σb<sub>i</sub>d<sub>i</sub> - Σb̄<sub>i</sub> + d̄<sub>0</sub><sup>2</sup> = -1/2</td></tr><tr><td></td><td>●</td><td>φ'Δf<sub>2</sub></td><td>Σb̂<sub>i</sub>d<sub>i</sub><sup>2</sup> - 2Σb̂<sub>i</sub>d<sub>i</sub> + Σb̂<sub>i</sub> - Σb̄<sub>i</sub> + d̄<sub>0</sub><sup>2</sup> = 0</td></tr></tbody></table>

We note that if $f_2 = J\nabla V$, we have $\frac{\sigma^2}{2} \stackrel{\circ}{\bullet} \sim \stackrel{\circ}{\bullet} - \stackrel{\circ}{\bullet}$. In this case, the order condition for $\stackrel{\circ}{\bullet}$ can be omitted and the two conditions of $\stackrel{\circ}{\bullet}$ and $\stackrel{\circ}{\bullet}$ are respectively replaced by

$$2 \sum \hat{b}_i d_i - \sum \hat{b}_i + \sum b_i \hat{c}_i - \sum \hat{b}_i d_i^2 + \sum \bar{b}_i - \bar{d}_0^2 = 0, \\ \sum \hat{b}_i c_i - \sum \hat{b}_i d_i^2 - \sum \bar{b}_i + \bar{d}_0^2 = 0.$$

**Remark 5.12.** In order for the method to satisfy $A_0 = \mathcal{L}$, the condition $\sum \hat{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 + h f_1(X_n + \frac{\sigma}{2}\sqrt{h\xi_n}) + \frac{5}{4}h f_2(X_n + \frac{\sigma}{2}\sqrt{h\xi_n}) \\ - \frac{1}{4}h f_2(X_n - 2h f_2(X_n + \frac{\sigma}{2}\sqrt{h\xi_n}) - \frac{\sigma}{2}\sqrt{h\xi_n}) + \sigma\sqrt{h\xi_n}, \quad (5.3) \\ \bar{X}_n = X_n + \frac{\sigma}{2}\sqrt{h\xi_n}.$$

If $f_2 = J\nabla V$, it needs two evaluations of $\nabla 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 [30], formulated in a different manner (See [47]).
**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.

## References

[1] A. Abdulle, I. Almuslimani, and G. Vilmart. Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations. *SIAM/ASA J. Uncertain. Quantif.*, 6(2):937-964, 2018.

[2] A. Abdulle, D. Cohen, G. Vilmart, and K. C. Zygalakis. High weak order methods for stochastic differential equations based on modified equations. *SIAM J. Sci. Comput.*, 34(3):A1800–A1823, 2012.

[3] A. Abdulle, G. Vilmart, and K. C. Zygalakis. High order numerical approximation of the invariant measure of ergodic SDEs. *SIAM J. Numer. Anal.*, 52(4):1600–1622, 2014.

[4] A. Abdulle, G. Vilmart, and K. C. Zygalakis. Long time accuracy of Lie-Trotter splitting methods for Langevin dynamics. *SIAM J. Numer. Anal.*, 53(1):1–16, 2015.

[5] A. Alamo and J. M. Sanz-Serna. A technique for studying strong and weak local errors of splitting stochastic integrators. *SIAM J. Numer. Anal.*, 54(6):3239–3257, 2016.

[6] S. Anmarkrud and A. Kværnø. Order conditions for stochastic Runge-Kutta methods preserving quadratic invariants of Stratonovich SDEs. *J. Comput. Appl. Math.*, 316:40–46, 2017.

[7] V. Bally and D. Talay. The law of the Euler scheme for stochastic differential equations. I. Convergence rate of the distribution function. *Probab. Theory Related Fields*, 104(1):43–60, 1996.

[8] G. Bogfjellmo. Algebraic structure of aromatic B-series. *ArXiv e-prints*, May 2015.

[9] N. Bou-Rabee and H. Owhadi. Long-run accuracy of variational integrators in the stochastic context. *SIAM J. Numer. Anal.*, 48(1):278–297, 2010.

[10] C.-E. Bréhier and G. Vilmart. High order integrator for sampling the invariant distribution of a class of parabolic stochastic PDEs with additive space-time noise. *SIAM J. Sci. Comput.*, 38(4):A2283–A2306, 2016.

[11] K. Burrage and P. M. Burrage. High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations. *Appl. Numer. Math.*, 22(1-3):81–101, 1996. Special issue celebrating the centenary of Runge-Kutta methods.

[12] K. Burrage and P. M. Burrage. Order conditions of stochastic Runge-Kutta methods by B-series. *SIAM J. Numer. Anal.*, 38(5):1626–1646, 2000.

[13] J. C. Butcher. The effective order of Runge-Kutta methods. In *Conf. on Numerical Solution of Differential Equations (Dundee, 1969)*, pages 133–139. Springer, Berlin, 1969.
[14] J. C. Butcher. An algebraic theory of integration methods. *Math. Comp.*, 26:79–106, 1972.

[15] J. C. Butcher. *Numerical methods for ordinary differential equations*. John Wiley & Sons, Ltd., Chichester, third edition, 2016.

[16] P. Chartier, E. Hairer, and G. Vilmart. Algebraic structures of B-series. *Found. Comput. Math.*, 10(4):407–427, 2010.

[17] P. Chartier and A. Murua. Preserving first integrals and volume forms of additively split systems. *IMA J. Numer. Anal.*, 27(2):381–405, 2007.

[18] K. Debrabant. Runge-Kutta methods for third order weak approximation of SDEs with multidimensional additive noise. *BIT*, 50(3):541–558, 2010.

[19] K. Debrabant and A. Kværnø. B-series analysis of stochastic Runge-Kutta methods that use an iterative scheme to compute their internal stage values. *SIAM J. Numer. Anal.*, 47(1):181–203, 2008/09.

[20] K. Debrabant and A. Kværnø. Composition of stochastic B-series with applications to implicit Taylor methods. *Appl. Numer. Math.*, 61(4):501–511, 2011.

[21] A. Debussche and E. Faou. Weak backward error analysis for SDEs. *SIAM J. Numer. Anal.*, 50(3):1735–1752, 2012.

[22] A. B. Duncan, T. Lelièvre, and G. A. Pavliotis. Variance reduction using nonreversible Langevin samplers. *J. Stat. Phys.*, 163(3):457–491, 2016.

[23] E. Hairer, C. Lubich, and G. Wanner. *Geometric numerical integration*, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. Structure-preserving algorithms for ordinary differential equations.

[24] E. Hairer and G. Wanner. On the Butcher group and general multi-value methods. *Computing (Arch. Elektron. Rechnen)*, 13(1):1–15, 1974.

[25] R. Z. Hasminskii. *Stochastic stability of differential equations*, volume 7 of *Monographs and Textbooks on Mechanics of Solids and Fluids: Mechanics and Analysis*. Sijthoff & Noordhoff, Alphen aan den Rijn—Germantown, Md., 1980. Translated from the Russian by D. Louvish.

[26] L. Isserlis. On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables. *Biometrika*, 12, 11 1918.

[27] Y. Komori, T. Mitsui, and H. Sugiura. Rooted tree analysis of the order conditions of ROW-type scheme for stochastic differential equations. *BIT*, 37(1):43–66, 1997.

[28] M. Kopec. Weak backward error analysis for Langevin process. *BIT*, 55(4):1057–1103, 2015.

[29] M. Kopec. Weak backward error analysis for overdamped Langevin processes. *IMA J. Numer. Anal.*, 35(2):583–614, 2015.
[30] B. Leimkuhler and C. Matthews. Rational construction of stochastic numerical methods for molecular sampling. *Appl. Math. Res. Express. AMRX*, (1):34-56, 2013.

[31] B. Leimkuhler, C. Matthews, and G. Stoltz. The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics. *IMA J. Numer. Anal.*, 36(1):13-79, 2016.

[32] T. Lelièvre, F. Nier, and G. A. Pavliotis. Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion. *J. Stat. Phys.*, 152(2):237-274, 2013.

[33] J. C. Mattingly, A. M. Stuart, and D. J. Higham. Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise. *Stochastic Process. Appl.*, 101(2):185-232, 2002.

[34] J. C. Mattingly, A. M. Stuart, and M. V. Tretyakov. Convergence of numerical time-averaging and stationary measures via Poisson equations. *SIAM J. Numer. Anal.*, 48(2):552-577, 2010.

[35] R. I. McLachlan, K. Modin, H. Munthe-Kaas, and O. Verdier. B-series methods are exactly the affine equivariant methods. *Numer. Math.*, 133(3):599-622, 2016.

[36] G. N. Milstein. Weak approximation of solutions of systems of stochastic differential equations. *Teor. Veroyatnost. i Primenen.*, 30(4):706-721, 1985.

[37] G. N. Milstein and M. V. Tretyakov. *Stochastic numerics for mathematical physics*. Scientific Computation. Springer-Verlag, Berlin, 2004.

[38] G. N. Milstein and M. V. Tretyakov. Computing ergodic limits for Langevin equations. *Phys. D*, 229(1):81-95, 2007.

[39] H. Munthe-Kaas and O. Verdier. Aromatic Butcher series. *Found. Comput. Math.*, 16(1):183-215, 2016.

[40] A. Rößler. Stochastic Taylor expansions for the expectation of functionals of diffusion processes. *Stochastic Anal. Appl.*, 22(6):1553-1576, 2004.

[41] A. Rößler. Rooted tree analysis for order conditions of stochastic Runge-Kutta methods for the weak approximation of stochastic differential equations. *Stoch. Anal. Appl.*, 24(1):97-134, 2006.

[42] A. Rößler. Runge-Kutta methods for Itô stochastic differential equations with scalar noise. *BIT*, 46(1):97-110, 2006.

[43] A. Rößler. Strong and weak approximation methods for stochastic differential equations—some recent developments. In *Recent developments in applied probability and statistics*, pages 127-153. Physica, Heidelberg, 2010.

[44] J. M. Sanz-Serna. Markov chain Monte Carlo and numerical differential equations. In *Current challenges in stability issues for numerical differential equations*, volume 2082 of *Lecture Notes in Math.*, pages 39-88. Springer, Cham, 2014.
[45] J. M. Sanz-Serna and A. Murua. Formal series and numerical integrators: some history and some new techniques. In *Proceedings of the 8th International Congress on Industrial and Applied Mathematics*, pages 311–331. Higher Ed. Press, Beijing, 2015.

[46] D. Talay and L. Tubaro. Expansion of the global error for numerical schemes solving stochastic differential equations. *Stochastic Anal. Appl.*, 8(4):483–509 (1991), 1990.

[47] G. Vilmart. Postprocessed integrators for the high order integration of ergodic SDEs. *SIAM J. Sci. Comput.*, 37(1):A201–A220, 2015.

[48] K. C. Zygalakis. On the existence and the applications of modified equations for stochastic differential equations. *SIAM J. Sci. Comput.*, 33(1):102–130, 2011.

## Appendix

Expression of $\tilde{\gamma}_2$ in Theorem 4.9.

$$ \begin{aligned} \tilde{\gamma}_2 = & \left( a(\bullet) - \frac{2}{\sigma^2} a(\bullet) - \frac{2}{\sigma^2} a(\bullet) + \frac{2}{\sigma^2} a(\bullet) + \frac{2}{\sigma^2} a(\bullet) - \frac{4}{\sigma^2} a(\bullet) - \frac{4}{\sigma^4} a(\bullet) \right. \\ & \left. + \frac{12}{\sigma^4} a(\bullet) - \frac{24}{\sigma^6} a(\bullet) \right) + \left( a(\bullet) - a(\bullet) - a(\bullet) + a(\bullet) + a(\bullet) - 2a(\bullet) \right. \\ & \left. - \frac{2}{\sigma^2} a(\bullet) + \frac{6}{\sigma^2} a(\bullet) - \frac{12}{\sigma^4} a(\bullet) \right) + \left( a(\bullet) - \frac{2}{\sigma^2} a(\bullet) + \frac{4}{\sigma^4} a(\bullet) + \frac{2}{\sigma^2} a(\bullet) \right. \\ & \left. - \frac{4}{\sigma^4} a(\bullet) + \frac{4}{\sigma^4} a(\bullet) - \frac{8}{\sigma^6} a(\bullet) \right) + \left( a(\bullet) - a(\bullet) - a(\bullet) - a(\bullet) + 2a(\bullet) \right) \\ & + \frac{2}{\sigma^2} a(\bullet) + a(\bullet) - 2a(\bullet) - \frac{4}{\sigma^2} a(\bullet) + \frac{8}{\sigma^2} a(\bullet) - \frac{16}{\sigma^4} a(\bullet) + \left( a(\bullet) \right. \\ & \left. - \frac{2}{\sigma^2} a(\bullet) - a(\bullet) + \frac{2}{\sigma^2} a(\bullet) + a(\bullet) + \frac{4}{\sigma^2} a(\bullet) - \frac{4}{\sigma^2} a(\bullet) + \frac{4}{\sigma^2} a(\bullet) - \frac{8}{\sigma^4} a(\bullet) \right) + \left( a(\bullet) \right. \\ & \left. + \left( a(\bullet) - a(\bullet) + a(\bullet) + a(\bullet) - a(\bullet) + a(\bullet) - \frac{2}{\sigma^2} a(\bullet) \right) + \left( a(\bullet) \right. \right. \\ & \left. \left. - \frac{2}{\sigma^2} a(\bullet) + \frac{4}{\sigma^4} a(\bullet) - \frac{8}{\sigma^6} a(\bullet) \right) + \left( a(\bullet) - a(\bullet) + 2a(\bullet) - \frac{2}{\sigma^2} a(\bullet) \right. \right. \\ & \left. \left. + \frac{2}{\sigma^2} a(\bullet) - \frac{6}{\sigma^2} a(\bullet) + \frac{12}{\sigma^4} a(\bullet) \right) + \left( a(\bullet) - \frac{2}{\sigma^2} a(\bullet) - \frac{2}{\sigma^2} a(\bullet) + \frac{4}{\sigma^4} a(\bullet) \right. \right. \\ & \left. \left. + \frac{4}{\sigma^2} a(\bullet) - \frac{12}{\sigma^4} a(\bullet) + \frac{24}{\sigma^6} a(\bullet) \right) \right) . \end{aligned} $$