# High order integration of stochastic dynamics on Riemannian manifolds with frozen flow methods

Eugen Bronasco, Adrien Busnot Laurent, Baptiste Nicolas Huguet

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

## To cite this version:

Eugen Bronasco, Adrien Busnot Laurent, Baptiste Nicolas Huguet. High order integration of stochastic dynamics on Riemannian manifolds with frozen flow methods. 2025. hal-05368742v1

HAL Id: hal-05368742

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

Preprint submitted on 17 Nov 2025 (v1), last revised 19 Dec 2025 (v2)

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
High order integration of stochastic dynamics on Riemannian manifolds with frozen flow methods

Eugen Bronasco<fnref n="1" />, Adrien Busnot Laurent<fnref n="2" /> and Baptiste Huguet<fnref n="2" />

July 15, 2025

Abstract

We present a new class of numerical methods for solving stochastic differential equations with additive noise on general Riemannian manifolds with high weak order of accuracy. In opposition to the popular approach with projection methods, the proposed methods are intrinsic: they only rely on geometric operations and avoid coordinates and embeddings. We provide a robust and general convergence analysis and an algebraic formalism of exotic planar Butcher series for the computation of order conditions at any high order. To illustrate the methodology, an explicit method of second weak order is introduced, and several numerical experiments confirm the theoretical findings and extend the approach for the sampling of the invariant measure of Riemannian Langevin dynamics.

Keywords: geometric numerical integration, stochastic differential equations, Riemannian manifolds, Riemannian Langevin, ergodicity, Lie-group methods, frozen flow, Butcher series, exotic series, post-Lie algebra, Hopf algebra, order conditions.

AMS subject classification (2020): 16T05, 41A58, 60H35, 37M25, 65L06, 70H45.

# 1 Introduction

Let ($\mathcal{M}$, $g$) be a smooth connected complete Riemannian manifold equipped with a metric $g$ and its Levi-Civita connection $\nabla$. Let $E_1, \dots, E_D$ be a frame, that is, $D$ smooth vector fields with $\dim(\mathcal{M}) \le D$ such that

$$ \mathrm{Span}_{\mathbb{R}}(E_1(y), \dots, E_D(y)) = T_y \mathcal{M}, \quad y \in \mathcal{M}. $$

Given a smooth vector field $F(x) = \sum_{d=1}^{D} f^d(x) E_d(x)$, we consider Stratonovich stochastic differential equations on $\mathcal{M}$ of the form,

$$ dX(t) = F(X(t))dt + \sqrt{2} \sum_{d=1}^{D} E_d(X(t)) \circ dW_d(t). \quad (1.1) $$

Equation (1.1) can be seen as a stochastic differential equation with additive noise w.r.t. the frame ($E_d$). Under reasonable assumptions, the flow of (1.1) is well-posed on $\mathcal{M}$. The aim of the present paper is to introduce a new class of intrinsic one-step integrators

$$ X_{n+1} = \Phi_h(X_n), \quad (1.2) $$

<footnote n="1">¹Université de Genève, Section de mathématiques, Genève, Switzerland. Eugen.Bronasco@unige.ch.</footnote>

<footnote n="2">²Univ Rennes, INRIA (Research team MINGuS), IRMAR (CNRS UMR 6625) and ENS Rennes, France. Adrien.Busnot-Laurent@inria.fr, Baptiste.Huguet@math.cnrs.fr.</footnote>
where $\Phi_h$ is a (random) local diffeomorphism of $\mathcal{M}$ for a timestep $h$ small enough, for solving dynamics of the form (1.1) with high order of accuracy in the weak sense. We take inspiration from Lie-group integrators [69, 76, 70, 71, 46, 72] and more precisely from Crouch-Grossman and commutator-free methods [28, 81, 24, 80]. In opposition to [63], we consider an approach with frames allowing us to define a new versatile class of methods that works on any smooth Riemannian manifold. An important feature of the new integrators is that they do not rely on the use of embeddings or on local coordinates. Their formulation, convergence analysis, order theory, and implementation are entirely intrinsic.

It is crucial that a numerical approximation respects the geometry of the problem and lies on $\mathcal{M}$, especially for the sampling of measures. In particular, our methods extend naturally for the sampling of measures interpreted as invariant measure of ergodic dynamics. Let us consider for instance the central example from molecular dynamics of Riemannian Langevin dynamics:

$$ 
\begin{aligned}
dX(t) &= -\sum_{d=1}^{D} (E_d[V]E_d + \nabla_{E_d}E_d)(X(t))dt + \sqrt{2}\sum_{d=1}^{D} E_d(X(t)) \circ dW_d(t) \quad (1.3) \\
&= -\nabla V(X(t))dt + \sqrt{2}dB_{\mathcal{M}}(t),
\end{aligned}
$$

where $V$ is a smooth potential, the frame $(E_d)$ is an orthonormal frame basis, and the vector field $F = -\nabla V - \sum \nabla_{E_d}E_d$ contains the Itô correction. Note that the choice $V = 0$ yields Brownian dynamics $B_{\mathcal{M}}(t)$ on $\mathcal{M}$. Riemannian Langevin dynamics (1.3) are naturally ergodic under mild assumptions, that is, the dynamics in long time behaves according to a deterministic probability measure, called the invariant measure, which takes the explicit form $d\mu_\infty = e^{-V}d\nu$. The generator of the SDE is $\mathcal{L} = -\nabla V \cdot \nabla + \Delta_g$, where $\Delta_g$ is the Laplace-Beltrami operator. Actually, an integrator of weak order $p$ automatically has at least order $p$ for the invariant measure [87]. We refer for instance to [32] and references therein for appropriate assumptions to obtain the ergodicity of a numerical scheme in the context of dynamics in $\mathbb{R}^d$. To the best of our knowledge, there does not exist any intrinsic method of order greater than one for sampling the invariant measure of Riemannian Langevin dynamics.

There exists a handful of extrinsic stochastic integrators on manifolds, that is, integrators that rely on an embedding of $\mathcal{M}$ in a Euclidean space of higher dimension. In the case of smooth embedded manifolds $\mathcal{M} = \{x \in \mathbb{R}^D, \zeta(x) = 0\}$ defined by constraints, the approach using projection methods is the most popular in the literature (see [59] and references therein). We cite in particular the constrained Euler method:

$$ X_{n+1} = X_n + hF(X_n) + \sqrt{2h}\xi_n + \lambda_n \nabla\zeta(X_n), \quad \zeta(X_{n+1}) = 0, $$

where $\lambda_n$ is a Lagrange multiplier and the $\xi_n$ are independent standard Gaussian random vectors. While the approach using projection methods is straightforward to implement and widely used in applications, these methods rely on an embedding of $\mathcal{M}$ into a vector space of much higher dimension and use non-intrinsic quantities, so that they often face severe timestep restrictions. Moreover, the algebraic structure associated to order theory is difficult and does not allow at the present time for the creation of methods of arbitrarily high order, nor for a systematic approach to stochastic backward error analysis [15]. Stochastic Lie group methods seldom appear in the literature [63, 30, 64, 3, 65, 27], but focus on low order approximations on specific manifolds (typically variants of Lie-Euler on the sphere), with an extrinsic convergence analysis, and no high order theory. In the recent work [10], a new approach for sampling (1.3) is developed using the following intrinsic integrator, called Riemannian Langevin method

$$ X_{n+1} = \exp^{\text{Riem}}(-h\nabla V(X_n) + \sqrt{2h}\xi_{n+1})X_n, $$
where exp<sup>Riem</sup> is the geodesic exponential for the Levi-Civita connection ∇. This method has order one for solving (1.3) in the weak sense and for the invariant measure. One could be tempted to say that such a method is exact for V = 0, but it still has only order one, in contrast to the Euclidean setting where it is exact in law.

We introduce the first intrinsic methods of high weak order for solving (1.1) using only evaluations of F and the frame. In particular, we highlight the following new explicit integrator of second weak order for solving (1.1), and using the minimal number of evaluations of F per step:

$$
\begin{align*}
H_n &= \exp \left( \sum_{d=1}^{D} \left( \frac{1}{2} h f^d(X_n) + \sqrt{h} \xi_n^{d,1} \right) E_d \right) X_n \\
X_{n+1} &= \exp \left( \sum_{d=1}^{D} \left( \left( \frac{\sqrt{2}}{2} - 1 \right) h f^d(X_n) + (2 - \sqrt{2}) h f^d(H_n) + (1 - \sqrt{2}) \sqrt{h} \xi_n^{d,1} + \sqrt{h} \xi_n^{d,2} \right) E_d \right) \\
&\quad \exp \left( \sum_{d=1}^{D} \left( \left( 1 - \frac{\sqrt{2}}{2} \right) h f^d(X_n) + (\sqrt{2} - 1) h f^d(H_n) + \sqrt{2h} \xi_n^{d,1} \right) E_d \right) X_n.
\end{align*}
$$

where the $\xi_n^{d,l}$ are independent standard Gaussian random variables (or discrete symmetric approximations with matching moments up to order 4) and we use frozen flows (see Section 2.3). An originality of our approach is that we present a general framework for the analysis of a new class of stochastic Lie-group methods at any order on any manifold in the weak sense using a new algebraic formalism similar to the one of Butcher series. The new methods outperform in accuracy, and versatility the previous attempts.

The calculations of the order conditions are intricate and require the use of appropriate
algebraic tools. We propose a new extension of exotic and Lie-Butcher series, called exotic Lie
series, for the systematic computation of order conditions of the new frozen flow methods at any
order in the weak sense. While Butcher trees and series [40, 21, 67, 22] represent naturally Taylor
expansions of deterministic flows in $\mathbb{R}^D$, planar trees and Lie-Butcher series were introduced for
the study of flows on manifolds [46, 80, 78] and later extended for the study of the connection
algebra [77, 34, 5, 39, 73]. On the other hand, the formalism of exotic series was introduced in
[54, 50, 14] for the creation of integrators for solving stochastic dynamics with additive noise
with high order in the weak sense and for the invariant measure. This formalism, combined
with the aromatic B-series [26, 47, 13, 52, 51], was extended in [55] into the exotic aromatic
series formalism for the numerical integration of SDEs on embedded manifolds with projection
methods. The geometric and algebraic properties of the exotic formalism of trees were later
studied in [15, 53] (see also [25, 66, 74, 12] in the deterministic setting). The modern approach
to such algebraic formalism relies extensively on Hopf algebras [25, 23, 62, 12, 14, 82, 15, 35, 20]
that we identify here in the context of stochastic numerics on manifolds.

The paper is organised the following way. We present in Section 2 the main numerical
results of this paper: a robust framework for the Riemannian convergence analysis of stochastic
integrators on manifolds, the new stochastic frozen flow methods, and the associated order
conditions for high weak order. Section 3 introduces the formalism of planar decorated and
exotic forests and translates the geometric operations in the connection algebra in terms of
algebraic operations on graphs. This approach is then applied for the explicit description at any
order of the weak Taylor expansion of the exact flow of (1.1), of the new frozen flow methods,
and of their associated order conditions. The Hoph algebra structure of exotic planar forests is
further detailed in Section 4. In particular, we generalise the Munthe-Kaas-Wright Hopf algebra
[78] to the stochastic context. Numerical experiments confirm the theoretical findings in Section 5 and outlooks and future works are presented in Section 6.

# 2 New frozen flow methods of high weak order

We generalise the standard Euclidean framework for weak convergence in the context of Riemannian manifolds, and we introduce new frozen flow integrators and their order conditions for high weak order. We emphasize that the following convergence analysis is new, even for the underlying deterministic differential equations (see also [29]).

## 2.1 Stochastic differential equations on manifolds

For $v \in TM$, we denote $|v| = g(v, v)^{1/2}$, its Riemannian norm. Let $o$ an arbitrary point on $M$. For $x \in M$, we define $r(x) = d(o, x)$ the Riemannian distance function. This function is 1-Lipschitz and the squared distance $r^2$ is smooth on $M\setminus Cut_o$. This allows us to define moments of a $M$-valued random variable.

Let us introduce the main classes of test functions. The space $C_c^\infty(M)$ denotes the class of smooth compactly supported functions. For $\phi \in C_c^\infty(M)$ we denote $d\phi$ its differential and $|d\phi|$ its Riemannian norm on $T^*M$, inherited from the norm on $TM$. We denote $v[\phi](x)$ the differential of a test function $\phi$ in a direction $v \in T_xM$. We shall use the wider class of test functions $C_p^\infty(M)$, which consists in smooth functions satisfying polynomial growth estimates of the form:

$$|E_{d_q}[\dots E_{d_1}[\phi]\dots] | (x) \le C(1+r(x)^K), \quad q = 0, 1, \dots$$

Given a vector field $F$ and a fixed decomposition in the frame $F = \sum f^d E_d$, we say that $F \in \mathcal{X}_P(M)$ if its components are Lipschitz continuous, $|df^d| \le C$, and satisfy: $f^d \in C_p^\infty(M)$, $d = 1, \dots, D$. We mention that considering test functions and vector fields with polynomial growth of their derivatives up to a given order $p+2$ would be sufficient for the creation of numerical schemes of order up to $p$.

**Assumption 2.1.** The vector fields $E_1, \dots, E_D$ are smooth and bounded. The vector field $F$ is in $\mathcal{X}_P(M)$. Moreover, there exist constants $\nu \ge 1$ and $\lambda \in \mathbb{R}$ such that on $M\setminus Cut_o$, we have

$$\mathcal{L}r^2 \le \nu + \lambda r^2, \qquad (2.1)$$

where $\mathcal{L}$ is the generator

$$\mathcal{L}\phi = F[\phi] + \sum_{d=1}^{D} E_d[E_d[\phi]].$$

Under the smoothness assumption on $F$, $E_1, \dots, E_D$, there exists a unique solution of (1.1), up to a stopping time $\tau$. Following [42], we recall that a stochastic process $X$, defined on a stochastic interval $[0, \tau]$, is a solution of (1.1) if and only if for all test functions $\phi \in C_c^\infty(M)$, we have

$$\phi(X(t)) = \phi(X_0) + \int_0^t F[\phi](X(t))dt + \sqrt{2} \sum_{d=1}^D \int_0^t E_d[\phi](X(t)) \circ dW_d(t), \quad 0 < t < \tau.$$

Unlike the Euclidean case, where Lipschitz conditions are sufficient to avoid finite time explosion, on a manifold, a blow-up can be induced from the geometry, even for a Brownian motion (see
[42] for examples). The Lyapunov-like assumption (2.1) ensures the completeness of the SDE. This means that the solution $X$ of (1.1) is defined on $\mathbb{R}_+$, i.e. $\tau = +\infty$ a.s. In addition with the completeness, it provides finite moments estimates (even finite exponential moments) for the solution of the SDE. This assumption is adapted from [88]. In the case where $(E_i)_i$ is a basis there exists a new metric $\tilde{g}$ such that $(E_i)_i$ is an orthonormal basis. Using this metric and the associated Levi-Civita connection $\tilde{\nabla}$, we find

$$ \mathcal{L}\phi = \Delta_{\tilde{g}}\phi + (F + \tilde{\nabla}_{E_i}E_i)[\phi], $$

which is the framework of [88]. The general case is treated likewise, by extracting local basis on a cover of $\mathcal{M}$.

There exist various simpler criteria which ensure that inequality (2.1) is satisfied. The first criterion is compactness, which can be used for classical spaces such as the sphere $S^n$ and the Lie group $SO_n(\mathbb{R})$. This condition being seldom satisfied, a second classical criterion is the Bakry-Émery criterion from [8]. This criterion applies to the case of equation (1.3), where the generator writes as $\mathcal{L} = \Delta - \nabla V \cdot \nabla$. One says that the potential $V$ satisfies the criterion if there exists $\kappa \in \mathbb{R}$ such that

$$ \mathrm{Ric} + \mathrm{Hess}(V) \geq \kappa. $$

The term Ric designs the Ricci tensor on $(\mathcal{M}, g)$ and is defined as the trace of the Riemann tensor. In this context, the symmetric operator $\mathrm{Ric} + \mathrm{Hess}(V)$ is called the Bakry-Émery curvature and the term $\mathrm{Hess}(V)$ can be interpreted as an additional curvature. Note that this handy criterion has been adapted in [6] for equation (1.1) with a lower bound of the operator $\mathrm{Ric} - \nabla F$. We also cite the moment conditions from [60] which imply stochastic completeness.

By a standard induction, we obtain the following estimates.

**Lemma 2.2.** Under Assumption 2.1, for all $q \ge 1$ we have $\mathcal{L}(r^{2q}) \le C(1 + r^{2q})$.

*Proof.* The proof comes from classical $\Gamma$ calculus for diffusions (see [9] for instance). For all $u \in C^2(\mathcal{M})$ and $\eta \in C^2(\mathbb{R})$, we have

$$ \mathcal{L}\phi(u) = \eta'(u)\mathcal{L}u + \eta''(u)\Gamma(u), $$

where

$$ \Gamma(u) = \sum_{i=1}^{D} |E_i[u]|^2. $$

Applying this formula to $u = r^2$ and $\eta: x \mapsto x^q$, we obtain the result. $\square$

A key tool for the weak expansion of the exact flow is the Kolmogorov equation. Under Assumption 2.1, the SDE (1.1) generates a semigroup, well-defined on $C_p^\infty(\mathcal{M})$, and Markovian. Namely, for all $\phi \in C_c^\infty(\mathcal{M})$, the function $u(t, x) = \mathbb{E}[\phi(X(t))|X(0) = x]$ satisfies the Kolmogorov equation

$$ \partial_t u = \mathcal{L}u, \quad u(0, x) = \phi(x). \qquad (2.2) $$

In order to compare an expansion of the exact solution an expansion of a numerical integrator, we need more regularity on the semigroup itself. Since the generator $\mathcal{L}$ is elliptic, the semigroup $u$ is $C^\infty$ on $]0, T[ \times \mathcal{M}$. We further assume that the semigroup preserves the $C_p^\infty$ regularity.

**Assumption 2.3 (Regularity).** For all test function $\phi \in C_p^\infty(\mathcal{M})$, the semigroup (2.2) satisfies $u \in C^\infty([0, T[, C_p^\infty(\mathcal{M}))$, that is, for all $k \ge 0$, there exist $C > 0$ and $\kappa \ge 0$, depending on $k$, $T$, and $\phi$ such that for all $x \in \mathcal{M}$

$$ \sup_{t \in ]0, T[} |\partial_t^k u(t, x)| \le C(1 + r^\kappa(x)). $$
Note that Assumption 2.3 is satisfied automatically in $\mathbb{R}^d$ with our regularity assumptions. For clarity, we leave the study of sufficient conditions to obtain polynomial growth estimates of the semigroup for future works.

We are now able to derive the Taylor expansion of the semigroup.

**Proposition 2.4.** Under Assumption 2.3, for all $\phi \in C_p^\infty(M)$ and all $h \le h_0$ small enough, $u$ has the Taylor expansion

$$u(h, x) = \phi(x) + \sum_{j=1}^{p} \frac{h^j}{j!} \mathcal{L}^j \phi(x) + h^{p+1} R_p^h(\phi, x), \quad x \in \mathcal{M}, \qquad (2.3)$$

with $|R_p^h(\phi, x)| \le C(1 + r(x)^K)$, where $C$ and $K$ do not depend on $x$ and $h$.

## 2.2 Weak order theory on Riemannian manifolds

Following the Euclidean works [87, 68], we present a general theory for achieving high order weak estimates for intrinsic methods. The approach does not rely on an embedding and does not assume compactness of $M$. In order to obtain estimates, we use the Riemannian structure of $M$. This structure is not needed to perform the calculations of order conditions of the new methods.

**Definition 2.5.** An intrinsic one-step integrator of the form (1.2) is of local weak order $p$ if for all $\phi \in C_p^\infty(M)$, there exists $C > 0$ such that the following estimate holds for all $x \in M$ and all $h \le h_0$ small enough,

$$|\mathbb{E}[\phi(X_1)|X_0 = x] - \mathbb{E}[\phi(X(h))|X_0 = x]| \le C(1 + r(x)^K)h^{p+1}.$$

Similarly, the integrator is of global weak order $p$ if for all $\phi \in C_p^\infty(M)$, $T > 0$, $Nh = T$, $h \le h_0$ small enough, and $x \in M$, the following estimate is satisfied:

$$\sup_{n=0,...,N} |\mathbb{E}[\phi(X_n)] - \mathbb{E}[\phi(X(nh))]| \le Ch^p.$$

Let us assume that the integrator satisfies a Taylor expansion, often called Talay-Tubaro expansion [86] in the Euclidean context, that is similar to (2.3).

**Assumption 2.6.** For all $\phi \in C_p^\infty(M)$ and all $h \le h_0$ small enough, the numerical integrator (1.2) has a weak Taylor expansion of the form

$$\mathbb{E}[\phi(X_1)|X_0 = x] = \phi(x) + \sum_{j=1}^{p} h^j \mathcal{A}_j \phi(x) + h^{p+1} R_p^h(\phi, x), \quad x \in \mathcal{M}, \qquad (2.4)$$

where the remainder satisfies $|R_p^h(\phi, x)| \le C(1 + r(x)^K)$ and the $\mathcal{A}_j$ are linear differential operators.

We adapt the standard assumption of bounded moments in the Riemannian context. We shall give sufficient assumptions on the integrator in Subsection 2.3 so that this assumption is satisfied.

**Assumption 2.7.** The numerical integrator (1.2) has bounded moments of any order $K \ge 0$,

$$\sup_{n=0,...,N} \mathbb{E}[r(X_n)^K] \le C(T), \quad Nh = T. \qquad (2.5)$$
Under the above assumptions, we obtain the following convergence result, in the spirit of the Euclidean theory [86, 87, 68].

**Theorem 2.8.** Consider a one-step integrator of the form (1.2) for solving (1.1). Under Assumptions 2.1, 2.3, 2.6, and 2.7, if the Taylor expansion (2.4) of the integrator satisfies

$$A_j = \frac{1}{j!} \mathcal{L}^j, \quad j = 1, \dots, p, \qquad (2.6)$$

then the integrator has global weak order p.

*Proof.* Let $X^x(t_n)$ be the exact solution of (1.1) at time $t_n = nh$ starting from $x \in \mathcal{M}$ and $X_n^x$ be the n-th step of the numerical solution. The global error satisfies, thanks to a telescopic sum argument,

$$\mathbb{E}[\phi(X_N) - \phi(X(T))] = \sum_{i=1}^{N} \mathbb{E}[\tilde{\phi}(X_1^{X_{N-i}}) - \tilde{\phi}(X^{X_{N-i}}(h))], \quad \tilde{\phi}(x) = \mathbb{E}[\phi(X^x(t_{i-1})].$$

Assumption 2.3 yields that $\tilde{\phi} \in C_p^\infty(\mathcal{M})$. Then, the local order p condition, together with Assumptions 2.6 and 2.7 give

$$\mathbb{E}[\phi(X_N) - \phi(X(T))] \le \sum_{i=1}^{N} Ch^{p+1} (1 + \mathbb{E}[r(X_{N-i}^x)^K]) \le CNh^{p+1} \le Ch^p.$$

Hence, the global weak order p.

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

## 2.3 Stochastic frozen flow integrators

The deterministic Lie-group methods are divided into two categories that correspond to the two kinds of coordinates: the Runge-Kutta-Munthe-Kaas (RKMK) methods (see [69, 76, 70, 71]) and the frozen flow methods. We choose the latter option in this paper for the sake of generality, as the RKMK approach relies on the existence of a transitive group action on $\mathcal{M}$. The frozen flow methods originally appeared as Crouch-Grossmann methods [28, 81]. The generalisation presented in the papers [24, 80] allowed the creation of more efficient methods, called commutator-free Lie-group methods.

The methods rely heavily on the concept of frozen vector fields and flows: given a decomposition of a (potentially random) vector field in the frame $G = \sum_d g^d E_d$, the frozen vector field $G_x \in \mathfrak{X}(\mathcal{M})$ at point $x \in \mathcal{M}$ is

$$G_x(p) = \sum_d g^d(x) E_d(p).$$

The solution of the ODE $y'(t) = G(y(t))$ with $y(0) = p$ is denoted $\exp(tG)p$. The frozen flow $\exp(tG_x)p$ is the exact flow driven by the frozen vector field $G_x$.

**Lemma 2.9.** If $G \in \mathfrak{X}_P(\mathcal{M})$ then, the frozen flow $\exp(tG_x)p$ is globally defined. Moreover, it satisfies

$$d(p, \exp(tG_x)p) \le C(1+r(x)) \qquad (2.7)$$

*Proof.* The coefficients of the ODE being globally Lipschitz continuous, there exists a global solution. Let us denote $\gamma: t \mapsto \exp(tG_x)p$. It is a smooth curve, joining $p$ and $\exp(G_x)p$. As $|E_d|$ are bounded, $\gamma$ has a bounded speed

$$|\dot{\gamma}_t|^2 \le D \sum_{d=1}^{D} |g^d(x)|^2 |E_d(\gamma_t)|^2 \le C(1+r^2(x)).$$
Therefore, we have

$$d(p, \exp(tG_x)p) \le \int_0^t |\dot{\gamma}_s| ds \le C(1+r(x))$$

Hence, the result.

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

**Remark 2.10.** The frozen vector fields depend on a choice of frame ($E_d$), but the definitions of the frozen vector field, and thus, of the frozen flow, do not depend on the choice of a decomposition of $G$ in this frame. Indeed, given two decomposition of $G$ in the frame ($E_d$)

$$G(p) = \sum_{d=1}^{D} g_d(p) E_d(p) = \sum_{d=1}^{D} \hat{g}_d(p) E_d(p),$$

we define two, presumably different, frozen vector fields $G_x$ and $\hat{G}_x$ associated to these decompositions. Now, let us extract a basis ($E_d$)$_{d \in I}$ of the frame valid in an open set $U \subset \mathcal{M}$. The vector field $G_x$ and $\hat{G}_x$ decompose uniquely in $U$ as

$$G_x(p) = \sum_{d \in I} \alpha_d E_d(p), \quad \hat{G}_x(p) = \sum_{d \in I} \hat{\alpha}_d E_d(p).$$

As $G(x) = G_x(x) = \hat{G}_x(x)$, we obtain $\alpha_d = \hat{\alpha}_d$ and $G_x = \hat{G}_x$. Therefore, the formulation of the integrators is independent of the choice of a decomposition in the frame.

The central assumption that motivates the use of these methods is that while computing the exact flow of a SDE is difficult, computing the flow associated to a frozen vector field often is a much simpler problem for a good choice of frame. For instance, in $\mathbb{R}^D$ with the standard frame basis $E_i = \partial_i$, the frozen flow is the Euler method: $\exp(G_x)p = p + G(x)$. A second example is given by homogeneous manifolds. In this context, a natural choice of frame is derived from a basis ($A_d$) of the Lie algebra $\mathfrak{g}$ by $E_d(y) = A_d \cdot y$, where we use for simplicity the same notation $\cdot$ for the Lie group and Lie algebra action on $\mathcal{M}$. On matrix homogeneous manifolds, the frozen flow is explicitly given by

$$\exp(G_x)p = \mathrm{Exp} \left( \sum_d g^d(x) A_d \right) p,$$

where $\mathrm{Exp}$ is the matrix exponential. We emphasize that a significant part of the manifolds used in numerical experiments and applications of stochastic dynamics on manifolds are homogeneous manifolds, such as Lie-groups, spheres, the Stiefel manifold, symmetric positive definite matrices, ... We further discuss the implementation of the methods in Section 5.1.

We consider the following new class of stochastic frozen flow methods, in the spirit of [81] for the deterministic part and [83, 31, 54] for the stochastic part,

$$ 
\begin{aligned}
H_n^i &= \exp\left(\sum_{d=1}^D \left(h \sum_{j=1}^s Z_{i,j,K}^0 f^d(H_n^j) + \sqrt{h} Z_{i,K}^d\right) E_d\right) \dots \\
&\quad \exp\left(\sum_{d=1}^D \left(h \sum_{j=1}^s Z_{i,j,1}^0 f^d(H_n^j) + \sqrt{h} Z_{i,1}^d\right) E_d\right) X_n, \\
X_{n+1} &= \exp\left(\sum_{d=1}^D \left(h \sum_{i=1}^s z_{i,K}^0 f^d(H_n^i) + \sqrt{h} z_K^d\right) E_d\right) \dots \\
&\quad \exp\left(\sum_{d=1}^D \left(h \sum_{i=1}^s z_{i,1}^0 f^d(H_n^i) + \sqrt{h} z_1^d\right) E_d\right) X_n,
\end{aligned}
\quad (2.8) $$
where the coefficients $z_{i,K}^0$, $Z_{i,j,k}^0$ are deterministic and $z_k^d$, $Z_{i,k}^d$ are of the form

$$z_k^d = \sum_{l=1}^{L} \hat{z}_{k,l} \xi_n^{d,l}, \quad Z_{i,k}^d = \sum_{l=1}^{L} \hat{Z}_{i,k,l} \xi_n^{d,l},$$

where the $\xi^{d,l}$ are independent standard Gaussian random variables or approximations of Gaussian random variables with the same $2p$ first moments, vanishing odd moments, and finite moments of all order. The dependence in $n$ of the Runge-Kutta coefficients $z_k^d$, $Z_{i,k}^d$ is omitted for clarity, and we enforce that the random coefficients between two steps in time are independent. The constant $K$ is the number of exponentials used per stage. It is well-known [46] that choosing $K = 1$ does not allow to reach arbitrarily high order already in the context of ODEs. In particular, the simplest method is the frozen flow Euler method

$$X_{n+1} = \exp \left( \sum_{d=1}^{D} \left( h f^d(X_n) + \sqrt{2h} \xi_n^d \right) E_d \right) X_n. \quad (2.9)$$

The assumptions of Subsection 2.1 are in general not sufficient to ensure bounded moments and bounds on the Talay-Tubaro remainder of the numerical integrators. Aside from the compactness assumption, the following assumption is sufficient.

**Proposition 2.11.** Assume that the frame $E_d$, $d = 1, \dots, D$, and the vector field $F$ satisfy Assumption 2.1 and

$$|E_{d_1}[E_{d_1}[f^d]]| \le C.$$

Assume further that for an infinite number of arbitrarily large $q$:

$$|E_{d_k}[\dots E_{d_1}[r^q]\dots]| \le C(1+r^{q-k}), \quad k \le 4.$$

Then, all consistent explicit frozen flow methods of the form (2.8) satisfy the bounded moments property and the Talay-Tubaro expansion of Assumptions 2.6 and 2.7.

*Proof.* Let us first show that the first remainder in the expansion (2.4) acts on the distance function by

$$|R_1^h(r^q, x)| \le C(1 + r(x)^q), \quad (2.10)$$

for an infinite number of arbitrarily large $q$. Without loss of generality, we consider the following method, where $G_y = \sum_{d=1}^D g_y^d E_d$ and $g_y^d = h f^d(y) + \sqrt{2h} \xi_n^d$,

$$X_{n+1} = \exp(G_{H_n}) X_n, \quad H_n = \exp(G_{X_n}) X_n.$$

A Taylor expansion around $x = X_0$ yields (see Section 3 for the detailed calculations),

$$f^d(H) = f^d(x) + G_x[\phi](x) + \int_0^1 (1-t) g_x^{d_1} g_x^{d_2} E_{d_2}[E_{d_1}[f^d]] (\exp(tG_x)x) dt,$$

and we obtain

$$|f^d(H)| \le C(\xi^d)(1+r(x)), \quad |\mathbb{E}[f^d(H) - f^d(x)]| \le Ch(1+r(x)),$$

where $C(\xi^d)$ is a polynomial in $|\xi^d|$. Then, a Taylor expansion for $\phi = r^q$ yields

$$\mathbb{E}[\phi(X_1)] = \phi(x) + h\mathcal{L}\phi(x) + h^2 R_1^h(\phi, x)$$
$$
\begin{align*}
R_1^h(\phi, x) ={}& \mathbb{E}\left[\frac{1}{h}(f^d(H) - f^d(x))E_d[\phi](x) + \frac{1}{2}f^{d_1}(H)f^{d_2}(H)E_{d_2}[E_{d_1}[\phi]](x)\right. \\
& \left. + \frac{1}{3}f^{d_1}(H)\xi^{d_2}\xi^{d_3}(E_{d_3}[E_{d_2}[E_{d_1}[\phi]]] + E_{d_3}[E_{d_1}[E_{d_2}[\phi]]] + E_{d_1}[E_{d_3}[E_{d_2}[\phi]]])\right. \\
& \left. + \frac{h}{6}f^{d_1}(H)f^{d_2}(H)f^{d_3}(H)E_{d_3}[E_{d_2}[E_{d_1}[\phi]]](x)\right. \\
& \left. + \frac{1}{h^2}\int_0^1 \frac{(1-t)^3}{3!}g_H^{d_1}g_H^{d_2}g_H^{d_3}g_H^{d_4}E_{d_4}[E_{d_3}[E_{d_2}[E_{d_1}[\phi]]]](\exp(tG_H)x)dt\right].
\end{align*}
$$

The estimate on $R_1^h(\phi, x)$ is straightforwardly obtained from (2.7). Indeed, one shows that the distance function satisfies

$$r(\exp(tG_H)x) \le C(\xi^d)(1 + r(x)),$$

where $C(\xi^d)$ is a polynomial in $|\xi^d|$.

The expansion (2.4) is computed explicitly in Theorem 3.16. The estimate on the remainder is obtained thanks to the regularity Assumption 2.3. For Assumption 2.7, we use the Talay-Tubaro expansion (2.4) applied to $\phi = r^q$:

$$\mathbb{E}[r(X_{n+1})^q] = \mathbb{E}[r(X_n)^q] + h\mathbb{E}[\mathcal{L}r(X_n)^q] + h^2\mathbb{E}[R_1^h(r^q, X_n)],$$

where $|\mathbb{E}[\mathcal{L}\phi(X_n)]| \le C(1 + \mathbb{E}[r(X_n)^q])$ by Proposition 2.2 and $|\mathbb{E}[R_1^h(r^q, X_n)]| \le C(1 + \mathbb{E}[r(X_n)^q])$ by equation (2.10) (for $q$ large enough). Thus, we find

$$\mathbb{E}[r(X_{n+1})^q] \le C(1+h)\mathbb{E}[r(X_n)^q] + Ch,$$

and a standard discrete Grönwall argument yields the bounded moments (2.5). $\square$

## 2.4 High-order frozen flow integrators

Let us present the order conditions of the new class of frozen flow methods (2.8) for solving the stochastic dynamics (1.1), as well as an example of second order frozen flow integrator.

We use the notation of factorial of multi-indices and sums for writing the order conditions.

**Definition 2.12.** If $k \in \{1, \dots, K\}^n$ is a multi-index that takes $r_1$ times the value $k_1, \dots, r_p$ times the value $k_p$, the factorial of $k$ is

$$k! = r_1! \dots r_p!.$$

The factorial sum is given for a finite set $S \subset \{1, \dots, K\}^n$ of multi-indices and a functional $f$ by $\sum_{k \in S} f(k) = \sum_{k \in S} \frac{1}{k!}f(k)$.

The weak order conditions for stochastic Runge-Kutta methods in $\mathbb{R}^d$ are presented in [83, 31, 54, 14] and require 4 order conditions for weak order two. On embedded manifolds, the approach with projection methods of [55] yields 28 order conditions, applies only in codimension one, and relies on an extrinsic approach with tedious calculations (see [55, App. D]). The approach proposed here is intrinsic, versatile, generalises the order theory of Lie-group methods to the stochastic setting, and yields a surprisingly low number of order conditions for weak order two.

**Theorem 2.13.** Let a frozen flow integrator of the form (2.8) that satisfies the order conditions in Table 1, then, under Assumption 2.3 and the assumptions of Proposition 2.11, the integrator has local weak order two for solving (1.1). The additional order conditions for weak order three are presented in Appendix A.
<table><caption>Table 1: Order conditions of frozen flow methods up to weak order 2. The sums on the indices that do not satisfy inequalities (except the $d_i$) are omitted for clarity. The order conditions do not depend on the dimension of the problem.</caption><thead><tr><th>Exotic forest π</th><th>Differential F(π)[φ]</th><th>Order condition a(π) = e(π)</th></tr></thead><tbody><tr><td>??</td><td>f<sup>i</sup>E<sub>i</sub>[φ]</td><td>z<sup>0</sup><sub>i,k<sub>1</sub></sub> = 1</td></tr><tr><td>??</td><td>E<sub>d<sub>1</sub></sub>[E<sub>d<sub>1</sub></sub>[φ]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>2</sub></sub> E[z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub>z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 1</td></tr><tr><td>??</td><td>f<sup>j</sup>E<sub>j</sub>[f<sup>i</sup>]E<sub>i</sub>[φ]</td><td>z<sup>0</sup><sub>i,k<sub>2</sub></sub> Z<sup>0</sup><sub>i,j,k<sub>1</sub></sub> = 1/2</td></tr><tr><td>??</td><td>E<sub>d<sub>1</sub></sub>[E<sub>d<sub>1</sub></sub>[f<sup>i</sup>]]E<sub>i</sub>[φ]</td><td>∑<sup>1</sup><sub>k<sub>2</sub>≥k<sub>3</sub></sub> E[Z<sup>d<sub>1</sub></sup><sub>i,k<sub>3</sub></sub> Z<sup>d<sub>1</sub></sup><sub>i,k<sub>2</sub></sub>]z<sup>0</sup><sub>i,k<sub>1</sub></sub> = 1/2</td></tr><tr><td>??</td><td>f<sup>j</sup>f<sup>i</sup>E<sub>j</sub>[E<sub>i</sub>[φ]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>2</sub></sub> z<sup>0</sup><sub>i,k<sub>2</sub></sub> E[Z<sup>d<sub>1</sub></sup><sub>i,k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 0</td></tr><tr><td>??</td><td>E<sub>d<sub>1</sub></sub>[f<sup>i</sup>]E<sub>d<sub>1</sub></sub>[E<sub>i</sub>[φ]]</td><td>∑<sup>1</sup><sub>k<sub>2</sub>≥k<sub>1</sub></sub> z<sup>0</sup><sub>i,k<sub>2</sub></sub> E[Z<sup>d<sub>1</sub></sup><sub>i,k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 1</td></tr><tr><td>??</td><td>f<sup>i</sup>E<sub>i</sub>[E<sub>d<sub>1</sub></sub>[E<sub>d<sub>1</sub></sub>[φ]]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub></sub> z<sup>0</sup><sub>i,k<sub>3</sub></sub> E[z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 1/2</td></tr><tr><td>??</td><td>f<sup>i</sup>E<sub>d<sub>1</sub></sub>[E<sub>i</sub>[E<sub>d<sub>1</sub></sub>[φ]]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>3</sub>≥k<sub>2</sub></sub> z<sup>0</sup><sub>i,k<sub>3</sub></sub> E[z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 0</td></tr><tr><td>??</td><td>f<sup>i</sup>E<sub>d<sub>1</sub></sub>[E<sub>d<sub>1</sub></sub>[E<sub>i</sub>[φ]]]</td><td>∑<sup>1</sup><sub>k<sub>3</sub>≥k<sub>1</sub>≥k<sub>2</sub></sub> z<sup>0</sup><sub>i,k<sub>3</sub></sub> E[z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 1/2</td></tr><tr><td>??</td><td>E<sub>d<sub>2</sub></sub>[E<sub>d<sub>2</sub></sub>[E<sub>d<sub>1</sub></sub>[E<sub>d<sub>1</sub></sub>[φ]]]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> E[z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub>] E[z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 1/2</td></tr><tr><td>??</td><td>E<sub>d<sub>2</sub></sub>[E<sub>d<sub>1</sub></sub>[E<sub>d<sub>2</sub></sub>[E<sub>d<sub>1</sub></sub>[φ]]]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> E[z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub>] E[z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 0</td></tr><tr><td>??</td><td>E<sub>d<sub>1</sub></sub>[E<sub>d<sub>2</sub></sub>[E<sub>d<sub>2</sub></sub>[E<sub>d<sub>1</sub></sub>[φ]]]]</td><td>∑<sup>1</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> E[z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub>] E[z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub>] = 0</td></tr></tbody></table>

Theorem 2.13 is a direct consequence of Theorem 3.14 and Theorem 3.16, postponed to Section 3. The approach relies on a new algebraic framework of exotic planar forests and series that gives the explicit expression of the operators $A_j$ in (2.6). It allows us to provide the order conditions for arbitrarily high weak order.

**Remark 2.14.** We shall prove in Section 3.5 that the coefficient maps $a$ and $e$ of the numerical and exact flow are characters of the shuffle algebra over exotic forests ($\mathcal{E}\mathcal{F}, \sqcup, \Delta$). Following [80], an important consequence is that the order conditions are not independent. Up to second order, the following "shuffle relations" are satisfied (and similarly for $e$):

$$ 
\begin{aligned}
a(??)^2 &= 2a(??), \\
a(??)a(??) &= a(??) + a(??) + a(??), \\
a(??)^2 &= 2a(??) + 2a(??) + 2a(??).
\end{aligned}
$$

The number of order conditions is not given by the number of entries in Table 1, as one has to take into account the shuffle relations. In particular, there are 2 conditions for order 1, 8 conditions for order 2, and 73 conditions for order 3 (against 2, 11, and 95 exotic forests of order 1,2, and 3). This is further discussed in Remark 3.17.

Let us now provide an example of a second order method. Let $\xi_n^{d,1}$, $\xi_n^{d,2}$, $d = 1, \dots, D$, be bounded independent random variables satisfying

$$ \mathbb{E}[(\xi_n^{d,l})^2] = 1, \quad \mathbb{E}[(\xi_n^{d,l})^4] = 3, \quad \mathbb{E}[(\xi_n^{d,l})^{2q+1}] = 0, \quad q = 0, 1, 2, \dots $$

For weak order two, one can use for instance

$$ \mathbb{P}(\theta = 0) = \frac{2}{3}, \quad \mathbb{P}(\theta = \pm\sqrt{3}) = \frac{1}{6} $$
so that the first four moments coincide with the ones of standard Gaussians.

The new method is the following. It is explicit, intrinsic, uses two evaluations of the vector field $F$ per step, and three frozen flow exponentials.

**Method SFF2** (New explicit frozen flow integrator of weak order two for equation (1.1))

$X_0 \in \mathcal{M}$
for $n = 0, \dots, N$ do

$$H_n^1 = X_n \\ H_n^2 = \exp \left( \sum_{d=1}^{D} \left( \frac{1}{2} h f^d(H_n^1) + \sqrt{h} \xi_n^{d,1} \right) E_d \right) X_n \quad (2.11)$$

$$X_{n+1} = \exp \left( \sum_{d=1}^{D} \left( \left( \frac{\sqrt{2}}{2} - 1 \right) h f^d(H_n^1) + \left( 2 - \sqrt{2} \right) h f^d(H_n^2) + \left( 1 - \sqrt{2} \right) \sqrt{h} \xi_n^{d,1} + \sqrt{h} \xi_n^{d,2} \right) E_d \right) \\ \exp \left( \sum_{d=1}^{D} \left( \left( 1 - \frac{\sqrt{2}}{2} \right) h f^d(H_n^1) + \left( \sqrt{2} - 1 \right) h f^d(H_n^2) + \sqrt{2h} \xi_n^{d,1} \right) E_d \right) X_n.$$

end for

**Example 2.15.** Brownian dynamics on $\mathcal{M}$ are defined as the solution of

$$dX(t) = -\frac{1}{2} \sum_{d=1}^{D} (\nabla_{E_d} E_d)(X(t)) dt + \sum_{d=1}^{D} E_d(X(t)) \circ dW_d(t),$$

where $E_d$ is an orthonormal frame of $\mathcal{M}$. In the specific case of a Lie group equipped with a bi-invariant metric, we find $\sum_d \nabla_{E_d} E_d = 0$ and a second order integrator is

$$X_{n+1} = \exp \left( \sum_{d=1}^{D} \left( \left( \frac{\sqrt{2}}{2} - 1 \right) \sqrt{h} \xi_n^{d,1} + \frac{\sqrt{2}}{2} \sqrt{h} \xi_n^{d,2} \right) E_d \right) \exp \left( \sum_{d=1}^{D} \sqrt{h} \xi_n^{d,1} E_d \right) X_n. \quad (2.12)$$

# 3 Order theory of stochastic frozen flow methods

This section is devoted to the formalisation of the order theory of stochastic frozen flow methods with the help of the new algebraic formalism of exotic Lie series governing the Taylor expansions of the exact flow and the frozen flow integrators.

## 3.1 Connection algebra of vector fields

The theory of the frozen flow methods relies on an additional geometric structure on $\mathcal{M}$. The idea is to define an affine connection $\triangleright$ on $TM$ whose geodesics are frozen flows. Such a connection is given by the choice of a right inverse to $(x, v) \in \mathcal{M} \times \mathbb{R}^D \mapsto v_d E_d(x) \in TM$, i.e. a choice of decomposing vector fields in the frame $(E_d)$. A convenient way to choose such an inverse is to extract a basis $(E_d)_{d \in I}$ from the frame. Then, the connection $\triangleright$ is uniquely defined by

$$\triangleright E_d = 0, \quad d \in I.$$

The extraction of a basis is licit only locally. However, for a given initial value, the new integrators (with bounded random variables) evolve in one step in an open neighbourhood $U$ of $\mathcal{M}$.
Taking the stepsize small enough, we can extract a frame basis from the frame $(E_d)$. Following Remark 2.10, we do not lose any generality in assuming that $(E_d)$ is a frame basis: the methods do not depend on the choice of a basis, but the Taylor expansion in the frame bundle does. We shall write our statements globally on $\mathcal{M}$ for the sake of simplicity, and all of them are valid in a neighbourhood of the starting point of the method.

The new connection can be seen as a product $\triangleright: \mathcal{X}(\mathcal{M}) \times \mathcal{X}(\mathcal{M}) \to \mathcal{X}(\mathcal{M})$ given by

$$F \triangleright G = \sum_{d=1}^{D} F[g^d]E_d, \quad G = \sum_{d=1}^{D} g^d E_d.$$

The product is well-posed as there is a unique decomposition of $G$ in the frame basis. The connection $\triangleright$ can be seen as the Weitzenböck connection w.r.t. the frame basis. We extend the notation to test functions: $F \triangleright \phi = F[\phi]$ (independent of the choice of connection).

**Remark 3.1.** Given a vector field $G$, the frozen flow $t \mapsto \exp(tG(p))p$ is the geodesic for the connection $\triangleright$ starting at $p$ with direction $G(p)$. Moreover, the parallel transport of $G(p)$ along any curve going from $p$ to $q$ is given by the frozen vector field $G_p(q)$.

This new connection $\triangleright$ defines a geometric structure far from the original structure associated to the metric and the Levi-Civita connection $\nabla$. Indeed, $\triangleright$ is neither compatible with the metric nor torsion-free. Let $[[F, G]]$ denote the Jacobi bracket of vector fields. The torsion $T$ of the connection $\triangleright$ is defined as

$$T(F, G) = F \triangleright G - G \triangleright F - [[F, G]].$$

This tensor is entirely determined by $T(E_i, E_j) = -[[E_i, E_j]]$, and does not vanish in general. An additional difference between $\triangleright$ and $\nabla$ concerns the curvature. As for a general connection, the curvature tensor $R$ is defined as

$$R(F, G)H = F \triangleright (G \triangleright H) - G \triangleright (F \triangleright H) - [[F, G]] \triangleright H,$$

and it vanishes in our context, i.e. $R = 0$. The flatness of $\triangleright$ explains the path-independence of parallel transport, raised in Remark 3.1.

**Remark 3.2.** In the deterministic context, the definition and analysis of Lie group methods were naturally tied to Lie group structures. For example, the connection $\triangleright$ naturally arise as the (-) Cartan connection associated to left/right-invariant vector fields. In this context, the torsion is constant and $-T(-, -)$ defines a Lie bracket, which naturally yields a post-Lie algebra structure [34, 39]. We emphasize that the frozen-flow methods, as well as their analysis, apply on any smooth manifold, not just on Lie groups.

Let us now consider the associative non-commutative concatenation product $\cdot$ as the differential operator

$$(F_n \cdots F_1) \triangleright \phi = \sum_{i_1, \dots, i_n} f_n^{i_1} \cdots f_1^{i_n} E_{i_n} [\dots E_{i_1}[\phi] \cdots],$$

defined on the tensor algebra $(T(\mathcal{X}(\mathcal{M})), \cdot)$. This product is called the frozen composition of differential operators.

We shall see that the Taylor expansions of the numerical and exact flows for the equation (1.1) typically write in $T(\mathcal{X}(\mathcal{M}))$. The number and complexity of the differential operators involved in the Taylor expansions of numerical and exact flows motivate the introduction of decorated and exotic planar forests.
## 3.2 Decorated planar forests

We introduce decorated planar forests and decorated Lie series.

**Definition 3.3.** A decorated planar tree $\tau$ is a connected, oriented graph in which every node has an outgoing edge, except for a node called the root. It is equipped with a map $\alpha: V(\tau) \to D$ from the set of nodes to the set of decorations. An ordered, possibly empty, collection of decorated planar trees is called a decorated planar forest. The sets of decorated planar trees and forests are denoted by $T_D$ and $F_D$, respectively.

The vector spaces spanned by $T_D$ and $F_D$ are denoted by $\mathcal{T}_D$ and $\mathcal{F}_D$. Decorated planar trees $\mathcal{T}_D$ are used as a convenient formalism for representing vector fields $\mathcal{X}(\mathcal{M})$ that appear in a Taylor expansion of a numerical or exact solution of a differential equation, while decorated planar forests $\mathcal{F}_D$ are used to represent differential operators $T(\mathcal{X}(\mathcal{M}))$ that appear in a Taylor expansion of a test function applied to a numerical or exact solution. We note that $\mathcal{F}_D = T(\mathcal{T}_D)$.

We use decorated planar forests $\mathcal{F}_D$ with $D = \{\bullet, \circ^1, \dots, \circ^L\}$ to represent the Taylor expansion of a test function applied to the one-step of a frozen flow method introduced in Section 2. We consider a map $\mathbb{F}^{\text{dec}}: \mathcal{F}_D \to T(\mathcal{X}(\mathcal{M}))$ defined as follows.

**Definition 3.4.** Let the map $\mathbb{F}^{\text{dec}}: \mathcal{F}_D \to T(\mathcal{X}(\mathcal{M}))$ be defined for $\pi \in \mathcal{F}_D$ and $\phi \in \mathcal{C}^{\infty}(\mathcal{M})$ as

$$ \mathbb{F}^{\text{dec}}(\pi) \triangleright \phi := \sum_{\substack{i_v=1 \\ v \in V(\pi)}}^{D} \left( \prod_{v \in V} E_{I_{\Pi}(v)}[\mathbb{F}^{\text{dec}}(v)^{i_v}] \right) E_{I_R}[\phi], $$

with $\mathbb{F}^{\text{dec}}(\bullet) = f$ and $\mathbb{F}^{\text{dec}}(\circ^l) = \xi^{d,l} E_d$ where $R$ is the set of roots of $\pi$, $V$ is the set of nodes of $\pi$, $\Pi(v)$ is the set of predecessors of $v$ ordered from right to left, $E_I[\phi] = E_{i_p}[\dots E_{i_1}[\phi]\dots]$ for $I = (i_p, \dots, i_1)$. The notation can also be frozen for $x, p \in \mathcal{M}$,

$$ (\mathbb{F}_x^{\text{dec}}(\pi) \triangleright \phi)(p) := \sum_{\substack{i_v=1 \\ v \in V(\pi)}}^{D} \left( \prod_{v \in V} E_{I_{\Pi}(v)}[\mathbb{F}^{\text{dec}}(v)^{i_v}] (x) \right) E_{I_R}[\phi](p), $$

For example, $\mathbb{F}^{\text{dec}}(??) \triangleright \phi = \xi^{1,i} E_i[f^j] f^k E_j[E_k[\phi]]$. Decorated Lie series are formal series of differential operators indexed by decorated planar forests. They do not converge in general, but their truncations at any order belong to $T(\mathcal{X}(\mathcal{M}))$.

**Definition 3.5.** Given a one form $a \in \mathcal{F}_D^*$, called a coefficient map, a decorated Lie S-series is the following formal power series in h in $T(\mathcal{X}(\mathcal{M}))$,

$$ S_h^{\text{dec}}(a) \triangleright \phi := \sum_{\pi \in F_D} h^{\lvert \pi \rvert} a(\pi) \mathbb{F}^{\text{dec}}(\pi) \triangleright \phi, \quad S_1^{\text{dec}}(a) = S_1^{\text{dec}}(a). $$

A decorated Lie series is only indexed on decorated trees and is the formal vector field,

$$ B_h^{\text{dec}}(a) := \sum_{\tau \in T_D} h^{\lvert \tau \rvert} a(\tau) \mathbb{F}^{\text{dec}}(\tau), \quad B_1^{\text{dec}}(a) = B_1^{\text{dec}}(a). $$

These series can also be frozen:

$$ S_{h,p}^{\text{dec}}(a) \triangleright \phi := \sum_{\pi \in F_D} h^{\lvert \pi \rvert} a(\pi) \mathbb{F}_p^{\text{dec}}(\pi) \triangleright \phi, \quad B_{h,p}^{\text{dec}}(a) := \sum_{\tau \in T_D} h^{\lvert \tau \rvert} a(\tau) \mathbb{F}_p^{\text{dec}}(\tau). $$
We endow the space $\mathcal{F}_D$ with the shuffle product $\sqcup$ defined inductively for $\tau, \gamma \in \mathcal{T}_D$ and $\pi, \eta \in \mathcal{F}_D$ as

$$ \tau\pi \sqcup \gamma\eta = \tau(\pi \sqcup \gamma\eta) + \gamma(\tau\pi \sqcup \eta). $$

For example,

$$ \begin{aligned} ?? \sqcup ?? &= ?? + ?? + ??, \\ ?? \sqcup ?? &= ?? + ?? + ?? + ?? + ?? + ?? + ?? \end{aligned} $$

A case of particular interest is when the coefficient map $a$ is a character over $(\mathcal{F}_D, \sqcup)$, that is, when $a(\pi_1 \sqcup \pi_2) = a(\pi_1)a(\pi_2)$. This condition is satisfied by the frozen flow, as expressed in Proposition 3.6.

**Proposition 3.6.** Let $\tilde{a} \in \mathcal{T}_D^*$, the frozen flow associated to the vector field $B^{\text{dec}}(\tilde{a})$ satisfies

$$ \phi(\exp(B_p^{\text{dec}}(\tilde{a}))) = \exp(B^{\text{dec}}(\tilde{a})) \triangleright \phi = S^{\text{dec}}(a) \triangleright \phi, $$

where $a \in \mathcal{F}_D^*$ extends $\tilde{a}$ to a character of $(\mathcal{F}_D, \sqcup)$, that is, $a(\pi_1 \sqcup \pi_2) = a(\pi_1)a(\pi_2)$. Moreover,

$$ a(\tau_1 \cdots \tau_n) = \frac{1}{n!} \tilde{a}(\tau_1) \cdots \tilde{a}(\tau_n). \qquad (3.1) $$

*Proof.* The proof is analogous to the proof of Lemma 2.4 from [80]. The fact that $a \in \mathcal{F}_D^*$ is a character with respect to the shuffle product follows from the formula (3.1) using the fact that there are $\binom{n+m}{n}$ terms in $\pi \sqcup \eta$ with $\pi$ and $\eta$ having $n$ and $m$ trees, respectively. $\square$

**Remark 3.7.** Representing Taylor expansions of vector fields with formal series indexed by tree-like structures is very natural from the geometric point of view. In [66, 74, 53], it is shown that exotic, aromatic and standard B-series are completely characterised by universal geometric properties of equivariance and locality. The extension of these properties to the planar case are open questions, even for standard Lie-Butcher series.

## 3.3 Exotic planar forests

This section adds planarity to the exotic forests from [54, 14, 15]. We use decorated planar forests to define exotic planar forests in Definition 3.8 which are the main object of interest for representing weak stochastic Taylor expansions. Exotic planar forests are decorated planar forests in which the decoration is used to denote the pairing between certain leaves. Let $S_N$ denote the group of permutations of natural numbers.

**Definition 3.8.** An exotic forest is a decorated forest $(\pi, \alpha) \in F_D$ with the decorations $D = \{\bullet\} \cup \mathbb{N}$, $\mathbb{N} = \{1, 2, 3, \dots\}$, that follows the following rules. If a natural number is used as a decoration, then it must decorate exactly two leaves, that is, $|\alpha^{-1}(n)| \in \{0, 2\}$ for all $n \in \mathbb{N}$. Two exotic forests $(\pi_1, \alpha_1)$ and $(\pi_2, \alpha_2)$ are considered to be identical if $\pi_1 = \pi_2$ and there exists a map $\varphi : D \to D$ with $\varphi|_{\mathbb{N}} \in S_N$ and $\varphi(\bullet) = \bullet$ such that $\alpha_1 = \varphi \circ \alpha_2$. The set of numbered leaves that correspond to the same number, that is $\alpha^{-1}(n)$ for $n \in \mathbb{N}$, is called a liana.

We gather the numbered nodes in $V_\circ$ and the black nodes in $V_\bullet$. The order of an exotic forest $\pi$ is then defined as

$$ |\pi| = |V_\bullet| + \frac{|V_\circ|}{2}. $$
The set of exotic forests is denoted by $EF$. The set of all exotic forests of order $p$ is $EF^p$, and the ones of order less than or equal to $p$ is $EF^{\le p}$. Exotic forests with one root are called exotic trees and form a set denoted $ET$. The corresponding vector spaces are denoted by $\mathcal{EF}$ and $\mathcal{ET}$, and similarly for the graded subspaces $\mathcal{EF}^p$ and $\mathcal{EF}^{\le p}$.

The exotic trees of order up to three are the following. Note that the order does not coincide with the number of black nodes in general.

??, ??, ??, ??, ??, ??, ??, ??, ??, ??, ??, ??, ??, ??.

We also note that the choice of the number used to denote the liana does not matter, for example, the following exotic forests are identical:

$$?? = ??.$$

Our main focus with exotic planar forests is to represent some specific elementary differentials in $\mathcal{T}(\mathcal{X}(M))$ and to translate the connection algebra to a natural algebraic structure on $\mathcal{EF}$.

The vector space $\mathcal{EF}$ of exotic forests forms an algebra $(\mathcal{EF}, \cdot)$ w.r.t. the non-commutative concatenation product. We emphasize that it does not coincide with the tensor algebra $T^*(\mathcal{ET})$ as for instance ?? cannot be expressed as the concatenation of exotic trees.

We introduce the new formalism of exotic Lie series for the study of order conditions of stochastic frozen flow methods. The formalism mixes the features of the exotic formalism [54, 50, 14, 15] for the stochastic part and the planar forests formalism [46, 80, 78] for Lie-group methods.

**Definition 3.9.** Let the map $\mathbb{F}: \mathcal{EF} \to T(\mathcal{X}(M))$ for $\pi \in EF$ and $\phi \in C^{\infty}(M)$ be defined as

$$\mathbb{F}(\pi) \triangleright \phi = \sum_{i_v=1, \dots, D} \delta_{i_L} \left( \prod_{v \in V_{\bullet}} E_{I_{\Pi(v)}} [f^{i_v}] \right) E_{I_R}[\phi],$$

where $R$ is the set of roots, $\Pi(v)$ is the set of predecessors of $v$ ordered from right to left, $E_I = E_{i_p} \dots E_{i_1}$ for $I = i_p, \dots, i_1$, and $\delta_{i_L} = 0$ if there exists $i_v \neq i_w$ with $\alpha(v) = \alpha(w)$ (that is, $(v, w)$ is a liana), and is 1 otherwise. The notation can also be frozen for $x, p \in M$:

$$(\mathbb{F}_x(\pi) \triangleright \phi)(p) = \sum_{i_v=1, \dots, D} \delta_{i_L} \left( \prod_{v \in V_{\bullet}} E_{I_{\Pi(v)}} [f^{i_v}](x) \right) E_{I_R}[\phi](p).$$

Exotic Lie S-series and exotic Lie-Butcher series are defined analogously to Definition 3.5: given a one form $a \in \mathcal{EF}^*$ (resp. $a \in \mathcal{ET}^*$), an exotic Lie S-series (resp. exotic Lie-Butcher series) is the following formal series

$$S_h(a) \triangleright \phi := \sum_{\pi \in EF} h^{\lvert \pi \rvert} a(\pi) \mathbb{F}(\pi) \triangleright \phi, \quad B_h(a) := \sum_{\tau \in ET} h^{\lvert \tau \rvert} a(\tau) \mathbb{F}(\tau),$$

and analogously for their frozen counterparts $S_{h,p}$ and $B_{h,p}$.

The exotic series arise naturally by considering the expectation of decorated series. Let $\Phi: \mathcal{EF} \to \mathcal{FD}$ be the map that forgets the pairings of the white leaves and sends each liana to a sum of paired decorated nodes $\sum_{l=1}^L \circ^{l_1} \circ^{l_2} \dots$. For example,

$$\Phi(??) = \sum_{i,j,k=1}^{L} ??, \quad \Phi(??) = \sum_{i,j,k=1}^{L} ??, \quad \Phi(??) = \sum_{i,j,k=1}^{L} ??.$$
**Proposition 3.10.** Let $a^{\text{dec}}$ be a character on $(\mathcal{F}_D, \mathbb{1})$ with $D = \{\bullet, \circ^1, \dots, \circ^L\}$, then

$$ \mathbb{E}[S^{\text{dec}}(a^{\text{dec}}) \triangleright \phi] = S(a) \triangleright \phi, $$

where $a = a^{\text{dec}} \circ \Phi$ is a character on $\mathcal{E}\mathcal{F}$ with respect to the shuffle product

For example, one finds

$$ \mathbb{E}[\mathbb{F}^{\text{dec}}(??)] = \mathbb{F}(?? + ?? + ??). $$

*Proof.* The proof is analogous to [14, Sect. 4.1] where the white leaves decorated by the same number are paired using the Isserlis formula [48]. That is,

$$ \mathbb{E}[S^{\text{dec}}(a^{\text{dec}}) \triangleright \phi] = \mathbb{E}\left[\sum_{\pi \in F_D} h^{\lvert \pi \rvert} a^{\text{dec}}(\pi) \mathbb{F}^{\text{dec}}(\pi) \triangleright \phi\right] = \sum_{\pi \in F_D} h^{\lvert \pi \rvert} a^{\text{dec}}(\pi) \mathbb{E}[\mathbb{F}^{\text{dec}}(\pi)] \triangleright \phi. $$

We use Isserlis formula to write

$$ \mathbb{E}[\mathbb{F}^{\text{dec}}(\pi, \alpha)] = \sum_{\alpha_e} \mathbb{F}(\pi, \alpha_e), $$

where $(\pi, \alpha)$ and $(\pi, \alpha_e)$ denote decorated and exotic planar forests as defined in Definitions 3.3 and 3.8 with the sum being over all $\alpha_e$ such that $\alpha_e^{-1}(\bullet) = \alpha^{-1}(\bullet)$ and $\alpha_e^{-1}(\mathbb{N})$ is a union of partitions into pairs of $\alpha^{-1}(\circ^l)$ for $l = 1, \dots, L$. We note that $\mathbb{E}[\mathbb{F}^{\text{dec}}(\pi)] = 0$ if $\pi$ contains an odd number of nodes of color $\circ^l$ for an $l \in \{1, \dots, L\}$. Therefore, we have,

$$ \mathbb{E}[S^{\text{dec}}(a^{\text{dec}}) \triangleright \phi] = \sum_{\pi \in EF} h^{\lvert \pi \rvert} a^{\text{dec}}(\Phi(\pi)) \mathbb{F}(\pi) \triangleright \phi = S(a^{\text{dec}} \circ \Phi) \triangleright \phi, $$

and the statement is proved. $\square$

**Example 3.11.** Consider the Euler frozen flow method (2.9). Using Proposition 3.6, its Taylor expansion in decorated forests is given by

$$ \begin{aligned} \phi(\exp(hF_p + \sqrt{2h}\xi^d E_d)p) &= \exp(hF + \sqrt{2h}\xi^d E_d) \triangleright \phi(p) = \mathbb{F}^{\text{dec}}(\exp(h?? + \sqrt{2h}??)) \triangleright \phi(p) \\ &= \left( \text{id} + h^{1/2}\sqrt{2}\xi^d E_d + h(F + \xi^{d_2}\xi^{d_1}E_{d_2} \cdot E_{d_1}) \right. \\ &\qquad \left. + h^{3/2}\left(\frac{1}{\sqrt{2}}\xi^d F \cdot E_d + \frac{1}{\sqrt{2}}\xi^d E_d \cdot F + \frac{\sqrt{2}}{3}\xi^{d_3}\xi^{d_2}\xi^{d_1}E_{d_3} \cdot E_{d_2} \cdot E_{d_1}\right) + \dots \right) \triangleright \phi \\ &= \mathbb{F}^{\text{dec}}\left( \mathbf{1} + h^{1/2}\sqrt{2}?? + h(?? + ??) + h^{3/2}\left(\frac{1}{\sqrt{2}}?? + \frac{1}{\sqrt{2}}?? + \frac{\sqrt{2}}{3}??\right) \right. \\ &\qquad \left. + h^2\left(\frac{1}{2}?? + \frac{1}{3}?? + \frac{1}{3}?? + \frac{1}{3}?? + \frac{1}{6}?? + \frac{1}{6}?? + \frac{\sqrt{2}}{6}??\right) + \dots \right) \triangleright \phi(p). \end{aligned} $$

Then, the expectation pairs the decorated nodes together as detailed in Proposition 3.10 and yields

$$ \begin{aligned} \mathbb{E}[\phi(\exp(hF_p + \sqrt{2h}\xi^d E_d)p)] &= \mathbb{F}\left(\mathbf{1} + h(?? + ??) \right. \\ &\qquad \left. + h^2\left(\frac{1}{2}?? + \frac{1}{3}?? + \frac{1}{3}?? + \frac{1}{3}?? + \frac{1}{6}?? + \frac{1}{6}?? + \frac{1}{6}??\right) + \dots \right) \triangleright \phi(p) \\ &= \left( \text{id} + h(F + E_d \cdot E_d) + h^2\left(\frac{1}{2}F \cdot F + \frac{1}{3}F \cdot E_d \cdot E_d + \frac{1}{3}E_d \cdot F \cdot E_d + \frac{1}{3}E_d \cdot E_d \cdot F \right. \right. \\ &\qquad \left. + \frac{1}{6}E_{d_2} \cdot E_{d_2} \cdot E_{d_1} \cdot E_{d_1} + \frac{1}{6}E_{d_2} \cdot E_{d_1} \cdot E_{d_2} \cdot E_{d_1} + \frac{1}{6}E_{d_1} \cdot E_{d_2} \cdot E_{d_2} \cdot E_{d_1} + \dots \right) \triangleright \phi(p). \end{aligned} $$
## 3.4 Conversion of geometric operations into exotic Lie series

We saw in Proposition 3.6 that the frozen flow rewrites as a tree series operation. The following result, proved analogously to [80, Lemma 2.3], provides an algebraic description of the Taylor expansion of a frozen vector field.

**Lemma 3.12.** Let a smooth map $\varphi: \mathcal{M} \to \mathcal{M}$ with expansion $\phi(\varphi(p)) = (S(a) \triangleright \phi)(p)$ with $a \in \mathcal{E}\mathcal{F}^*$, then the frozen vector field $F_\varphi: p \to F_{\varphi(p)}(p)$ satisfies

$$F_{\varphi} \triangleright \phi = (S(a) \triangleright hF) \triangleright \phi = B(\tilde{a}) \triangleright \phi, \quad \tilde{a}(\tau) = a(B^{-}(\tau)), \quad \tilde{a} \in \mathcal{E}\mathcal{T}^*.$$

Let two smooth maps $\varphi^1, \varphi^2: \mathcal{M} \times \mathcal{M} \to \mathcal{M}$, typically $\varphi_x^k(p) = \exp(F_{\psi^k(x)}) \cdot p$ with $\psi^k: \mathcal{M} \to \mathcal{M}$. We use two possible ways to compose the $\varphi^k$. The standard composition of the maps would read

$$\varphi^1 \circ \varphi^2(p) = \varphi_{\varphi_p^2(p)}^1(\varphi_p^2(p)).$$

On the other hand, the frozen flow methods (2.8) rely on the frozen composition:

$$\varphi^1 \cdot \varphi^2(p) = \varphi_p^1 \cdot \varphi_p^2(p) = \varphi_p^1(\varphi_p^2(p)).$$

The decorated and exotic Lie series allow to represent naturally all the geometric operations needed for the analysis of the frozen-flow methods (2.8). We present the following result in the context of exotic series (with $h = 1$ for simplicity), and we mention that decorated series satisfy analogous identities. Let the de concatenation product $\Delta: \mathcal{E}\mathcal{F} \to \mathcal{E}\mathcal{F} \otimes \mathcal{E}\mathcal{F}$ be defined as

$$\Delta(\pi_{1,...,n}) = \sum_{k=0}^{n} \pi_{1,...,k} \otimes \pi_{k+1,...,n}, \quad \text{with } \pi_i \in \text{Irr},$$

where $\pi_{1,...,n} = \pi_1 \cdots \pi_n$, $\pi_\emptyset = \mathbf{1}$, and Irr is the space of irreducible exotic forests introduced in Section 4.3. The composition and the Munthe-Kaas-Wright coproduct are further discussed in Section 4.

**Theorem 3.13.** Let two smooth maps $\varphi^1, \varphi^2: \mathcal{M} \times \mathcal{M} \to \mathcal{M}$ with expansions $\phi(\varphi_x^k(p)) = (S_x(a^k) \triangleright \phi)(p)$ with $a \in \mathcal{E}\mathcal{F}^*$. The frozen composition is given by

$$\phi(\varphi^1 \cdot \varphi^2) = (S(a^2) \cdot S(a^1)) \triangleright \phi = S(a^2 \cdot a^1) \triangleright \phi, \quad a^2 \cdot a^1 = \mu \circ (a^2 \otimes a^1) \circ \Delta.$$

Let two smooth maps $\varphi^1, \varphi^2: \mathcal{M} \to \mathcal{M}$ with expansions $\phi(\varphi^k(p)) = (S(a^k) \triangleright \phi)(p)$ with $a^k \in \mathcal{E}\mathcal{F}^*$. The composition of exotic series as differential operators is given by the composition law, dual of the MKW coproduct,

$$\phi(\varphi^1 \circ \varphi^2) = S(a^2) \triangleright (S(a^1) \triangleright \phi) = S(a^2 * a^1) \triangleright \phi, \quad a^2 * a^1 = \mu \circ (a^2 \otimes a^1) \circ \Delta_{MKW}.$$

*Proof.* The expression for the frozen composition is derived similarly to [80, Lemma 2.2]. The composition of flows derives from Theorem 4.6. $\square$

## 3.5 Exotic Lie series of the exact and numerical flows

The exotic Lie series formalism allows us to give explicit expressions of the Taylor expansions of both exact and numerical flows. The generator of equation (1.1) rewrites as

$$\mathcal{L}[\phi] = \sum_{d=1}^{D} f^d E_d[\phi] + \sum_{d=1}^{D} E_d[E_d[\phi]] = \mathbb{F}(?? + ??)[\phi] = S(l)[\phi], \quad l = \delta?? + \delta??.$$

The formal expansion (2.3) then rewrites straightforwardly in terms of forests. The proof is a generalisation of the works [81, 84, 14].
**Theorem 3.14.** The Taylor expansion of the exact flow of (1.1) is given by the exotic Lie series

$$u(h,p) = \exp(h\mathcal{L})[\phi](p) = S_h(e) \triangleright \phi(p), \quad e = \exp^*(l) = \sum_{n=0}^{\infty} \frac{l^{*n}}{n!}.$$

In addition, let $\alpha(\pi)$ be the number of different ways the exotic planar forest $\pi$ can be obtained by successive applications of the Grossman-Larson product of a black node, or of a pair of decorated nodes. Then, the coefficient $e$ is given by

$$e(\pi) = \frac{\alpha(\pi)}{|\pi|!}.$$

The first terms are the following

$$\exp(t\mathcal{L}) = \mathbb{F}(1 + t(?? + ??) + \frac{t^2}{2}(?? + ?? + ?? + 2?? + ?? + ?? + ??) + \dots),$$

where the second order term is $\frac{t^2}{2}\mathcal{L}^2 = \frac{t^2}{2}S_t(l) * S_t(l) = \frac{t^2}{2}S_t(l * l)$.

Let us now present the expansion in exotic series of the frozen flow integrators. For the sake of clarity, the random Runge-Kutta coefficients $z_k^d, z_{i,k}^d$ are chosen as centered Gaussians. They can be replaced with bounded approximations of Gaussians with the same first moments up to order $2p$ and the expansion of the integrator will be an exotic series up to order $p$.

**Definition 3.15.** For $K > 0$, a labelling in $Lab_K(\pi)$ of a decorated forest $\pi$ is a multi-index map $k: V \to \{1, \dots, K\}$ such that $k_v \ge k_w$ if $v$ and $w$ are in the same tree, with same height and $v$ is on the right of $w$.

For instance, a possible labelling in $Lab_8(??)$ is $??$. Then, the Talay-Tubaro expansion of the new methods (2.8) is also given by an exotic Lie series.

**Theorem 3.16.** The Talay-Tubaro expansion of a frozen flow method of the form (2.8) is an exotic Lie series:

$$\mathbb{E}[\phi(X_1)|X_0 = p] = S_h(a) \triangleright \phi(p),$$

where a is given by

$$a(\pi) = \mathbb{E} \left[ \sum_{i_v=1, \dots, s} \sum_{k \in Lab_K(\pi)} \sum_{v \in V_o} \delta_{d_L} \prod_{(v,w) \in E} Z_{i_w, i_v, k_v}^0 \prod_{v \in V_o} Z_{i_v, k_v}^{d_v} \prod_{r \in R \cap V_o} z_{i_r, k_r}^0 \prod_{r \in R \cap V_o} z_{k_r}^{d_r} \right],$$

where $\delta_{d_L} = 0$ if there exists $d_v \ne d_w$ with $(v, w) \in L$, and is 1 else. Moreover, the coefficient map a is a character on $(\mathcal{E}\mathcal{F}, \sqcup)$, that is, $a(\pi_1 \sqcup \pi_2) = a(\pi_1)a(\pi_2)$.

The proof is a stochastic generalisation of [80].

*Proof.* Let us consider one step of the algorithm, with one Gaussian for the sake of simplicity: $Z_{i,k}^d = \hat{Z}_{i,k}\xi^d, z_k^d = \hat{z}_k\xi^d$. Let the intermediate steps

$$g_k^i = \exp \left( h \sum_{j=1}^{s} Z_{i,j,k}^0 f(H^j) + \sqrt{h} \hat{Z}_{i,k} \sum_{d} \xi_d E_d \right) \circ g_{k-1}^i, \quad g_0^i = \text{id}.$$
Theorem 3.13 yields $\phi(g_k^i) = S_h^{\text{dec}}(\alpha_k^i) \triangleright \phi$ and $\phi(H^i) = S_h^{\text{dec}}(\alpha_K^i) \triangleright \phi$ with $\alpha_0^i = \delta_1$, $\alpha_k^i(1) = 1$, and

$$ \alpha_k^i(\tau_1 \dots \tau_n) = \sum_{p=0}^{n} \frac{1}{(n-p)!} \alpha_{k-1}^i(\tau_1 \dots \tau_p) \left( \sum_{j_{p+1}=1}^{s} Z_{i,j_{p+1},k}^0 \alpha_K^{j_{p+1}}(B^{-\tau_{p+1}}) \mathbf{1}_{\tau_{p+1} \neq ??} + \hat{Z}_{i,k} \mathbf{1}_{\tau_{p+1} = ??} \right) \dots \left( \sum_{j_n=1}^{s} Z_{i,j_n,k}^0 \alpha_K^{j_n}(B^{-\tau_n}) \mathbf{1}_{\tau_n \neq ??} + \hat{Z}_{i,k} \mathbf{1}_{\tau_n = ??} \right). $$

An induction then yields

$$ \alpha_K^i(\tau_1 \dots \tau_n) = \sum_{k_1 \le \dots \le k_n} \frac{1}{k!} \left( \sum_{j_1=1}^s Z_{i,j_1,k_1}^0 \alpha_K^{j_1}(B^{-\tau_1}) \mathbf{1}_{\tau_1 \ne ??} + \hat{Z}_{i,k_1} \mathbf{1}_{\tau_1 = ??} \right) \dots \left( \sum_{j_n=1}^s Z_{i,j_n,k_n}^0 \alpha_K^{j_n}(B^{-\tau_n}) \mathbf{1}_{\tau_n \ne ??} + \hat{Z}_{i,k_n} \mathbf{1}_{\tau_n = ??} \right). $$

Similar reasoning yields the desired expansion for the last stage of the method as $S_h^{\text{dec}}(a^{\text{dec}})$, where $a^{\text{dec}}$ is a character by construction. Proposition 3.10 yields the expected formula for the coefficient map $a$, as well as the character property. $\square$

Examples of computations of the coefficient map $a$ associated to frozen flow methods are presented in Table 1 and Appendix A. A more involved example of order 7 is the following:

$\pi = ??$

$$ a(\pi) = \mathbb{E} \sum_{i_1, \dots, i_9=1}^{s} \sum_{\substack{k_1 \ge k_2 \ge k_3 \\ k_4 \ge k_5, k_6 \ge k_7 \ge k_8}}^{!} \sum_{d_1, d_2=1}^{D} Z_{i_8, i_9, k_9}^{d_2} Z_{i_3, i_8, k_8}^{0} Z_{i_3, i_7, k_7}^{d_1} Z_{i_3, i_6, k_6}^{d_1} Z_{i_1, i_5, k_5}^{0} Z_{i_1, i_4, k_4}^{0} z_{i_3, k_3}^{0} z_{k_2}^{d_2} z_{i_1, k_1}^{0} $$

**Remark 3.17.** Order conditions for order $p$ are obtained by requiring $(a-e)(\pi) = 0$ for all $\pi \in EF\setminus\{1\}$ of order $|\pi| \le p$. This results in the set of order conditions being indexed by the elements of $EF\setminus\{1\}$. The fact that the coefficient maps $a$ and $e$ are characters of $(\mathcal{E}\mathcal{F}, \sqcup)$ induces relations between order conditions as detailed in Remark 2.14. In particular, $(a-e)(\pi \sqcup \eta) = 0$ if $(a-e)(\pi) = 0$ and $(a-e)(\eta) = 0$ using,

$$ (a-e)(\pi \sqcup \eta) = a(\pi)a(\eta) - e(\pi)e(\eta) = 0. $$

To reduce the number of superfluous order conditions, we consider the quotient of $\mathcal{E}\mathcal{F}$ by the ideal $\mathcal{I} := \langle \pi \sqcup \eta : \pi, \eta \in \mathcal{E}\mathcal{F}\setminus\{1\} \rangle$. Using the fact that the order conditions corresponding to $\pi \sqcup \eta$ are satisfied automatically, we take $a-e$ to be a functional over $\mathcal{E}\mathcal{F}/\mathcal{I}$, that is $a-e \in (\mathcal{E}\mathcal{F}/\mathcal{I})^*$. Dualizing $\mathcal{E}\mathcal{F}/\mathcal{I}$ yields the space of primitive elements $\text{Prim}$ of $\mathcal{E}\mathcal{F}$ with respect to the deshuffle coproduct $\Delta_{\sqcup}$ defined as

$$ \text{Prim} := \text{Span}\{\pi \in \mathcal{E}\mathcal{F} : \Delta_{\sqcup}(\pi) = \pi \otimes \mathbf{1} + \mathbf{1} \otimes \pi\}. $$

Therefore, using the identification $(a-e)(\pi) = \langle a-e, \pi \rangle$, we have $a-e \in \text{Prim}$. The detailed description of $\text{Prim}$ is left for future work.
# 4 Exotic planar forests and their algebraic structure

In this section, we focus on the algebraic structures of exotic planar forests. We first recall the well-known D-algebra structure of the connection algebra of vector fields in Section 4.1 and decorated planar forests in Section 4.2. We then explore how these structures change in the context of exotic planar forests in Section 4.3. In particular, we show that exotic planar forests form both the Grossman-Larson and Munthe-Kaas-Wright Hopf algebras. We also observe that the D-algebra structure is preserved. However, exotic planar forests, unlike decorated planar forests, do not form a post-Hopf algebra, since a connected exotic forest can have multiple roots, which violates the Leibniz rule for the grafting product.

## 4.1 D-algebra of vector fields

We recall the connection algebra $\mathcal{X}(\mathcal{M})$ introduced in Section 3.1 with the connection $\triangleright$. The deshuffle coproduct $\Delta_{\triangleright}: T(\mathcal{X}(\mathcal{M})) \to T(\mathcal{X}(\mathcal{M})) \otimes T(\mathcal{X}(\mathcal{M}))$ is, for $A, B \in T(\mathcal{X}(\mathcal{M}))$, $F \in \mathcal{X}(\mathcal{M})$,

$$ \Delta_{\triangleright}(AB) = \Delta_{\triangleright}A \cdot \Delta_{\triangleright}B, \quad \Delta_{\triangleright}(F) = \mathbf{1} \otimes F + F \otimes \mathbf{1}, $$

where we recall $(A \otimes B) \cdot (C \otimes D) = (A \cdot C) \otimes (B \cdot D)$. Following the extension [34] of the so-called Guin-Oudom process [79], the product $\triangleright$ extends to $T(\mathcal{X}(\mathcal{M}))$, for $A, B, C \in T(\mathcal{X}(\mathcal{M}))$, $F, G \in \mathcal{X}(\mathcal{M})$, by

$$ \begin{align*} \mathbf{1} \triangleright A &= A, & A \triangleright \mathbf{1} &= 0, \\ (F \cdot A) \triangleright G &= F \triangleright (A \triangleright G) - (F \triangleright A) \triangleright G, \\ A \triangleright (B \cdot C) &= (A_{(1)} \triangleright B) \cdot (A_{(2)} \triangleright C), \end{align*} $$

where we use the Sweedler notation w.r.t. $\Delta_{\triangleright}$. We note that the extension of the Guin-Oudom process does not require $(\mathcal{X}, -T, \triangleright)$ to have a post-Lie structure. The triple $(T(\mathcal{X}(\mathcal{M})), \cdot, \triangleright)$ is called a D-algebra [78]. That is, $T(\mathcal{X}(\mathcal{M}))$ is a graded associative algebra and for all $F \in \mathcal{X}(\mathcal{M})$ and $A, B \in T(\mathcal{X}(\mathcal{M}))$, we have, $A \triangleright F \in \mathcal{X}(\mathcal{M})$, as well as,

$$ \begin{align*} F \triangleright (A \cdot B) &= (F \triangleright A) \cdot B + A \cdot (F \triangleright B), \\ F \triangleright (A \triangleright B) &= (F \triangleright A) \triangleright B + (F \cdot A) \triangleright B. \end{align*} $$

The tensor algebra is equipped with the associative non-commutative Grossman-Larson product

$$ A * B = A_{(1)} \cdot (A_{(2)} \triangleright B), \quad (A * B) \triangleright C = A \triangleright (B \triangleright C). $$

Then, $(T(\mathcal{X}(\mathcal{M})), \cdot, \Delta_{\triangleright})$ and $(T(\mathcal{X}(\mathcal{M})), \cdot, \Delta_{\triangleright})$ are Hopf algebras. By the Cartier-Quillen-Milnor-Moore theorem, the Hopf algebra $(T(\mathcal{X}(\mathcal{M})), \cdot, \Delta_{\triangleright})$ coincides with the universal enveloping algebra $U(\mathcal{X}(\mathcal{M}))$ of $(\mathcal{X}(\mathcal{M}), \llbracket - , - \rrbracket)$. The Grossman-Larson product satisfies, in particular, for $F, G \in \mathcal{X}(\mathcal{M})$:

$$ \llbracket F, G \rrbracket \triangleright \phi = (F * G - G * F) \triangleright \phi. $$

## 4.2 Decorated planar forests

We recall the discussion of the decorated planar forests in Section 3.2. We present an alternative definition of the decorated planar forests in Definition 4.1.
**Definition 4.1.** Let $\tau_1, \dots, \tau_n \in T_D$ be decorated planar trees, then, a decorated planar forest $\pi \in F_D$ is a concatenation $\tau_1 \cdots \tau_n$. Let $\pi \in F_D$ be a decorated planar forest, then, a decorated planar tree $\tau \in T_D$ is defined as $B_d^+(\pi)$ where the map $B_d^+ : F_D \to T_D$ attaches all roots of $\pi$ to a new root decorated by $d$.

We define the left grafting $\rightsquigarrow: \mathcal{T}_D \otimes \mathcal{T}_D \to \mathcal{T}_D$, for $\tau, \gamma_1, \dots, \gamma_n \in T_D$, as

$$ \tau \rightsquigarrow B_d^+(\mathbf{1}) := B_d^+(\tau), \quad \tau \rightsquigarrow B_d^+(\gamma_1 \cdots \gamma_n) := B_d^+(\tau \cdot \gamma_1 \cdots \gamma_n + \sum_{i=1}^n \gamma_1 \cdots (\tau \rightsquigarrow \gamma_i) \cdots \gamma_n). $$

It is extended to $\rightsquigarrow: \mathcal{F}_D \otimes \mathcal{F}_D \to \mathcal{F}_D$, for $\tau, \gamma \in T_D$, $\pi, \eta, \mu \in F_D$, by Guin-Oudom process,

$$ \begin{align*} \mathbf{1} \rightsquigarrow \pi &:= \pi, & \pi \rightsquigarrow \mathbf{1} &:= 0, \\ (\tau \cdot \pi) \rightsquigarrow \gamma &:= \tau \rightsquigarrow (\pi \rightsquigarrow \gamma) - (\tau \rightsquigarrow \pi) \rightsquigarrow \gamma, \\ \pi \rightsquigarrow (\eta \cdot \mu) &:= (\pi_{(1)} \rightsquigarrow \eta) \cdot (\pi_{(2)} \rightsquigarrow \mu). \end{align*} $$

We note that $(\mathcal{F}_D, \cdot, \rightsquigarrow)$ forms the free D-algebra over the set of decorations $D$, [78]. The universal property of $(\mathcal{F}_D, \cdot, \rightsquigarrow)$ allows us to define a D-algebra homomorphism $\mathbb{F}^{\text{dec}} : \mathcal{F}_D \to \mathcal{T}(\mathfrak{X}(\mathcal{M}))$ by $\mathbb{F}^{\text{dec}}(\mathbb{F}^{\text{dec}}(\mathcal{T})) = f^d$ satisfying

$$ \mathbb{F}^{\text{dec}}(\pi_1 \cdot \pi_2) = \mathbb{F}^{\text{dec}}(\pi_1)\mathbb{F}^{\text{dec}}(\pi_2), \quad \mathbb{F}^{\text{dec}}(\pi_1 \rightsquigarrow \pi_2) = \mathbb{F}^{\text{dec}}(\pi_1) \triangleright \mathbb{F}^{\text{dec}}(\pi_2). $$

**Remark 4.2.** The space of decorated planar forests forms the tensor algebra $T(\mathcal{T}_D)$ of the space of decorated planar trees. It is the universal enveloping Lie algebra $\mathcal{F}_D = \mathcal{U}(\text{Lie}(\mathcal{T}_D))$ of the post-Lie algebra $\text{Lie}(\mathcal{T}_D)$ with the Lie bracket $[x, y] = x \cdot y - y \cdot x$ for $x, y \in \text{Lie}(\mathcal{T}_D)$. Moreover, $(\mathcal{F}_D, \cdot, \rightsquigarrow, \Delta_\wedge, \rightsquigarrow)$ is a post-Hopf algebra [61] and $\text{Lie}(\mathcal{T}_D)$ is the space of primitive elements with respect to the deshuffle coproduct $\Delta_\wedge$, that is,

$$ \text{Lie}(\mathcal{T}_D) = \{\pi \in \mathcal{F}_D : \Delta_\wedge(\pi) = \pi \otimes \mathbf{1} + \mathbf{1} \otimes \pi\}. $$

We emphasize that $\mathbb{F}^{\text{dec}} : \mathcal{F}_D \to T(\mathfrak{X}(\mathcal{M}))$ does not restrict to a Lie algebra morphism $\text{Lie}(\mathcal{T}_D) \to \mathfrak{X}(\mathcal{M})$.

The Grossman-Larson product over decorated planar forests is defined as

$$ \pi \diamond \eta := \pi_{(1)} \cdot (\pi_{(2)} \rightsquigarrow \eta), \quad \text{for } \pi, \eta \in \mathcal{F}_D, $$

and $\mathbb{F}^{\text{dec}}(\pi \diamond \eta) = \mathbb{F}^{\text{dec}}(\pi) * \mathbb{F}^{\text{dec}}(\eta)$ due to the fact that $\mathbb{F}^{\text{dec}}$ is a D-algebra homomorphism. We generalize the concatenation and Grossman-Larson products by specifying the order in which the trees of the operands are concatenated.

Let $\pi, \eta \in \mathcal{F}_D$ and let $\omega$ denote an admissible ordering of the roots of $\pi \cdot \eta$. An ordering $\omega$ is *admissible* if the order of roots in $\pi$ and $\eta$ is preserved, that is, if $r_1$ and $r_2$ are roots of $\pi$ or $\eta$ and $r_1 \le r_2$, then $r_1 \le r_2$ according to $\omega$. We define $\cdot_\omega : \mathcal{F}_D \times \mathcal{F}_D \to \mathcal{F}_D$ to be the concatenation product with roots arranged according to the ordering $\omega$. The Grossman-Larson product $\diamond_\omega : \mathcal{F}_D \times \mathcal{F}_D \to \mathcal{F}_D$ is defined as follows,

$$ \pi \diamond_{\omega} \eta := \pi_{(1)} \cdot_{\omega} (\pi_{(2)} \rightsquigarrow \eta), $$

where the ordering $\omega$ is defined over the roots of $\pi \cdot \eta$. For example

$$ \text{??} \cdot_{\omega} \text{??} = \text{??}, \quad \text{for } \omega : c \le a \le b \le d, $$
$$ \begin{aligned} \text{?? } \diamond_{\omega} \text{??} = &+ \text{??} + (\text{??} \leftrightharpoons \text{??})\text{??} + \text{??}(\text{??} \leftrightharpoons \text{??}) \\ &+ (\text{??} \leftrightharpoons \text{??})\text{??} + \text{??}(\text{??} \leftrightharpoons \text{??}) + \text{??} \leftrightharpoons \text{??}. \end{aligned} $$

We note that $\pi \cdot_{\omega} \eta = \pi \cdot \eta$ and $\pi \diamond_{\omega} \eta = \pi \diamond \eta$ for $\omega$ which takes the roots of $\pi$ to be smaller than the roots of $\eta$. If the admissible ordering $\omega$ in $\pi \cdot_{\omega} \eta$ does not specify the order between some roots in $\pi \cdot \eta$, then, the roots of $\pi$ are taken to be smaller than the roots of $\eta$. The generalised concatenation and Grossman-Larson products are used in the context of planar exotic forests to describe the structure of the Grossman-Larson Hopf algebra.

## 4.3 Exotic planar forests

The algebraic structures described in this section are applied in Section 3 for the description of the composition law and the order conditions of stochastic frozen flow methods.

We construct the space $\mathcal{EF}$ from the space of decorated planar forest $\mathcal{F}_D$ and inherit a number of algebraic structures from it. There are two possible approaches: with an appropriate subspace or with a quotient. We use the first approach to define the algebras and the second approach to define coalgebras.

Let us consider decorated forests $\mathcal{F}_D$ with $D = \{??, ?? : d \in \mathbb{N}\}$ and let $\mathcal{F}_D^L$ denote its quotient by the subspace $\mathcal{K}$ spanned by the forests in which a numbered node is found at a position which is not a leaf. It can be checked that $\mathcal{K}$ is an ideal with respect to concatenation, shuffle, grafting, and Grossman-Larson products, and, therefore, $\mathcal{F}_D^L$ forms an algebra when endowed with one of these products. Let $\overline{\mathcal{F}}_D^L$ denote its completion with respect to the grading given by the number of nodes, that is, $\overline{\mathcal{F}}_D^L$ is a space of formal sums over $F_D^L$. We define the subspace $\mathcal{F}_D^\mathcal{E}$ of $\overline{\mathcal{F}}_D^L$, corresponding to exotic forests, spanned by

$$ \varphi(\pi, \alpha_e) := \sum_{\alpha \in P(\alpha_e)} (\pi, \alpha), $$

with $\varphi : \mathcal{EF} \to \mathcal{F}_D^\mathcal{E}$ an isomorphism. $P(\alpha_e)$ is the set of decorations $\alpha$ with $\alpha^{-1}(\bullet) = \alpha_e^{-1}(\bullet)$ and

$$ \alpha(v_1) = \alpha(v_2) \quad \text{if} \quad \alpha_e(v_1) = \alpha_e(v_2), \quad \text{for } v_1, v_2 \in V(\pi). $$

Some examples of the values of $\varphi$ are

$$ \varphi(??) = \sum_{i=1}^{\infty} ??, \quad \varphi(??) = \sum_{i,j=1}^{\infty} ?? $$

We take $\varphi : \mathcal{EF} \to \mathcal{F}_D^\mathcal{E}$ to be an algebra homomorphism with respect to the generalised concatenation, shuffle, grafting, and generalised Grossman-Larson products which induces the corresponding algebraic structures over $\mathcal{EF}$. For instance, the shuffle product $\sqcup : \mathcal{EF} \times \mathcal{EF} \to \mathcal{EF}$ permutes the roots of forests. Let $\pi \cdot_{\omega} \eta$ be the generalised concatenation obtained by concatenating the forests $\pi$ and $\eta$ with the roots permuted according to the admissible ordering $\omega$ of the roots of $\pi \cdot \eta$. The shuffle product is defined as

$$ \pi \sqcup \eta = \sum_{\omega} \pi \cdot_{\omega} \eta, $$

where the sum is over admissible orderings. For example, we find

$$ \text{?? } \sqcup \text{??} = \text{??} + \text{??} + \text{??}, $$
$$?? \sqcup ?? = ?? + ?? + ?? + ?? + ?? + ??,$$

with **1** being the neutral element for $\sqcup$.

Let us now consider the subspace of $\mathcal{F}_D^L$ orthogonal to $\mathcal{F}_D^\mathcal{E}$ with respect to the inner product with which $F_D^L$ forms an orthonormal basis of $\mathcal{F}_D^L$. We denote it by $\mathcal{I}$ and we can see that it is a coideal in the coalgebra $(\mathcal{F}_D^L, \Delta_\omega)$ using the fact that the deshuffle coproduct is adjoint to the shuffle product. Taking a quotient of $(\mathcal{F}_D^L, \Delta_\omega)$ by $\mathcal{I}$ gives a definition of the coalgebra of exotic forests $(\mathcal{E}\mathcal{F}, \Delta_\omega)$ through the isomorphism

$$\psi(\pi, \alpha_e) := (\pi, \alpha_e) + \mathcal{I} \in \mathcal{F}_D^L/\mathcal{I}. \qquad (4.1)$$

We refer to [15] for a precise analogous definition with non-planar exotic forests.

An exotic forests $\pi$ is irreducible if it cannot be written as a concatenation of two exotic forests, that is, $\pi \neq \gamma \cdot \eta$ for $\gamma, \eta \in EF\{\mathbf{1}\}$. The space of irreducible exotic forests is denoted by $\text{Irr} \subset \mathcal{E}\mathcal{F}$. An exotic forests $\pi$ is connected if it cannot be written as $\pi = \gamma \cdot_\omega \eta$ for any $\gamma, \eta \in EF\{\mathbf{1}\}$ and any ordering $\omega$ of the roots of $\gamma \cdot \eta$. The set of connected exotic forests is denoted by $EF_C$ and the corresponding space by $\text{Conn} := \text{Span } EF_C$. We note that $\text{Conn} \subset \text{Irr}$.

**Remark 4.3.** Analogously to Remark 4.2, the space of planar exotic forests forms the tensor algebra $T(\text{Irr})$ of $\text{Irr}$, and, thus, a universal enveloping Lie algebra $\mathcal{E}\mathcal{F} = \mathcal{U}(\text{Lie}(\text{Irr}))$ with the Lie bracket $[x, y] = x \cdot y - y \cdot x$ on $\mathcal{E}\mathcal{F}$. The first terms are

$$\text{Lie}(\text{Irr}) = \text{Span}\{??, ??, ??, ??, ??, ??, ??, ??, ??, [??, ??], \dots\}.$$

In order to show that the space of exotic planar forests forms a Hopf algebra, we need to prove the compatibility between the algebraic and coalgebraic structures. We consider the product $\tilde{\otimes}_\omega : \mathcal{E}\mathcal{F} \otimes \mathcal{E}\mathcal{F} \to \mathcal{E}\mathcal{F}$ defined as

$$\pi \tilde{\otimes}_\omega \eta = \pi \diamond_\omega \eta - \pi \cdot_\omega \eta, \quad \text{for } \pi \in EF_C, \eta \in EF, \\ (\pi \cdot_\omega \gamma) \tilde{\otimes}_\omega \eta = \pi \tilde{\otimes}_\omega (\gamma \tilde{\otimes}_\omega \eta) - (\pi \tilde{\otimes}_\omega \gamma) \tilde{\otimes}_\omega \eta, \quad \text{for } \gamma \in EF.$$

Intuitively, the product $\tilde{\otimes}_\omega$ is a modification of the generalised Grossman-Larson product $\diamond_\omega$ in which we require each element of $EF_C$ in the left operand to attach at least one root to a node of the right operand. For example, let $\omega : a \leq c \leq b$, then,

$$?? \tilde{\otimes}_\omega ?? = ?? + ?? + ??.$$

We note the following property

$$\pi \diamond_\omega \eta = \pi_{(1)} \cdot_\omega (\pi_{(2)} \tilde{\otimes}_\omega \eta), \quad \text{for } \pi, \eta \in \mathcal{E}\mathcal{F},$$

where $\Delta_\omega(\pi) = \pi_{(1)} \otimes \pi_{(2)}$ is the deshuffle coproduct defined on $\mathcal{E}\mathcal{F}$. We note that $\pi \diamond \eta = \pi \diamond_\omega \eta$ with $\omega$ being the order for which $\pi \cdot_\omega \eta = \pi \cdot \eta$. We emphasize that the product $\diamond_\omega$ on $\mathcal{E}\mathcal{F}$ is associative, since $(\mathcal{E}\mathcal{F}, \diamond_\omega)$ is a subalgebra of $(\overline{\mathcal{F}}_D^L, \diamond_\omega)$, where $\diamond_\omega$ is associative. This follows from the fact that $\diamond_\omega$ can be defined as $\diamond$, with the roots arranged according to $\omega$ in the final step.

**Theorem 4.4.** The space of exotic forests $\mathcal{E}\mathcal{F}$ together with Grossman-Larson product and deshuffle coproduct forms a Hopf algebra.
Proof. We note that the algebra structure is obtained by considering a subalgebra of $(\overline{\mathcal{F}}_D^L, \diamond)$, while the coalgebra structure is obtained by taking a quotient of the coalgebra $(\overline{\mathcal{F}}_D^L, \Delta_\omega)$. This means that the algebra and coalgebra structures are not automatically compatible. We prove the compatibility condition of the generalised Grossman-Larson product for any ordering $\omega$,

$$ \Delta_{\omega}(\pi \diamond_{\omega} \eta) = \Delta_{\omega}(\pi) \diamond_{\omega} \Delta_{\omega}(\eta), $$

by induction on the number of connected components in $\pi$. Assume the compatibility is proven for all forests with the number of components less or equal to the number of components in $\pi$. Consider $\tilde{\pi} \cdot_{\omega} \pi$ with $\tilde{\pi} \in EF_C$ and ordering $\omega$ of roots of $\tilde{\pi} \cdot \pi \cdot \eta$ with an arbitrary order between the roots of $\tilde{\pi} \cdot \pi$ and $\eta$, then,

$$ \begin{aligned} \Delta_{\omega}((\tilde{\pi} \cdot_{\omega} \pi) \diamond_{\omega} \eta) &= \Delta_{\omega}(\tilde{\pi} \diamond_{\omega} (\pi \diamond_{\omega} \eta) - (\tilde{\pi} \tilde{\cdot}_{\omega} \pi) \diamond_{\omega} \eta) \\ &= \Delta_{\omega}(\tilde{\pi}) \diamond_{\omega} \Delta_{\omega}(\pi) \diamond_{\omega} \Delta_{\omega}(\eta) - \Delta_{\omega}(\tilde{\pi} \tilde{\cdot}_{\omega} \pi) \diamond_{\omega} \Delta_{\omega}(\eta) \\ &= \Delta_{\omega}(\tilde{\pi} \diamond_{\omega} \pi - \tilde{\pi} \tilde{\cdot}_{\omega} \pi) \diamond_{\omega} \Delta_{\omega}(\eta) = \Delta_{\omega}(\tilde{\pi} \cdot_{\omega} \pi) \diamond_{\omega} \Delta_{\omega}(\eta), \end{aligned} $$

where we use the associativity of $\diamond_{\omega}$ and coassociativity of $\Delta_{\omega}$. It remains to show that

$$ \Delta_{\omega}(\tilde{\pi} \diamond_{\omega} \eta) = \Delta_{\omega}(\tilde{\pi}) \diamond_{\omega} \Delta_{\omega}(\eta), \quad \text{for } \tilde{\pi} \in \text{Conn}, \eta \in \mathcal{E}\mathcal{F}. $$

We note that $\tilde{\pi}$ has multiple roots which, when grafted onto different connected components of $\eta$, connect them into a single component. Therefore, we have

$$ \Delta_{\omega}(\tilde{\pi} \diamond_{\omega} \eta) = (\tilde{\pi} \diamond_{\omega} \eta_{(1)}) \otimes \eta_{(2)} + \eta_{(1)} \otimes (\tilde{\pi} \diamond_{\omega} \eta_{(2)}) = \Delta_{\omega}(\tilde{\pi}) \diamond_{\omega} \Delta_{\omega}(\eta). $$

This proves the compatibility between the Grossman-Larson product and deshuffle coproduct. Therefore, $\mathcal{E}\mathcal{F}$ forms a graded connected bialgebra, that is, a Hopf algebra. $\square$

**Remark 4.5.** The space $(\mathcal{E}\mathcal{F}, \cdot, \Delta_{\omega}, \tilde{\cdot})$ is not a post-Hopf algebra, as defined in [61], since in general

$$ \pi \tilde{\cdot} (\pi_1 \cdot \pi_2) \neq (\pi_{(1)} \tilde{\cdot} \pi_1) \cdot (\pi_{(2)} \tilde{\cdot} \pi_2) \text{ in } \mathcal{E}\mathcal{F}. $$

Moreover, $(\mathcal{E}\mathcal{F}, \cdot, \Delta_{\omega}, \tilde{\cdot})$ is not a post-Hopf algebra for the same reason. This is due to the fact that a connected planar exotic forest can have multiple roots. However, we see that $(\mathcal{E}\mathcal{F}, \cdot, \tilde{\cdot})$ is a D-subalgebra of $(\overline{\mathcal{F}}_D^L, \cdot, \tilde{\cdot})$, as defined in [78], with the grading of the associative algebra $(\mathcal{E}\mathcal{F}, \cdot)$ given by the number of roots. A straightforward extension of the arguments in [79] fails to show that the spaces $(\mathcal{E}\mathcal{F}, \cdot)$ and $(\mathcal{E}\mathcal{F}, \tilde{\cdot})$ are either brace or symmetric brace algebras for the same reason.

The Connes-Kreimer Hopf algebra structure is the dual of the Grossman-Larson one in a Euclidean setting. On homogeneous spaces, it generalises into the Munthe-Kaas-Wright Hopf algebra [78, 35].

Let an admissible cut $c$ of an exotic planar tree $\tau$ be a subset of edges of $\tau$ such that it contains at most one edge from each path in $\tau$ from a leaf to the root. If $e = (u, v)$ is in $c$, then all incoming edges of $v$ which are to the left of $e$ are also in $c$. Moreover, given a liana with nodes $v_1$ and $v_2$, if any edge in the path from $v_1$ to the root is in $c$, then the path from $v_2$ to the root must also contain an edge in $c$.

Given an admissible cut $c$, $\tau \backslash c$ is an exotic planar forest. Let $R^c(\tau)$ be the component of $\tau \backslash c$ which contains the root of $\tau$. Let $P^c(\tau)$ be the linear combination of exotic planar forests obtained by taking the components of $\tau \backslash c$ that do not contain the root of $\tau$ and shuffling the roots while preserving the order of the roots cut off from the same node.
**Theorem 4.6.** Let the Munthe-Kaas-Wright coproduct on $\mathcal{E}\mathcal{F}$ be, for $\tau \in ET$ and $\pi \in EF$, defined as

$$ \Delta_{MKW}(\tau) := \sum_{\text{adm. cut } c} P^c(\tau) \otimes R^c(\tau), \\ \Delta_{MKW}(\pi) := (\text{id} \otimes B^-) \Delta_{MKW}(B^+(\pi)), $$

where the sum is over admissible cuts $c$. Then, $(\mathcal{E}\mathcal{F}, \sqcup, \Delta_{MKW})$ is a Hopf algebra dual to the planar Grossman-Larson Hopf algebra $(\mathcal{E}\mathcal{F}, \diamond, \Delta_{\sqcup})$. Its convolution product,

$$ a * b = (a \otimes b) \circ \Delta_{MKW}, $$

is called the composition law.

*Proof.* We recall that the dual of a Hopf algebra is a Hopf algebra. It is clear that $\sqcup$ and $\Delta_{\sqcup}$ are dual. It remains to show that the Munthe-Kaas-Wright coproduct $\Delta_{MKW}$ is dual to the Grossman-Larson product $\diamond$. We recall that Munthe-Kaas-Wright coproduct and Grossman-Larson product are dual over the planar decorated forests. The Munthe-Kaas-Wright coproduct over exotic forests is obtained by taking the quotient over $\mathcal{I}$ defined in (4.1) for the construction of deshuffle coproduct. This is the dual approach to the way the Grossman-Larson product over exotic forests is obtained. $\square$

Some examples of the values of $\Delta_{MKW}$ are

$$ \begin{aligned} \Delta_{MKW}(??) &= ?? \otimes \mathbf{1} + \mathbf{1} \otimes ??, \\ \Delta_{MKW}(??) &= ?? \otimes \mathbf{1} + ?? \otimes ?? + \mathbf{1} \otimes ??, \\ \Delta_{MKW}(??) &= ?? \otimes \mathbf{1} + 2?? \otimes ?? + (?? + ??) \otimes ?? + (?? + ??) \otimes ?? \\ &\quad + (?? + ??) \otimes ?? + 2?? \otimes ?? + \mathbf{1} \otimes ?? \end{aligned} $$

The Hopf algebra structures on $\mathcal{E}\mathcal{F}$ translate to the structures on $T(\mathcal{X}(M))$ via the use of the elementary differential map $\mathbb{F}$ defined in Section 3.3. We recall that decorated planar forest form a free D-algebra with the concatenation and grafting products [78] which is generated by the set of decorations $D$. The elementary differential map over decorated planar forests $\mathbb{F}^{\text{dec}} : \mathcal{F}_D \to T(\mathcal{X}(M))$ is a D-algebra homomorphism that is fully characterized by its values on the nodes decorated by $D = \{\bullet, \circ^l : 1 \le l \le L\}$,

$$ \mathbb{F}^{\text{dec}}(\bullet) = f, \quad \mathbb{F}^{\text{dec}}(\circ^l) = \xi^{d,l} E_d, $$

where the $\xi^{d,l}$ are independent standard Gaussian random variables. The map $\mathbb{F}^{\text{dec}}$ thus defined is equivalent to the map $\mathbb{F}^{\text{dec}}$ from Definition 3.4.

The space of exotic forests $\mathcal{E}\mathcal{F}$ is isomorphic to a D-subalgebra $\mathcal{F}_D^\mathcal{E}$ of the D-algebra of decorated forests with $D = \{??, ?? : d \in \mathbb{N}\}$ through the isomorphism $\varphi : \mathcal{E}\mathcal{F} \to \mathcal{F}_D^\mathcal{E}$ as is described in the beginning of Section 4.3, and, thus, inherits the elementary differential homomorphism as $\mathbb{F} := \tilde{\mathbb{F}}^{\text{dec}} \circ \varphi$ with $\tilde{\mathbb{F}}^{\text{dec}} : \mathcal{F}_D^\mathcal{E} \to T(\mathcal{X}(M))$ defined as

$$ \tilde{\mathbb{F}}^{\text{dec}}(??) = f, \quad \tilde{\mathbb{F}}^{\text{dec}}(??) = E_d. $$

This definition is equivalent to the one given in Definition 3.9.
**Proposition 4.7.** The elementary differential $\mathcal{F}: \mathcal{E}\mathcal{F} \to T(\mathcal{X}(\mathcal{M}))$ is a Hopf algebra morphism:

$$ \mathcal{F}(\pi_1 \rightsquigarrow \pi_2) = \mathcal{F}(\pi_1) \triangleright \mathcal{F}(\pi_2), \quad \mathcal{F}(\pi_1 \cdot \pi_2) = \mathcal{F}(\pi_1) \cdot \mathcal{F}(\pi_2), \quad \mathcal{F}(\pi_1 \diamond \pi_2) = \mathcal{F}(\pi_1) * \mathcal{F}(\pi_2). $$

In particular, the frozen composition is represented by concatenation and the composition of differential operators is represented by the Grossman-Larson product:

$$ \begin{aligned} (\mathcal{F}(\pi_1) \cdot \mathcal{F}(\pi_2)) \triangleright \phi(p) &= \mathcal{F}_p(\pi_1) \triangleright (\mathcal{F}_p(\pi_2) \triangleright \phi)(p) = \mathcal{F}(\pi_1 \cdot \pi_2) \triangleright \phi(p), \\ (\mathcal{F}(\pi_1) * \mathcal{F}(\pi_2)) \triangleright \phi(p) &= \mathcal{F}(\pi_1) \triangleright (\mathcal{F}(\pi_2) \triangleright \phi)(p) = \mathcal{F}(\pi_1 \diamond \pi_2) \triangleright \phi(p). \end{aligned} $$

# 5 Numerical experiments

We provide several numerical experiments to confirm the theoretical findings and the numerical properties of the new schemes. In particular, we present a practical guide for the implementation of frozen flow methods, convergence curves for the weak error on the orthogonal group, a comparison with the sampling geometric methods from the literature on the sphere, and an ergodicity test for the simulation of generalised Cauchy measures in a context where the assumptions of our convergence analysis are not satisfied.

## 5.1 Practical implementation of frozen flow methods

The implementation of frozen flow methods relies on the ability to find local frames and to compute the solution of the frozen equation

$$ y'(t) = \sum_{d} \alpha_{d} E_{d}(y(t)), \quad y(0) = y_{0}. $$

On homogeneous manifolds, we recall that global frames can be derived from any basis $(A_d)$ of the Lie algebra $\mathfrak{g}$ with the derivative of Lie-group action: $E_d(p) = A_d \cdot p$. On matrix homogeneous manifolds, the solution of the frozen flow equation coincides with the matrix exponential:

$$ y(t) = \mathrm{Exp}(\sum_{d} \alpha_{d} A_{d}) \cdot y_{0}. $$

For sampling the invariant measure of Riemannian Langevin dynamics (1.3), the algorithms require at each step a local orthonormal frame basis and the solution of the frozen flow equation. While one can change the orthonormal local frame between different steps, it is crucial to rely only on one fixed frame per step. On Lie groups, choosing an orthonormal basis $(A_d)$ of the Lie algebra $\mathfrak{g}$ for the definition of the $E_d$ naturally yields an orthonormal frame. On general Riemannian manifolds, one can provide an orthonormal basis of $T_{X_0}\mathcal{M}$ and parallel transport it around $X_0$ to obtain a local orthonormal frame basis.

In general, solving the frozen flow equation is not straightforward. It is equally difficult to solving the geodesic equation, as done in [10]. Similar to the idea of retractions (see, for instance, the textbook [4]), the frozen flow could be replaced by any approximation with sufficiently high order of accuracy and that lies on the manifold. A standard example of such an approximation is the Cayley map [46]. In our numerical experiments, we rely on exact expressions of the frozen flows using either Lie geometry or coordinates for the sake of simplicity.

We emphasize that the geometric operations needed for stochastic intrinsic methods are available on a wide collection on manifolds and in a variety of programming language (see, for instance, [36, 75, 7]).
## 5.2 Brownian dynamics on the special orthogonal group

Let $\mathcal{M} = \text{SO}_p(\mathbb{R})$ be the compact Lie group of special orthogonal matrices of size $p \times p$. Let $(A_d)$ be an orthonormal basis of its Lie algebra $\mathfrak{so}_d$ of antisymmetric matrices and for the bi-invariant metric $g(A, B) = \text{Tr}(A^T B)$. Let $E_d(y) = A_d y$ be the associated right-invariant orthonormal frame basis. As the manifold is smooth and compact, the geometric assumptions are all satisfied automatically and our analysis applies.

Following Example 2.15, Brownian dynamics on $\text{SO}_p(\mathbb{R})$ are given by

$$dX(t) = \sum_{d=1}^{D} E_d(X(t)) \circ dW_d(t), \quad D = \dim(\mathfrak{so}_p) = \frac{p(p-1)}{2}.$$

In this context, the implementation of the new frozen flow methods is very natural as it relies only on fast linear algebra operations. In particular, the new second order method (2.12) becomes (under a time rescaling)

$$X_{n+1} = \text{Exp} \left( \sum_{d=1}^{D} \left( \left( \frac{\sqrt{2}}{2} - 1 \right) \sqrt{h} \xi_n^{d,1} + \frac{\sqrt{2}}{2} \sqrt{h} \xi_n^{d,2} \right) A_d \right) \text{Exp} \left( \sum_{d=1}^{D} \sqrt{h} \xi_n^{d,1} A_d \right) X_n.$$

Note that there is no timestep restriction for the scheme to be well-posed. This is in striking contrast to projection methods, which would rely on a tedious non-geometric implementation, on the embedding of $\text{SO}_p(\mathbb{R})$ into $\mathbb{R}^{p \times p}$ (of more than double dimension), and would require a small enough timestep for the projection map to be well-defined, even in this simple context of a compact manifold.

We consider the weak approximation of Brownian dynamics on $\text{SO}_3(\mathbb{R})$ with the Euler frozen flow method (2.9) and the new second order method (2.12) (with discrete bounded random variables with correct first moments - see Section 2.3) for the test function

$$\phi(x) = \exp\left(-\frac{\text{Tr}((x - I_d)^T (x - I_d))}{2d}\right),$$

the initial condition $X_0 = I_d$, final time $T = 1$, and different timesteps to observe the order of convergence. We emphasize that the explicit expression of Brownian dynamics $X(t)$ on $\text{SO}_3(\mathbb{R})$ are not given by $\text{Exp}(\sum A_d W_d(t))X_0$, but rather by a stochastic Magnus expansion (see, for instance, [19, 18, 46, 91, 49]):

$$X(t) = \text{Exp}(\Omega(t))X_0, \quad d\Omega(t) = d\text{Exp}_{\Omega(t)}^{-1}\left(\sum_d A_d \circ dW_d(t)\right), \quad \Omega(0) = 0,$$

whose first terms are

$$\begin{aligned} \Omega(t) = & \sum_{d_1} A_{d_1} W_{d_1}(t) - \frac{1}{2} \int_0^t \sum_{d_1, d_2} [A_{d_2}, A_{d_1}] W_{d_2}(s) \circ dW_{d_1}(s) \\ & + \frac{1}{12} \int_0^t \sum_{d_1, d_2, d_3} [A_{d_3}, [A_{d_2}, A_{d_1}]] W_{d_3}(s) W_{d_2}(s) \circ dW_{d_1}(s) \\ & + \frac{1}{4} \int_0^t \int_0^s \sum_{d_1, d_2, d_3} [[A_{d_3}, A_{d_2}], A_{d_1}] W_{d_3}(u) \circ dW_{d_2}(u) \circ dW_{d_1}(s) + \dots \end{aligned}$$

For simplicity, the reference solution for $\mathbb{E}[\phi(X(T))]$ is chosen as the output of the new second order scheme with the timestep $h_{\text{ref}} = 2^{-12}$. We observe on Figure 1 the correct order of convergence of the two new methods. The second order integrator reaches the Monte-Carlo error threshold almost instantly, even for relatively large timestep.
<figure><img src="image_5.png" /><figcaption>Figure 1: Brownian excursion at time $T = 1$ on $SO_3(\mathbb{R})$ represented by the trajectories of the three columns on the sphere $S^2$ (left) and error curve for the weak approximation of Brownian dynamics on $SO_3(\mathbb{R})$ with $\phi(x) = \exp(-\text{Tr}((x - I_d)^T(x - I_d))/2d)$, $T = 1$ and $M = 10^8$ trajectories (right).</figcaption></figure>

## 5.3 Ergodic dynamics on the sphere

Our second test case focuses on the sampling of Riemannian Langevin dynamics on the sphere $S^2$, following the experiment of [10]. While the approximation for the invariant measure is not the focus of the new methods implemented in this paper, this example allows us to conveniently compare the performance of the new frozen flow methods to the alternatives from the literature. The sphere $S^d$ can be seen as a homogeneous manifold with respect to the Lie group $G = SO_{d+1}(\mathbb{R})$, so that we could consider a natural frame derived from a basis of the Lie algebra $so_{d+1}$ (which would typically fit the simulation of a stochastic rigid body). However, we need an orthonormal frame to sample from the Riemannian Langevin dynamics (1.3), so that we use a different choice of frame relying on the standard spherical coordinates. We consider the coordinates on $S^2$ minus the north/south poles:

$$y_{\theta, \varphi} = \begin{pmatrix} \cos(\theta) \cos(\varphi) \\ \cos(\theta) \sin(\varphi) \\ \sin(\theta) \end{pmatrix} \in S^2, \quad \theta \in [-\frac{\pi}{2}, \frac{\pi}{2}], \quad \varphi \in [0, 2\pi[$

with the associated orthonormal frame basis:

$$E_1(y_{\theta, \varphi}) = \begin{pmatrix} -\sin(\theta) \cos(\varphi) \\ -\sin(\theta) \sin(\varphi) \\ \cos(\theta) \end{pmatrix}, \quad E_2(y_{\theta, \varphi}) = \begin{pmatrix} -\sin(\varphi) \\ \cos(\varphi) \\ 0 \end{pmatrix},$$

and a second frame associated to the spherical coordinates on $S^2$ minus $[\pm 1, 0, 0]$. Using two frames guarantees a bounded Lipschitz constant for the chosen frame and does not impact the sampled measure. As the manifold is compact, the assumptions of our analysis are satisfied automatically. The frozen flow is computed explicitly in coordinates. In order to obtain a fair comparison between the methods, we use Gaussian random variables.

We now compare the numerical behaviour of the frozen flow Euler method (2.9), the new second order method (2.11), the Riemannian Langevin method [10], and the extrinsic projection methods of order one [59] for sampling the invariant measure of Riemannian Langevin dynamics (1.3) on $S^2$. We consider the potential $V(y_{\theta, \varphi}) = -\sin(\theta)$, the associated vector field (containing the Ito correction)

$$F = -\nabla V - \nabla_{E_1} E_1 - \nabla_{E_2} E_2,$$
and the test function $\phi(x) = x_3^2$. In this context, the generator is $\mathcal{L}\phi = -\phi'(\nabla V) + \frac{1}{2}\Delta_{S^2}$, with the spherical Laplacian $\Delta_{S^2}$. Assuming ergodicity (which is observed numerically for the methods studied here), the numerical schemes approximate the quantity

$$ \frac{\int_{S^2} \phi(x) e^{-V(x)} d\mathrm{vol}(x)}{\int_{S^2} e^{-V(x)} d\mathrm{vol}(x)} = 3 - \frac{2}{\tanh(1)}. $$

We observe the weak error curves in Figure 2, as well as the expected orders of convergence. The new second-order frozen flow method outperforms the other integrators and reaches the Monte-Carlo threshold with much larger timestep.

<figure><img src="image_6.png" /><figcaption>Figure 2: A trajectory of the new second order method (2.12) (left) and the convergence curve for the weak approximation of stochastic rigid body equation on the sphere (right) with $\phi(x) = x_3^2$, $T = 10$ and $M = 10^8$ trajectories.</figcaption></figure>

## 5.4 Generalised Cauchy measures

The generalised Cauchy measures are probability measures $\mu_\beta$ defined on $\mathbb{R}^n$ by

$$ \mu_{\beta}(dx) = (1 + |x|^2)^{-\beta} \frac{dx}{Z_{\beta}}, \quad Z_{\beta} = \frac{\Gamma(\beta - n/2)}{\Gamma(\beta)}, $$

with $\beta > n/2$. They can be interpreted as finite dimensional approximations of Gaussian distributions, as in [11], and are related to the fast diffusion equation for which they play the role of the heat kernel and are known as Barenblatt' solutions (see [89]). A generalised Cauchy measure is the ergodic distribution of the dynamics

$$ dX_t = -2(\beta - 1)X_t dt + \sqrt{2(1 + |X_t|^2)}dB(t), \quad (5.1) $$

where $B(t)$ is a 2-dimensional Brownian motion. A direct numerical simulation of this SDE presents numerous difficulties as the multiplicative non-Lipschitz noise brings tedious stability problems, so that most standard Euclidean numerical methods fail. However, the noise can be interpreted as additive at the cost of a change of metric $g_x = (1 + |x|^2)^{-1}\langle \cdot, \cdot \rangle$. This change is not simply a transport of difficulties from noise to geometry, as it reinterprets the problem in a convenient space. The same approach has been fruitful in the study of the associated spectral gap, i.e. the speed of convergence of the semigroup to the ergodic measure, in [44]. This context
calls for methods on manifolds, as Euclidean methods are not adapted. We emphasize that the process $(X_t)_t$ does not have bounded moments at all orders. Thus, the convergence results of the standard Euclidean analysis, and of our analysis as well, do not hold. This makes the simulation of generalised Cauchy measures a convenient toy problem for testing our methods.

We are concerned with the simulation of generalised Cauchy measure on $\mathbb{R}^2$. As explained in [43], the adapted manifold can be interpreted as the revolution surface

$$S = \{(\tanh(r) \cos(\theta), \tanh(r) \sin(\theta), h(r)); (r, \theta) \in \mathbb{R}_+ \times \mathbb{S}^1\}, \quad (5.2)$$

with $h(r) = \operatorname{argsinh}(\cosh(r)) - \cosh^{-1}(r)\sqrt{\cosh^2(r) + 1}$. correction of the function h However, it is more handy to see it through a global map, as $\mathbb{R}_+ \times \mathbb{S}^1$, with coordinate $(r, \theta)$, endowed with the metric

$$g = \begin{pmatrix} 1 & 0 \\ 0 & \tanh^2(r) \end{pmatrix}$$

in the basis $(\partial_r, \partial_\theta)$. This map comes from fitting polar coordinates on $\mathbb{R}^2$,

$$(x, y) = (\sinh(r) \cos(\theta), \sinh(r) \sin(\theta)),$$

which trivialise the metric $g$. We define the following orthonormal frame basis on $\mathbb{R}_+^* \times \mathbb{S}^1$

$$E_1(r, \theta) = \partial_r, \quad E_2(r, \theta) = \frac{1}{\tanh(r)} \partial_\theta.$$

In this manifold, generalised Cauchy measures correspond to Gibbs measures associated to the potential $V_\beta(r, \theta) = 2(\beta - 1)\log(\cosh(r))$ and the generator

$$\mathcal{L} = E_1^2 + E_2^2 + \left( \frac{1}{\tanh(r)} - (2\beta - 1) \tanh(r) \right) E_1.$$

Note that the frozen flow is explicit on $S$. For all $\alpha, \beta \in \mathbb{R}$ and $(r_0, \theta_0) \in \mathbb{R}_+^* \times \mathbb{S}^1$, we find

$$\exp(t(\alpha E_1 + \beta E_2)) \cdot (r_0, \theta_0) = \left( r_0 + \alpha t, \theta_0 + \frac{\beta}{\alpha} \log \left( \frac{\sinh(r_0 + \alpha t)}{\sinh(r_0)} \right) \right),$$

for all $t \ge 0$ such that $r_t$ stays positive. This brings an additional difficulty in this experiment. Although the potential $f^1$ blows up at the origin, the coefficient

$$\alpha = h f^1(r) + \sqrt{2h\xi}$$

with $\xi$ a Gaussian random variable, can be negative and the frozen flow not defined. The use of bounded random variable is instrumental in avoiding this instability. In particular, the method (2.9) with Rademacher random variables is well-posed.

Our numerical experiment aims to illustrate the ergodicity of the simplest frozen flow method (2.9), which, to the best of our knowledge, has not yet been established in the literature. We use the first eigenfunctions of the generator $\phi^1(r, \theta) = \sinh(r)^2 - (\beta - 2)^{-1}$, $\phi^2(r, \theta) = \sinh(r) \cos(\theta)$, which both satisfy $\int \phi^k d\mu_\beta = 0$. As a consequence of the Poincaré inequality, the variance of the semigroup associated to $\mathcal{L}$ converges at exponential speed to 0. The coefficient of the exponential is given by the spectral gap (see [45] for the spectral gap expression). The test functions $\phi^k$ are extremal for the Poincaré inequality. Hence, for these particular test functions, the exponential convergence can be observed not only for the variance but also pointwise. In Figure 3, we
illustrate numerically the exponential convergence to equilibrium of the numerical methods, in the spirit of [32] in the Euclidean setting. In particular, we observe a convergence in $Ce^{-4(\beta-2)t}$ and $Ce^{-2(\beta-1)t}$ for our specific choice of test functions, confirming numerically the results in [45]. For the first eigenfunction $\phi^1$, the error for the invariant measure reaches the discretisation bias with exponential speed. We observe that the bias threshold evolves linearly in $h$ (up to Monte-Carlo error), as expected for a method of order one. For $\phi^2$, we do not observe the bias of the Euler frozen flow method (2.9), hinting that for this specific test function, the method might be of higher order of accuracy for the invariant measure. For the sake of comparison, we mention that the standard Euclidean integrators applied to (5.1) face severe instability issues and thus fail to converge, let alone be ergodic.

<figure><img src="image_7.png" /><figcaption>Figure 3: Evolution in time of the error for the invariant measure of the Euler frozen flow method (2.9) for the simulation of generalised Cauchy measures on the manifold (5.2) with $\beta = 4$, $M = 10^8$ trajectories, and the test functions $\phi^1(r, \theta) = \sinh(r)^2 - (\beta - 2)^{-1}$ (left) and $\phi^2(r, \theta) = \sinh(r) \cos(\theta)$ (right).</figcaption></figure>

# 6 Conclusion

In this paper, we provided a brand-new class of intrinsic discretisations for the approximation of stochastic dynamics on Riemannian manifolds with high weak order of accuracy. We presented a new robust convergence analysis, that naturally generalises the Euclidean analysis, and the new algebraic formalism of planar exotic trees, forests, and series for the study of the order conditions. The analysis applies in particular for the creation of high-order sampling methods for Riemannian Langevin dynamics. The new methods were tested on a variety of different manifolds, where we observe that the new second order frozen flow method outperforms the other approaches from the literature in terms of accuracy.

The tools developed in the present paper open several avenues for future research. Similar to the Euclidean setting, we will extend the analysis to general SDEs with multiplicative noise and to the creation of integrators with high-order of accuracy for sampling the invariant measure, and low weak order (in the spirit of [56, 90]). This requires generalisations of the exotic forests formalism [20], the challenging extension of the integration by parts of trees in the planar case [54, 55, 14, 15], and the extension of stochastic backward error analysis for the invariant measure [85, 92, 32, 2, 15] on manifolds. While the focus of this paper is on discretisations and improvements in accuracy only, our approach could naturally be combined with popular sampling techniques, such as MLMC [37, 38], perturbations [58, 33, 1],... Then, following the work [10], one could wonder whether a purely Riemannian approach of algebraic order theory could be used. Our approach follows the one of deterministic Lie-group methods and thus relies
on a curvature-free connection, unrelated to the natural Levi-Civita connection on Riemannian manifolds. The study of such new methods is already open in the deterministic setting and relies on the challenging general understanding of the connection algebra [5]. Finally, the new algebraic formalism of planar exotic series is interesting in itself and could be studied for its universal combinatorial, algebraic, and geometric properties, but also for its potential applications in different fields, in the spirit of the use of Butcher series and their extensions in rough paths [41, 57], renormalisation theory [17], variational calculus [52, 51] or approximation of PDEs [16]. These projects will be studied in upcoming works.

**Acknowledgements.** The authors would like to thank R. Bergmann, E. Celledoni, H. Munthe-Kaas, and B. Owren for helpful discussions and advice on manipulating the package Manifolds.jl [7]. The authors acknowledge the support of the Swiss National Science Foundation, projects No 200020_214819, and No. 200020_192129, of the Research Council of Norway through project 302831 "Computational Dynamics and Stochastics on Manifolds" (CODYSMA), and of the French program ANR-11-LABX-0020-0 (Labex Lebesgue). The computations were performed on the Baobab cluster of the University of Geneva.

## References

[1] A. Abdulle, G. A. Pavliotis, and G. Vilmart. Accelerated convergence to equilibrium and reduced asymptotic variance for Langevin dynamics using Stratonovich perturbations. *C. R. Math. Acad. Sci. Paris*, 357(4):349-354, 2019.

[2] 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.

[3] M. Ableidinger and E. Buckwar. Weak stochastic Runge-Kutta Munthe-Kaas methods for finite spin ensembles. *Applied Numerical Mathematics*, 118:50-63, 2017.

[4] P.-A. Absil, R. Mahony, and R. Sepulchre. *Optimization algorithms on matrix manifolds*. Princeton University Press, Princeton, NJ, 2008. With a foreword by Paul Van Dooren.

[5] M. J. H. Al-Kaabi, K. Ebrahimi-Fard, D. Manchon, and H. Z. Munthe-Kaas. Algebraic aspects of connections: from torsion, curvature, and post-lie algebras to Gavrilov's double exponential and special polynomials. *Submitted*, 2022.

[6] A. Antonyuk and A. Antonyuk. Nonexplosion and solvability of nonlinear diffusion equations on noncompact manifolds. *Ukr. Math. J.*, 59:1632-1652, 2007.

[7] S. D. Axen, M. Baran, R. Bergmann, and K. Rzecki. Manifolds.jl: An extensible julia framework for data analysis on manifolds. *ACM Transactions on Mathematical Software*, 49(4), dec 2023.

[8] D. Bakry. Un critère de non-explosion pour certaines diffusions sur une variété riemannienne complète. *C.R. Acad. Sc. Paris*, 303(1):23-26, 1986.

[9] D. Bakry, I. Gentil, and M. Ledoux. *Analysis and geometry of Markov diffusion operators*, volume 348 of *Grundlehren Math. Wiss*. Cham: Springer, 2014.

[10] K. Bharath, A. Lewis, A. Sharma, and M. V. Tretyakov. Sampling and estimation on manifolds using the Langevin diffusion. *arXiv:2312.14882*, 2023.

[11] S. G. Bobkov and M. Ledoux. Weighted Poincaré-type inequalities for Cauchy and other convex measures. *Ann. Probab.*, 37(2):403-427, 2009.

[12] G. Bogfjellmo. Algebraic structure of aromatic B-series. *J. Comput. Dyn.*, 6(2):199-222, 2019.
[13] G. Bogfjellmo, E. Celledoni, R. I. McLachlan, B. Owren, and G. R. W. Quispel. Using aromas to search for preserved measures and integrals in Kahan's method. *Math. Comp.*, 93(348):1633–1653, 2024.

[14] E. Bronasco. Exotic B-series and S-series: algebraic structures and order conditions for invariant measure sampling. *Found. Comput. Math.*, pages 1–31, 2024.

[15] E. Bronasco and A. Laurent. Hopf algebra structures for the backward error analysis of ergodic stochastic differential equations. *arXiv:2407.07451*, 2024.

[16] Y. A. Bronsard, Y. Bruned, and K. Schratz. Approximations of dispersive PDEs in the presence of low-regularity randomness. *arXiv preprint arXiv:2205.02156*, 2022.

[17] C. Brouder. Runge-Kutta methods and renormalization. *Eur. Phys. J. C*, 12(3):521–534, 2000.

[18] K. Burrage and P. M. Burrage. High strong order methods for non-commutative stochastic ordinary differential equation systems and the Magnus formula. *Phys. D*, 133(1-4):34–48, 1999. Predictability: quantifying uncertainty in models of complex phenomena (Los Alamos, NM, 1998).

[19] P. M. Burrage. *Runge-Kutta methods for stochastic differential equations*. PhD thesis, Dept. Maths., Univ. Queensland, Australia, 1999.

[20] A. Busnot Laurent, K. Debrabant, and A. Kværnø. Optimal stochastic Runge-Kutta methods for the weak integration of stochastic dynamics. *Preprint*, 2025.

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

[22] J. C. Butcher. *B-series: algebraic analysis of numerical methods*. Springer, 2021.

[23] D. Calaque, K. Ebrahimi-Fard, and D. Manchon. Two interacting Hopf algebras of trees: a Hopf-algebraic approach to composition and substitution of B-series. *Adv. in Appl. Math.*, 47(2):282–308, 2011.

[24] E. Celledoni, A. Marthinsen, and B. Owren. Commutator-free Lie group methods. *Future Generation Computer Systems*, 19(3):341–352, 2003.

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

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

[27] X. Cheng, J. Zhang, and S. Sra. Efficient sampling on Riemannian manifolds via Langevin MCMC. *Advances in Neural Information Processing Systems*, 35:5995–6006, 2022.

[28] P. E. Crouch and R. Grossman. Numerical integration of ordinary differential equations on manifolds. *Journal of Nonlinear Science*, 3:1–33, 1993.

[29] C. Curry and A. Schmeding. Convergence of lie group integrators. *Numerische Mathematik*, 144:357–373, 2020.

[30] R. L. Davidchack, T. Ouldridge, and M. Tretyakov. New Langevin and gradient thermostats for rigid body dynamics. *The Journal of chemical physics*, 142(14), 2015.

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

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

[33] A. B. Duncan, T. Lelièvre, and G. A. Pavliotis. Variance reduction using nonreversible Langevin samplers. *J. Stat. Phys.*, 163(3):457–491, 2016.
[34] K. Ebrahimi-Fard, A. Lundervold, and H. Z. Munthe-Kaas. On the Lie enveloping algebra of a post-Lie algebra. *J. Lie Theory*, 25(4):1139–1165, 2015.

[35] K. Ebrahimi-Fard and L. Rahm. A survey on the Munthe-Kaas–Wright Hopf algebra. *Journal of Computational Dynamics*, 2024.

[36] K. Engø, A. Marthinsen, and H. Z. Munthe-Kaas. Diffman: An object-oriented matlab toolbox for solving differential equations on manifolds. *Applied numerical mathematics*, 39(3-4):323–347, 2001.

[37] M. B. Giles. Multilevel Monte Carlo path simulation. *Oper. Res.*, 56(3):607–617, 2008.

[38] M. B. Giles. Multilevel Monte Carlo methods. *Acta Numer.*, 24:259–328, 2015.

[39] E. Grong, H. Z. Munthe-Kaas, and J. Stava. Post-Lie algebra structure of manifolds with constant curvature and torsion. *arXiv preprint arXiv:2305.02688*, 2023.

[40] 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.

[41] M. Hairer and D. Kelly. Geometric versus non-geometric rough paths. *Ann. Inst. Henri Poincaré Probab. Stat.*, 51(1):207–251, 2015.

[42] E. P. Hsu. *Stochastic analysis on manifolds*, volume 38 of *Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI, 2002.

[43] B. Huguet. *Stochastic calculus on manifold and application to functional inequalities*. PhD thesis, Université de Bordeaux, 2020.

[44] B. Huguet. Intertwining relations for diffusions in manifolds and applications to functional inequalities. *Stoch. Proc. and their App.*, 145:1–28, 2022.

[45] B. Huguet. Poincaré inequalities and integrated curvature-dimension criterion for generalised cauchy and convex measures. *Bernoulli*, 30(3):2207–2227, 2024.

[46] A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna. Lie-group methods. In *Acta numerica, 2000*, volume 9 of *Acta Numer.*, pages 215–365. Cambridge Univ. Press, Cambridge, 2000.

[47] A. Iserles, G. R. W. Quispel, and P. S. P. Tse. B-series methods cannot be volume-preserving. *BIT Numer. Math.*, 47(2):351–378, 2007.

[48] 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.

[49] K. Kamm, S. Pagliarani, and A. Pascucci. On the stochastic Magnus expansion and its application to SPDEs. *J. Sci. Comput.*, 89(3):Paper No. 56, 31, 2021.

[50] A. Laurent. *Algebraic Tools and Multiscale Methods for the Numerical Integration of Stochastic Evolutionary Problems*. PhD thesis, University of Geneva, 2021.

[51] A. Laurent. The Lie derivative and Noether’s theorem on the aromatic bicomplex for the study of volume-preserving numerical integrators. *J. Comput. Dyn.*, 11(1):10–22, 2024.

[52] A. Laurent, R. I. McLachlan, H. Z. Munthe-Kaas, and O. Verdier. The aromatic bicomplex for the description of divergence-free aromatic forms and volume-preserving integrators. *Forum Math. Sigma*, 11:Paper No. e69, 2023.

[53] A. Laurent and H. Z. Munthe-Kaas. The universal equivariance properties of exotic aromatic B-series. *Foundations of Computational Mathematics*, pages 1–32, 2024.

[54] A. Laurent and G. Vilmart. Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs. *Math. Comp.*, 89(321):169–202, 2020.

[55] A. Laurent and G. Vilmart. Order conditions for sampling the invariant measure of ergodic stochastic differential equations on manifolds. *Found. Comput. Math.*, 22(3):649–695, 2022.
[56] B. Leimkuhler and C. Matthews. Rational construction of stochastic numerical methods for molecular sampling. *Appl. Math. Res. Express. AMRX*, 2013(1):34–56, 2013.

[57] A. Lejay. Constructing general rough differential equations through flow approximations. *Electron. J. Probab.*, 27:Paper No. 7, 24, 2022.

[58] 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.

[59] T. Lelièvre, M. Rousset, and G. Stoltz. *Free energy computations*. Imperial College Press, London, 2010. A mathematical perspective.

[60] X.-M. Li. Stochastic differential equations on noncompact manifolds: moment stability and its topological consequences. *Probab. Theory Relat. Fields*, 100:417–428, 1994.

[61] Y. Li, Y. Sheng, and R. Tang. Post-hopf algebras, relative rota-baxter operators and solutions to the yang-baxter equation. *Journal of Noncommutative Geometry*, 18(2):605–630, 2023.

[62] A. Lundervold and H. Munthe-Kaas. Hopf algebras of formal diffeomorphisms and numerical integration on manifolds. In *Combinatorics and physics*, volume 539 of *Contemp. Math.*, pages 295–324. Amer. Math. Soc., Providence, RI, 2011.

[63] S. J. A. Malham and A. Wiese. Stochastic Lie group integrators. *SIAM J. Sci. Comput.*, 30(2):597–617, 2008.

[64] G. Marjanovic, M. J. Piggott, and V. Solo. A simple approach to numerical methods for stochastic differential equations in Lie groups. In *2015 54th IEEE Conference on Decision and Control (CDC)*, pages 7143–7150. IEEE, 2015.

[65] G. Marjanovic and V. Solo. Numerical methods for stochastic differential equations in matrix Lie groups made simple. *IEEE Trans. Automat. Control*, 63(12):4035–4050, 2018.

[66] 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.

[67] R. I. McLachlan, K. Modin, H. Munthe-Kaas, and O. Verdier. Butcher series: a story of rooted trees and numerical methods for evolution equations. *Asia Pac. Math. Newsl.*, 7(1):1–11, 2017.

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

[69] H. Munthe-Kaas. Lie-Butcher theory for Runge-Kutta methods. *BIT Numer. Math.*, 35:572–587, 1995.

[70] H. Munthe-Kaas. Runge-Kutta methods on Lie groups. *BIT Numer. Math.*, 38:92–111, 1998.

[71] H. Munthe-Kaas. High order Runge-Kutta methods on manifolds. *Applied Numerical Mathematics*, 29(1):115–127, 1999.

[72] H. Munthe-Kaas. Geometric integration on symmetric spaces. *J. Comput. Dyn.*, 11(1):43–58, 2024.

[73] H. Munthe-Kaas and J. Stava. Lie admissible triple algebras: The connection algebra of symmetric spaces. *Submitted*, 2023.

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

[75] H. Munthe-Kaas and O. Verdier. Integrators on homogeneous spaces: Isotropy choice and connections. *Foundations of Computational Mathematics*, 16:899–939, 2016.

[76] H. Munthe-Kaas and A. Zanna. Numerical integration of differential equations on homogeneous manifolds. In *Foundations of Computational Mathematics: Selected Papers of a Conference Held at Rio de Janeiro, January 1997*, pages 305–315. Springer, 1997.
[77] H. Z. Munthe-Kaas and A. Lundervold. On post-Lie algebras, Lie-Butcher series and moving frames. *Foundations of Computational Mathematics*, 13:583-613, 2013.

[78] H. Z. Munthe-Kaas and W. M. Wright. On the Hopf algebraic structure of Lie group integrators. *Found. Comput. Math.*, 8(2):227-257, 2008.

[79] J.-M. Oudom and D. Guin. On the Lie enveloping algebra of a pre-Lie algebra. *J. K-Theory*, 2(1):147-167, 2008.

[80] B. Owren. Order conditions for commutator-free Lie group methods. *Journal of Physics A: Mathematical and General*, 39(19):5585, 2006.

[81] B. Owren and A. Marthinsen. Runge-Kutta methods adapted to manifolds and based on rigid frames. *BIT Numer. Math.*, 39(1):116-142, 1999.

[82] L. Rahm. An operadic approach to substitution in Lie-Butcher series. *Forum Math. Sigma*, 10:Paper No. e20, 29, 2022.

[83] A. Rößler. Runge-Kutta methods for Stratonovich stochastic differential equation systems with commutative noise. In *Proceedings of the 10th International Congress on Computational and Applied Mathematics (ICCAM-2002)*, volume 164/165, pages 613–627, 2004.

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

[85] T. Shardlow. Modified equations for stochastic differential equations. *BIT Numer. Math.*, 46(1):111-125, 2006.

[86] D. Talay. Second-order discretization schemes of stochastic differential systems for the computation of the invariant law. *Stochastics: An International Journal of Probability and Stochastic Processes*, 29(1):13–36, 1990.

[87] 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.

[88] J. Thompson. Brownian motion and the distance to a submanifold. *Potential Anal.*, 45(3):485-508, 2016.

[89] J. L. Vázquez. *The porous medium equation. Mathematical theory*. Oxford Math. Monogr. Oxford: Oxford University Press, 2007.

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

[91] Z. Wang, Q. Ma, Z. Yao, and X. Ding. The Magnus expansion for stochastic differential equations. *J. Nonlinear Sci.*, 30(1):419-447, 2020.

[92] 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.
# Appendices

## A Order conditions for third weak order

For the sake of completion, we detail in Table 2 the 95 additional order conditions for weak order 3 with stochastic frozen flow methods (2.8). Recall that these conditions are not independent and satisfy shuffle relations (see Remark 2.14).

<table><thead><tr><th>Exotic forest π</th><th>a(π)</th><th>e(π)</th></tr></thead><tbody><tr><td>??</td><td>z<sup>0</sup><sub>i,k3</sub> Z<sup>0</sup><sub>i,j,k2</sub> Z<sup>0</sup><sub>j,k,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> z<sup>0</sup><sub>i,k3</sub> Z<sup>0</sup><sub>i,j,k2</sub> Z<sup>0</sup><sub>i,k,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> z<sup>0</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> Z<sup>d1</sup><sub>j,k2</sub> Z<sup>d1</sup><sub>j,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k3</sub> z<sup>0</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> Z<sup>d1</sup><sub>j,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k3≥k1</sub> z<sup>0</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> Z<sup>d1</sup><sub>j,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2≥k3</sub> z<sup>0</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> Z<sup>d1</sup><sub>i,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k3≥k2</sub> z<sup>0</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> Z<sup>d1</sup><sub>i,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k3≥k1≥k2</sub> z<sup>0</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> Z<sup>d1</sup><sub>i,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2≥k3≥k4</sub> z<sup>0</sup><sub>i,k5</sub> Z<sup>d2</sup><sub>i,k4</sub> Z<sup>d2</sup><sub>i,k3</sub> Z<sup>d1</sup><sub>i,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2≥k3≥k4</sub> z<sup>0</sup><sub>i,k5</sub> Z<sup>d2</sup><sub>i,k4</sub> Z<sup>d1</sup><sub>i,k3</sub> Z<sup>d2</sup><sub>i,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2≥k3≥k4</sub> z<sup>0</sup><sub>i,k5</sub> Z<sup>d1</sup><sub>i,k4</sub> Z<sup>d2</sup><sub>i,k3</sub> Z<sup>d2</sup><sub>i,k2</sub> Z<sup>d1</sup><sub>i,k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> Z<sup>d1</sup><sub>j,k4</sub> Z<sup>0</sup><sub>i,k3</sub> z<sup>0</sup><sub>j,k2</sub> z<sup>0</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k2≥k1</sub> Z<sup>0</sup><sub>j,k,k3</sub> z<sup>0</sup><sub>j,k2</sub> z<sup>0</sup><sub>i,k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2,k3≥k4</sub> Z<sup>d1</sup><sub>j,k4</sub> Z<sup>0</sup><sub>i,k3</sub> z<sup>0</sup><sub>j,k2</sub> z<sup>0</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k2≥k1,k3≥k4</sub> Z<sup>d1</sup><sub>j,k4</sub> Z<sup>0</sup><sub>i,k3</sub> z<sup>0</sup><sub>j,k2</sub> z<sup>0</sup><sub>i,k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> Z<sup>d1</sup><sub>j,k4</sub> Z<sup>0</sup><sub>i,k3</sub> z<sup>0</sup><sub>j,k2</sub> z<sup>0</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> Z<sup>d1</sup><sub>j,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k2≥k1</sub> Z<sup>d1</sup><sub>j,k4</sub> Z<sup>0</sup><sub>i,j,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2,k3≥k4</sub> Z<sup>0</sup><sub>i,j,k4</sub> Z<sup>d1</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> Z<sup>d1</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> Z<sup>d1</sup><sub>i,k4</sub> Z<sup>0</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k2≥k1,k3≥k4</sub> Z<sup>0</sup><sub>i,j,k4</sub> Z<sup>d1</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2,k3≥k4≥k5</sub> Z<sup>d2</sup><sub>i,k5</sub> Z<sup>d2</sup><sub>i,k4</sub> Z<sup>d1</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2,k3≥k4≥k5</sub> Z<sup>d2</sup><sub>i,k5</sub> Z<sup>d1</sup><sub>i,k4</sub> Z<sup>d2</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2</sub> Z<sup>d2</sup><sub>i,k5</sub> Z<sup>d2</sup><sub>i,k4</sub> Z<sup>d1</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k2≥k1,k3≥k4≥k5</sub> Z<sup>d2</sup><sub>i,k5</sub> Z<sup>d1</sup><sub>i,k4</sub> Z<sup>d2</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k2≥k1,k3≥k4≥k5</sub> Z<sup>d1</sup><sub>i,k5</sub> Z<sup>d2</sup><sub>i,k4</sub> Z<sup>d2</sup><sub>i,k3</sub> z<sup>0</sup><sub>i,k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>1/3</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2≥k3≥k4</sub> z<sup>0</sup><sub>k1,k3</sub> z<sup>0</sup><sub>j,k2</sub> z<sup>0</sup><sub>i,k1</sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>1</sup><sub>k1≥k2≥k3≥k4</sub> Z<sup>0</sup><sub>i,j,k4</sub> z<sup>0</sup><sub>i,k3</sub> z<sup>d1</sup><sub>k2</sub> z<sup>d1</sup><sub>k1</sub></td><td>1/6</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂</sub> Z<sup>0</sup><sub>i,j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₁≥k₂</sub> Z<sup>0</sup><sub>i,j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/6</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃</sub> Z<sup>d₁</sup><sub>j,k₄</sub> z<sup>0</sup><sub>j,k₃</sub> z<sup>0</sup><sub>i,k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₂≥k₁≥k₃</sub> Z<sup>d₁</sup><sub>j,k₄</sub> z<sup>0</sup><sub>j,k₃</sub> z<sup>0</sup><sub>i,k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂</sub> Z<sup>d₁</sup><sub>j,k₄</sub> z<sup>0</sup><sub>j,k₃</sub> z<sup>0</sup><sub>i,k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₁≥k₂</sub> Z<sup>d₁</sup><sub>j,k₄</sub> z<sup>0</sup><sub>j,k₃</sub> z<sup>0</sup><sub>i,k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₂≥k₃≥k₁</sub> Z<sup>d₁</sup><sub>j,k₄</sub> z<sup>0</sup><sub>j,k₃</sub> z<sup>0</sup><sub>i,k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₂≥k₁</sub> Z<sup>d₁</sup><sub>j,k₄</sub> z<sup>0</sup><sub>j,k₃</sub> z<sup>0</sup><sub>i,k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₂</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/6</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₂</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₁≥k₂,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₂</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₁</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₁</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₁</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₁</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₁≥k₂,k₄≥k₅</sub> Z<sup>d₂</sup><sub>i,k₅</sub> Z<sup>d₁</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>2/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₁≥k₂,k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> Z<sup>d₂</sup><sub>i,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>0</sup><sub>j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/6</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂≥k₄≥k₅</sub> Z<sup>0</sup><sub>j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₃≥k₁≥k₂≥k₄≥k₅</sub> Z<sup>0</sup><sub>j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/6</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂≥k₄≥k₅</sub> Z<sup>0</sup><sub>j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₃≥k₂≥k₄≥k₅</sub> Z<sup>0</sup><sub>j,k₄</sub> z<sup>0</sup><sub>i,k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/6</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₄≥k₂≥k₃≥k₄</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₄≥k₂≥k₃≥k₄</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/3</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>d₁</sup><sub>i,k₅</sub> z<sup>0</sup><sub>j,k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₂</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>1/6</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>0</sup><sub>i,k₅</sub> z<sup>d₂</sup><sub>k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
    <tr>
      <td>??</td>
      <td>Σ<sub>k₁≥k₂≥k₃≥k₄≥k₅</sub> Z<sup>0</sup><sub>i,k₅</sub> z<sup>d₂</sup><sub>k₄</sub> z<sup>d₂</sup><sub>k₃</sub> z<sup>d₁</sup><sub>k₂</sub> z<sup>d₁</sup><sub>k₁</sub></td>
      <td>0</td>
    </tr>
  </tbody>
</table>
<table><caption>Table 2: Order conditions of frozen flow methods for weak order 3. The expectations and sums that do not satisfy inequalities are omitted for conciseness.</caption><tbody><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>5</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>5</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>5</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>5</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>5</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>5</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>5</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>5</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>5</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>5</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>5</sub>≥k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>5</sub>≥k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>5</sub>≥k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub></sub> z<sup>0</sup><sub>i,k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>1/6</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr><tr><td>??</td><td>Σ<sup>!</sup><sub>k<sub>1</sub>≥k<sub>2</sub>≥k<sub>3</sub>≥k<sub>4</sub>≥k<sub>5</sub>≥k<sub>6</sub></sub> z<sup>d<sub>3</sub></sup><sub>k<sub>6</sub></sub> z<sup>d<sub>2</sub></sup><sub>k<sub>5</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>4</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>3</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>2</sub></sub> z<sup>d<sub>1</sub></sup><sub>k<sub>1</sub></sub></td><td>0</td></tr></tbody></table>