3]http://www.math.kit.edu/ianm3/ jahnke/

# Splitting integrators for stochastic Lie–Poisson systems

Charles-Edouard Bréhier  Charles-Edouard Bréhier
Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille
Jordan, 43 blvd. du 11 novembre 1918, F–69622 Villeurbanne cedex, France
<EMAIL_ADDRESS>http://math.univ-lyon1.fr/ brehier/ ,  David Cohen
David Cohen
Mathematical Sciences, Chalmers University of Technology and University of
Gothenburg, SE–41296 Gothenburg, Sweden<EMAIL_ADDRESS>http://www.math.chalmers.se/ cohend/ and  Tobias Jahnke  Tobias Jahnke
Department of Mathematics, Institute for Applied and Numerical Mathematics,
Karlsruhe Institute of Technology, Englerstr. 2, D–76131 Karlsruhe, Germany
<EMAIL_ADDRESS>[

###### Abstract.

We study stochastic Poisson integrators for a class of stochastic Poisson
systems driven by Stratonovich noise. Such geometric integrators preserve
Casimir functions and the Poisson map property. For this purpose, we propose
explicit stochastic Poisson integrators based on a splitting strategy, and
analyse their qualitative and quantitative properties: preservation of Casimir
functions, existence of almost sure or moment bounds, asymptotic preserving
property, and strong and weak rates of convergence. The construction of the
schemes and the theoretical results are illustrated through extensive
numerical experiments for three examples of stochastic Lie–Poisson systems,
namely: stochastically perturbed Maxwell–Bloch, rigid body and sine–Euler
equations.

###### Key words and phrases:

stochastic Poisson systems; splitting schemes; Poisson integrators; strong and
weak rates of convergence; asymptotic preserving schemes; Maxwell–Bloch
equations; rigid body equations; sine–Euler equations

##  1\. Introduction

Hamiltonian ordinary differential equations and their generalisation, Poisson
systems, are extensively used as mathematical models to describe the dynamical
evolution of various physical systems in science and engineering. The flow of
a Hamiltonian system is known to be a symplectic map, whereas the flow of a
Poisson system is a Poisson map. This is a key property that the flow of a
numerical method should also have. The recent years have thus witnessed a
large amount of research activities in the design and numerical analysis of
symplectic numerical schemes, resp. Poisson integrators, for deterministic
(non-canonical) Hamiltonian systems, see for instance the classical monographs
[46, 30, 24, 9] and references therein.

This research has naturally come to the realm of stochastic Hamiltonian
systems. Without being too exhaustive, we mention the works [37, 36, 53, 10,
41, 12, 33, 51, 5, 20, 47, 54, 26, 25, 31, 15] on the numerical analysis of
symplectic methods for stochastic Hamiltonian systems.

Since symplectic methods for stochastic Hamiltonian systems offer advantages
compared to standard numerical methods, as observed in the above list of
references, it is natural to ask if one can derive numerical integrators
respecting the structure of stochastic Poisson systems of the Stratonovich
form

$\left\\{\begin{aligned} \,\mathrm{d}y(t)&=B(y(t))\nabla
H(y(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}B(y(t))\nabla\widehat{H}_{k}(y(t))\circ\,\,\mathrm{d}W_{k}(t)\\\
y(0)&=y_{0},\end{aligned}\right.$

with Hamiltonian functions
$H,\widehat{H}_{1},\ldots,\widehat{H}_{m}\colon\mathbb{R}^{d}\to\mathbb{R}$, a
structure matrix $B\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$, and
independent standard real-valued Wiener processes $W_{1},\ldots,W_{m}$, see
Section 2 for details on the notation.

Stochastic Poisson systems are popular models to describe diverse random
phenomena, see below and [8, 39, 40, 32, 28, 29, 6, 50, 1, 17, 52] for
instance. However, to the best of our knowledge, there has been no general
study of integrators for stochastic Poisson systems which respect their
geometric properties in the literature so far except the recent work [27]. In
this manuscript, we intend to fill this gap and we study the notion of
stochastic Poisson integrators (see Definition 3 and [27, Theorem 3.1]): such
integrators need to be Poisson maps (see Definition 2) and to preserve the
Casimir functions of the system. Imposing these conditions is natural: indeed
we prove that the flow of the stochastic Poisson system is a Poisson map (see
Theorem 1) and also preserves Casimir functions. In addition, the present
notion of stochastic Poisson integrators is a natural generalisation of the
notion of Poisson integrators for deterministic Poisson systems.

The main contribution of this manuscript is the analysis of a class of
explicit stochastic Poisson integrators, see equation (3.1.2), based on a
splitting strategy. The splitting strategy is often applicable for stochastic
Lie–Poisson systems, which have a structure matrix $B(y)$ which depends
linearly on $y$. The construction of the scheme is illustrated for stochastic
perturbations of three systems which have been studied extensively in the
deterministic case: Maxwell–Bloch, rigid body and sine–Euler equations. Note
that these examples give stochastic differential equations (SDEs) with non-
globally Lipschitz drift and diffusion coefficients, thus standard explicit
schemes such as the Euler–Maruyama method are not expected to converge
(strongly or weakly) or even to satisfy moment bounds. Instead, under
appropriate assumptions, we prove that the proposed integrators converge
strongly and weakly, with rates $1/2$ and $1$ respectively, see Theorem 2.
Indeed, if one assumes that the system admits a Casimir function with compact
level sets (which is the case for the rigid body and the sine–Euler
equations), both the exact and the numerical solutions of the stochastic
Poisson systems remain bounded almost surely, uniformly with respect to the
time step size. Our main convergence result, Theorem 2, is illustrated with
extensive numerical experiments.

On top of that, we study the properties of the stochastic Poisson systems (see
Subsection 2.5) and stochastic Poisson integrators (see Subsection 3.3) in a
multiscale regime, namely when the Wiener processes are approximated by a
smooth noise. The proposed splitting schemes are asymptotic preserving in this
diffusion approximation regime, in the sense of the notion recently introduced
in [11]. This property, which is not satisfied by standard integrators, is
illustrated with numerical experiments.

Let us now compare our contributions with existing works. As already
mentioned, the notion of stochastic Poisson integrators is a natural
generalisation of the notion of Poisson integrators for deterministic systems.
In the stochastic case, we are only aware of the recent work [27], where
techniques which differ from ours are employed. First, in [27], the proof that
the flow is a Poisson map consists in using the Darboux–Lie theorem to rewrite
the stochastic Poisson system into a canonical form, i. e. as a stochastic
Hamiltonian system, for which it is already known that the flow is a
symplectic map. On the contrary, our approach to prove Theorem 1 below is more
direct and extends the approach considered in [24, Chapter VII] for the
deterministic case. Second, the authors of [27] design stochastic Poisson
integrators by starting from a stochastic symplectic scheme for the canonical
version, and then by coming back to the original variables. Note that the
transformations between the non canonical and canonical variables are often
found by solving partial differential equations, and that symplectic schemes
are usually implicit. Our approach is more direct and leads to explicit
splitting schemes. In particular, for the stochastic rigid body system, the
scheme proposed in [27] is based on the midpoint rule and is thus implicit,
whereas the scheme proposed in our work is explicit, and we are able to prove
strong and weak convergence results.

The paper is organized as follows. Section 2 is devoted to the setting and to
the description of the main properties of stochastic Poisson systems, namely
the preservation of Casimir functions and the Poisson map property (Theorem 1,
see Subsection 2.4). The three main examples of stochastic Lie–Poisson systems
(Maxwell–Bloch, rigid body and sine–Euler equations) are introduced in
Subsection 2.3. The diffusion approximation regime is presented in Subsection
2.5. Section 3 presents the main theoretical contributions of this work: we
introduce the notion of stochastic Poisson integrators (Definition 3) and we
propose a class of such integrators using a splitting technique. The main
convergence result (Theorem 2) is stated and proved in Subsection 3.2 (using
an auxiliary result proved in Appendix A): under appropriate assumptions, the
proposed explicit splitting stochastic Poisson integrators converge in strong
and weak senses, with orders $1/2$ and $1$ respectively. The asymptotic
preserving property in the diffusion approximation regime is studied in
Section 3.3. Finally, Section 4 presents numerical experiments using the
proposed splitting stochastic Poisson integrators and variants for the three
examples of stochastic Lie–Poisson systems (Maxwell–Bloch, rigid body and
sine–Euler equations). We illustrate various qualitative and quantitative
properties, which show the superiority of the proposed schemes compared with
existing methods.

##  2\. Stochastic Poisson and Lie–Poisson systems

In this section, we set notation and introduce the stochastic differential
equations studied in this article, namely _stochastic (Lie–)Poisson systems_.
We then state the main properties of such systems and give several examples
for which stochastic Poisson integrators are designed and analysed in Sections
3 and 4. We conclude this section with a diffusion approximation result
justifying why considering stochastic Poisson systems with a Stratonovich
interpretation of the noise is relevant.

###  2.1. Setting and stochastic Poisson dynamics

Let $d,m$ be positive integers: $d$ is the dimension of the considered system
and $m$ is the dimension of the stochastic perturbation. We study _stochastic
Poisson systems_ of the type

(1) $\left\\{\begin{aligned} \,\mathrm{d}y(t)&=B(y(t))\nabla
H(y(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}B(y(t))\nabla\widehat{H}_{k}(y(t))\circ\,\,\mathrm{d}W_{k}(t),\\\
y(0)&=y_{0},\end{aligned}\right.$

with _Hamiltonian functions_
$H,\widehat{H}_{1},\ldots,\widehat{H}_{m}\colon\mathbb{R}^{d}\to\mathbb{R}$,
with _structure matrix_ $B\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$, and
with independent standard real-valued Wiener processes $W_{1},\ldots,W_{m}$
defined on a probability space $(\Omega,\mathcal{F},\mathbb{P})$. The noise in
the SDE (1) is understood in the Stratonovich sense. The initial value $y_{0}$
is assumed to be non-random for ease of presentation, but the results of this
paper can be extended to the case of random $y_{0}$ (independent of
$W_{1},\ldots,W_{m}$ and satisfying appropriate moment bounds).

Henceforth we assume at least that $H\in\mathcal{C}^{2}$,
$\widehat{H}_{1},\ldots,\widehat{H}_{m}\in\mathcal{C}^{3}$, and that
$B\in\mathcal{C}^{2}$. The gradient is denoted by $\nabla$, e.g. $\nabla
H(y)=\bigl{(}\frac{\partial H(y)}{\partial y_{1}},\ldots,\frac{\partial
H(y)}{\partial y_{d}})\in\mathbb{R}^{d}$. The structure matrix $B$ is assumed
to satisfy the following properties.

  * •

Skew-symmetry: for every $y\in\mathbb{R}^{d}$ and for all
$i,j\in\\{1,\ldots,d\\}$, one has

$B_{ij}(y)=-B_{ji}(y);$

  * •

Jacobi identity: for every $y\in\mathbb{R}^{d}$ and for all
$i,j,k\in\\{1,\ldots,d\\}$, one has

$\sum_{\ell=1}^{d}\left(\frac{\partial B_{ij}(y)}{\partial
y_{\ell}}B_{lk}(y)+\frac{\partial B_{jk}(y)}{\partial
y_{\ell}}B_{li}(y)+\frac{\partial B_{ki}(y)}{\partial
y_{\ell}}B_{lj}(y)\right)=0.$

Sometimes the structure matrix $B$ is referred to as the Poisson matrix. In
many applications the structure matrix $B$ depends linearly on $y$: if there
is a family of real numbers $\bigl{(}b_{ij}^{k}\bigr{)}_{1\leq i,j,k\leq d}$
such that

(2) $\displaystyle B_{ij}(y)=\sum_{k=1}^{d}b_{ji}^{k}y_{k}$

for all $y\in\mathbb{R}^{d}$ and $i,j=1,\ldots,d$, then the system (1) is
called a stochastic _Lie–Poisson system_. Examples are provided below in
Section 2.3. A stochastic _Hamiltonian_ system is obtained if $d$ is even and
$B(y)=J^{-1}$ for all $y\in\mathbb{R}^{d}$, where

$J=\begin{pmatrix}0&Id\\\ -Id&0\end{pmatrix}.$

If $\widehat{H}_{k}=0$ for all $k=1,\ldots,m$, then the SDE (1) reduces to a
classical deterministic (Lie–)Poisson or Hamiltonian system; cf. [24].

Properties and numerical approximations of stochastic Hamiltonian systems and
of deterministic Poisson systems have been extensively studied in the
literature, see the references in the introduction. The results presented in
this work are generalisations to the above stochastic Poisson case, with a
special focus on stochastic Lie–Poisson systems.

Under the previous regularity assumptions, the drift coefficient $y\mapsto
B(y)\nabla H(y)$ is of class $\mathcal{C}^{1}$, and, for all $k=1,\ldots,m$,
the diffusion coefficient $y\mapsto B(y)\nabla\widehat{H}_{k}(y)$ is of class
$\mathcal{C}^{2}$. As a consequence, the stochastic differential equation (1)
is locally well-posed: for any deterministic initial condition
$y_{0}\in\mathbb{R}^{d}$, there exists a random time $\tau$, which is almost
surely positive, such that (1) admits a unique solution $t\in[0,\tau)\mapsto
y(t)$ with $y(0)=y_{0}$. Below we will present a criterion to ensure global
well-posedness ($\tau=\infty$ almost surely for any initial condition
$y_{0}$). This criterion is applied to study the examples presented below.

###  2.2. Properties of stochastic Poisson systems

Deterministic and stochastic Poisson systems have several geometric properties
which we discuss in this section. Let
$\mathcal{H}\colon\mathbb{R}^{d}\to\mathbb{R}$ be a mapping of class
$\mathcal{C}^{2}$. The evolution of $\mathcal{H}(y)$ along a solution $y(t)$
of the stochastic Poisson system (1) is described by

(3)
$\displaystyle\,\mathrm{d}\mathcal{H}(y(t))=\\{\mathcal{H},H\\}(y(t))\,\mathrm{d}t+\sum_{k=1}^{m}\\{\mathcal{H},\widehat{H}_{k}\\}(y(t))\circ\,\,\mathrm{d}W_{k}(t),$

where the Poisson bracket $\\{\cdot,\cdot\\}$ associated with the structure
matrix $B$ is defined by

(4) $\displaystyle\\{F,G\\}(y)=\displaystyle\sum_{i,j=1}^{d}\frac{\partial
F(y)}{\partial y_{i}}B_{ij}(y)\frac{\partial G(y)}{\partial y_{j}}=\nabla
F(y)^{T}B(y)\nabla G(y).$

The identity (3) is proved using the chain rule for solutions of SDEs written
in the Stratonovich formulation. The fact that the same structure matrix $B$
appears in both the deterministic and stochastic parts of the system (1) is
important to express (3) using the Poisson bracket defined by (4), which
depends only on $B$ but not on the Hamiltonian functions
$H,\widehat{H}_{1},\ldots,\widehat{H}_{m}$. This assumption on the system (1)
is the key to study the geometric properties of such a system and its
numerical discretisation, such as the preservation of Casimir functions or the
Poisson map property. Most of the properties stated below would not hold if
different structure matrices were considering in the stochastic terms. If
$\\{\mathcal{H},H\\}=0$ and if the system is deterministic (i.e. if
$\widehat{H}_{k}=0$ for all $k=1,\ldots,m$), then the equality (3) implies
that $t\mapsto\mathcal{H}(y(t))$ is constant, i.e. that $\mathcal{H}$ is
preserved by the flow of the deterministic Poisson system. Since every smooth
Hamiltonian has the property that $\\{H,H\\}=0$ this means, in particular,
that the flow of a deterministic Poisson systems $\dot{y}=B(y)\nabla H(y)$
preserves the Hamiltonian $H$ (see for instance [24, Sect. IV.1 and VII.2]).
In the stochastic case, however, the Hamiltonian is in general not preserved.
Precisely, Equation (3) yields the following sufficient condition

$\\{\mathcal{H},H\\}=\\{\mathcal{H},\widehat{H}_{1}\\}=\ldots=\\{\mathcal{H},\widehat{H}_{m}\\}=0$

to obtain $\,\mathrm{d}\mathcal{H}(y(t))=0$ and hence preservation of
$\mathcal{H}$ by the flow of the stochastic Poisson system (1).

In addition, deterministic and stochastic Poisson systems may have conserved
quantities called Casimir functions.

######  Definition 1.

A function $C\colon\mathbb{R}^{d}\to\mathbb{R}$ of class $\mathcal{C}^{2}$ is
called a _Casimir_ function of the stochastic Poisson system (1) if for all
$y\in\mathbb{R}^{d}$ one has

$\nabla C(y)^{T}B(y)=0.$

Observe that the definition of a Casimir function for stochastic and
deterministic Poisson systems only depends on the structure matrix $B$, but
not on the Hamiltonian functions $H,\widehat{H}_{1},\ldots,\widehat{H}_{m}$. A
Casimir function $C$ satisfies

$\\{C,H\\}=\\{C,\widehat{H}_{1}\\}=\ldots=\\{C,\widehat{H}_{m}\\}=0,$

see the definition of the Poisson bracket in equation (4). As a consequence,
owing to (3), any Casimir function $C$ is preserved by the flow of the
stochastic Poisson system (1), i.e. $C(y(t))=C(y_{0})$ for all $t\in[0,\tau]$,
independently of the choice of the Hamiltonian functions
$H,\widehat{H}_{1},\ldots,\widehat{H}_{m}$ (since the same structure matrix
$B$ appears in both the deterministic and stochastic parts of (1) in order to
have preservation of Casimir functions). The preservation of Casimir functions
is a desirable feature for a numerical method when applied to the problem (1).

A criterion to ensure global well-posedness of the dynamics can be stated
based on the preservation of Casimir functions by solutions of stochastic
Poisson systems: it suffices to assume the existence of a Casimir function
with compact level sets.

######  Proposition 1.

Assume that the stochastic Poisson system (1) admits a Casimir function $C$
such that for all $c\in\mathbb{R}$ the level sets
$\\{y\in\mathbb{R}^{d}\,\colon\,C(y)=c\\}$ are compact. Then for any initial
condition $y_{0}\in\mathbb{R}^{d}$ the SDE (1) admits a unique global solution
$\bigl{(}y(t)\bigr{)}_{t\geq 0}$, with $y(0)=y_{0}$, such that almost surely
one has, for all $t\geq 0$, $C(y(t))=C(y_{0})$ and

$\|y(t)\|\leq R(y_{0})=\max_{y\in\mathbb{R}^{d},C(y)=C(y_{0})}\|y\|.$

###### Proof.

The proof of Proposition 1 follows from a straightforward truncation argument:
let $R=R(y_{0})+1$, and introduce mappings $H^{R}$ of class $\mathcal{C}^{2}$
and
$\widehat{H}_{1}^{R},\ldots,\widehat{H}_{m}^{R}:\mathbb{R}^{d}\to\mathbb{R}$,
of class $\mathcal{C}^{3}$, with compact support included in the ball
$\\{y\in\mathbb{R}^{d}\,\colon\,\|y\|\leq R\\}$, and such that
$H^{R}(y)=H(y)$, $\widehat{H}_{k}^{R}(y)=\widehat{H}_{k}(y)$ for all $y$ with
$\|y\|\leq R(y_{0})$ and $k=1,\ldots,m$. Set $f^{R}(y)=B(y)\nabla H^{R}(y)$
and $\hat{f}_{k}^{R}(y)=B(y)\widehat{H}_{k}^{R}(y)$ for all
$y\in\mathbb{R}^{d}$ and $k=1,\ldots,m$. Then $f^{R}$ is globally Lipschitz
continuous and, for all $k=1,\ldots,m$, $\hat{f}_{k}^{R}$ is of class
$\mathcal{C}^{2}$, bounded and with bounded derivatives. By the standard well-
posedness result for SDEs with globally Lipschitz continuous nonlinearities
(when written in Itô form), the SDE

$\,\mathrm{d}y^{R}(t)=f^{R}(y^{R}(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}\hat{f}_{k}^{R}(y^{R}(t))\circ\,\,\mathrm{d}W_{k}(t)$

admits a unique global solution $\bigl{(}y^{R}(t)\bigr{)}_{t\geq 0}$ with
$y^{R}(0)=y_{0}$. Due to the discussion above, this solution preserves the
Casimir function $C$: $C(y^{R}(t))=C(y^{R}(0))=C(y_{0})$ for all $t\geq 0$.
Since the level sets of the Casimir function $C$ are assumed to be compact, by
the definition of $R(y_{0})$, one has $\|y^{R}(t)\|\leq R(y_{0})<R$ for all
$t\geq 0$. This yields the equalities $f^{R}(y^{R}(t))=B(y^{R}(t))\nabla
H(y^{R}(t))$ and
$\hat{f}_{k}^{R}(y^{R}(t))=B(y^{R}(t))\widehat{H}_{k}(y^{R}(t))$ for all
$t\geq 0$ and $k=1,\ldots,m$. Thus $\bigl{(}y^{R}(t)\bigr{)}_{t\geq 0}$ is in
fact a global solution of the SDE (1) (without truncation parameter $R$). This
concludes the proof of the existence of a global solution. Since the
uniqueness is a consequence of the local well-posedness of (1) (by local
Lipschitz continuity of the nonlinearities), the sketch of proof of
Proposition 1 is completed. ∎

To conclude this subsection, we would like to remark that the above analysis
of the preservation properties of stochastic Poisson systems is only valid
when considering stochastic differential equations with multiplicative noise
interpreted in the Stratonovich sense. Other behaviours are observed for Itô
SDEs, see for instance [14, 18], where Hamiltonian and Poisson Itô SDEs and
their numerical discretisations are studied. Indeed, if the noise is
interpreted in the Itô sense, the Hamiltonian function $H$ or the Casimir
functions $C$ are not preserved. Instead, one observes a linear drift in the
expectation of these quantities, which is due to the second-order contribution
appearing when using Itô’s formula instead of the chain rule. In the sequel,
we only study stochastic Poisson systems (1) with noise interpreted in the
Stratonovich sense. We refer to Subsection 2.5 below, where the relevance of
considering the Stratonovich interpretation is justified by a diffusion
approximation result.

###  2.3. Examples of stochastic Poisson systems

In this subsection, we first give an example of a stochastic Poisson systems
which is not a stochastic Lie–Poisson system. We then provide three examples
of stochastic Lie–Poisson systems for which stochastic Poisson integrators
based on a splitting strategy are designed and studied in Sections 3 and 4.

######  Example 1.

Stochastic Lotka–Volterra system. A two-dimensional stochastic Lotka–Volterra
system of the form

(5) $\,\mathrm{d}\begin{pmatrix}y_{1}\\\
y_{2}\end{pmatrix}=\begin{pmatrix}0&y_{1}y_{2}\\\
-y_{1}y_{2}&0\end{pmatrix}\left(\nabla
H(y)\,\,\mathrm{d}t+\sum_{k=1}^{m}\nabla\widehat{H}_{k}(y)\circ\,\,\mathrm{d}W_{k}\right)$

with the structure matrix $B(y_{1},y_{2})=\begin{pmatrix}0&y_{1}y_{2}\\\
-y_{1}y_{2}&0\end{pmatrix}$ and the Hamiltonian function
$H(y_{1},y_{2})=y_{1}-\ln(y_{1})+y_{2}-2\ln(y_{2})$ gives a stochastic Poisson
system (with arbitrary Hamiltonian functions
$\widehat{H}_{1},\ldots,\widehat{H}_{m}$). The stochastic prey-predator model
studied in [45] is obtained taking $m=1$ and
$\widehat{H}_{1}(y)=-\ln(y_{1})+\ln(y_{2})$.

Observe that the stochastic Lotka-Volterra system (5) does not admit Casimir
functions and that it is not a stochastic Lie–Poisson system (since the
mapping $y\mapsto B(y)$ is quadratic).

Random perturbations of higher dimensional Lotka–Volterra systems, see [24]
for deterministic problems and [22] for Itô SDEs, could also fit in the
general framework presented above.

Let us now present the three examples of stochastic Lie–Poisson systems for
which numerical integrators are built and studied in this article.

######  Example 2.

Stochastic Maxwell–Bloch system. Let $d=3$. The deterministic Maxwell–Bloch
equations from laser-matter dynamics read

$\left\\{\begin{aligned} \dot{y}_{1}&=y_{2}\\\ \dot{y}_{2}&=y_{1}y_{3}\\\
\dot{y}_{3}&=-y_{1}y_{2}.\end{aligned}\right.$

This system is relevant to study phenomena of self-induced transparency in
lasers, see for instance [19, 44]. This system is a deterministic Lie–Poisson
system with Poisson matrix, Hamiltonian and Casimir functions given by

$B(y)=\begin{pmatrix}0&-y_{3}&y_{2}\\\ y_{3}&0&0\\\
-y_{2}&0&0\end{pmatrix},\quad H(y)=\frac{1}{2}y_{1}^{2}+y_{3},\quad
C(y)=\frac{1}{2}(y_{2}^{2}+y_{3}^{2}),$

respectively, for all $y=(y_{1},y_{2},y_{3})\in\mathbb{R}^{3}$.

In this article, we consider the following stochastic version of the Maxwell-
Bloch system:

(6) $\,\mathrm{d}y=B(y)\left(\nabla
H(y)\,\,\mathrm{d}t+\sigma_{1}\nabla\widehat{H}_{1}(y)\circ\,\,\mathrm{d}W_{1}(t)+\sigma_{3}\nabla\widehat{H}_{3}(y)\circ\,\,\mathrm{d}W_{3}(t)\right),$

where $\widehat{H}_{1}(y)=\frac{1}{2}y_{1}^{2}$ and
$\widehat{H}_{3}(y)=y_{3}$, $\sigma_{1},\sigma_{3}\geq 0$, driven by two
independent Wiener processes $W_{1}$ and $W_{3}$. Observe that the Casimir
function $C$ does not have compact level sets in this example. The criterion
given in Proposition 1 thus cannot be applied to ensure well-posedness of the
stochastic differential equation (6). In addition, the theoretical strong and
weak convergence results stated below cannot be applied to this example.
However, it is legitimate to introduce a stochastic Poisson integrator and
investigate its behaviour with numerical experiments; this will be presented
in Section 4 below.

######  Example 3.

Stochastic rigid body system. Let $d=3$. The equations governing the
deterministic rigid body motion read

$\left\\{\begin{aligned} \dot{y}_{1}&=(I_{3}^{-1}-I_{2}^{-1})y_{3}y_{2}\\\
\dot{y}_{2}&=(I_{1}^{-1}-I_{3}^{-1})y_{1}y_{3}\\\
\dot{y}_{3}&=(I_{2}^{-1}-I_{1}^{-1})y_{2}y_{1},\end{aligned}\right.$

where the unknown $y=(y_{1},y_{2},y_{3})\in\mathbb{R}^{3}$ represents the
angular momentum in the body frame and the positive distinct real numbers
$I_{1},I_{2},I_{3}$ are referred to as the principal moments of inertia, see
for instance [24, Sect. VII.5].

The system above is a deterministic Lie–Poisson system, with $d=3$, where the
Poisson matrix is given by

$B(y)=\begin{pmatrix}0&-y_{3}&y_{2}\\\ y_{3}&0&-y_{1}\\\
-y_{2}&y_{1}&0\end{pmatrix}$

and the Hamiltonian function is given by

$H(y)=\frac{1}{2}\left(\frac{y_{1}^{2}}{I_{1}}+\frac{y_{2}^{2}}{I_{2}}+\frac{y_{3}^{2}}{I_{3}}\right),$

for all $y\in\mathbb{R}^{3}$. The system admits the quadratic Casimir function
given by

$C(y)=y_{1}^{2}+y_{2}^{2}+y_{3}^{2},$

for all $y\in\mathbb{R}^{3}$.

In this article, we consider the following stochastic version of the rigid
body system

$\displaystyle\,\mathrm{d}\begin{pmatrix}y_{1}\\\ y_{2}\\\ y_{3}\end{pmatrix}$
$\displaystyle=B(y)\left(\nabla
H(y)\,\,\mathrm{d}t+\nabla\widehat{H}_{1}(y)\circ\,\,\mathrm{d}W_{1}(t)+\nabla\widehat{H}_{2}(y)\circ\,\,\mathrm{d}W_{2}(t)\right.$
(7)
$\displaystyle\quad+\left.\nabla\widehat{H}_{3}(y)\circ\,\,\mathrm{d}W_{3}(t)\right),$

where $\widehat{H}_{k}(y)=\frac{y_{k}^{2}}{\widehat{I}_{k}}$, for $k=1,2,3$,
with $\widehat{I}_{1},\widehat{I}_{2},\widehat{I}_{3}$ positive and pairwise
distinct real numbers. The system (3) is a stochastic Lie–Poisson system.
Owing to Proposition 1, it is globally well-posed, since the Casimir function
$C$ has compact level sets.

Observe that the system (3) is a generalisation of the stochastic rigid body
motion equations studied in [17, 6, 4, 16], where $m=1$. It would be
straightforward to adapt the results presented below concerning the
discretisation of (3) to the case $m=1$, this is left to the interested
reader. In the sequel, we only consider the case $m=3$, i. e. the system (3).

######  Example 4.

Stochastic sine–Euler system. The sine–Euler equations consist of a finite-
dimensional truncation of the two-dimensional Euler equations in fluid
dynamics. These equations were first proposed in [55] and yield the
deterministic Lie–Poisson system

$\dot{\omega}_{\bf m}=\sum_{n_{1},n_{2}=-M,{\bf n}\neq{\bf
0}}^{M}\frac{\sin(\frac{2\pi}{N}{\bf m}\times{\bf n})}{|{\bf
n}|^{2}}\omega_{{\bf m}+{\bf n}}\omega_{-\bf n},$

where $M$ is an arbitrary positive integer. In the current example, all
indices are understood modulo $N=2M+1$, and one sets ${\bf m}\times{\bf
n}=m_{1}n_{2}-m_{2}n_{1}$ for ${\bf n}=(n_{1},n_{2})$ and ${\bf
m}=(m_{1},m_{2})$.

The unknown is $\omega=(\omega_{\bf n})_{n_{1},n_{2}=-M}^{M}$, where the
complex numbers $\omega_{\bf n}$ satisfy the Hermitian symmetry property
$\omega_{-\bf n}=\omega_{\bf n}^{\star}$ (where $w^{\star}$ denotes the
complex conjugate of a complex number $w$), and $\omega_{(0,0)}=0$. The
dimension of the system is $d=N^{2}-1=(2M+1)^{2}-1$, but under the Hermitian
symmetry property there are only $d/2$ independent complex-valued components.

In this article, we consider the case $M=1$, thus $N=3$ and $d=8$. The unknown
is written as

$\omega=(\omega_{(1,0)},\omega_{(1,1)},\omega_{(0,1)},\omega_{(-1,1)},\omega_{(1,0)}^{\star},\omega_{(1,1)}^{\star},\omega_{(0,1)}^{\star},\omega_{(-1,1)}^{\star}).$

Introduce the fundamental cell

$K=\\{(0,1),(0,-1),(1,1),(-1,-1),(1,0),(-1,0),(-1,1),(1,-1)\\}.$

Then the deterministic sine–Euler system above can be written as a Lie–Poisson
system (see for instance [21, 34])

(8) $\dot{\omega}=B(\omega)\nabla H(\omega),$

where the Poisson matrix is given by

$B(\omega)=\frac{\sqrt{3}}{2}\begin{pmatrix}0&\omega_{(-1,1)}&\omega_{(1,1)}&\omega_{(0,1)}&0&-\omega_{(0,1)}^{*}&-\omega_{(-1,1)}^{*}&-\omega_{(1,1)}^{*}\\\
-\omega_{(-1,1)}&0&\omega_{(-1,1)}^{*}&-\omega_{(0,1)}^{*}&\omega_{(0,1)}&0&-\omega_{(1,0)}&\omega_{(1,0)}^{*}\\\
-\omega_{(1,1)}&-\omega_{(-1,1)}^{*}&0&\omega_{(1,1)}^{*}&\omega_{(-1,1)}&\omega_{(1,0)}^{*}&0&-\omega_{(1,0)}\\\
-\omega_{(0,1)}&\omega_{(0,1)}^{*}&-\omega_{(1,1)}^{*}&0&\omega_{(1,1)}&-\omega_{(1,0)}&\omega_{(1,0)}^{*}&0\\\
0&-\omega_{(0,1)}&-\omega_{(-1,1)}&-\omega_{(1,1)}&0&\omega_{(-1,1)}^{*}&\omega_{(1,1)}^{*}&\omega_{(0,1)}^{*}\\\
\omega_{(0,1)}^{*}&0&-\omega_{(1,0)}^{*}&\omega_{(1,0)}&-\omega_{(-1,1)}^{*}&0&\omega_{(-1,1)}&-\omega_{(0,1)}\\\
\omega_{(-1,1)}^{*}&\omega_{(1,0)}&0&-\omega_{(1,0)}^{*}&-\omega_{(1,1)}^{*}&-\omega_{(-1,1)}&0&\omega_{(1,1)}\\\
\omega_{(1,1)}^{*}&-\omega_{(1,0)}^{*}&\omega_{(1,0)}&0&-\omega_{(0,1)}^{*}&\omega_{(0,1)}&-\omega_{(1,1)}&0\end{pmatrix}$

and with the Hamiltonian function given by

$\displaystyle H(\omega)=\frac{1}{2}\sum_{{\bf n}\in K}\frac{\omega_{\bf
n}\omega_{-\bf n}}{|{\bf n}|^{2}}$
$\displaystyle=\omega_{(1,0)}\omega_{(1,0)}^{*}+\frac{1}{2}\omega_{(1,1)}\omega_{(1,1)}^{*}+\omega_{(0,1)}\omega_{(0,1)}^{*}+\frac{1}{2}\omega_{(-1,1)}\omega_{(-1,1)}^{*}$
$\displaystyle=H_{(1,0)}(\omega)+H_{(1,1)}(\omega)+H_{(0,1)}(\omega)+H_{(-1,1)}(\omega).$

The prefactor $\sqrt{3}/2$ in $B(\omega)$ originates from the equality
$\sin(\pm 2\pi/3)=\sin(\mp 4\pi/3)=\pm\sqrt{3}/2$. Recall that the
deterministic Poisson system preserves the Hamiltonian $H$. In addition, it
admits two Casimir functions, defined by

$C_{1}(\omega)=\frac{1}{2}\sum_{{\bf n}\in K}\omega_{\bf n}\omega_{\bf
n}^{\star}=\omega_{(1,0)}\omega_{(1,0)}^{*}+\omega_{(1,1)}\omega_{(1,1)}^{*}+\omega_{(0,1)}\omega_{(0,1)}^{*}+\omega_{(-1,1)}\omega_{(-1,1)}^{*}$

and

$C_{2}(\omega)=\sum_{{\bf n},{\bf m}\in K}\cos{\left(\frac{2\pi}{3}({\bf
n}\times{\bf m})\right)}\omega_{\bf n}\omega_{\bf m}\omega_{-\bf n-\bf m}.$

In this article, we consider the following stochastic version of the
sine–Euler equations

$\displaystyle\,\mathrm{d}\omega$ $\displaystyle=B(\omega)\left(\nabla
H(\omega)\,\,\mathrm{d}t+\sigma_{(1,0)}\nabla\widehat{H}_{(1,0)}(\omega)\circ\,\,\mathrm{d}W_{(1,0)}(t)+\sigma_{(1,1)}\nabla\widehat{H}_{(1,1)}(\omega)\circ\,\,\mathrm{d}W_{(1,1)}(t)\right.$
(9)
$\displaystyle\quad\left.+\sigma_{(0,1)}\nabla\widehat{H}_{(0,1)}(\omega)\circ\,\,\mathrm{d}W_{(0,1)}(t)+\sigma_{(-1,1)}\nabla\widehat{H}_{(-1,1)}(\omega)\circ\,\,\mathrm{d}W_{(-1,1)}(t)\right),$

where, for $\bf k\in\\{(1,0),(1,1),(0,1),(-1,1)\\}$, the Hamiltonian functions
are $\widehat{H}_{\bf k}(\omega)=H_{\bf k}(\omega)=\frac{\omega_{\bf
k}\omega_{\bf k}^{\star}}{2|\bf k|^{2}}$, $\sigma_{\bf k}\geq 0$ are
nonnegative real numbers, and $W_{\bf k}$ are independent standard Wiener
process.

The SDE (4) is a stochastic Lie–Poisson system. Owing to Proposition 1, it is
globally well-posed, since the Casimir function $C_{1}$ has compact level
sets.

Let us provide a possible physical interpretation of the stochastic sine–Euler
system (4): the Hamiltonian function can be decomposed as
$H(\omega)=H_{(1,0)}(\omega)+H_{(1,1)}(\omega)+H_{(0,1)}(\omega)+H_{(-1,1)}(\omega)$,
where $H_{\bf k}$ can be interpreted as the kinetic energy of the modes, or
periodic waves, which are parallel to the direction $\bf
k\in\\{(1,0),(1,1),(0,1),(-1,1)\\}$. In the stochastic version introduced
above, noise acts independently on each of these modes.

To conclude this subsection on examples of stochastic (Lie–)Poisson systems,
let us mention other systems which may be treated using the techniques
developed in this work: appropriate spatial discretisations of (stochastic)
Vlasov–Poisson equations [49], random perturbations of the full rigid body [1]
(with a rotation matrix giving the orientation of the body in a fixed frame),
stochastic models of fluid dynamics [23], or reduced models based on a
Gaussian wavepacket of the time-dependent $N$-body Schrödinger equation [24,
Chap VII.6.1-VII.6.4]. These examples are not considered in this article and
may be studied in future works.

###  2.4. The Poisson map property

We proceed with showing that the flow of stochastic Poisson systems (1)
satisfies a property which is a generalisation of the symplecticity property
for deterministic, respectively stochastic, Hamiltonian systems, see e.g [24],
respectively [38].

Let $D_{y}$ denote the Jacobian operator, i.e.
$D_{y}f(y)=(\partial_{j}f_{i}(y))_{1\leq i,j\leq d}$ for a smooth function
$f\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$. The transpose of a matrix $M$ is
denoted by $M^{T}$.

######  Definition 2.

Let $U\subset\mathbb{R}^{d}$ be an open set. A transformation $\varphi\colon
U\to\mathbb{R}^{d}$ is called a _Poisson map_ for the problem (1), if one has,
almost surely, for all $y\in\mathbb{R}^{d}$,

$D_{y}\varphi(y)B(y)D_{y}\varphi(y)^{T}=B(\varphi(y)).$

Observe that a composition of Poisson maps is a Poisson map. This property
will be used in the design of stochastic Poisson splitting integrators in the
next sections.

The main result of this section is stated below.

######  Theorem 1.

Introduce the flow $(t,y)\mapsto\varphi_{t}(y)$ of the stochastic Poisson
system (1) with coefficients of class $\mathcal{C}^{3}$. Assume that the flow
is globally well defined and of class $\mathcal{C}^{1}$ with respect to the
variable $y$. Then, for all $t\geq 0$, $\varphi_{t}$ is a Poisson map: almost
surely, for all $y\in\mathbb{R}^{d}$, one has

$D_{y}\varphi_{t}(y)B(y)D_{y}\varphi_{t}(y)^{T}=B(\varphi_{t}(y)).$

This result can be compared with the similar statement in [27, Th. 2.1]. To
prove their result, the authors of [27] use the Darboux–Lie theorem to perform
a change of coordinates, and replace the considered stochastic Poisson system
by a stochastic canonical Hamiltonian system. Our strategy to prove Theorem 1
is more direct and mimics the analysis of the deterministic case, as proposed
in [24, Ex. 6.VII.4].

###### Proof.

For all $t\geq 0$ and all $y\in\mathbb{R}^{d}$, let
$\Phi_{t}(y)=D_{y}\varphi_{t}(y)$. To simplify notation, $\varphi_{t}$ stands
for $\varphi_{t}(y)$ in the computations below. For convenience, we set
$\widehat{H}_{0}=H$ in this proof.

First, it is well-known that $t\mapsto\Phi_{t}(y)$ is the solution of the
variational equation

$\,\mathrm{d}\Phi_{t}=D_{y}\bigl{(}B(\varphi_{t})\nabla
H(\varphi_{t})\bigr{)}\Phi_{t}\,\mathrm{d}t+\sum_{k=1}^{m}D_{y}\bigl{(}B(\varphi_{t})\nabla\widehat{H}_{k}(\varphi_{t})\bigr{)}\Phi_{t}\circ\,\mathrm{d}W_{k}(t)$

with $\Phi_{0}=Id$, see for instance [7, Theorem 2.3.32]. We claim that, in
the case of the stochastic Poisson system (1), the variational equation above
may be rewritten as

(10)
$\,\mathrm{d}\Phi_{t}=\mathcal{S}_{0}(\varphi_{t})\Phi_{t}\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\Phi_{t}\circ\,\mathrm{d}W_{k}(t),$

where, for all $k=0,\ldots,m$, one has

$\displaystyle\mathcal{S}_{k}(\varphi_{t})$
$\displaystyle=\mathcal{M}_{k}(\varphi_{t})+\mathcal{L}_{k}(\varphi_{t})$
$\displaystyle\mathcal{M}_{k}(\varphi_{t})$
$\displaystyle=\left[\big{(}\partial_{1}B(\varphi_{t})\big{)}\nabla\widehat{H}_{k}(\varphi_{t})\mid\ldots\mid\big{(}\partial_{d}B(\varphi_{t})\big{)}\nabla\widehat{H}_{k}(\varphi_{t})\right]$
$\displaystyle\mathcal{L}_{k}(\varphi_{t})$
$\displaystyle=B(\varphi_{t})\nabla^{2}\widehat{H}_{k}(\varphi_{t}).$

The identification of the matrices $\mathcal{S}_{k}(\varphi_{t})$ above is
based on the following computation: if $[M]_{i,j}$ denotes the $(i,j)$-th
entry of a matrix $M$, then applying the chain rule (recall that the SDE is
interpreted in the Stratonovich sense) yields

$\displaystyle\left[D_{y}\big{(}B(\varphi_{t})\nabla\widehat{H}_{k}(\varphi_{t})\big{)}\right]_{i,j}=\sum_{l=1}^{d}\frac{\partial}{\partial
y_{j}}\big{(}B_{il}(\varphi_{t})\partial_{l}\widehat{H}_{k}(\varphi_{t})\big{)}$
$\displaystyle=\sum_{l=1}^{d}\left(\sum_{n=1}^{d}\big{(}\partial_{n}B_{il}(\varphi_{t})\big{)}\left[\frac{\partial\varphi_{t}}{\partial
y_{j}}\right]_{n}\right)\partial_{l}\widehat{H}_{k}(\varphi_{t})+\sum_{l=1}^{d}B_{il}(\varphi_{t})\sum_{n=1}^{d}\big{(}\partial_{l}\partial_{n}\widehat{H}_{k}(\varphi_{t})\big{)}\left[\frac{\partial\varphi_{t}}{\partial
y_{j}}\right]_{n}$
$\displaystyle=\sum_{n=1}^{d}\left(\sum_{l=1}^{d}\big{(}\partial_{n}B_{il}(\varphi_{t})\big{)}\partial_{l}\widehat{H}_{k}(\varphi_{t})\right)\left[\frac{\partial\varphi_{t}}{\partial
y_{j}}\right]_{n}+\sum_{n=1}^{d}\left(\sum_{l=1}^{d}B_{il}(\varphi_{t})\partial_{l}\partial_{n}\widehat{H}_{k}(\varphi_{t})\right)\left[\frac{\partial\varphi_{t}}{\partial
y_{j}}\right]_{n}$
$\displaystyle=\sum_{n=1}^{d}\Big{[}\big{(}\partial_{n}B(\varphi_{t})\big{)}\nabla\widehat{H}_{k}(\varphi_{t})\Big{]}_{i}[\Phi_{t}]_{n,j}+\sum_{n=1}^{d}\Big{[}B(\varphi_{t})\nabla^{2}\widehat{H}_{k}(\varphi_{t})\Big{]}_{i,n}[\Phi_{t}]_{n,j}.$

Let us now define

(11) $\displaystyle\delta(t):=\Phi_{t}B(y)\Phi_{t}^{T}-B(\varphi_{t}).$

Since $\varphi_{0}(y)=y$ for all $y\in\mathbb{R}^{d}$ and $\Phi_{0}=Id$, one
has $\delta(0)=0$. To prove that $\varphi_{t}$ is a Poisson map for all $t\geq
0$ almost surely, it suffices to check that $\delta(t)=0$ for all $t\geq 0$.

We will show that $t\mapsto\delta(t)$ is a solution of the linear equation

(12) $\displaystyle\,\mathrm{d}\delta(t)$
$\displaystyle=\bigg{(}\mathcal{S}_{0}(\varphi_{t})\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\circ\,\mathrm{d}W_{k}(t)\bigg{)}\delta(t)$
$\displaystyle\qquad+\delta(t)\bigg{(}\mathcal{S}_{0}(\varphi_{t})\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\circ\,\mathrm{d}W_{k}(t)\bigg{)}^{T}.$

Since the initial value is $\delta(0)=0$, by uniqueness of the solution we
obtain $\delta(t)=0$ for all $t\geq 0$, thus $\varphi_{t}$ is a Poisson map
for all $t\geq 0$.

It remains to prove that $\delta$ is indeed a solution of (12). On the one
hand, from the variational equation (10) for $\Phi_{t}$ above, applying the
product rule yields

(13) $\displaystyle\,\mathrm{d}\delta(t)$
$\displaystyle=\bigg{(}\mathcal{S}_{0}(\varphi_{t})\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\circ\,\mathrm{d}W_{k}(t)\bigg{)}\Phi_{t}B(y)\Phi_{t}^{T}$
$\displaystyle\qquad+\Phi_{t}B(y)\Phi_{t}^{T}\bigg{(}\mathcal{S}_{0}(\varphi_{t})\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\circ\,\mathrm{d}W_{k}(t)\bigg{)}^{T}-\,\mathrm{d}B(\varphi_{t}).$

On the other hand, one has the identity

(14) $\displaystyle\,\mathrm{d}B(\varphi_{t})$
$\displaystyle=\bigg{(}\mathcal{S}_{0}(\varphi_{t})\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\circ\,\mathrm{d}W_{k}(t)\bigg{)}B(\varphi_{t})$
$\displaystyle\qquad+B(\varphi_{t})\bigg{(}\mathcal{S}_{0}(\varphi_{t})\,\mathrm{d}t+\sum_{k=1}^{m}\mathcal{S}_{k}(\varphi_{t})\circ\,\mathrm{d}W_{k}(t)\bigg{)}^{T}.$

Combining (13) and (14) provides the claim that $t\mapsto\delta(t)$ is
solution of (12). The proof of the identity (14) requires to exploit the
assumptions on the structure matrix $B$ as follows. First, note that the
$(i,j)$-th entry of the left-hand side of (14) satisfies

$\displaystyle\big{[}\,\mathrm{d}B(\varphi_{t})\big{]}_{i,j}$
$\displaystyle=\nabla B_{ij}^{T}(\varphi_{t}){\circ\,}\,\mathrm{d}\varphi_{t}$
$\displaystyle=\nabla B_{ij}^{T}(\varphi_{t})\left(B(\varphi_{t})\nabla
H(\varphi_{t})\,\,\mathrm{d}t+\sum_{k=1}^{m}B(\varphi_{t})\nabla\widehat{H}_{k}(\varphi_{t})\circ\,\,\mathrm{d}W_{k}(t)\right),$

using the chain rule. Recall that
$\mathcal{S}_{k}(\varphi_{t})=\mathcal{M}_{k}(\varphi_{t})+\mathcal{L}_{k}(\varphi_{t})$.

On the one hand, all the terms involving
$\mathcal{L}_{k}(\varphi_{t})=B(\varphi_{t})\nabla^{2}\widehat{H}_{k}(\varphi_{t})$
appearing on the right-hand side of (14) vanish: this follows from the
symmetry of the Hessians $\nabla^{2}\widehat{H}_{k}$ and the skew-symmetry of
$B$ which yields (omitting the argument “$(\varphi_{t})$” everywhere)

$\displaystyle\mathcal{L}_{k}B+B\mathcal{L}_{k}^{T}$
$\displaystyle=B\nabla^{2}\widehat{H}_{k}B+B\nabla^{2}\widehat{H}_{k}B^{T}=0$

for $k=0,1,\ldots,m$.

On the other hand, using the skew-symmetry of $B$ and the Jacobi identity (and
omitting “$(\varphi_{t})$” again), one has

$\displaystyle\Big{[}\mathcal{S}_{k}B+B\mathcal{S}_{k}^{T}\Big{]}_{ij}$
$\displaystyle=\Big{[}\mathcal{M}_{k}B+B\mathcal{M}_{k}^{T}\Big{]}_{ij}$
$\displaystyle=\sum_{l=1}^{d}\big{[}\mathcal{M}_{k}\big{]}_{il}B_{lj}+\sum_{l=1}^{d}B_{il}\big{[}\mathcal{M}_{k}\big{]}_{jl}$
$\displaystyle=\sum_{l=1}^{d}\sum_{n=1}^{d}\big{(}\partial_{l}B_{in}\big{)}\partial_{n}\widehat{H}_{k}B_{lj}+\sum_{l=1}^{d}B_{il}\sum_{n=1}^{d}\big{(}\partial_{l}B_{jn}\big{)}\partial_{n}\widehat{H}_{k}$
$\displaystyle=\sum_{n=1}^{d}\sum_{l=1}^{d}\Big{(}\big{(}\partial_{l}B_{in}\big{)}B_{lj}+\big{(}\partial_{l}B_{jn}\big{)}B_{il}\Big{)}\partial_{n}\widehat{H}_{k}$
$\displaystyle=-\sum_{n=1}^{d}\sum_{l=1}^{d}\Big{(}\big{(}\partial_{l}B_{ni}\big{)}B_{lj}+\big{(}\partial_{l}B_{jn}\big{)}B_{li}\Big{)}\partial_{n}\widehat{H}_{k}$
$\displaystyle=\sum_{n=1}^{d}\sum_{l=1}^{d}\Big{(}\big{(}\partial_{l}B_{ij}\big{)}B_{ln}\Big{)}\partial_{n}\widehat{H}_{k}$
$\displaystyle=\nabla B_{ij}^{T}B\nabla\widehat{H}_{k}$

for all $k=0,\ldots,m$.

This concludes the proof of (14), which as already explained above gives the
identity $\delta(t)=0$ for all $t\geq 0$. In the end, this concludes the proof
that $\varphi_{t}$ is a Poisson map for all $t\geq 0$, almost surely.

∎

One of the objectives of this work is to design and study integrators for the
stochastic Poisson system (1), which preserve its geometric structure, namely
the Poisson map property (Definition 2 and Theorem 1), and the preservation of
the Casimir functions. This leads to define and analyze so-called stochastic
Poisson integrators, see Section 3.

###  2.5. Stochastic Poisson systems obtained by diffusion approximation

The goal of this subsection is to describe a class of multiscale stochastic
systems, depending on a parameter $\epsilon\in(0,1)$, such that the stochastic
Poisson system (1) is obtained as a limit when $\epsilon\to 0$. In particular,
this approximation result justifies why considering stochastic Poisson systems
with a Stratonovich interpretation of the noise is relevant.

For all $\epsilon\in(0,1)$, introduce the multiscale system

(15) $\left\\{\begin{aligned}
\,\mathrm{d}y^{\epsilon}(t)&=B(y^{\epsilon}(t))\nabla
H(y^{\epsilon}(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}B(y^{\epsilon}(t))\nabla\widehat{H}_{k}(y^{\epsilon}(t))\frac{\xi_{k}^{\epsilon}(t)}{\epsilon}\,\mathrm{d}t,\\\
\,\mathrm{d}\xi_{k}^{\epsilon}(t)&=-\frac{\xi_{k}^{\epsilon}(t)}{\epsilon^{2}}\,\mathrm{d}t+\frac{1}{\epsilon}\,\mathrm{d}W_{k}(t),\quad
k=1,\ldots,m,\end{aligned}\right.$

with initial values $y^{\epsilon}(0)=0$ and $\xi_{k}^{\epsilon}(0)=0$, for all
$k=1,\ldots,m$. Note that $\xi_{k}^{\epsilon}$ is an Ornstein–Uhlenbeck
process, for all $k=1,\ldots,m$. Compared with the stochastic Poisson system
(1), $y^{\epsilon}$ solves a random ordinary differential equation, since the
noise $\circ\,\mathrm{d}W_{k}(t)$ is replaced by
$\frac{\xi_{k}^{\epsilon}(t)}{\epsilon}\,\mathrm{d}t$. Observe that if $C$ is
a Casimir function of the stochastic Poisson system (1), then one has
$dC(y^{\epsilon}(t))=0$ owing to the chain rule, thus
$C(y^{\epsilon}(t))=C(y^{\epsilon}(0))$ for all $t\geq 0$. Assuming that the
Casimir function $C$ has compact level sets, like in Proposition 1, ensures
that the system (15) is globally well-posed, for all $\epsilon\in(0,1)$, and
that

$\underset{\epsilon\in(0,1)}{\sup}~{}\underset{t\geq
0}{\sup}~{}\|y^{\epsilon}(t)\|\leq R(y_{0}).$

When $\epsilon\to 0$, one has the following result.

######  Proposition 2.

Assume that the stochastic Poisson system (1) admits a Casimir function $C$
such that for all $c\in\mathbb{R}$, the level sets
$\\{y\in\mathbb{R}^{d}\,\colon\,C(y)=c\\}$ are compact.

For all $t\geq 0$,

$y^{\epsilon}(t)\underset{\epsilon\to 0}{\to}y(t),$

where the convergence holds in distribution ($y^{\epsilon}(t)$ and $y(t)$ are
random variables). In addition, for all $y_{0}\in\mathbb{R}^{d}$, all
$T\in(0,\infty)$ and any function $\phi$ of class $\mathcal{C}^{3}$, there
exists a real number $c(T,y_{0},\phi)\in(0,\infty)$ such that for all
$\epsilon\in(0,1)$

$\underset{0\leq t\leq
T}{\sup}~{}\big{|}\mathbb{E}[\phi(y^{\epsilon}(t))]-\mathbb{E}[\phi(y(t))]\big{|}\leq
c(T,y_{0},\phi)\epsilon.$

Proposition 2 is a diffusion approximation result: roughly, a diffusion
process is the solution of a SDE driven by a Wiener process. Here, the
diffusion process $y$ is approximated by the solutions $y^{\epsilon}$ of ODEs.

The proof of Proposition 2 is omitted. Indeed, this is a standard result in
the literature: we refer for instance to the monograph [43, Chapter 11 and 18]
for a presentation of homogenization techniques, and references therein for a
historical perspective. Proposition 2 fits in the class of Wong–Zakai
approximation results (where a smooth approximation of a Wiener noise leads to
a SDE driven by Stratonovich noise). We also refer to [11, Proposition 2.6]
for a similar statement (with $m=1$), and references therein for ideas of
proof. Note that the weak error estimate can be obtained using a variant of
the proof of [11, Proposition 2.4], decomposing the error in terms of
solutions of Kolmogorov and Poisson equations. The details of the proof are
left to the interested readers.

Let us provide an heuristic argument which justifies the diffusion
approximation result: for all $k=1,\ldots,m$, one has the identity

$\frac{\xi_{k}^{\epsilon}(t)}{\epsilon}\,\mathrm{d}t=\,\mathrm{d}W_{k}(t)-\epsilon\,\mathrm{d}\xi_{k}^{\epsilon}(t).$

Then, the contribution of $\epsilon\,\mathrm{d}\xi_{k}^{\epsilon}(t)$ vanishes
in the limit $\epsilon\to 0$, and only the contribution $\,\mathrm{d}W_{k}(t)$
remains at the limit. There may be different interpretations of the noise at
the limit: at least, Itô and Stratonovich interpretations are possible
candidates. However, recall that Casimir functions $C$ of the stochastic
Poisson system (1) are preserved by the solution $y^{\epsilon}$ of (15), for
all $\epsilon>0$. The Itô interpretation of the noise is not consistent with
this preservation property, whereas the Stratonovich one is, due to the chain
rule. As a consequence, the Stratonovich interpretation is the natural
candidate for the diffusion approximation limit. Checking rigorously that
indeed $y^{\epsilon}(t)\to y(t)$ in distribution requires additional arguments
which are omitted in this work.

######  Remark 1.

Let
$(t,y,\xi_{1},\ldots,\xi_{m})\mapsto\varphi^{\epsilon}(t,y,\xi_{1},\ldots,\xi_{m})$
define the flow map associated with (15). Then for all $t\geq 1$,
$\epsilon\in(0,1)$ and all $\xi_{1},\ldots,\xi_{m}\in\mathbb{R}$, the mapping

$y\in\mathbb{R}^{d}\mapsto\varphi^{\epsilon}(t,y,\xi_{1},\ldots,\xi_{m})$

is a Poisson map in the sense of Definition 2. This may be proved by
modifications of the proof of Theorem 1, using the chain rule. The details are
left to the reader.

The multiscale system (15) has components evolving at different time scales:
the component $y^{\epsilon}$ evolves at a time scale of order ${\rm O}(1)$,
whereas the Ornstein–Uhlenbeck processes
$\xi_{1}^{\epsilon},\ldots,\xi_{m}^{\epsilon}$ evolve at a time scale of order
${\rm O}(\epsilon^{-2})$. The definition of effective integrators for the
multiscale system (15), which avoid prohibitive time step size restrictions of
the type $h={\rm O}(\epsilon^{2})$, and which lead to consistent
discretisation of $y(t)$ when $\epsilon\to 0$, is a crucial and challenging
problem. This question is briefly studied in Section 3.3 below: we define so-
called asymptotic preserving schemes (in the spirit of [11]), employing the
preservation of the geometric structure satisfied by the stochastic Poisson
integrators introduced in the next section.

##  3\. Stochastic Poisson integrators

In Section 2, we have proved that the flow of a stochastic Poisson system of
the type (1) satisfies two key properties: it is a Poisson map (see Theorem 1)
and it preserves Casimir functions which are associated with the structure
matrix $B$. Having the methodology of geometric numerical integration in mind,
this motivates us to introduce the concept of a stochastic Poisson integrator
for the stochastic Poisson system (1), see Definition 3. We then present and
analyse a general strategy to derive efficient stochastic Poisson integrators,
based on a splitting technique, which can be implemented easily for some
stochastic Lie–Poisson systems. We then proceed with a convergence analysis of
the proposed splitting integrators: Theorem 2 states that, under appropriate
conditions, the scheme has in general strong and weak convergence rates equal
to $1/2$ and $1$, respectively. First, the analysis is performed for an
auxiliary problem (20) with globally Lipschitz continuous nonlinearities, see
Lemma 1. Second, if the system admits a Casimir function with compact level
sets, the auxiliary convergence result is applied to get strong and weak error
estimates for the SDE (1). Finally, we show that the proposed stochastic
Poisson integrators based on a splitting technique satisfy an asymptotic
preserving property when considering the multiscale SDE (15) in the diffusion
approximation regime.

###  3.1. Definition and splitting integrators for stochastic (Lie–)Poisson
systems

Let us recall that symplectic, respectively Poisson, integrators preserve the
key features of deterministic and stochastic Hamiltonian systems, respectively
deterministic Poisson systems. Such geometric numerical integrators offer
various benefits over classical time integrators in the deterministic setting,
see for instance [24, 30, 9]. We shall now state the definition and study the
properties of stochastic Poisson integrators for stochastic Poisson systems
(1). On the one hand, this extends the definition and analysis of
deterministic Poisson integrators (see [27, Th. 3.1] for another approach). On
the other hand, this extends the definition and analysis of stochastic
symplectic integrators for stochastic Hamiltonian systems.

We first consider general stochastic Poisson integrators, and then focus the
discussion on a class of splitting integrators.

####  3.1.1. Stochastic Poisson integrators

The following notation is used below. The time step size is denoted by $h>0$.
A numerical scheme is defined as follows: for all $n\geq 1$,

(16) $y^{[n]}=\Phi_{h}(y^{[n-1]},\Delta_{n}W_{1},\ldots,\Delta_{n}W_{m}),$

with Wiener increments $\Delta_{n}W_{k}=W_{k}(nh)-W_{k}((n-1)h)$,
$k=1,\ldots,m$. The Wiener increments are independent centered real-valued
Gaussian random variables with variance $h$. The mapping $\Phi_{h}$ is
referred to as the integrator.

######  Definition 3.

A numerical scheme (16) for the stochastic Poisson system (1) is called a
_stochastic Poisson integrator_ if

  * •

for all $h>0$ and all $\Delta w_{1},\ldots,\Delta w_{m}\in\mathbb{R}$, the
mapping

$y\mapsto\Phi_{h}(y,\Delta w_{1},\ldots,\Delta w_{m})$

is a Poisson map (in the sense of Definition 2),

  * •

if $C$ is a Casimir of the stochastic Poisson system (1), then $\Phi_{h}$
preserves $C$, precisely

$C(\Phi_{h}(y,\Delta w_{1},\ldots,\Delta w_{m}))=C(y)$

for all $y\in\mathbb{R}^{d}$, $h>0$ and $\Delta w_{1},\ldots,\Delta
w_{m}\in\mathbb{R}$.

As in the deterministic case, it is seen that standard integrators like the
Euler–Maruyama scheme are not (stochastic) Poisson integrators. In addition,
it is a difficult task to construct Poisson integrators for the general
Poisson systems, see [24, Chapter VII.4.2] for deterministic problems.
Therefore, the design of stochastic Poisson integrators requires to exploit
the special structure for each considered problem. In this article, we focus
on constructing and analyzing stochastic Poisson integrators for stochastic
Lie–Poisson systems. More precisely, we propose explicit Poisson integrators
for a large class of stochastic Lie–Poisson systems using a splitting
strategy. In Section 4, we will exemplify this strategy for three models
introduced in Section 2: the stochastic Maxwell–Bloch equations (Example 2 and
Subsection 4.1), the stochastic free rigid body equations (3) (Example 3 and
Subsection 4.2), as well as the stochastic sine–Euler equations (Example 4 and
Subsection 4.3).

####  3.1.2. Splitting integrators for stochastic Poisson systems

We first propose an abstract splitting integrator for general stochastic
Poisson systems (1). We then focus on stochastic Lie–Poisson systems (17) and
propose implementable stochastic Poisson integrators for this class of SDEs,
which includes the three examples mentioned above.

The key observation made in [34, p.3044] is that a wide class of deterministic
Lie–Poisson systems can be split into subsystems which are all linear. This
was used in [34] for the construction of very efficient geometric integrators
for deterministic Lie–Poisson systems. Inspired by [34], we propose and
analyse efficient explicit Poisson integrators for stochastic Lie–Poisson
systems. On an abstract level, our splitting approach is not restricted to
Lie–Poisson systems and could also be applied to general stochastic Poisson
systems (1).

Let us consider a stochastic Poisson system of the type (1), and assume that
the Hamiltonian $H$ can be split as follows:

$H=\sum_{k=1}^{p}H_{k}.$

for some $p\geq 1$, where the Hamiltonian functions $H_{1},\ldots,H_{p}$ have
the same regularity as $H$.

To define the abstract splitting schemes for (1), it is convenient to define
the flows associated to the subsystems:

  * •

for each $k=1,\ldots,p$, let
$(t,y)\in\mathbb{R}^{+}\times\mathbb{R}^{d}\mapsto\varphi_{k}(t,y)$ be the
flow associated with the ordinary differential equation
$\dot{y}_{k}=B(y_{k})\nabla H_{k}(y_{k})$;

  * •

for each $k=1,\ldots,m$, let
$(t,y)\in\mathbb{R}\times\mathbb{R}^{d}\mapsto\widehat{\varphi}_{k}(t,y)$ be
the flow associated with the ordinary differential equation
$\dot{y}_{k}=B(y_{k})\nabla\widehat{H}_{k}(y_{k})$.

Note that it is sufficient to consider
$\varphi_{1}(t,\cdot),\ldots,\varphi_{p}(t,\cdot)$ for $t\geq 0$, however the
mappings
$\widehat{\varphi}_{1}(t,\cdot),\ldots,\widehat{\varphi}_{k}(t,\cdot)$ need to
be considered for $t\in\mathbb{R}$.

Below, we shall also use the notation $\exp(hY_{H_{k}})=\varphi_{k}(h,\cdot)$
and $\exp(hY_{\widehat{H}_{k}})=\widehat{\varphi}_{k}(h,\cdot)$, where
$Y_{H_{k}}=B\nabla H_{k}$, resp. $Y_{\widehat{H}_{k}}=B\nabla\widehat{H}_{k}$,
to denote the vector fields of the corresponding differential equations. For
the definition of the splitting integrators below, it is essential to note
that the exact solution of the Stratonovich stochastic differential equation
$\,\mathrm{d}y_{k}=B(y_{k})\nabla\widehat{H}_{k}(y_{k})\circ\,\mathrm{d}W_{k}(t)$
is given by $y_{k}(t)=\widehat{\varphi}_{k}(W_{k}(t),y_{k}(0))$.

As explained above, closed-form expressions for the flows $\varphi_{k}$ and
$\widehat{\varphi}_{k}$ are unknown in general but can be obtained for a wide
class of stochastic Lie–Poisson systems

(17) $\left\\{\begin{aligned} &\,\mathrm{d}y(t)=B(y(t))\nabla
H(y(t))\,\,\mathrm{d}t+B(y(t))\sum_{k=1}^{m}\nabla\widehat{H}_{k}(y(t))\circ\,\,\mathrm{d}W_{k}(t),\\\
&B_{ij}(y)=\sum_{k=1}^{d}b_{ji}^{k}y_{k}\quad\text{for}\quad
i,j=1,\ldots,d,\end{aligned}\right.$

where the structure matrix $B(y)$ depends linearly on $y$. For the examples of
stochastic Lie–Poisson systems introduced in Section 2.3, below we design
explicit splitting schemes which can be easily implemented by a splitting
strategy. In the sequel, we analyse the geometric and convergence properties
of splitting integrators in an abstract framework, where it is not assumed
that the flows $\varphi_{k}$ and $\hat{\varphi}_{k}$ can be computed exactly.
In particular, the assumption that the structure matrix $B$ depends linearly
on $y$ is not required in the analysis. Note also that expressions of the
flows may also be known for some stochastic Poisson systems which are not
Lie–Poisson problems, in which case the abstract analysis would also be
applicable.

We are now in position to define splitting integrators for the stochastic
Poisson system (1), which will be exemplified in the case of stochastic
Lie–Poisson systems (17). This general splitting integrator is given by

$\displaystyle\Phi_{h}(\cdot)$ $\displaystyle=\Phi_{h}(\cdot,\Delta
W_{1},\ldots,\Delta
W_{m})=\exp(hY_{H_{p}})\circ\exp(hY_{H_{p-1}})\circ\ldots\circ\exp(hY_{H_{1}})$
(18) $\displaystyle\circ\exp(\Delta W_{m}Y_{\widehat{H}_{m}})\circ\exp(\Delta
W_{m-1}Y_{\widehat{H}_{m-1}})\circ\ldots\circ\exp(\Delta
W_{1}Y_{\widehat{H}_{1}}).$

It is immediate to check the following fundamental result.

######  Proposition 3.

The splitting integrator (3.1.2) is a stochastic Poisson integrator, in the
sense of Definition 3, for the stochastic Poisson system (1).

###### Proof.

Observe that for any $h>0$ and any real numbers $\Delta w_{1},\ldots,\Delta
w_{m}$, the mapping $\Phi_{h}(\cdot,\Delta w_{1},\ldots,\Delta w_{m})$ is a
composition of flow maps $\varphi_{k}(h,\cdot)$, $k=1,\ldots,p$ and
$\widehat{\varphi}_{k}(\Delta w_{k},\cdot)$, $k=1,\ldots,m$. Owing to Theorem
1, all of these flow maps are Poisson maps (since they are flow maps of either
deterministic or stochastic Poisson systems).

In addition, if $C$ is a Casimir function of the stochastic Poisson system
(1), then $C$ is preserved by each of the flow maps $\varphi_{k}(h,\cdot)$,
$k=1,\ldots,p$ and $\widehat{\varphi}_{k}(\Delta w_{k},\cdot)$,
$k=1,\ldots,m$. Indeed, recall that the definition of a Casimir only depends
on the structure matrix $B$, and not on the Hamiltonian functions, and all the
associated vectors fields are of the type $Y_{H_{k}}=B\nabla H_{k}$ and
$Y_{\widehat{H}_{k}}=B\nabla\widehat{H}_{k}$: the associated flow maps thus
preserve $C$. As a consequence, the general splitting integrator
$\Phi_{h}(\cdot,\Delta w_{1},\ldots,\Delta w_{m})$ also preserves the Casimir
functions $C$ of the stochastic Poisson system (1). This concludes the proof
that the splitting scheme (3.1.2) is a stochastic Poisson integrator. ∎

Before proceeding to the convergence analysis for the splitting integrators
(3.1.2), it is worth exploiting the fact that they are stochastic Poisson
integrators to state that the numerical solution remains bounded if the
considered stochastic Poisson system (1) admits a Casimir function $C$ with
compact level sets. We refer to Proposition 1 for the statement of a similar
result for the solution of the stochastic Poisson system (1), in particular
the assumption on compact level sets.

######  Proposition 4.

Assume that the stochastic Poisson system (1) admits a Casimir function $C$
which has compact level sets. Consider the stochastic Poisson integrator
$y^{[n+1]}=\Phi_{h}(y^{[n]})$ given by (3.1.2). Then, for any initial
condition $y^{[0]}=y_{0}\in\mathbb{R}^{d}$, for all $t\geq 0$, almost surely
one has the following bound for the numerical solution

$\underset{h>0}{\sup}~{}\underset{n\geq 0}{\sup}~{}\|y^{[n]}\|\leq
R(y^{[0]})=\max_{y\in\mathbb{R}^{d},C(y)=C(y^{[0]})}\|y\|.$

###### Proof.

The splitting scheme (3.1.2) is a stochastic Poisson integrator (owing to
Proposition 3), thus it preserves the Casimir function $C$: therefore for all
$n\geq 1$,

$C(y^{[n]})=C(y^{[n-1]})=\ldots=C(y^{[0]}).$

Note that $R(y^{[0]})<\infty$, since by assumption the Casimir function $C$
has compact level sets. Therefore one obtains

$\|y^{[n]}\|\leq R(y^{[0]})$

for all $n\geq 0$ by the definition of $R(y^{[0]})$. This concludes the proof.
∎

######  Remark 2.

The stochastic Poisson integrator (3.1.2) employs a Lie–Trotter splitting
strategy. Changing the orders of integration of the deterministic and
stochastic parts yields the following alternative to (3.1.2)

$\displaystyle\Phi_{h}(\cdot)$ $\displaystyle=\Phi_{h}(\cdot,\Delta
W_{1},\ldots,\Delta W_{m})$ $\displaystyle=$ $\displaystyle\exp(\Delta
W_{m}Y_{\widehat{H}_{m}})\circ\exp(\Delta
W_{m-1}Y_{\widehat{H}_{m-1}})\circ\ldots\circ\exp(\Delta
W_{1}Y_{\widehat{H}_{1}})$
$\displaystyle\circ\exp(hY_{H_{p}})\circ\exp(hY_{H_{p-1}})\circ\ldots\circ\exp(hY_{H_{1}}).$

This alternative scheme is also a stochastic Poisson integrator, which
satisfies Propositions 3 and 4. The theoretical analysis of that scheme and
associated numerical experiments are not reported in the present article.

######  Remark 3.

A numerical method of weak order $2$ can be designed using the strategy
developed in [2]. The integrator is a combination of three mappings and
depends on an additional random variable $\gamma_{n}$, uniformly distributed
in $\\{-1,1\\}$:

(19)
$y^{[n]}=\Phi_{h,\gamma_{n}}(y^{[n-1]})=\Phi_{h/2}^{det,S}\circ\Phi_{h,\gamma_{n}}^{sto}(\cdot,\Delta_{n}W_{1},\ldots,\Delta_{n}W_{m})\circ\Phi_{h/2}^{det,S}(y^{[n-1]}),$

where

$\Phi_{h/2}^{det,S}=\exp(\frac{h}{4}Y_{H_{1}})\circ\ldots\circ\exp(\frac{h}{4}Y_{H_{p-1}})\circ\exp(\frac{h}{2}Y_{H_{p}})\circ\exp(\frac{h}{4}Y_{H_{p-1}})\circ\ldots\circ\exp(\frac{h}{4}Y_{H_{1}})$

is obtained using a Strang splitting integrator with time step size $h/2$ for
the deterministic part of the equation, and

$\Phi_{h,\gamma_{n}}^{sto}(\cdot,\Delta_{n}W_{1},\ldots,\Delta_{n}W_{m})=\begin{cases}\exp(\Delta
W_{m}Y_{\widehat{H}_{m}})\circ\ldots\circ\exp(\Delta
W_{1}Y_{\widehat{H}_{1}}),\quad\gamma_{n}=1\\\ \exp(\Delta
W_{1}Y_{\widehat{H}_{1}})\circ\ldots\circ\exp(\Delta
W_{m}Y_{\widehat{H}_{m}}),\quad\gamma_{n}=-1\end{cases}$

is obtained using a Lie–Trotter splitting integrator
$\Phi_{h,\gamma_{n}}^{sto}$ with time step size $h$ applied to the stochastic
part of the equation, where the order of the integration depends on
$\gamma_{n}$.

It is straightforward to check that the numerical scheme (19) is a stochastic
Poisson integrator, using the same arguments as in the proof of Proposition 3.
Numerical experiments which illustrate the behaviour of this scheme and weak
convergence with order $2$ will be reported below in Section 4. However, we do
not give details concerning the theoretical analysis of the scheme (19).

We also refer to [3, 42] for other possible constructions of higher order
splitting methods for SDEs. Finally, another possible strategy to design
higher order integrators would be to use modified equations, like in [1].

###  3.2. Convergence analysis

The objective of this section is to prove a general strong and weak
convergence result for stochastic Poisson integrators (3.1.2) defined by the
splitting strategy. Note that we assume that the stochastic Poisson system (1)
admits a Casimir function with compact level sets: as explained above, this
condition ensures global well-posedness for the continuous problem, and
provides almost sure bounds for the exact and numerical solutions
(Propositions 1 and 4). As a consequence, the general convergence result can
be applied to get strong and weak convergence rates for the proposed explicit
stochastic Poisson integrator (3.1.2), when applied to the stochastic rigid
body system (Example 3) and to the stochastic sine–Euler system (Example 4),
see Theorems 6 and 7 below respectively. Note that these two SDEs do not have
globally Lipschitz continuous coefficients, so for those examples standard
explicit schemes such as the Euler–Maruyama method may fail to converge
strongly. The fact that the proposed scheme is a stochastic Poisson integrator
is essential to perform the convergence analysis. However, the general
convergence result below cannot be applied to the stochastic Maxwell–Bloch
system – the generalisation of the result to that example is not treated in
the present work.

######  Theorem 2.

Assume that the stochastic Poisson system (1) admits a Casimir function with
compact level sets.

  * Strong convergence

Assume that $B$ is of class $\mathcal{C}^{2}$, that the mappings
$H_{1},\ldots,H_{p}$ are of class $\mathcal{C}^{2}$, and that the mappings
$\widehat{H}_{1},\ldots,\widehat{H}_{m}$ are of class $\mathcal{C}^{3}$. Then
the stochastic Poisson integrator (3.1.2) has strong order of convergence
equal to $1/2$: for all $T\in(0,\infty)$ and all $y_{0}\in\mathbb{R}^{d}$,
there exists a real number $c(T,y_{0})\in(0,\infty)$ such that

$\underset{0\leq n\leq N}{\sup}~{}\left(\mathbb{E}\left[\left\lVert
y\left(nh\right)-y^{[n]}\right\rVert^{2}\right]\right)^{1/2}\leq
c(T,y_{0})h^{\frac{1}{2}},$

with time step size $h=T/N$, and $y^{[0]}=y_{0}=y(0)$.

If $m=1$, then the strong order of convergence is equal to $1$.

  * Weak convergence

Assume that $B$ is of class $\mathcal{C}^{5}$, that the mappings
$H_{1},\ldots,H_{p}$ are of class $\mathcal{C}^{5}$, and that the mappings
$\widehat{H}_{1},\ldots,\widehat{H}_{m}$ are of class $\mathcal{C}^{6}$. Then
the stochastic Poisson integrator (3.1.2) has weak order of convergence equal
to $1$: for all $T\in(0,\infty)$ and all $y_{0}\in\mathbb{R}^{d}$, and any
test function $\phi\colon\mathbb{R}^{d}\to\mathbb{R}$ of class
$\mathcal{C}^{4}$ with bounded derivatives, there exists a real number
$c(T,y_{0},\phi)\in(0,\infty)$ such that

$\underset{0\leq n\leq
N}{\sup}~{}\left|\mathbb{E}\left[\phi\left(y\left(nh\right)\right)\right]-\mathbb{E}\left[\phi\left(y^{[n]}\right)\right]\right|\leq
c(T,y_{0},\phi)h.$

The convergence theorem stated above concerning the strong and weak rates of
convergence of the stochastic Poisson integrator (3.1.2) applied to the
stochastic Poisson system (1) is an immediate consequence of the following
auxiliary result, which is stated for a general SDE of the type

(20)
$\,\mathrm{d}z(t)=\sum_{k=1}^{p}f_{k}(z(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}\widehat{f}_{k}(z(t))\circ\,\,\mathrm{d}W_{k}(t),$

with functions $f_{k}$ and $\widehat{f}_{k}$ which are globally Lipschitz
continuous.

######  Lemma 1.

Consider the auxiliary splitting scheme

(21)
$z^{[n]}=\varphi_{p}(h,\cdot)\circ\ldots\circ\varphi_{1}(h,\cdot)\circ\widehat{\varphi}_{m}(\Delta
W_{m}^{n},\cdot)\circ\ldots\circ\widehat{\varphi}_{1}(\Delta
W_{1}^{n},\cdot)(z^{[n-1]}),$

with $z^{[0]}=z_{0}=z(0)$, associated with the auxiliary SDE (20), where
$\varphi_{k}$ is the flow associated with the ODE $\dot{z}_{k}=f_{k}(z_{k})$,
$k=1,\ldots,p$, and $\widehat{\varphi}_{k}$ is the flow associated with the
ODE $\dot{z}_{k}=\widehat{f}_{k}(z_{k})$.

  * Strong convergence

Assume that the mappings $f_{1},\ldots,f_{p}$ are of class $\mathcal{C}^{1}$
with bounded derivatives, and that the mappings
$\widehat{f}_{1},\ldots,\widehat{f}_{m}$ are bounded and of class
$\mathcal{C}^{2}$ with bounded first and second order derivatives. Then the
auxiliary scheme (21) has strong order of convergence equal to $1/2$: for all
$T\in(0,\infty)$ and all $z_{0}\in\mathbb{R}^{d}$, there exists a real number
$c(T,z_{0})\in(0,\infty)$ such that

(22) $\underset{0\leq n\leq N}{\sup}~{}\left(\mathbb{E}\left[\left\lVert
z\left(nh\right)-z^{[n]}\right\rVert^{2}\right]\right)^{1/2}\leq
c(T,z_{0})h^{\frac{1}{2}}.$

In the commutative noise case, i.e. if
$\widehat{f}_{k}^{\prime}(z)\widehat{f}_{\ell}(z)=\widehat{f}_{\ell}^{\prime}(z)\widehat{f}_{k}(z)$
for all $k,\ell=1,\ldots,m$, the strong order of convergence is equal to $1$.

  * Weak convergence

Assume that the mappings $f_{1},\ldots,f_{p}$ are of class $\mathcal{C}^{4}$
with bounded derivatives, and that the mappings
$\widehat{f}_{1},\ldots,\widehat{f}_{m}$ are bounded and of class
$\mathcal{C}^{5}$ with bounded first and second order derivatives. Then the
auxiliary scheme (21) has weak order of convergence equal to $1$: for all
$T\in(0,\infty)$ and all $z_{0}\in\mathbb{R}^{d}$, and any test function
$\phi:\mathbb{R}^{d}\to\mathbb{R}$ of class $\mathcal{C}^{4}$, there exists a
real number $c(T,z_{0},\phi)\in(0,\infty)$ such that

(23) $\underset{0\leq n\leq
N}{\sup}~{}\left|\mathbb{E}\left[\phi\left(z\left(nh\right)\right)\right]-\mathbb{E}\left[\phi\left(z^{[n]}\right)\right]\right|\leq
c(T,z_{0},\phi)h.$

The proof of Lemma 1 is postponed to Appendix A. Let us now check how Theorem
2 is a straightforward corollary of Lemma 1. Note that if $m=1$, the
commutative noise case condition is satisfied.

###### Proof of Theorem 2.

Owing to Propositions 1 and 4, the exact and numerical solutions of the SDE
(1), resp. scheme (3.1.2), satisfy the almost sure bounds

$\underset{t\in[0,T]}{\sup}~{}\left\lVert y(t)\right\rVert\leq
R(y_{0}),\quad\underset{N\geq 1}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left\lVert y^{[n]}\right\rVert\leq R(y_{0}),$

where $R(y_{0})=\max_{y\in\mathbb{R}^{d},C(y)=C(y_{0})}\|y\|$, and
$R(y_{0})<\infty$ since $C$ has compact level sets by assumption.

Using the same construction as in the proof of Proposition 1, one can define
compactly supported functions $f_{k}$ and $\widehat{f}_{k}$, such that
$f_{k}(y)=B(y)\nabla H_{k}(y)$ and
$\widehat{f}_{k}(y)=B(y)\nabla\widehat{H}_{k}(y)$ for all $y\in\mathbb{R}^{d}$
such that $\|y\|\leq R(y_{0})$. In addition, $f_{k}$ is at least of class
$\mathcal{C}^{1}$ and $\widehat{f}_{k}$ is at least of class
$\mathcal{C}^{2}$.

Note that with this choice, $y(t)=z(t)$ and $y^{[n]}=z^{[n]}$ for all
$t\in[0,T]$ and all $n\in\\{0,\ldots,N\\}$, where $\bigl{(}z(t)\bigr{)}_{t\geq
0}$ is the solution of the auxiliary SDE (20) and
$\bigl{(}z^{[n]}\bigr{)}_{n\geq 0}$ is obtained by the auxiliary scheme (21).
It remains to apply Lemma 1 to conclude. Note also that it is not necessary to
assume that the functions $\phi$, $B$, $H_{1},\ldots,H_{p}$,
$\widehat{H}_{1},\ldots,\widehat{H}_{m}$ and their derivatives are bounded.
This is due to the boundedness of the exact and numerical solutions provided
by the preservation of the Casimir function $C$ and the compact level sets
assumption. ∎

######  Remark 4.

If one considers the following variant of the stochastic Poisson system (1)

$\,\mathrm{d}y(t)=B(y(t))\nabla
H(y(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}B(y(t))\nabla\widehat{H}_{k}(y(t))\circ\,\,\mathrm{d}W(t)$

driven by a single Wiener process $W$ (that is $W_{1}=\ldots=W_{m}=W$), the
associated variant of the proposed stochastic Poisson integrator (3.1.2) reads

$\displaystyle\Phi_{h}$
$\displaystyle=\exp(hY_{H_{p}})\circ\exp(hY_{H_{p-1}})\circ\ldots\circ\exp(hY_{H_{1}})$
$\displaystyle\circ\exp(\Delta WY_{\widehat{H}_{m}})\circ\exp(\Delta
WY_{\widehat{H}_{m-1}})\circ\ldots\circ\exp(\Delta WY_{\widehat{H}_{1}}).$

This scheme is not consistent when $m\geq 2$. In the proof of the convergence
result Theorem 2, more precisely in the proof of Lemma 1, the independence of
the Wiener processes $W_{1},\ldots,W_{m}$ plays a crucial role.

###  3.3. Asymptotic preserving schemes in the diffusion approximation regime

The objective of this section is to use the proposed splitting stochastic
Poisson integrators (3.1.2) in order to define effective numerical schemes for
the discretisation of the multiscale system (15) described in Section 2.5. The
challenge is to obtain a good behaviour of the numerical scheme when
$\epsilon\to 0$. On the one hand, one needs to avoid time step size
restrictions of the type $h={\rm O}(\epsilon)$ or $h={\rm O}(\epsilon^{2})$,
which would be prohibitive when $\epsilon$ is small. On the other hand, it
would be desirable to have a convergence (in distribution) of the type
$y^{\epsilon,[n]}\underset{\epsilon\to 0}{\to}y^{[n]}$, for all fixed $h>0$
and $n\geq 1$, to reproduce the diffusion approximation result Proposition 2
at the discrete time level. Indeed, if the two requirements above are
satisfied, the integrator can be used to approximate both (1) and (15),
without the need to adapt the time step size $h$ when $\epsilon$ vanishes.

The class of numerical methods which satisfy the two requirements above is
known as _asymptotic preserving_ schemes. We refer to the recent work [11]
where asymptotic preserving schemes were introduced for a class of stochastic
differential equations of the type (15). Note that a standard Euler–Maruyama
scheme does not satisfy the asymptotic preserving property. Recall that for
this notion of asymptotic preserving schemes, the convergence is understood in
the sense of convergence in distribution of random variables. Using the
splitting strategy allows us to design other examples of asymptotic preserving
schemes for (15), such that the corresponding limit scheme (obtained when
$\epsilon\to 0$ with fixed time step size $h>0$) is the splitting stochastic
Poisson integrator (3.1.2).

We propose the following integrator for the multiscale system (15): for any
$\epsilon\in(0,1)$ and any time step size $h>0$, for all $n\geq 1$, set

$\displaystyle y^{\epsilon,[n]}$
$\displaystyle=\exp(hY_{H_{p}})\circ\ldots\circ\exp(hY_{H_{1}})$ (24)
$\displaystyle\circ\exp\left(\frac{h\xi_{m}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{m}}\right)\circ\ldots\circ\exp\left(\frac{h\xi_{1}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{1}}\right)(y^{\epsilon,[n-1]}),$

where, for each $k=1,\ldots,m$, the Ornstein–Uhlenbeck process
$\xi_{k}^{\epsilon}$ is discretised using the linear implicit Euler scheme

$\xi_{k}^{\epsilon,[n]}=\xi_{k}^{\epsilon,[n-1]}-\frac{h}{\epsilon^{2}}\xi_{k}^{\epsilon,[n]}+\frac{\Delta_{n}W_{k}}{\epsilon}=\frac{1}{1+\frac{h}{\epsilon^{2}}}\Bigl{(}\xi_{k}^{\epsilon,[n-1]}+\frac{\Delta_{n}W_{k}}{\epsilon}\Bigr{)}.$

Note that $C(y^{\epsilon,[n]})=C(y^{\epsilon,[0]})$ for all $n\geq 0$, if $C$
is a Casimir function of the stochastic Poisson system (1). If $C$ has compact
level sets, this yields the following variant of the bound of Proposition 4,

$\underset{\epsilon\in(0,1)}{\sup}~{}\underset{h>0}{\sup}~{}\underset{n\geq
0}{\sup}~{}\|y^{[n]}\|\leq
R(y^{[0]})=\max_{y\in\mathbb{R}^{d},C(y)=C(y^{[0]})}\|y\|,$

which is uniform over $\epsilon$.

Observe that for all $n\geq 1$ and $h>0$, one has

$\frac{h\xi_{k}^{\epsilon,[n]}}{\epsilon}=\Delta_{n}W_{k}+\epsilon\bigl{(}\xi_{k}^{\epsilon,[n-1]}-\xi_{k}^{\epsilon,[n]}\bigr{)}\underset{\epsilon\to
0}{\to}\Delta_{n}W_{k}.$

By a recursion argument, it is then straightforward to check that

$y^{\epsilon,[n]}\underset{\epsilon\to 0}{\to}y^{[n]},$

for all $n\geq 0$ and for all fixed $h>0$, where $y^{[n]}$ is given by the
splitting scheme (3.1.2). As a consequence, the scheme (3.3) is an asymptotic
preserving scheme, in the sense of [11]: the following diagram commutes

$\begin{CD}y^{\epsilon,[N]}@>{N\to\infty}>{}>y^{\epsilon}(T)\\\
@V{}V{\epsilon\to 0}V@V{}V{\epsilon\to 0}V\\\ y^{[N]}@>{\Delta t\to
0}>{}>y(T)\end{CD}$

when the time step size is given by $h=T/N$. In other words:

  * •

for each fixed $\epsilon>0$, the scheme (3.3) is consistent with (15) when
$h\to 0$,

  * •

for each fixed $h>0$, the proposed scheme (3.3) converges to a limiting scheme
when $\epsilon\to 0$, which is given here by the abstract splitting scheme
(3.1.2),

  * •

the limiting scheme (3.1.2) is consistent when $h\to 0$ with (1), which is the
limit when $\epsilon\to 0$ of (15).

We refer to the recent work [11] for a general analysis of asymptotic
preserving schemes for stochastic differential equations. As explained in
[11], the construction of asymptotic preserving schemes for SDEs, in
particular to obtain equations interpreted in the Stratonovich sense, may be
subtle. Here we do not employ a predictor-corrector strategy as in [11] (which
is used to get the Stratonovich interpretation instead of the Itô one), since
we directly use exact flows of the appropriate subsystems in the splitting
procedure: in the present paper, the Stratonovich interpretation is obtained
in a natural way.

The property of being asymptotic preserving is a qualitative property of a
numerical scheme. Let us now briefly discuss the behaviour of weak error
estimates of the asymptotic preserving scheme (3.3) when $\epsilon$ is small.
For each fixed $\epsilon>0$, it is expected that the proposed asymptotic
preserving scheme (3.3) has a weak order of convergence equal to $1$ in
general: for test functions $\phi\colon\mathbb{R}^{d}\to\mathbb{R}$ of class
$\mathcal{C}^{4}$, one has

$\left|\mathbb{E}[\varphi(y^{\epsilon,[N]})]-\mathbb{E}[\varphi(y^{\epsilon}(T))]\right|\leq
c_{\epsilon}(T,\varphi)h,$

where $h=\frac{T}{N}$ and the real number $c_{\epsilon}(T,\varphi)$ may depend
on $\epsilon$ and diverge when $\epsilon\to 0$. In order to have a
computational cost independent of the parameter $\epsilon$, it would be
desirable to establish that the proposed scheme is _uniformly accurate_ : one
would need to prove error estimates of the type

$\underset{\epsilon\in(0,\epsilon_{0})}{\sup}~{}\big{|}\mathbb{E}[\varphi(y^{\epsilon,[N]})]-\mathbb{E}[\varphi(y^{\epsilon}(T))]\big{|}\leq
c(T,\varphi)h^{\alpha},$

which are uniform with respect to $\epsilon\in(0,\epsilon_{0})$ (with
arbitrary $\epsilon_{0}>0$), in other words $c(T,\varphi)$ is independent of
$\epsilon$. Observe that a reduction of the order of convergence, namely
$\alpha<1$, may happen. Proving the uniform accuracy property of the scheme
(3.3) is beyond the scope of this work. However, in the numerical experiments
reported below, we investigate whether such uniform weak error estimates hold
for the considered problems.

######  Remark 5.

It is possible to define a variant of the asymptotic preserving scheme (3.3),
using a midpoint approximation for the the Ornstein–Ulenbeck components:

$\xi_{k}^{\epsilon,[n]}=\xi_{k}^{\epsilon,[n-1]}-\frac{h}{2\epsilon^{2}}\left(\xi_{k}^{\epsilon,[n-1]}+\xi_{k}^{\epsilon,[n]}\right)+\frac{\Delta_{n}W_{k}}{\epsilon},$

in which case the definition of $y^{\epsilon,[n]}$ needs to be modified as
follows:

$\displaystyle y^{\epsilon,[n]}$
$\displaystyle=\exp(hY_{H_{p}})\circ\exp(hY_{H_{p-1}})\circ\ldots\circ\exp(hY_{H_{1}})$
$\displaystyle\circ\exp\left(\frac{h(\xi_{k}^{\epsilon,[n-1]}+\xi_{k}^{\epsilon,[n]})}{2\epsilon}Y_{\widehat{H}_{m}}\right)\circ\ldots\circ\exp\left(\frac{h(\xi_{k}^{\epsilon,[n-1]}+\xi_{k}^{\epsilon,[n]})}{2\epsilon}Y_{\widehat{H}_{1}})(y^{\epsilon,[n-1]}\right).$

That scheme is also asymptotic preserving.

##  4\. Numerical experiments

In this section, we illustrate the behaviour of the stochastic Poisson
integrators which have been proposed and analysed in Section 3. We choose to
present numerical experiments for the three examples of stochastic Lie–Poisson
systems introduced in Section 2. On the one hand, we illustrate the
qualitative properties of the proposed splitting stochastic Poisson
integrators, compared with standard methods, by considering the temporal
evolution of Casimir functions. On the other hand, we investigate and state
strong and weak orders of convergence (which are consequences of Theorem 2),
and we illustrate the quantitative error estimates obtained above. In
addition, we illustrate the asymptotic preserving property for the multiscale
versions of the considered systems. Note that, in general, the theoretical
convergence results cannot be applied to the stochastic Maxwell–Bloch system
(Example 2), since no Casimir functions with compact level sets is known for
that example. However, the theoretical results can be applied to the
stochastic rigid body system (Example 3) and to the stochastic sine–Euler
system (Example 4).

###  4.1. Explicit stochastic Poisson integrators for stochastic Maxwell–Bloch
equations

This subsection presents explicit stochastic Poisson integrators for the
stochastic Maxwell–Bloch system (6) (Example 2). We first give a detailed
construction of the splitting scheme, which gives a stochastic Poisson
integrator. We then illustrate its qualitative properties (preservation of the
Casimir function) and strong and weak error estimates of the proposed scheme
by numerical experiments. Finally, we illustrate the asymptotic preserving
property (Section 3.3) for a multiscale version of the system.

####  4.1.1. Presentation of the splitting scheme for the stochastic
Maxwell–Bloch system

Recall that the stochastic Maxwell–Bloch system (6) introduced in Example 2 is
of the type

$\,\mathrm{d}y=B(y)\left(\nabla
H(y)\,\,\mathrm{d}t+\sigma_{1}\nabla\widehat{H}_{1}(y)\circ\,\,\mathrm{d}W_{1}(t)+\sigma_{3}\nabla\widehat{H}_{3}(y)\circ\,\,\mathrm{d}W_{3}(t)\right).$

To apply the strategy described in Section 3.1 and construct explicit
stochastic Poisson integrators, we follow the approach from [44] for the
deterministic Maxwell–Bloch system. The Hamiltonian function $H$ is split as
$H=H_{1}+H_{3}$, with $H_{1}(y)=\widehat{H}_{1}(y)=\frac{1}{2}y_{1}^{2}$ and
$H_{3}(y)=\widehat{H}_{3}(y)=y_{3}$. The two associated deterministic
subsystems can be solved exactly as follows. On the one hand, the
deterministic subsystem corresponding with the vector field $Y_{H_{1}}=B\nabla
H_{1}$ is given by

$\left\\{\begin{aligned} \dot{y}_{1}&=0\\\ \dot{y}_{2}&=y_{3}y_{1}\\\
\dot{y}_{3}&=-y_{2}y_{1}.\end{aligned}\right.$

Observe that $y_{1}$ may be considered as a constant and thus $(y_{2},y_{3})$
is solution of a linear ordinary differential equation (it is the standard
harmonic oscillator): the exact solution of the first subsystem is thus given
by

$\exp(tY_{H_{1}})y(0)=\begin{pmatrix}1&0&0\\\
0&\cos(y_{1}(0)t)&\sin(y_{1}(0)t)\\\
0&-\sin(y_{1}(0)t)&\cos(y_{1}(0)t)\end{pmatrix}y(0)$

for all $t\in\mathbb{R}$ and $y(0)\in\mathbb{R}^{3}$.

On the other hand, the deterministic subsystem corresponding with the vector
field $Y_{H_{3}}=B\nabla H_{3}$ is given by

$\left\\{\begin{aligned} \dot{y}_{1}&=y_{2}\\\ \dot{y}_{2}&=0\\\
\dot{y}_{3}&=0.\end{aligned}\right.$

The exact solution of the second subsystem is thus given by

$\exp(tY_{H_{3}})y(0)=\begin{pmatrix}1&t&0\\\ 0&1&0\\\ 0&0&1\end{pmatrix}y(0)$

for all $t\in\mathbb{R}$ and $y(0)\in\mathbb{R}^{3}$.

In the case of the stochastic Maxwell–Bloch system (6), the splitting
integrator (3.1.2) then reads

(25) $\Phi_{h}=\exp(hY_{H_{3}})\circ\exp(hY_{H_{1}})\circ\exp(\sigma_{3}\Delta
W_{3}Y_{\widehat{H}_{3}})\circ\exp(\sigma_{1}\Delta
W_{1}Y_{\widehat{H}_{1}}),$

where for all $y\in\mathbb{R}^{3}$ one has

$\exp(\sigma_{1}\Delta W_{1}Y_{\widehat{H}_{1}})y=\begin{pmatrix}1&0&0\\\
0&\cos(y_{1}\sigma_{1}\Delta W_{1})&\sin(y_{1}\sigma_{1}\Delta W_{1})\\\
0&-\sin(y_{1}\sigma_{1}\Delta W_{1})&\cos(y_{1}\sigma_{1}\Delta
W_{1})\end{pmatrix}y$

and

$\exp(\sigma_{3}\Delta
W_{3}Y_{\widehat{H}_{3}})y=\begin{pmatrix}1&\sigma_{3}\Delta W_{3}&0\\\
0&1&0\\\ 0&0&1\end{pmatrix}y.$

Owing to Proposition 3, the explicit splitting scheme (25) is a stochastic
Poisson integrator.

####  4.1.2. Preservation of the Casimir of the stochastic Maxwell–Bloch
system

Let us first illustrate the qualitative behaviour of the stochastic Poisson
integrator (25) introduced above. Figure 1 illustrates the preservation of the
Casimir $C(y)=\frac{1}{2}(y_{2}^{2}+y_{3}^{2})$ by the stochastic Poisson
integrator (25). In this numerical experiment, the initial value is
$y(0)=(1,2,3)$ and the final time is $T=1$. We consider the two cases where
the system (6) is driven by a single Wiener process:
$(\sigma_{1},\sigma_{3})=(1,0)$ and $(\sigma_{1},\sigma_{3})=(0,1)$. Similar
results would be obtained if the system was driven by two independent Wiener
processes ($\sigma_{1}=\sigma_{3}=1$ for instance). In Figure 1, we compare
the numerical solutions given by the classical Euler–Maruyama scheme (applied
to the Itô formulation of the system), the stochastic midpoint scheme from
[36], and the explicit splitting scheme (25). The time step size is equal to
$h=0.01$. To implement the implicit stochastic midpoint scheme, a truncation
of the noise with threshold $A=\sqrt{4|\log(h)|}$ is applied (see [36] for
details). To be able to compare the results for different schemes, we use this
truncation in all experiments where the implicit stochastic midpoint scheme is
involved. As shown in Proposition 3, we observe that the Casimir function
$C(y)=\frac{1}{2}(y_{2}^{2}+y_{3}^{2})$ is preserved when using the stochastic
Poisson integrator (25). The Casimir function is also preserved when using the
stochastic midpoint scheme: indeed, this integrator is known to preserve
quadratic invariants, see [17].

Figure 1. Stochastic Maxwell–Bloch system: preservation of the Casimir by the
Euler–Maruyama scheme ($\times$), the midpoint scheme ($\circ$), and the
explicit stochastic Poisson integrator ($\diamond$). Left:
$(\sigma_{1},\sigma_{3})=(1,0)$. Right: $(\sigma_{1},\sigma_{3})=(0,1)$.

####  4.1.3. Strong and weak convergence of the explicit stochastic Poisson
integrator for the stochastic Maxwell–Bloch system

The preservation of the Casimir $C$ is not sufficient to ensure almost sure
boundedness of the numerical solution, which is instrumental to deduce Theorem
2 from Lemma 1. As a consequence, we are not able to state a convergence
result for the stochastic Poisson integrator (25) in general. However, when
$\sigma_{3}=0$, it is possible to show the following result.

######  Proposition 5.

Consider a numerical discretisation of the stochastic Maxwell–Bloch system (6)
by the stochastic Poisson integrator (25). Assume that $\sigma_{3}=0$. Then,
the strong order of convergence and the weak order of convergence of this
scheme are equal to $1$.

###### Proof.

Let us prove that, when $\sigma_{3}=0$, the following bounds are satisfied
almost surely:

(26) $\left\\{\begin{aligned} \|y^{[n]}\|&\leq(1+h)^{n}\|y^{[0]}\|\\\
\|y(t)\|&\leq e^{t}\|y(0)\|\end{aligned}\right.$

for all $n\geq 0$, $h\in(0,h_{0})$ and $t\geq 0$.

The proof of the bounds above is straightforward. On the one hand, for all
$y\in\mathbb{R}^{3}$, all $h>0$ and $t\in\mathbb{R}$, one has

$\|e^{tY_{H_{1}}}y\|=\|e^{tY_{\widehat{H}_{1}}}y\|=\|y\|$

and

$\displaystyle\|e^{hY_{H_{3}}}y\|^{2}$
$\displaystyle=(y_{1}+hy_{2})^{2}+y_{2}^{2}+y_{3}^{2}=\|y\|^{2}+2hy_{1}y_{2}+h^{2}y_{2}^{2}$
$\displaystyle\leq(1+h)^{2}\|y\|^{2}.$

Therefore

$\|y^{[n]}\|\leq\|e^{hY_{H_{3}}}\circ e^{hY_{H_{1}}}\circ
e^{\sigma_{1}\Delta_{n}W_{1}Y_{\widehat{H}_{1}}}(y^{[n-1]})\|\leq(1+h)\|y^{[n-1]}\|$

thus $\|y^{[n]}\|\leq(1+h)^{n}\|y^{[0]}\|$.

On the other hand, let $\mathcal{H}(y)=\|y\|^{2}$. Then one has
$\\{\mathcal{H},H_{1}\\}=\\{\mathcal{H},\widehat{H}_{1}\\}=0$ and
$\\{\mathcal{H},H_{3}\\}(y)=2y_{1}y_{2}\leq\mathcal{H}(y)$ for all
$y\in\mathbb{R}^{3}$ (recall that the Poisson bracket is defined by (4)).
Using (3), one thus obtains
$\,\mathrm{d}\mathcal{H}(y(t))\leq\mathcal{H}(y(t))$ for all $t\geq 0$ and the
bound for the exact solution follows from Gronwall’s lemma.

Let $T\in(0,\infty)$, one can then repeat the arguments used in the proof of
Theorem 2 as a corollary of Lemma 1, using the almost sure bounds

$\underset{t\in[0,T]}{\sup}~{}\left\lVert y(t)\right\rVert\leq
R(y_{0},T),\quad\underset{N\geq 1}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left\lVert y^{[n]}\right\rVert\leq R(y_{0},T)$

with $R(y_{0},T)=e^{T}\|y_{0}\|$. The details are omitted. Note that the
strong order of convergence is equal to $1$ since the system is driven by a
single Wiener process ($m=1$, the commutative noise case condition is
satisfied). This concludes the proof of Proposition 5. ∎

When $\sigma_{3}>0$, one can prove the following moment bound for the
numerical solution:

$\operatorname{E}[\|y^{[n]}\|^{2}]\leq
e^{(1+\frac{\sigma_{3}^{2}}{2})nh}\operatorname{E}[\|y^{[0]}\|^{2}].$

This follows from the inequality

$\displaystyle\operatorname{E}[\|e^{\sigma_{3}\Delta W_{3}Y_{H_{3}}}y\|^{2}]$
$\displaystyle=\operatorname{E}\bigl{[}(y_{1}+\sigma_{3}\Delta
W_{3}y_{2})^{2}+y_{2}^{2}+y_{3}^{2}\bigr{]}=\|y\|^{2}+\sigma_{3}^{2}h|y_{2}|^{2}$
$\displaystyle\leq(1+\sigma_{3}^{2}h)\|y\|^{2}$

and a recursion argument. A similar moment bound holds for the exact solution,
however these moment bounds are not sufficient to prove a strong convergence
result.

The objectives of this subsection are first to illustrate Proposition 5, and
second to investigate the behaviour of the strong and weak errors when the
condition $\sigma_{3}=0$ is removed. Whether it is possible to prove strong
and weak convergence estimates, or convergence in probability results, for
this problem in the general case is left open for future works.

We first illustrate the convergence of the strong error. In this numerical
experiment, the initial value is $y(0)=(1,2,3)$ and the final time is $T=1$.
We consider the two cases where the system (6) is driven by a single Wiener
process: $(\sigma_{1},\sigma_{3})=(1,0)$ and $(\sigma_{1},\sigma_{3})=(0,1)$.
Similar results would be obtained if the system was driven by two independent
Wiener processes ($\sigma_{1}=\sigma_{3}=1$ for instance). The reference
solution is computed using each scheme with time step size
$h_{\text{ref}}=2^{-16}$, and the schemes are applied with the range of time
step sizes $h=2^{-5},\ldots,2^{-13}$. The expectation is approximated
averaging the error over $M_{s}=500$ independent Monte Carlo samples.

Like in Figure 1, we compare the splitting integrator (25) with the standard
Euler–Maruyama scheme, and the stochastic midpoint scheme from [36]. A
truncation of the noise is used, see above for details. To be able to compare
the results for different schemes, we use this truncation in all experiments
where the implicit stochastic midpoint scheme is involved.

For the cases where the system is driven by a single Wiener process
($\sigma_{3}=0$ or $\sigma_{1}=0$), the results of the numerical experiment
are presented in Figure 2: we observe a strong order of convergence equal to
$1$ for the proposed explicit stochastic Poisson integrator (25). This
confirms the result of Proposition 5 when $\sigma_{3}=0$. We also conjecture
that the stochastic Poisson integrator (25) has strong order of convergence
equal to $1$ when $\sigma_{1}=0$.

Figure 2. Stochastic Maxwell–Bloch system with a single Wiener process: Strong
errors of the Euler–Maruyama scheme ($\times$), the midpoint scheme ($\circ$),
and the explicit stochastic Poisson integrator ($\diamond$). Left:
$(\sigma_{1},\sigma_{3})=(1,0)$. Right: $(\sigma_{1},\sigma_{3})=(0,1)$.

For the case where the system is driven by two Wiener processes, the results
of the numerical experiment are presented in Figure 3, with
$\sigma_{1}=\sigma_{3}=1$. Based on the observed convergence behaviour, we
conjecture that the strong order of the proposed integrator is $1/2$. This
result is not covered by the theoretical analysis performed in this article.

Figure 3. Stochastic Maxwell–Bloch system with
$(\sigma_{1},\sigma_{3})=(1,1)$: Strong errors of the Euler–Maruyama scheme
($\times$), the midpoint scheme ($\circ$), and the explicit stochastic Poisson
integrator ($\diamond$).

We now illustrate the weak convergence of the stochastic Poisson integrator
(25). For this numerical experiment, we set $\sigma_{1}=\sigma_{3}=1$, the
initial value is $y(0)=(1,2,3)$ and the final time is $T=1$. The reference
solution is computed using each scheme with time step size
$h_{\text{ref}}=2^{-16}$, and the schemes are applied with the range of time
step sizes $h=2^{-6},\ldots,2^{-12}$. The expectation is approximated
averaging the error over $M_{s}=10^{9}$ independent Monte Carlo samples.
Finally, the test function is given by $\phi(y)=\sin(2\pi y_{1})+\sin(2\pi
y_{2})+\sin(2\pi y_{3})$.

The results are presented in Figure 4. According to the observed rate of
convergence, we conjecture that the weak order of the proposed integrator is
$1$, but this result is not covered by the theoretical analysis performed in
this article.

Figure 4. Stochastic Maxwell–Bloch system: Weak errors of the explicit
stochastic Poisson integrator.

Numerical experiments illustrating the behaviour of the weak error when the
system is driven by a single noise are not reported: indeed, as seen in Figure
2, the stochastic Poisson integrator has strong order of convergence equal to
$1$ if $\sigma_{3}=0$ (rigorous result, Proposition 5) or if $\sigma_{1}=0$
(conjecture). In those cases, the weak error behaves like the strong error and
the rate of convergence is $1$.

To conclude this subsection, let us provide a numerical experiment using the
scheme (19) of weak order $2$ presented in Remark 3. For this numerical
experiment, all the values of the parameters are the same as for Figure 4,
except $\sigma_{1}=\sigma_{3}=10^{-3}$. The results are presented on Figure 5.
We observe that the weak convergence seems to be of order $2$ for the scheme
(19), but for small values of $h$ the error saturates due to the Monte Carlo
approximation. This is illustrated on the right figure, which gives results
for different values of the Monte Carlo sample size.

Figure 5. Stochastic Maxwell–Bloch system: Weak errors for the weak order $2$
scheme (19). Left: $10^{9}$ Monte Carlo samples. Right: $10^{6}$ to $10^{9}$
Monte Carlo samples.

####  4.1.4. Asymptotic preserving splitting scheme for the stochastic
Maxwell–Bloch system

In this subsection, we consider the multiscale version (15), parametrized by
$\epsilon$, of the stochastic Maxwell–Bloch system. Based on the expression
(25) for the stochastic Poisson integrator, applying the general asymptotic
preserving scheme (3.3) introduced in Section 3.3 gives the scheme

(27) $\left\\{\begin{aligned}
y^{\epsilon,[n]}&=\exp(hY_{H_{3}})\circ\exp(hY_{H_{1}})\circ\exp(\frac{\sigma_{3}h\xi_{3}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{3}})\circ\exp(\frac{\sigma_{1}h\xi_{1}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{1}})(y^{\epsilon,[n-1]}),\\\
\xi_{k}^{\epsilon,[n]}&=\xi_{k}^{\epsilon,[n-1]}-\frac{h}{\epsilon^{2}}\xi_{k}^{\epsilon,[n]}+\frac{\Delta_{n}W_{k}}{\epsilon}=\frac{1}{1+\frac{h}{\epsilon^{2}}}\Bigl{(}\xi_{k}^{\epsilon,[n-1]}+\frac{\Delta_{n}W_{k}}{\epsilon}\Bigr{)},\quad
k=1,2.\end{aligned}\right.$

The initial values are $y^{\epsilon,[0]}=y^{[0]}=y(0)$ and
$\xi_{k}^{\epsilon,[0]}=0$, $k=1,2$.

First, let us illustrate the qualitative behaviour of the scheme, for
different values of $\epsilon$. For this numerical experiment,
$\sigma_{1}=\sigma_{3}=0.1$, the initial value is $y(0)=(1,2,3)$ and the final
time is $T=1$. The time step size is equal to $h=10^{-3}$. In Figure 6, we
illustrate the preservation of the Casimir, up to an error of size
$O(10^{-14})$, for the asymptotic preserving scheme applied with
$\epsilon=1,0.1,0.001$ (left) and for the stochastic Poisson integrator
(4.2.1), formally $\epsilon=0$ (right). In Figure 7, we plot the evolution of
the approximation of the trajectory $t_{n}\mapsto
y(t_{n})=(y_{1}(t_{n}),y_{2}(t_{n}),y_{3}(t_{n}))$, for different values of
$\epsilon=1,0.1,0.001,0$. We observe that the trajectories are more regular
when $\epsilon$ is large and converge to the solution of the stochastic
Poisson integrator (4.2.1) as $\epsilon$ tends to $0$.

Figure 6. Stochastic Maxwell–Bloch with a single noise: evolution of the
Casimir using the asymptotic preserving scheme (27). Left:
$\epsilon=1,0.1,0.001$. Right: $\epsilon=0$.

Figure 7. Stochastic Maxwell–Bloch system: trajectories of the numerical
solution using the asymptotic preserving scheme (27). Top: left $\epsilon=1$,
right $\epsilon=0.1$. Bottom: left $\epsilon=0.001$, right $\epsilon=0$.

Finally, the last experiment of this subsection illustrates the uniform
accuracy property of the splitting scheme (27) with respect to the parameter
$\epsilon$ in the weak sense. For this numerical experiment,
$\sigma_{1}=\sigma_{3}=0.1$, the initial value is $y(0)=(1,2,3)$ and the final
time is $T=1$. The reference solution is computed using each scheme with time
step size $h_{\text{ref}}=2^{-16}$, and the schemes are applied with the range
of time step sizes $h=2^{-6},\ldots,2^{-12}$. The expectation is approximated
averaging the error over $M_{s}=10^{9}$ independent Monte Carlo samples.
Finally, the test function is given by $\phi(y)=\sin(2\pi y_{1})+\sin(2\pi
y_{2})+\sin(2\pi y_{3})$. The parameter $\epsilon$ takes the following values:
$\epsilon=10^{-2},10^{-3},10^{-4},10^{-5}$. The results are seen in Figure 8.
We observe that the weak error seems to be bounded uniformly with respect to
$\epsilon$, with an order of convergence $1$. For a standard method such as
the Euler–Maruyama scheme, the behaviour would be totally different: for fixed
time step size $h$, the error is expected to be bounded away from $0$ when
$\epsilon$ goes to $0$. Based on this numerical experiment, we conjecture that
the asymptotic preserving scheme is uniformly accurate.

Figure 8. Stochastic Maxwell–Bloch: Weak errors of the asymptotic preserving
scheme (27) for $\epsilon=10^{-2},10^{-3},10^{-4},10^{-5}$.

###  4.2. Explicit stochastic Poisson integrators for the stochastic rigid
body system

This subsection presents explicit stochastic Poisson integrators for the
stochastic rigid body system (3) (Example 3). We first give a detailed
construction of the splitting scheme, which gives a stochastic Poisson
integrator. We then illustrate its qualitative properties (preservation of the
Casimir function) and strong and weak error estimates of the proposed scheme
by numerical experiments. Finally, we illustrate the asymptotic preserving
property (see Section 3.3) for a multiscale version of the system.

Below we state and prove Proposition 6, which gives strong and weak rates of
convergence of the proposed explicit scheme. This is a non-trivial result
since the coefficients of the considered stochastic differential equation are
not globally Lipschitz continuous.

####  4.2.1. Presentation of the splitting scheme for the stochastic rigid
body system

To apply the strategy described in Section 3.1 and construct explicit
stochastic Poisson integrators, we first follow the approach from [34, 35] for
the deterministic rigid body system. The Hamiltonian function $H$ is split as
$H=H_{1}+H_{2}+H_{3}$, with $H_{j}=\frac{1}{2}\frac{y_{j}^{2}}{I_{j}}$ for
$j=1,2,3$. Recall also that the Hamiltonian functions appearing in the
stochastic part of the dynamics are given by
$\widehat{H}_{j}=\frac{1}{2}\frac{y_{j}^{2}}{\widehat{I}_{j}}$, for $j=1,2,3$.

The application of the general splitting integrator (3.1.2) for the stochastic
rigid body system (3) requires to compute the exact solutions of the
deterministic subsystems

$\dot{y}_{j}=B(y_{j})\nabla H_{j}(y_{j})$

and of the stochastic subsystems

$\,\mathrm{d}y_{j}=B(y_{j})\nabla\widehat{H}_{j}(y_{j})\circ\,\mathrm{d}W_{j}(t).$

As explained for instance in [34, 35] for the deterministic system, it is
straightforward to solve such subsystems. We only provide the details when
$j=1$. In that case, the deterministic subsystem is of the type

$\begin{cases}\dot{y}_{1}&=0\\\ \dot{y}_{2}&=y_{1}y_{3}/I_{1}\\\
\dot{y}_{3}&=-y_{1}y_{2}/I_{1}.\end{cases}$

The first equation yields that $y_{1}$ is constant, i. e. $y_{1}(t)=y_{1}(0)$
for all $t\geq 0$. As a consequence, $(y_{2},y_{3})$ is a solution of a linear
ordinary differential equation.

The deterministic subsystem when $j=1$ thus admits the exact solution

$\begin{pmatrix}y_{1}(t)\\\ y_{2}(t)\\\
y_{3}(t)\end{pmatrix}=\begin{pmatrix}1&0&0\\\ 0&\cos(\theta t)&\sin(\theta
t)\\\ 0&-\sin(\theta t)&\cos(\theta t)\end{pmatrix}y(0),$

where $\theta=\frac{y_{1}(0)}{I_{1}}$. Similarly, the stochastic subsystem
when $j=1$ is written as

$\begin{cases}\,\mathrm{d}y_{1}&=0\\\
\,\mathrm{d}y_{2}&=y_{1}y_{3}/\widehat{I}_{1}\circ\,\,\mathrm{d}W_{1}\\\
\,\mathrm{d}y_{3}&=-y_{1}y_{2}/\widehat{I}_{1}\circ\,\,\mathrm{d}W_{1}\end{cases}$

and it admits the exact solution

$\begin{pmatrix}y_{1}(t)\\\ y_{2}(t)\\\
y_{3}(t)\end{pmatrix}=\begin{pmatrix}1&0&0\\\ 0&\cos(\theta
W_{1}(t))&\sin(\theta W_{1}(t))\\\ 0&-\sin(\theta W_{1}(t))&\cos(\theta
W_{1}(t))\end{pmatrix}y(0).$

where $\theta=\frac{y_{1}(0)}{\widehat{I}_{1}}$.

The solutions of the deterministic and stochastic subsystems when $j=2,3$ have
similar expressions, which are not written here for brevity.

Finally, setting
$\Phi_{h}^{\text{det}}=\exp(hY_{H_{3}})\circ\exp(hY_{H_{2}})\circ\exp(hY_{H_{1}})$
and $\Phi_{\Delta W}^{\text{stoch}}=\exp(\Delta
W_{3}Y_{\widehat{H}_{3}})\circ\exp(\Delta
W_{2}Y_{\widehat{H}_{2}})\circ\exp(\Delta W_{1}Y_{\widehat{H}_{1}})$, the
general splitting integrator (3.1.2) applied to the stochastic rigid body
system (3) gives

$\displaystyle\Phi_{h}=\Phi_{h}^{\text{det}}\circ\Phi_{\Delta
W}^{\text{stoch}}$
$\displaystyle=\exp(hY_{H_{3}})\circ\exp(hY_{H_{2}})\circ\exp(hY_{H_{1}})$
(28) $\displaystyle\quad\circ\exp(\Delta
W_{3}Y_{\widehat{H}_{3}})\circ\exp(\Delta
W_{2}Y_{\widehat{H}_{2}})\circ\exp(\Delta W_{1}Y_{\widehat{H}_{1}}).$

Owing to Proposition 3, the explicit splitting scheme (4.2.1) is a stochastic
Poisson integrator. In particular, it preserves the Casimir function
$C(y)=y_{1}^{2}+y_{2}^{2}+y_{3}^{2}$, which has compact level sets.

####  4.2.2. Preservation of the Casimir of the stochastic rigid body system

Let us first illustrate the qualitative behaviour of the stochastic Poisson
integrator (4.2.1) introduced above. In this numerical experiment, the moments
of inertia are $I=(2,1,2/3)$, $\widehat{I}=(1,2,3)$, the initial value is
$y(0)=(\cos(1.1),0,\sin(1.1))$ and the final time is $T=20$. In Figure 9, we
compare the numerical solutions given by the classical Euler–Maruyama scheme
(applied to the Itô formulation of the system), the stochastic midpoint scheme
from [36], and the explicit splitting scheme (4.2.1). The time step size is
equal to $h=0.2$ ($T/h=100$). A truncation of the noise is used for this
experiment, see above for details. As proved in Proposition 3, we observe that
the Casimir function $C(y)=y_{1}^{2}+y_{2}^{2}+y_{3}^{2}$ is preserved when
using the stochastic Poisson integrator (4.2.1). The Casimir function is also
preserved when using the stochastic midpoint scheme: indeed, this integrator
is known to preserve quadratic invariants, see [17].

In addition, a plot of the evolution of the Hamiltonian for the three schemes
(middle figure), and of the trajectory on the sphere of the proposed splitting
scheme (right figure) are presented in Figure 9.

Figure 9. Stochastic rigid body system: Qualitative behaviour of the
Euler–Maruyama scheme ($\times$), the midpoint scheme ($\circ$), and the
explicit stochastic Poisson integrator ($\diamond$). Left: preservation of the
Casimir. Middle: evolution of the Hamiltonian. Right: trajectory on the sphere
for the scheme (4.2.1).

####  4.2.3. Strong and weak convergence of the explicit stochastic Poisson
integrator for the stochastic rigid body system

The general Theorem 2 is applicable in the case of stochastic rigid body
system since the Casimir function $C$ has compact level sets.

######  Proposition 6.

Consider a numerical discretisation of the stochastic rigid body system (3) by
the stochastic Poisson integrator (4.2.1). Then, the strong order of
convergence of this scheme is $1/2$ and the weak order of convergence is $1$.

###### Proof.

The stochastic Poisson system (3) admits the Casimir function $y\mapsto
C(y)=y_{1}^{2}+y_{2}^{2}+y_{3}^{2}$, which has compact level sets. It then
suffices to apply the general convergence result, Theorem 2, which concludes
the proof of Proposition 6. ∎

Let us first illustrate the strong convergence result. We compare the
behaviours of the three integrators introduced above: the Euler–Maruyama
scheme, the stochastic midpoint scheme, and the explicit stochastic Poisson
integrator (4.2.1). Note that Proposition 6 is valid only for the splitting
scheme (4.2.1). For this numerical experiment, the moments of inertia are
$I=(2,1,2/3)$, $\widehat{I}=(1,2,3)$, the initial value is
$y(0)=(\cos(1.1),0,\sin(1.1))$ and the final time is $T=1$. The reference
solution is computed using each scheme with time step size
$h_{\text{ref}}=2^{-16}$, and the schemes are applied with the range of time
step sizes $h=2^{-5},\ldots,2^{-13}$. The expectation is approximated
averaging the error over $M_{s}=500$ independent Monte Carlo samples.

The results are presented in Figure 10: we observe a strong order of
convergence equal to $1/2$ for the proposed explicit stochastic Poisson
integrator (4.2.1), which confirms the result of Proposition 6.

Figure 10. Stochastic rigid body system: Strong errors for the Euler–Maruyama
scheme ($\times$), the midpoint scheme ($\circ$), and the explicit stochastic
Poisson integrator ($\diamond$).

We now illustrate the weak convergence of the stochastic Poisson integrator
(4.2.1). For this numerical experiment, the moments of inertia are
$I=\hat{I}=(2,1,2/3)$, the initial value is $y(0)=(\cos(1.1),0,\sin(1.1))$ and
the final time is $T=1$. The reference solution is computed using each scheme
with time step size $h_{\text{ref}}=2^{-16}$, and the schemes are applied with
the range of time step sizes $h=2^{-6},\ldots,2^{-12}$. The expectation is
approximated averaging the error over $M_{s}=10^{9}$ independent Monte Carlo
samples. Finally, the test function is given by $\phi(y)=\sin(2\pi
y_{1})+\sin(2\pi y_{2})+\sin(2\pi y_{3})$. The results are presented in Figure
11. We observe a weak order $1$ for the proposed explicit stochastic Poisson
integrator (4.2.1), which confirms the result of Proposition 6.

Figure 11. Stochastic rigid body system: Weak error for the explicit
stochastic Poisson integrator (4.2.1).

To conclude this subsection, let us provide a numerical experiment using the
scheme of weak order $2$ presented in Remark 3. To do so, we consider the
following variant of the stochastic rigid body system (3):

$\displaystyle\,\mathrm{d}\begin{pmatrix}y_{1}\\\ y_{2}\\\ y_{3}\end{pmatrix}$
$\displaystyle=B(y)\left(\nabla
H(y)\,\,\mathrm{d}t+\sigma_{1}\nabla\widehat{H}_{1}(y)\circ\,\,\mathrm{d}W_{1}(t)+\sigma_{2}\nabla\widehat{H}_{2}(y)\circ\,\,\mathrm{d}W_{2}(t)\right.$
$\displaystyle\quad+\left.\sigma_{3}\nabla\widehat{H}_{3}(y)\circ\,\,\mathrm{d}W_{3}(t)\right),$

with nonnegative real numbers $\sigma_{1},\sigma_{2},\sigma_{3}$. For this
numerical experiment, all the values of the parameters are the same as for
Figure 11, except the values of the additional parameters
$\sigma_{1},\sigma_{2},\sigma_{3}$: one has either three Wiener processes with
$(\sigma_{1},\sigma_{2},\sigma_{3})=(10^{-3},10^{-3},10^{-3})$ (left figure),
or a single Wiener process, with three possible choices
$(\sigma_{1},\sigma_{2},\sigma_{3})=(10^{-3},0,0)$,
$(\sigma_{1},\sigma_{2},\sigma_{3})=(0,10^{-3},0)$ and
$(\sigma_{1},\sigma_{2},\sigma_{3})=(0,0,10^{-3})$ (right figure). The results
are presented in Figure 12. We observe that the weak convergence seems to be
of order $2$ for the scheme (19), however for small values of $h$ the error
saturates due to the Monte Carlo approximation.

Figure 12. Stochastic rigid body system: Weak error for the integrator (19).
Left: $(\sigma_{1},\sigma_{2},\sigma_{3})=(10^{-3},10^{-3},10^{-3})$. Right:
$(\sigma_{1},\sigma_{2},\sigma_{3})=(10^{-3},0,0)$,
$(\sigma_{1},\sigma_{2},\sigma_{3})=(0,10^{-3},0)$ and
$(\sigma_{1},\sigma_{2},\sigma_{3})=(0,0,10^{-3})$.

####  4.2.4. Asymptotic preserving splitting scheme for the stochastic rigid
body system

In this subsection, we consider the multiscale version (15), parametrized by
$\epsilon$, of the stochastic rigid body system. Based on the expression
(4.2.1) for the stochastic Poisson integrator, applying the general asymptotic
preserving scheme (3.3) introduced in Section 3.3 gives the scheme

(29) $\left\\{\begin{aligned}
y^{\epsilon,n}&=\exp(hY_{H_{3}})\circ\exp(hY_{H_{2}})\circ\exp(hY_{H_{1}})\\\
&\quad\circ\exp(\frac{h\xi_{3}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{3}})\circ\exp(\frac{h\xi_{2}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{2}})\circ\exp(\frac{h\xi_{1}^{\epsilon,[n]}}{\epsilon}Y_{\widehat{H}_{1}})\\\
\xi_{k}^{\epsilon,[n]}&=\xi_{k}^{\epsilon,[n-1]}-\frac{h}{\epsilon^{2}}\xi_{k}^{\epsilon,[n]}+\frac{\Delta_{n}W_{k}}{\epsilon}=\frac{1}{1+\frac{h}{\epsilon^{2}}}\Bigl{(}\xi_{k}^{\epsilon,[n-1]}+\frac{\Delta_{n}W_{k}}{\epsilon}\Bigr{)},\quad
k=1,2,3.\end{aligned}\right.$

The initial values are $y^{\epsilon,[0]}=y^{[0]}=y(0)$ and
$\xi_{k}^{\epsilon,[0]}=0$, $k=1,2,3$.

First, let us illustrate the qualitative behaviour of the scheme (29), for
different values of $\epsilon$. For this numerical experiment, the moments of
inertia are $I=\hat{I}=(2,1,2/3)$, the initial value is
$y(0)=(\cos(1.1),0,\sin(1.1))$ and the final time is $T=1$. The time step size
is equal to $h=10^{-4}$. In Figure 13, we plot the evolution of the
Hamiltonian (top) and the Casimir (bottom) for the asymptotic preserving
scheme (29) applied with $\epsilon=1,0.1,0.001$ (left) and for the stochastic
Poisson integrator (4.2.1), formally $\epsilon=0$ (right). We observe the
preservation of the Casimir function. In Figure 14, we plot the evolution of
the approximation of the trajectory $t_{n}\mapsto
y(t_{n})=(y_{1}(t_{n}),y_{2}(t_{n}),y_{3}(t_{n}))$, for different values of
$\epsilon=1,0.1,0.001,0$. We observe that the trajectories are more regular
when $\epsilon$ is large and converge to the solution of the stochastic
Poisson integrator (4.2.1) as $\epsilon$ tends to $0$.

Figure 13. Stochastic rigid body system: evolution of the Hamiltonian (top)
and the Casimir (bottom) of the numerical solution using the asymptotic
preserving scheme (29). Left: $\epsilon=1,0.1,0.001$. Right: $\epsilon=0$.

Figure 14. Stochastic rigid body system: Trajectory of the numerical solutions
using the asymptotic preserving scheme (29). Top: $\epsilon=1,0.1$. Bottom:
$\epsilon=0.001,0$.

Finally, the last experiment of this subsection illustrates the uniform
accuracy property of the splitting scheme (29) with respect to the parameter
$\epsilon$ in the weak sense. For this numerical experiment, the moments of
inertia are $I=(2,1,2/3)$ and $\hat{I}=(20,10,20/3)$, the initial value is
$y(0)=(\cos(1.1),0,\sin(1.1))$ and the final time is $T=1$. The reference
solution is computed using each scheme with time step size
$h_{\text{ref}}=2^{-16}$, and the schemes are applied with the range of time
step sizes $h=2^{-6},\ldots,2^{-12}$. The expectation is approximated
averaging the error over $M_{s}=10^{8}$ independent Monte Carlo samples. The
test function is given by $\phi(y)=\sin(2\pi y_{1})+\sin(2\pi y_{2})+\sin(2\pi
y_{3})$. The parameter $\epsilon$ takes the following values:
$\epsilon=10^{-2},10^{-3},10^{-4},10^{-5}$. The results are seen in Figure 15.
We observe that the weak error seems to be bounded uniformly with respect to
$\epsilon$, with an order of convergence $1$. For a standard method such as
the Euler–Maruyama scheme, the behaviour would be totally different: for fixed
time step size $h$, the error is expected to be bounded away from $0$ when
$\epsilon$ goes to $0$. Based on this numerical experiment, we conjecture that
the asymptotic preserving scheme is uniformly accurate.

Figure 15. Stochastic rigid body system: Weak errors of the asymptotic
preserving scheme (29) for $\epsilon=10^{-2},10^{-3},10^{-4},10^{-5}$.

###  4.3. Explicit stochastic Poisson integrators for the stochastic
sine–Euler system

The last example of a stochastic Lie–Poisson system studied in this work is
the stochastic version of the sine–Euler equations (4) introduced in Example
4. Like for the previous examples, an explicit stochastic Poisson integrator
is designed using a splitting strategy. We both illustrate the qualitative
behaviour of the proposed integrator (preservation of Casimir functions) and
strong error estimates. Note that numerical experiments which would illustrate
weak error estimates or the asymptotic preserving property are not reported
for this example: indeed, the results would be similar to those presented for
the two other examples above.

Recall that the stochastic sine–Euler system (4) is of the type

$\displaystyle\,\mathrm{d}\omega$ $\displaystyle=B(\omega)\left(\nabla
H(\omega)\,\,\mathrm{d}t+\sigma_{(1,0)}\nabla\widehat{H}_{(1,0)}(\omega)\circ\,\,\mathrm{d}W_{(1,0)}(t)+\sigma_{(1,1)}\nabla\widehat{H}_{(1,1)}(\omega)\circ\,\,\mathrm{d}W_{(1,1)}(t)\right.$
$\displaystyle\quad\left.+\sigma_{(0,1)}\nabla\widehat{H}_{(0,1)}(\omega)\circ\,\,\mathrm{d}W_{(0,1)}(t)+\sigma_{(-1,1)}\nabla\widehat{H}_{(-1,1)}(\omega)\circ\,\,\mathrm{d}W_{(-1,1)}(t)\right),$

This is a stochastic Lie–Poisson system. In order to design a splitting
integrator, note that the Hamiltonian function $H$ can be split as
$H=H_{(1,0)}+H_{(1,1)}+H_{(0,1)}+H_{(-1,1)}$, with $H_{\bf
k}(\omega)=\widehat{H}_{\bf k}(\omega)=\frac{\omega_{\bf k}\omega_{\bf
k}^{\star}}{|\bf k|^{2}}$.

Like for the other examples, the deterministic subsystems

$\dot{\omega}_{\bf k}=B(\omega_{\bf k})\nabla H_{\bf k}(\omega_{\bf k})$

and the stochastic subsystems

$\,\mathrm{d}\omega_{\bf k}=B(\omega_{\bf k})\nabla\widehat{H}_{\bf
k}(\omega_{\bf k})\circ\,\mathrm{d}W_{\bf k}(t)$

can be solved exactly: indeed, for each subsystem, the variable $\omega_{\bf
k}$ is preserved and the three other variables evolve following a linear
differential equation. We refer to [21, 34] for the explanation of this idea
for the deterministic subsystems. The treatment of the stochastic subsystems
is straightforward using the exact solution of the subsystems
$\dot{\omega}_{\bf k}=B(\omega_{\bf k})\nabla\widehat{H}_{\bf k}(\omega_{\bf
k})$. The splitting scheme for the stochastic sine–Euler SDE (4) then reads

$\displaystyle\Phi_{h}=\Phi_{h}^{\text{det}}\circ\Phi_{\Delta
W}^{\text{stoch}}$
$\displaystyle=\exp(h\Omega_{H_{(-1,1)}})\circ\exp(h\Omega_{H_{(0,1)}})\circ\exp(h\Omega_{H_{(1,1)}})\circ\exp(h\Omega_{H_{(1,0)}})$
$\displaystyle\quad\circ\exp(\sigma_{(-1,1)}\Delta
W_{(-1,1)}\Omega_{\widehat{H}_{(-1,1)}})\circ\exp(\sigma_{(0,1)}\Delta
W_{(0,1)}\Omega_{\widehat{H}_{(0,1)}})$ (30)
$\displaystyle\quad\circ\exp(\sigma_{(1,1)}\Delta
W_{(1,1)}\Omega_{\widehat{H}_{(1,1)}})\circ\exp(\sigma_{(1,0)}\Delta
W_{(1,0)}\Omega_{\widehat{H}_{(1,0)}}),$

where we denote by $\Omega_{H_{\bf k}}$, resp. by $\Omega_{\widehat{H}_{\bf
k}}$, the exact flow of the ODE subsystem, resp. SDE subsystem, with index
${\bf k}\in\\{(1,0),(1,1),(0,1),(-1,1)\\}$.

####  4.3.1. Preservation of Casimir functions for the stochastic sine–Euler
system

Recall that the stochastic sine–Euler system admits two Casimir functions
$C_{1}$ and $C_{2}$ given in Example 4. Owing to Proposition 3, the numerical
scheme (4.3) is a stochastic Poisson integrator, in particular it preserves
the two Casimir functions $C_{1}$ and $C_{2}$. We numerically illustrate this
property in Figure 16, where one sample of the numerical solution is computed
with the time step size $h=0.02$. The initial value is
$\omega(0)=(\omega_{(1,0)}(0),\omega_{(1,1)}(0),\omega_{(0,1)}(0),\omega_{(-1,1)}(0))=(0.1+0.3i,0.2+0.3i,0.3+0.2i,0.4+0.1i)$.
In addition, $\sigma_{k}=1$ for all $k$ in this experiment. Figure 16 confirms
that the two Casimir functions are preserved by the proposed integrator (4.3).

Figure 16. Stochastic sine–Euler system: Preservation of the two Casimir
functions by the explicit stochastic Poisson integrator (4.3).

####  4.3.2. Strong convergence of the explicit stochastic Poisson integrator
for the stochastic sine–Euler system

Like for the stochastic rigid body system, we have the following convergence
result due to the preservation of the Casimir function $C_{1}$.

######  Proposition 7.

Consider a numerical discretisation of the stochastic sine–Euler system (4) by
the explicit stochastic Poisson integrator (4.2.1). Then, the strong order of
convergence of this scheme is $1/2$, and the weak order of convergence is $1$.
If the system is driven by a single Wiener process, the strong order of
convergence is $1$.

As already explained, the convergence result above is not trivial since the
scheme is explicit: since the coefficients of the equations are not globally
Lipschitz continuous, the explicit Euler–Maruyama scheme does not converge in
the strong sense.

###### Proof.

The stochastic Poisson system (4) admits the Casimir function $\omega\mapsto
C_{1}(\omega)=|\omega_{1}|^{2}+\ldots+|\omega_{4}|^{2}$, which has compact
level sets. It then suffices to apply the general convergence result, Theorem
2, which concludes the proof of Proposition 7. ∎

We now numerically illustrate the strong rate of convergence of the proposed
integrator (4.3) when applied to the SDE (4). The final time is $T=1$ and the
initial value is

$\omega(0)=(\omega_{(1,0)}(0),\omega_{(1,1)}(0),\omega_{(0,1)}(0),\omega_{(-1,1)}(0))=(0.1+0.3i,0.2+0.3i,0.3+0.2i,0.4+0.1i)$.
The reference solution is computed using the proposed scheme with time step
size $h_{\text{ref}}=2^{-15}$, and the scheme is applied with the range of
time step sizes $h=2^{-5},\ldots,2^{-13}$. The expectation is approximated
averaging the error over $M_{s}=500$ independent Monte Carlo samples. The
results are presented in Figure 17. On the left, $\sigma_{1}=1$ and
$\sigma_{j}=0$, for $j=2,\ldots,4$: we observe an order of convergence equal
to $1$, which confirms the result in Proposition 7 when the system is driven
by a single Wiener process. On the right, $\sigma_{j}=1$, for $j=1,\ldots,4$,
which means that the system is driven by four independent Wiener processes. We
observe an order of convergence equal to $1/2$, which confirms the result in
Proposition 7.

Figure 17. Stochastic sine–Euler system: Strong errors of the explicit
stochastic Poisson integrator (4.3). Left: single noise. Right: multiple
noise.

##  Appendix A Proof of Lemma 1

The goal of this section is to present the proof of Lemma 1. The arguments are
standard and not specific to the present stochastic Poisson systems and
integrators, and they are provided for the convenience of the reader.

Recall that the auxiliary SDE (20) is given by

$\,\mathrm{d}z(t)=\sum_{k=1}^{p}f_{k}(z(t))\,\,\mathrm{d}t+\sum_{k=1}^{m}\widehat{f}_{k}(z(t))\circ\,\,\mathrm{d}W_{k}(t),$

and that the associated splitting integrator (21) is given by

(31)
$z^{[n]}=\varphi_{p}(h,\cdot)\circ\ldots\circ\varphi_{1}(h,\cdot)\circ\widehat{\varphi}_{m}(\Delta
W_{m}^{n},\cdot)\circ\ldots\circ\widehat{\varphi}_{1}(\Delta
W_{1}^{n},\cdot)(z^{[n-1]}).$

###### Proof of Lemma 1.

To establish the strong convergence estimate (22), we apply the fundamental
theorem on the strong order of convergence, see [38, Theorem 1.1]. Note that
the Itô formulation of the SDE (20) is an equation with globally Lipschitz
nonlinearities, due to the assumptions on $f_{1},\ldots,f_{p}$ and
$\widehat{f}_{1},\ldots,\widehat{f}_{m}$. Since the coefficients $f_{k}$ and
$\widehat{f}_{k}$ do not depend on time, it is sufficient to prove the
following local error estimates:

$\displaystyle\left(\mathbb{E}\left[\left\lVert
z(h)-z^{[1]}\right\rVert^{2}\right]\right)^{1/2}\leq c\bigl{(}1+\left\lVert
z_{0}\right\rVert^{2}\bigr{)}^{1/2}h\quad\text{and}$ (32)
$\displaystyle\quad\left\lVert\mathbb{E}\left[z(h)-z^{[1]}\right]\right\rVert\leq
c\bigl{(}1+\left\lVert z_{0}\right\rVert^{2}\bigr{)}^{1/2}h^{3/2},$

for some real number $c\in(0,\infty)$ which does not depend on the time step
size $h$. In the proof, the value of $c$ may change from line to line.

The proof is based on the comparison of Stratonovich–Taylor expansions of the
exact and numerical solutions. On the one hand, the exact solution of the
auxiliary SDE (20) satisfies the following Stratonovich–Taylor expansion
formula: for all $t\geq 0$,

$\displaystyle z(t)-z(0)$
$\displaystyle=\sum_{k=1}^{p}\int_{0}^{t}f_{k}(z(s))\,\mathrm{d}s+\sum_{k=1}^{m}\int_{0}^{t}\widehat{f}_{k}(z(s))\circ\,\mathrm{d}W_{k}(s)$
$\displaystyle=\sum_{k=1}^{p}f_{k}(z(0))t+\sum_{k=1}^{m}\widehat{f}_{k}(z(0))W_{k}(t)$
$\displaystyle\quad+\sum_{k,\ell=1}^{m}\widehat{f}_{k}^{\prime}(z(0))\widehat{f}_{\ell}(z(0))\int_{0}^{t}\int_{0}^{s}\,\mathrm{d}W_{\ell}(r)\circ\,\mathrm{d}W_{k}(s)$
$\displaystyle\quad+R_{e}(t,z_{0}),$

where the random variable $R_{e}(t,z_{0})$ satisfies

$\left(\mathbb{E}\left[\left\lVert
R_{e}(t,z_{0})\right\rVert^{2}\right]\right)^{1/2}\leq c(1+\left\lVert
z_{0}\right\rVert^{2})^{1/2}t^{3/2}.$

On the other hand, we claim that the numerical solution satisfies the
following local error expansion:

$\displaystyle z^{[1]}-z^{[0]}$
$\displaystyle=\sum_{k=1}^{p}f_{k}(z^{[0]})h+\sum_{k=1}^{m}\widehat{f}_{k}(z^{[0]})W_{k}(h)$
$\displaystyle\quad+\frac{1}{2}\sum_{k=1}^{m}\widehat{f}_{k}^{\prime}(z^{[0]})\widehat{f}_{k}(z^{[0]})W_{k}(h)^{2}$
$\displaystyle\quad+\sum_{1\leq\ell<k\leq
m}\widehat{f}_{k}^{\prime}(z^{[0]})\widehat{f}_{\ell}(z^{[0]})W_{k}(h)W_{\ell}(h)$
$\displaystyle\quad+R_{s}(h,z_{0}),$

where the random variable $R_{s}(h,z_{0})$ satisfies

$\left(\mathbb{E}\left[\left\lVert
R_{s}(h,z_{0})\right\rVert^{2}\right]\right)^{1/2}\leq c(1+\left\lVert
z_{0}\right\rVert^{2})^{1/2}h^{3/2}.$

Finally, recall that for all $h>0$ and all $k=1,\ldots,m$,

  * •

$\int_{0}^{h}\int_{0}^{s}\,\mathrm{d}W_{k}(r)\circ\,\mathrm{d}W_{k}(s)=\frac{1}{2}W_{k}(h)^{2}$,

  * •

$\mathbb{E}[\int_{0}^{h}\int_{0}^{s}\,\mathrm{d}W_{\ell}(r)\circ\,\mathrm{d}W_{k}(s)]=0$
if $\ell\neq k$,

  * •

$\mathbb{E}\left[\left\lVert\int_{0}^{h}\int_{0}^{s}\,\mathrm{d}W_{\ell}(r)\circ\,\mathrm{d}W_{k}(s)\right\rVert^{2}\right]={\rm
O}(h^{2})$.

Comparing the Stratonovich–Taylor expansions of the exact and numerical
solutions then provides the local error estimates (A). As already explained,
the local error estimates (A) imply the strong error estimate (22). Thus, it
remains to prove the claim above for the local error expansion of the
numerical solution.

Note that the arguments below illustrate why the independence of
$W_{1},\ldots,W_{m}$ is essential for the consistency of the scheme (see
Remark 4). To simplify the notation, we write ${\rm O}(h^{\alpha})$ for random
variables $r$ which satisfy $\left(\mathbb{E}\left[\left\lVert
r\right\rVert^{2}\right]\right)^{1/2}\leq c(1+|z^{[0]}|^{2})^{1/2}h^{\alpha}$
for some real number $c\in(0,\infty)$ which does not depend on $h$.

The local error $z^{[1]}-z^{[0]}$ is decomposed as

$z^{[1]}-z^{[0]}=(Z_{p}-Z_{p-1})+\ldots+(Z_{1}-\widehat{Z}_{m})+(\widehat{Z}_{m}-\widehat{Z}_{m-1})+\ldots+(\widehat{Z}_{1}-\widehat{Z_{0}}),$

where $\widehat{Z}_{0}=z^{[0]}$,
$\widehat{Z}_{k}=\widehat{\varphi}_{k}(W_{k}(h),\widehat{Z}_{k-1})$ for all
$k=1,\ldots,m$, $Z_{0}=\widehat{Z}_{m}$ and $Z_{k}=\varphi_{k}(h,Z_{k-1})$.
Observe that $z^{[1]}=Z_{p}$.

For all $k\in\\{1,\ldots,m\\}$, the random variables $\widehat{Z}_{k-1}$ and
$W_{k}(h)$ are independent (since by construction $\widehat{Z}_{k-1}$ only
depends on $W_{1}(h),\ldots,W_{k-1}(h)$). Since $\widehat{\varphi}_{k}$ is the
flow associated with the ODE $\dot{z}_{k}=\widehat{f}_{k}(z_{k})$, a
Stratonovich–Taylor expansion for the solution
$t\in[0,h]\mapsto\widehat{\varphi}_{k}(W_{k}(t),\widehat{Z}_{k-1})$ of the
Stratonovich SDE
$\,\mathrm{d}\widehat{z}_{k}=\widehat{f}_{k}(\widehat{z}_{k})\circ\,\mathrm{d}W_{k}(t)$
yields

$\displaystyle\widehat{Z}_{k}-\widehat{Z}_{k-1}$
$\displaystyle=\widehat{\varphi}_{k}(W_{k}(h),\widehat{Z}_{k-1})-\widehat{Z}_{k-1}$
$\displaystyle=W_{k}(h)\widehat{f}_{k}(\widehat{Z}_{k-1})+\frac{1}{2}W_{k}(h)^{2}\widehat{f}_{k}^{\prime}(\widehat{Z}_{k-1})\widehat{f}_{k}(\widehat{Z}_{k-1})+{\rm
O}(h^{3/2}).$

Using this result successively for $k=1,\ldots,m$, one gets the more precise
results

$\displaystyle\widehat{Z}_{k}-z^{[0]}$
$\displaystyle=\sum_{\ell=1}^{k}W_{\ell}(h)\widehat{f}_{\ell}(\widehat{Z}_{\ell-1})+{\rm
O}(h)$ $\displaystyle={\rm O}(h^{1/2})$
$\displaystyle=\sum_{\ell=1}^{k}W_{{\ell}}(h)\widehat{f}_{\ell}(z^{[0]})+{\rm
O}(h)$

using that $\widehat{Z}_{\ell-1}=z^{[0]}+O(h^{1/2})$ in the last step.
Finally, one has

$\displaystyle\widehat{Z}_{k}-\widehat{Z}_{k-1}$
$\displaystyle=W_{k}(h)\widehat{f}_{k}(z^{[0]})+\frac{1}{2}W_{k}(h)^{2}\widehat{f}_{k}^{\prime}(z^{[0]}){\widehat{f}}_{k}(z^{[0]})$
$\displaystyle\quad+\sum_{\ell=1}^{k-1}W_{k}(h)W_{\ell}(h)\widehat{f}_{k}^{\prime}(z^{[0]}){\widehat{f}}_{\ell}(z^{[0]})+{\rm
O}(h^{3/2}).$

Similarly, for all $k\in\\{1,\ldots,p\\}$, one obtains in a first step

$Z_{k}-Z_{k-1}=\varphi_{k}(h,Z_{k-1})-Z_{k-1}=hf_{k}(Z_{k-1})+{\rm O}(h^{2}),$

and in a second step (recall that $Z_{0}=\widehat{Z}_{m}=z^{[0]}+{\rm
O}(h^{1/2})$)

$Z_{k}-Z_{k-1}=hf_{k}(Z_{0})+{\rm O}(h^{2})=hf_{k}(z^{[0]})+{\rm O}(h^{3/2}).$

Summing all the local error terms then concludes the proof of the claim, i.e.
of the Taylor–Stratonovich expansion formula for the numerical solution.

Note that in the commutative noise case (in particular, when $m=1$), one has
the identity

$\sum_{1\leq\ell<k\leq
m}\widehat{f}_{k}^{\prime}(z^{[0]})\widehat{f}_{\ell}(z^{[0]})W_{k}(h)W_{\ell}(h)=\sum_{1\leq
k\neq\ell\leq
m}\widehat{f}_{k}^{\prime}(z(0))\widehat{f}_{\ell}(z(0))\int_{0}^{h}\int_{0}^{s}\,\mathrm{d}W_{\ell}(r)\circ\,\mathrm{d}W_{k}(s).$

Using the two Taylor–Stratonovich expansion formulas written above, one
obtains the following version of the local error estimates (A) in the
commutative noise case:

$\displaystyle\left(\mathbb{E}\left[\left\lVert
z(h)-z^{[1]}\right\rVert^{2}\right]\right)^{1/2}\leq c\bigl{(}1+\left\lVert
z_{0}\right\rVert^{2}\bigr{)}^{1/2}h^{3/2}\quad\text{and}$
$\displaystyle\quad\left\lVert\mathbb{E}\left[z(h)-z^{[1]}\right]\right\rVert\leq
c\bigl{(}1+\left\lVert z_{0}\right\rVert^{2}\bigr{)}^{1/2}h^{3/2}.$

Owing to the fundamental theorem on the strong order of convergence, see [38,
Theorem 1.1], the strong order of convergence is thus equal to $1$ in the
commutative noise case.

This concludes the proof of the strong error estimate (22), in the general and
in the commutative noise cases.

Let us now prove the weak error estimate (23). To simplify the notation, we
only prove the weak error estimate for $n=N$. Indeed, the extension to the
case $n=0,\ldots,N$ with a uniform upper bound is straightforward and details
are omitted.

The proof is a variant of the Talay–Tubaro argument, see [48]: the weak error
is written in terms of the solution of a backward Kolmogorov equation. First,
the infinitesimal generator $\mathcal{L}$ of the auxiliary SDE (20) is written
as

$\mathcal{L}=\sum_{k=1}^{p}\mathcal{L}_{k}+\sum_{k=1}^{m}\widehat{\mathcal{L}}_{k},$

where, for all $k\in\\{1,\ldots,p\\}$, $\mathcal{L}_{k}=f_{k}\cdot\nabla$ is
associated with the ODE $\dot{z}_{k}=f_{k}(z_{k})$, and for all
$k\in\\{1,\ldots,m\\}$,
$\widehat{\mathcal{L}}_{k}=\frac{1}{2}f_{k}^{\prime}f_{k}\cdot\nabla+\frac{1}{2}f_{k}f_{k}^{T}:\nabla^{2}$
is associated with the SDE
$\,\mathrm{d}z_{k}=\widehat{f}_{k}(z_{k})\circ\,\mathrm{d}W_{k}(t)$ rewritten
in Itô form as
$\,\mathrm{d}z_{k}=\frac{1}{2}f_{k}^{\prime}(z_{k})f_{k}(z_{k})+\widehat{f}_{k}(z_{k})\,\mathrm{d}W_{k}(t)$.
The formula for $\mathcal{L}$ above employs the independence of the Wiener
processes $W_{1},\ldots,W_{m}$.

Introduce the auxiliary function $u$ which is the solution of the backward
Kolmogorov equation

$\partial_{t}u(t,z)=\mathcal{L}u(t,z)$

for all $t\in[0,T]$ and $z\in\mathbb{R}^{d}$, with the initial condition
$u(0,\cdot)=\phi$, assumed to be of class $\mathcal{C}^{4}$ with bounded
derivatives. Since the mappings $f_{1},\ldots,f_{p}$ are of class
$\mathcal{C}^{4}$ and the mappings $\widehat{f}_{1},\ldots,\widehat{f}_{m}$
are bounded and of class $\mathcal{C}^{5}$, with bounded derivatives of any
order, the function $u$ is of class $\mathcal{C}_{b}^{2,4}$, with bounded
derivatives of any order. This means that $t\mapsto u(t,z)$ is of class
$\mathcal{C}^{2}$ for all $z\in\mathbb{R}^{d}$, and that $z\mapsto u(t,z)$ is
of class $\mathcal{C}^{4}$, with bounded derivatives. The proof is omitted,
see for instance [13] for such results (after transformation to the Itô
formulation).

The weak error at time $T=Nh$ is expressed in terms of the solution $u$ of the
Kolmogorov equation as follows (using a standard telescoping sum argument):

$\displaystyle\mathbb{E}[\phi(z(Nh))]-\mathbb{E}[\phi(z^{[N]})]$
$\displaystyle=\mathbb{E}[u(Nh,z_{0})]-\mathbb{E}[u(0,z^{[N]})]$
$\displaystyle=\sum_{n=0}^{N-1}\bigl{(}\mathbb{E}[u(T-t_{n},z^{[n]})]-\mathbb{E}[u(T-t_{n+1},z^{[n+1]})]\bigr{)}$
$\displaystyle=\sum_{n=0}^{N-1}\bigl{(}\mathbb{E}[u(T-t_{n},z^{[n]})]-\mathbb{E}[u(T-t_{n+1},z^{[n]})]\bigr{)}$
$\displaystyle\quad+\sum_{n=0}^{N-1}\bigl{(}\mathbb{E}[u(T-t_{n+1},z^{[n]})]-\mathbb{E}[u(T-t_{n+1},z^{[n+1]})]\bigr{)},$