# A flexible split-step scheme for MV-SDEs

Xingyuan Chena
<EMAIL_ADDRESS>Gonçalo dos Reisa,b,***G.d.R. acknowledges support
from the _Fundaç $\tilde{\text{a}}$o para a Ci$\hat{e}$ncia e a Tecnologia_
(Portuguese Foundation for Science and Technology) through the project
UIDB/00297/2020 (Centro de Matemática e Aplicaç$\tilde{\text{o}}$es
CMA/FCT/UNL).
<EMAIL_ADDRESS>
(a School of Mathematics, University of Edinburgh, The King’s Buildings,
Edinburgh, UK
b Centro de Matemática e Aplicaç$\tilde{\text{o}}$es (CMA), FCT, UNL, Portugal
() )

###### Abstract

We present an implicit Split-Step explicit Euler type Method (dubbed SSM) for
the simulation of McKean–Vlasov Stochastic Differential Equations (MV-SDEs)
with drifts of super-linear growth in space, Lipschitz in measure and non-
constant Lipschitz diffusion coefficient. The scheme is designed to leverage
the structure induced by the interacting particle approximation system,
including parallel implementation and the solvability of the implicit
equation.

The scheme attains the classical $1/2$ root mean square error (rMSE)
convergence rate in stepsize and closes the gap left by [18] regarding
efficient implicit methods and their convergence rate for this class of
McKean-Vlasov SDEs. A sufficient condition for mean-square contractivity of
the scheme is presented. Several numerical examples are presented including a
comparative analysis to other known algorithms for this class (taming and
adaptive time-stepping) across parallel and non-parallel implementations.

Keywords: McKean-Vlasov equations, split-step methods, interacting particle
systems, superlinear growth 2000 MSC: 65C05 , 65C30 , 65C35

###### Contents

  1. 1 Introduction
  2. 2 The split-step methods for MV-SDEs
    1. 2.1 Notation and Spaces
    2. 2.2 Framework
    3. 2.3 The split-step method (SSM) for MV-SDEs: convergence and stability
  3. 3 Examples of interest
    1. 3.1 Example: the stochastic Ginzburg Landau equation
    2. 3.2 Example: the mean field FitzHugh-Nagumo model from Neuroscience
    3. 3.3 Example: Polynomial drift (non-Lipschitz measure dependency)
    4. 3.4 Stability of the SSM: linear, non-linear and the Cucker-Smale model case
  4. 4 Proof of the convergence result for the split-step method (SSM)
    1. 4.1 Proof of the main convergence result
    2. 4.2 Proof of the stability theorem
  5. A Short description of the Taming and time-adaptive method
    1. A.1 Explicit Euler particle system for the super-linear case: taming
    2. A.2 Time adaptive schemes

##  1 Introduction

The aim of this paper is to develop a numerical scheme for simulating McKean-
Vlasov Stochastic Differential Equations (MV-SDEs) with drifts of super-linear
growth and Lipschitz in the measure component, and Lipschitz diffusion
coefficients. MV-SDEs differ from standard SDEs by means of the presence of
the law of the solution process in the coefficients and their dynamics are of
the following type

$\mathrm{d}X_{t}=\Big{(}v(t,X_{t})+b(t,X_{t},\mu_{t}^{X})\Big{)}\mathrm{d}t+\sigma(t,X_{t},\mu_{t}^{X})\mathrm{d}W_{t},\quad
X_{0}\in L_{0}^{m}(\mathbb{R}^{d}),$ (1.1)

where $\mu_{t}^{X}$ denotes the law of the solution process $X$ at time $t$,
$W$ is a multidimensional Brownian motion, $v,b,\sigma$ are measurable maps
along with a sufficiently integrable initial condition $X_{0}$ (in $L_{0}^{m}$
for $m\geq 2$). These equations were introduced by McKean in the sixties [42]
and have been the target of much research since. There is a rich literature on
well-posedness [47, 6, 44, 46] and we point to the recent summary work [29]
highlighting recent developments in regularity estimates, exponential
ergodicity, long time large deviations (also [19, 1]), comparison theorems and
the Vlasov-Fokker-Plank equations associated to MV-SDEs. One particular
element of interest is the so-called propagation of chaos (PoC) introduced by
Kac [32] and further studied in the MV-SDE literature [47, 43, 35, 36]. The
PoC phenomena states that a MV-SDE is the limit of a certain weakly
interacting particle systems (of standard SDEs) as the system’s size increases
to infinity. Namely, (1.1) is the limit, as $N\to\infty$, of the
$N$-dimensional system of $\mathbb{R}^{d}$-valued interacting particles
$X^{i,N}$ ,

$\displaystyle\mathrm{d}{X}_{t}^{i,N}=\Big{(}v\big{(}t,X_{t}^{i,N}\big{)}+b\big{(}t,{X}_{t}^{i,N},\mu^{X,N}_{t}\big{)}\Big{)}\mathrm{d}t+\sigma\Big{(}t,{X}_{t}^{i,N},\mu^{X,N}_{t}\Big{)}\mathrm{d}W_{t}^{i},\qquad
X^{i,N}_{0}=X_{0}^{i},$ (1.2)

with $\mu^{X,N}_{t}$ being the empirical measure given as
$\mu^{X,N}_{t}(\mathrm{d}x):=\frac{1}{N}\sum_{j=1}^{N}\delta_{X_{t}^{j,N}}(\mathrm{d}x)$
and $\delta_{{X}_{t}^{j,N}}$ is the Dirac measure at point ${X}_{t}^{j,N}$,
$W^{i}$ are independent Brownian motions and with independent and identically
distributed initial conditions $X_{0}^{i}$ across $i=1,\cdots,N$. In essence,
the law $\mu^{X}_{t}$ in $\mathbb{R}^{d}$ of (1.1) is approximated by the
empirical average $\mu^{X,N}_{t}$ generated by the
$(\mathbb{R}^{d})^{N}$-system (a high-dimensional system). This methodology,
appealing to the interacting particle system, includes a well-known quantified
speed of convergence result (summarised in Proposition 2.2 below) providing a
path for a numerical approximation: from $X$ to $\\{X^{i,N}\\}_{i}$ to
$\\{X^{i,N,\pi}\\}_{i}$ with $X^{\cdot,N,\pi}$ the numerical approximation of
the each particle is the SDE system.

It is important to highlight that although $\\{X^{i,N}\\}_{i}$ is a high-
dimensional SDE and any existing numerical method a priori applies
straightforwardly, all rates and results obtained in this straightforwardly
way depend on the system’s dimension $Nd$. The constants then explode as $N$
increases. Part of the difficulty of the method is to be able to show that
such constants, rates and results are independent of $N$ albeit depending on
$d$. This step has its own non-trivial difficulties which we discuss further
below by drawing on prior work [18].

In this work, we focus on the class of MV-SDE with drifts of super-linear
growth in their spatial components [19, 18, 45, 4, 33], encapsulated in the
function $v$, where $b,\sigma$ are uniformly Lipschitz (in space and measure),
and $\sigma$ is non-constant. This class of MV-SDEs appears in several
practical models in science, for example, in neuroscience [3, 12] introduce
the mean-field FitzHugh-Nagumo model for a neurons networks in the brain; [11]
discuss individual-based and swarming Cucker-Smale interaction models; and
models of battery electrodes [21, 26]. These equation do not have explicit or
closed form solutions and numerical approximations are needed. Moreover,
standard explicit numerical methods suitable for the Lipschitz case fail to
converge on the super-linear growth setting. This is exemplified by the
‘particle corruption’ phenomena [18, Section 4.1] for MV-SDE particle system
numerics. This phenomena is akin to the divergence of SDE schemes in the
superlinear growth setting as highlighted in the seminal work [30].

The numerical approximation of McKean–Vlasov equations in the continuous case
was initiated in [13] and has been investigated further in a number of more
recent works. We briefly mention a few on numerical schemes for MV-SDEs under
the superlinear setting and in the Brownian framework. Tamed Euler schemes
appeared first [18], shortly followed by tamed Milstein schemes [5, 33]
(appeared simultaneously) and Milstein schemes for delay MV-SDEs [4]. In [34]
it is proposed a tamed scheme (and a new wellposedness result) for MV-SDE
featuring common noise in the particle system (which (1.2) does not) and where
the diffusion is also allowed to grow super-linearly. Adaptive time-stepping
methods come as an alternative to taming and in the context of MV-SDE they are
proposed in [45].

Outside the superlinear setting, [31] discuss computational complexity of MV-
SDE algorithms (uniformly Lipschitz drift and constant diffusion coefficient)
and we point the reader there for an in-depth overview of the state of the art
in that regard. Numerical approximations for MV-SDEs with non-Lipschitz
conditions in measure and space exist [17] but impose a linear growth
condition on the coefficients. The case of simulating MV-SDEs with
discontinuous coefficients has been addressed [37]. An alternative to the
empirical measure approximation of (1.2) is to use projection-type estimation
of the marginal densities [7] where the error analysis requires
differentiability of the coefficients. Variance reduction technique have been
analysed for the class of MV-SDE, namely, importance sampling [20], antithetic
multilevel Monte Carlo sampling [4] and antithetic sampling [8]. There also
recent progress in the jump-diffusion setting [2, 10].

Motivation. For the superlinear growth case described above, we are motivated
by an open question left in [18] regarding implicit-type numerical methods for
MV-SDEs. All the methods described above are of explicit time-stepping nature
and these are known to lose some of the geometric properties of the original
system. For instance, taming destroys the strict dissipativity of the drift
map which then raises questions regarding the stability of the scheme’s output
across long time horizons (see [49]). Implicit methods for MV-SDEs are largely
unexplored. The notable exceptions are [39, 18] which are also starting point
for this work. In [39] the authors study MV-SDEs and associated particle
systems with drifts of (symmetric) convolution kernel-type (as in [1]) and
constant diffusion coefficient. The theoretical (Bacry-Emery) machinery
employed there yields a critical logarithmic Sobolev inequality estimate that
is used to establish a concentration inequality for their implicit Euler
scheme. Assumption-wise, their setting and the setting of this work do not
cover each other but agree over a small class. It is unclear how to extend
their methodology for non-constant possibly degenerate diffusion coefficients.

More recently, a classic implicit method to deal with the super-linear growth
was proposed [18] (for general diffusion coefficients and without a concavity
assumption) but only convergence was shown and under stronger restrictions
than expected: the measure component of the drift is Lipschitz in
Wasserstein-1 metric not just in Wasserstein-2 (plus uniformly bounded measure
dependency). At the core of these difficulties was the use of stopping times
arguments which revealed themselves difficult to handle with the measure
dependency. The critical point is the proof of Lemma 5.12 (p41) in [18] and
the calculations executed in p.48-49. Lastly, upon inspection of that proof,
we emphasise that general $\theta$-methods for SDEs [28] will face the same
difficulties.

Our contributions. In this work, we revisit the framework of [18] and propose
a split-step numerical scheme inspired by the earlier work [28]. Our
contributions can be summarised as follows,

  1. (I)

Main results. We proposed a split-step method (SSM), see (2.3)-(2.4) below,
for this class of MV-SDEs. We prove its convergence and recover the
$1/2$-convergence rate in rMSE under the same general assumptions as the tamed
[18] or the adaptive time-stepping method [45]. No differentability or non-
degeneracy assumptions are imposed and stopping-time arguments are fully
avoided.

We provide the stability analysis of Mean-square contractivity for the SSM
(non-constant diffusion coefficient). To the best of our knowledge this has
not yet been discussed for MV-SDE schemes in general. The stability of the SSM
provides a theoretical foundation for carrying out simulation with larger
timestep and we point to positive results by way of numerical simulation with
the Cucker-Smale flocking model [22] where the SSM outperforms both taming
[18] and adaptive time-stepping [45] algorithms.

In regards to known findings, the SSM here overcomes the limitations of the
implicit method proposed in [18, Section 3.2] and extends its scope of
application. In the context of standard SDE simulation (en the coefficients
are independent of the measure), our results lift the differentiability
restriction of [28, Assumption 3.1], allow for time-dependence and can take
advantage of (the possible) concavity of the map $v$ (in (1.1)) – this is a
mild improvement of known SDE results.

  2. (II)

Computational gains. The scheme is designed so that the superlinear term can
be split from the main equation in a way that optimises/minimizes the
computational cost of the inversion method (from the implicit component). This
flexibility is understood in the following way: given an MV-SDE it is left to
the user to choose which terms form $v$ and which terms form $b$ in (1.1)
(within restrictions). Moreover, since there is a split $v$ Vs $b$ the user
may decide to add & subtract convenient terms to the drift (see Section 3.4.2
below; also [14, Eqs. (37) and (38)]).

The SSM, within the structure of (1.1), means that the split-step part of the
SSM can be solved independently of the measure component. In other words, the
SSM is amenable to a parallel implementation (e.g., [41]). Namely, one
parallelizes the task of solving $N$-times an $\mathbb{R}^{d}$-system in
opposition to the non-parallelizable task of solving the
$\mathbb{R}^{dN}$-system. The computational gains of the parallel
implementations for the SMM are shown to be on par of those for the explicit
taming [18] and adaptive time-stepping [45] algorithms.

  3. (III)

Comparative study against known literature. We provide a comparative analysis
against the taming [18] and adaptive time-stepping [45] methods, across
parallel and non-parallel implementations. The numerical study covers four
examples of interest with each highlighting a different flavour of these 3
algorithms including stability experiments.

We close this introduction with a two further comments. The general setting of
[19, 18, 45, 4, 33] is based on drifts maps
$(t,x,\mu)\mapsto\widehat{b}(t,x,\mu)$ and it is this map that satisfies a
one-sided Lipschitz condition in space and a uniform Lipschitz condition in
measure. In (1.1), we specify $\widehat{b}$ to be
$(t,x,\mu)\mapsto\widehat{b}(t,x,\mu)=v(t,x)+b(t,x,\mu)$ where the polynomial
growth is fully captured by $v$ while $b$ remains uniformly Lipschitz in its
variables. We have not found any example of relevance in the general framework
of $\widehat{b}$ that allows for spatial-measure cross-dependency satisfying
the assumptions, in other words, separating $\widehat{b}$ into $v$ and $b$ is
natural and within the scope of known models ($v$ would take the role of a
confining potential in scope of particle interaction models [27, 1]).

Two settings outside the scope of this work that will be explored in future
research are the case of non-Lipschitz measure dependency from [39, 1] and the
super-linear diffusion case of [34] (also [45, Section 3.1]).

Organization of the paper. In section 2 we set notation and review a few
concepts necessary for the main constructions. In Section 2.3 we state the
SSM, the main theorem of convergence and the stability results for the scheme.
Section 3 provides several numerical examples for comparison with other
methods. Examples covering the stability analysis in the super-linear case,
notably the Cucker-Smale model, are also given in Section 3. All proofs are
postponed to Section 4. For convenience of the reader a short description of
the Taming and time-adaptive algorithms is provided in A.

##  2 The split-step methods for MV-SDEs

###  2.1 Notation and Spaces

Let $\mathbb{N}$ be the set of natural numbers starting at $1$, $\mathbb{R}$
denotes the real numbers. For $a,b\in\mathbb{N}$ with $a\leq b$, we define
$\llbracket a,b\rrbracket:=[a,b]\cap\mathbb{N}=\\{a,\cdots,b\\}$. For
$x,y\in\mathbb{R}^{d}$ denote the scalar product of vectors by $\langle
x,y\rangle$; and $|x|=(\sum_{j=1}^{d}x_{j}^{2})^{1/2}$ the Euclidean distance.
Let $\mathbbm{1}_{A}$ be the indicator function of set
$A\subset\mathbb{R}^{d}$. For a matrix $A\in\mathbb{R}^{d\times n}$ we denote
by $A^{\intercal}$ its transpose and its Frobenius norm by
$|A|=\textrm{Trace}\\{AA^{\intercal}\\}^{1/2}$. Let
$I_{d}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ be the identity map. For collections
of vectors, let the upper indices denote the distinct vectors, whereas the
lower index is a vector component, i.e. $x^{l}_{j}$ denote the $j$-th
component of $l$-th vector.

We introduce over $\mathbb{R}^{d}$ the space of probability measures
$\mathcal{P}(\mathbb{R}^{d})$ and its subset $\mathcal{P}_{2}(\mathbb{R}^{d})$
of those with finite second moment. The space
$\mathcal{P}_{2}(\mathbb{R}^{d})$ is Polish under the Wasserstein distance

$\displaystyle
W^{(2)}(\mu,\nu)=\inf_{\pi\in\Pi(\mu,\nu)}\Big{(}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}|x-y|^{2}\pi(dx,dy)\Big{)}^{\frac{1}{2}},\quad\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{d}),$

where $\Pi(\mu,\nu)$ is the set of couplings for $\mu$ and $\nu$ such that
$\pi\in\Pi(\mu,\nu)$ is a probability measure on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ such that
$\pi(\cdot\times\mathbb{R}^{d})=\mu$ and $\pi(\mathbb{R}^{d}\times\cdot)=\nu$.

Let our probability space be a completion of
$(\Omega,\mathbb{F},\mathcal{F},\mathbb{P})$ with
$\mathbb{F}=(\mathcal{F}_{t})_{t\geq 0}$ carrying a Brownian motion
$W=(W^{1},\cdots,W^{l})$ with $l$-dimensions and generating the probability
space’s filtration, augmented by all $\mathbb{P}$-null sets, and with an
additionally sufficiently rich sub $\sigma$-algebra $\mathcal{F}_{0}$
independent of $W$. We denote by
$\mathbb{E}[\cdot]=\mathbb{E}^{\mathbb{P}}[\cdot]$ the usual expectation
operator with respect to $\mathbb{P}$. Let our probability space be an
atomless Polish space. As usual in probability when discussing random
variables, “i.i.d.” stands for “independent and identically distributed”.

We fix a finite terminal time $T<\infty$ and use the following notation for
spaces, which are standard in the (McKean-Vlasov) literature [18]. Let $p\geq
1$. Let $L_{t}^{p}\left(\mathbb{R}^{d}\right)$ define the space of
$\mathbb{R}^{d}$-valued, $\mathcal{F}_{t}$-measurable random variables $X$,
that satisfy $\mathbb{E}\left[|X|^{p}\right]^{1/p}<\infty$. Define
$\mathbb{S}^{p}$ as the space of $\mathbb{R}^{d}$ -valued,
$\mathcal{F}_{\cdot}$-adapted processes $Z$, that satisfy
$\mathbb{E}\big{[}\sup_{0\leqslant t\leqslant
T}|Z(t)|^{p}\big{]}^{1/p}<\infty$. Given two processes $(X_{t})_{t\geq 0}$ and
$(Y_{t})_{t\geq 0}$ let $\langle X_{t},Y_{t}\rangle$ denote their cross-
variation up to time $t$.

Throughout the text $C$ denotes a generic constant positive real number that
may change from line to line, may depend on the problem’s data but is always
independent of the constants $h,M,N$ (associated with the numerical scheme and
specified below).

###  2.2 Framework

Let $W$ be an $l$-dimensional Brownian motion and take the measurable maps
$v:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$,
$b:[0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}^{d}$
and
$\sigma:[0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}^{d\times
l}$. The MV-SDE of interest for this work is Equation (1.1) (for some $m\geq
1$),

where $\mu_{t}^{X}$ denotes the law of the process $X$ at time $t$, i.e.,
$\mu_{t}^{X}=\mathbb{P}\circ X_{t}^{-1}$. We make the following assumption on
the coefficients throughout.

######  Assumption 2.1.

Assume that $v,b$ and $\sigma$ are $1/2$-Hölder continuous in time, uniformly
in $x\in\mathbb{R}^{d}$ and $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Assume
that $b,\sigma$ are uniformly Lipschitz in the sense that there exists
$L_{b},L_{\sigma}\geq 0$ such that for all $t\in[0,T]$ all
$x,x^{\prime}\in\mathbb{R}^{d}$ and all
$\mu,\mu^{\prime}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ we have

$\displaystyle|b(t,x,\mu)-b(t,x^{\prime},\mu^{\prime})|^{2}$
$\displaystyle\leq L_{b}(|x-x^{\prime}|^{2}+W^{(2)}(\mu,\mu^{\prime})^{2})$
$\displaystyle|\sigma(t,x,\mu)-\sigma(t,x^{\prime},\mu^{\prime})|^{2}$
$\displaystyle\leq
L_{\sigma}(|x-x^{\prime}|^{2}+W^{(2)}(\mu,\mu^{\prime})^{2}).$

For $v$, there exist $L_{v}\in\mathbb{R}$, $L_{\hat{v}}>0$, $q\in\mathbb{N}$
and $q>1$ such that for all $t\in[0,T]$, $x,x^{\prime}\in\mathbb{R}^{d}$, we
have

$\displaystyle\langle x-x^{\prime},v(t,x)-v(t,x^{\prime})\rangle$
$\displaystyle\leq L_{v}|x-x^{\prime}|^{2}$ $\displaystyle\textrm{(One-sided
Lipschitz)},$ $\displaystyle|v(t,x)-v(t,x^{\prime})|$ $\displaystyle\leq
L_{\hat{v}}(1+|x|^{q}+|x^{\prime}|^{q})|x-x^{\prime}|$
$\displaystyle\textrm{(Locally Lipschitz)}.$

Regarding the structural choice the drift being $v+b$ with only $v$ containing
the superlinear growth component, as opposed to a single drift map
$\widehat{b}$ in the style of [18, 45], we point the reader to the comment
given at the end of this work’s introduction. We consider the loss of
generality to be negligible.

Immediate well-known properties can be derived from this assumption.

######  Remark 2.2 (Implied properties).

Under Assumption 2.1, define $\widehat{L}_{v}={1}/{2}+L_{v}$,
$C_{T}=\sup_{t\in[0,T]}|v(t,0)|^{2}/2$, and let $C>0$. Then for all
$t\in[0,T]$, $x\in\mathbb{R}^{d}$ and $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$
one has

$\displaystyle\langle x,v(t,x)\rangle=\langle x-0,v(t,x)-v(t,0)\rangle+\langle
x,v(t,0)\rangle\leq L_{v}|x|^{2}+|x||v(t,0)|\leq
C_{T}+\widehat{L}_{v}|x|^{2}.$

Moreover, for $\psi\in\\{b,\sigma\\}$ one has

$\displaystyle\langle x,\psi(t,x,\mu)\rangle\leq
C(1+|x|^{2}+W^{(2)}(\mu,\delta_{0})^{2})\quad\textrm{and}\quad|\psi(t,x,\mu)|^{2}\leq
C(1+|x|^{2}+W^{(2)}(\mu,\delta_{0})^{2}).$

These assumptions allow for a far larger class of models than the standard
globally Lipschitz drift assumption. Under the one-sided Lipschitz drift
assumption, [19, Theorem 3.3] provides a result for wellposedness. We use a
particularised version of this result as we also require $1/2$-Hölder time-
continuity to control the growth in time of the coefficient. Any globally
space-measure Lipschitz drift such as the Kuramoto model [20, 18], can of
course be simulated using our techniques. In Section 3 we numerically
illustrate our findings and present several examples with one-sided Lipschitz
drifts that satisfy these assumptions including the adjusted Ginzburg Landau
equation. There are also many other classes of models where one-sided
Lipschitz drifts appear such as Kinetic models e.g. in [25] and Self-
stabilizing diffusions [11], [39] and also [40]; the Czirók model for
collective motion [24]; meta-stability of interacting particle systems [27].

######  Theorem 2.3 (Theorem 3.3 in [19]).

Let Assumption 2.1 hold. Further, assume for some $m\geq 2$, $X_{0}\in
L_{0}^{m}(\mathbb{R}^{d})$. Then, there exists a unique solution
$X\in\mathbb{S}^{m}([0,T])$ to the MV-SDE (1.1). For some positive constant
$C$ we have

$\displaystyle\mathbb{E}\big{[}\sup_{t\in[0,T]}|X_{t}|^{m}\big{]}\leq
C\left(1+\mathbb{E}[|X_{0}|^{m}]\right)e^{CT}.$

If the law $\mu^{X}$ is known beforehand, then the MV-SDE reduces to a
classical SDE with added time-dependency. Typically this is not the case and
usually the MV-SDE is approximated by a particle system.

The interacting particle system approximation. We approximate $X$ of MV-SDE
(1.1) (driven by the Brownian motion $W$), using the $N$-dimensional system
(1.2) of $\mathbb{R}^{d}$-valued interacting particles $\\{X^{i,N}\\}_{i}$
each particle starting from its i.i.d. random initial condition
${X}_{0}^{i,N}=X_{0}^{i}$ for ${i\in\llbracket 1,N\rrbracket}$ and driven by
its own independent Brownian motion $W^{i}$ (also independent of the BM $W$
appearing in (1.1); with a slight abuse of notation to avoid re-defining the
probability space’s filtration). In (1.2) we approximate the law $\mu_{t}^{X}$
of $X_{t}$ by the empirical measure generated by the particle system, namely,
$\mu^{X,N}_{t}(\mathrm{d}x):=\frac{1}{N}\sum_{j=1}^{N}\delta_{X_{t}^{j,N}}(\mathrm{d}x)$
and $\delta_{x}$ is the Dirac measure at point $x\in\mathbb{R}^{d}$.

Propagation of chaos. In order to show that the particle approximation (1.2)
is of effective use, one requires a PoC result. Although different types
exist, in this work, the convergence rate of the scheme is shown in root Mean
Square Error (rMSE) and thus a pathwise PoC result is needed. For this end, we
introduce the system of non interacting particles (compare with (1.2))

$\displaystyle\mathrm{d}X_{t}^{i}=\Big{(}v(t,X_{t}^{i})+b(t,X_{t}^{i},\mu^{X^{i}}_{t})\Big{)}\mathrm{d}t+\sigma(t,X_{t}^{i},\mu^{X^{i}}_{t})\mathrm{d}W_{t}^{i},\quad
X_{0}^{i}=X_{0}^{i}\,,\quad t\in[0,T]\,,$ (2.1)

which are of course just $N$ MV-SDEs (each in $\mathbb{R}^{d}$) and since the
$X^{i}$s are independent, $\mu^{X^{i}}_{t}=\mu^{X}_{t}$ for all $i$ (and
$\mu^{X}_{t}$ the law of the solution to (1.1)). Under global Lipschitz
conditions, one can then prove the following convergence result (see e.g.,
[47, 43, 35, 15])

$\displaystyle\lim_{N\rightarrow\infty}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\big{[}\sup_{0\leq t\leq
T}|X_{t}^{i,N}-X_{t}^{i}|^{2}\,\big{]}=0\,.$

All SDEs appearing below have i.i.d. initial condition $X_{0}^{i}$,
${i\in\llbracket 1,N\rrbracket}$, and we work on the interval $[0,T]$. Recall
that we want to associate a particle system to the MV-SDE and show its
_propagation of chaos_ , [18, Theorem 3.1] (and [18, Proposition 5.1] for the
well-posedness of the interacting particle system (1.2)).

######  Proposition 2.4 (Propagation of chaos).

Let the assumptions of Theorem 2.3 hold for some $m>4$. Let (for
${i\in\llbracket 1,N\rrbracket}$) $X^{i}\in\mathbb{S}^{m}$ be the solution to
(2.1).

Then, there exists a unique solution $\\{X^{i,N}\\}_{i}$ to (1.2) and for any
$1\leq p\leq m$ there exists $C>0$ independent of $N$ such that

$\displaystyle\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}[\sup_{t\in[0,T]}|X^{i,N}_{t}|^{p}]\leq
C_{p}(1+\mathbb{E}[\,|X_{0}^{\cdot}|^{2}]).$ (2.2)

Moreover, we have the following convergence result

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}[\sup_{0\leq t\leq
T}|X_{t}^{i}-X_{t}^{i,N}|^{2}]\leq C\begin{cases}N^{-1/2}&,\text{if }d<4,\\\
N^{-1/2}\log(N)\quad&,\text{if }d=4,\\\ N^{-2/d}&,\text{if }d>4.\end{cases}$

This result shows that the particle scheme will converge to the MV-SDE with a
given rate. Therefore, to show convergence between the numerical scheme
proposed below and the MV-SDE, we only need to show that the numerical version
of the particle scheme converges to the “true” particle scheme.

######  Remark 2.5 (Optimising the PoC convergence rate).

We highlight a result from [16] and later reviewed in [45] in the context of
numerical methods for MV-SDEs. The PoC rate can be improved for the case
$d=4$, namely the $\log(N)$ term can be omitted (under restrictions). This
holds under the additional constraint of a constant diffusion coefficient
$\sigma$, and a bounded drift with bounded derivatives. It does not cover the
super-linear growth case here, nonetheless, it is reasonable to expect that
the result can be lifted to match the drift condition in this work and a
diffusion coefficient that is bounded (and sufficiently smooth).

###  2.3 The split-step method (SSM) for MV-SDEs: convergence and stability

The numerical scheme proposed in this work, and dubbed Split-Step Method
(SSM), improves strongly on the implicit numerical scheme proposed in [18]. It
is an enhanced variant of the split-step backward Euler scheme for standard
SDEs [28, Eq. (3.8)-(3.9)] and here further optimized for the MV-SDE setting.

Define the uniform partition of $[0,T]$ as $\pi:=\\{t_{n}:=nh:n\in\llbracket
0,M\rrbracket,h:=T/M\\}$ for a prescribed $M\in\mathbb{N}$. Define recursively
the split-step method to approximate (1.2) as follows: for $i\in\llbracket
1,N\rrbracket$ set $\hat{X}_{0}^{i,N}=X^{i}_{0}$, then, iteratively over
$n\in\llbracket 0,M-1\rrbracket$ for all $i\in\llbracket 1,N\rrbracket$,

$\displaystyle Y_{n}^{i,\star,N}$
$\displaystyle=\hat{X}_{n}^{i,N}+hv\left(t_{n},Y_{n}^{i,\star,N}\right),\quad\quad\hat{\mu}^{Y,N}_{n}(dx):=\frac{1}{N}\sum_{j=1}^{N}\delta_{Y_{n}^{j,\star,N}}(dx)$
(2.3) $\displaystyle\hat{X}_{n+1}^{i,N}$
$\displaystyle=Y_{n}^{i,\star,N}+b(t_{n},Y_{n}^{i,\star,N},\hat{\mu}^{Y,N}_{n})h+\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}^{Y,N}_{n})\Delta
W_{n}^{i}.$ (2.4)

We state immediately the main convergence result between the continuous time
extension of the scheme (2.3)-(2.4) and the solution to the MV-SDE (1.1). The
proof is postponed to Section 4.1.

######  Theorem 2.6.

Let the assumptions of Proposition 2.4 hold. Let $p\geq 1$ for sufficiently
large $m$ (depending on $p$). Take the collection
$\\{{\hat{X}}_{n}^{i,N}\\}_{i,n}$ for $n\in\llbracket 0,M\rrbracket$,
$i\in\llbracket 1,N\rrbracket$ generated through the scheme (2.3)-(2.4) under
the timestep constraint expressed through the one-sided Lipschitz constant
$L_{v}$ of $v$ (Assumption 2.1)

$\displaystyle\begin{cases}h>0\ \textrm{ and }\
h\leq\frac{1}{1+2L_{v}}&,\text{if }L_{v}>-\frac{1}{2},\\\ h>0&,\text{if
}L_{v}\leq-\frac{1}{2}.\end{cases}$

Then, the following assertions hold. There exists a continuous-time extension,
$(\hat{X}_{t}^{i,N})_{t\in[0,T]}$ to the SSM (2.3)-(2.4) (and given in (4.8)
below), satisfying:

  1. 1.

Uniformly bounded $p$ moments. Given $p\geq 1$ there exist constant $C>0$ such
that

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}$
$\displaystyle\mathbb{E}\Big{[}\sup_{0\leq t\leq
T}|\hat{X}_{t}^{i,N}|^{2p}\Big{]}<C\Big{(}1+\mathbb{E}\big{[}\,|\hat{X}_{0}|^{2p}\big{]}\Big{)}<\infty.$

  2. 2.

Take $(X^{i,N})_{i}$ as the solution of the interacting particle system (1.2).
Then, scheme (2.3)-(2.4) converges to $X^{i,N}$ with a strong global
convergence rate of ${1}/{2}$ in root mean square error (rMSE) over $[0,T]$,
namely,

$\displaystyle\textrm{rMSE}=\sqrt{\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Bigg{[}\sup_{0\leq t\leq
T}|X_{t}^{i,N}-\hat{X}_{t}^{i,N}|^{2}\Bigg{]}}=\mathcal{O}(h^{\frac{1}{2}}).$
(2.5)

  3. 3.

Let $(X^{i})_{i}$ be the solution of the non-interacting particle system
(2.1).

Then, we have the following MSE convergence

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq T}|X_{t}^{i}-\hat{X}_{t}^{i,N}|^{2}\Big{]}\leq
C\begin{cases}h+N^{-1/2}&\text{if }d<4,\\\ h+N^{-1/2}\log(N)\quad&\text{if
}d=4,\\\ h+N^{-2/d}&\text{if }d>4.\end{cases}$

We remind the reader about Remark 2.5 concerning the PoC’s convergence rate.
The proof of the result is postponed to Section 4.1. The constraint of
$L_{v}<-1/2$ is not sharp. In fact, it can be replaced by $L_{v}<-\varepsilon$
for some $\varepsilon\in(0,1)$ at the expense of another constant growing
proportionally to $1/\varepsilon$. We choose for simplicity $\varepsilon=1/2$,
see Remark 4.3 and the definition of $\widehat{L}_{v}$ in Remark 2.2 for
further details.

Theorem 2.6 shows the SSM to have a strong convergence rate of $1/2$ (rMSE)
which is the same as the taming Euler-Maruyama scheme [18] and the time
adaptive Euler-Maruyama scheme [45]. Also, the complexity of the particle
system is $\mathcal{O}(N^{2})$ in the worst situation, however, in several
examples of section 3 the complexity for the calculation of the interaction
term is of $\mathcal{O}(N)$.

After the convergence study of Theorem 2.6 we introduce the notion of mean-
square contractivity and study the stability of the SSM (2.3)-(2.4).

######  Definition 2.7 (Mean-square contractivity).

Suppose that $X_{0}\in L_{0}^{m}(\mathbb{R}^{d})$ and $Z_{0}\in
L_{0}^{m}(\mathbb{R}^{d})$ for some sufficiently large $m\in\mathbb{N}$. Take
two numerical solutions of a same numerical scheme $\hat{X}_{n}^{i,N}$ and
$\hat{Z}_{n}^{i,N}$ of (1.1) with $\hat{X}_{0}^{i}$ and $\hat{Z}_{0}^{i}$
being i.i.d. copies of $X_{0}$ and $Z_{0}$ respectively. The scheme is called
mean-square contractive if we have

$\displaystyle\lim_{n\rightarrow\infty}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\,|\hat{X}_{n}^{i,N}-\hat{Z}_{n}^{i,N}|^{2}\Big{]}=0.$

The next result shows when the SSM is mean-square contractive.

######  Theorem 2.8.

Let Assumption 2.1 hold and we assume a further particular form of the
Lipschitz condition of $b,\sigma$. Namely let
$L_{b},L_{\hat{b}},L_{\sigma},L_{\hat{\sigma}}\geq 0$ be such that

$\displaystyle|b(t,x,\mu)-b(t,x^{\prime},\mu^{\prime})|^{2}$
$\displaystyle\leq
L_{b}|x-x^{\prime}|^{2}+L_{\hat{b}}W^{(2)}(\mu,\mu^{\prime})^{2},$
$\displaystyle|\sigma(t,x,\mu)-\sigma(t,x^{\prime},\mu^{\prime})|^{2}$
$\displaystyle\leq
L_{\sigma}|x-x^{\prime}|^{2}+L_{\hat{\sigma}}W^{(2)}(\mu,\mu^{\prime})^{2},$

for all $t\in[0,T]$, $x,x^{\prime}\in\mathbb{R}^{d}$ and
$\mu,\mu^{\prime}\in\mathcal{P}_{2}(\mathbb{R}^{d})$. Suppose that $X_{0}\in
L_{0}^{m}(\mathbb{R}^{d})$ and $Z_{0}\in L_{0}^{m}(\mathbb{R}^{d})$ for a
sufficiently large $m\in\mathbb{N}$, and let $\hat{X}_{0}^{i}$ and
$\hat{Z}_{0}^{i}$ be i.i.d. copies of $X_{0}$ and $Z_{0}$ respectively.

Set $h>0$. Define two families $\\{(X_{n}^{i,N},Y^{i,\star,N}_{n})\\}_{i,n}$
and $\\{(Z_{n}^{i,N},G^{i,\star,N}_{n})\\}_{i,n}$ as the output of the SSM
(2.3)-(2.4) with input initial conditions $\\{X_{0}^{i,N}\\}_{i}$ and
$\\{Z_{0}^{i,N}\\}_{i}$ respectively. Then, for any $n\in\mathbb{N}$,

$\displaystyle\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\,|\hat{X}_{n}^{i,N}-\hat{Z}_{n}^{i,N}|^{2}\Big{]}\leq(1+\beta
h)^{n}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}^{i,N}_{0}-\hat{Z}^{i,N}_{0}|^{2}\Big{]},$

where

$\displaystyle\beta=\frac{2(L_{v}-\alpha)+h(L_{b}+L_{\hat{b}})}{1-2hL_{v}}\quad\textrm{and}\quad\alpha=-\frac{1}{2}(1+L_{\sigma}+L_{\hat{\sigma}}+L_{b}+L_{\hat{b}}).$
(2.6)

Under the choice of $h$ of Theorem 2.6, the quantity $1+\beta h$ is always
positive. If $L_{v}<\alpha\leq-\frac{1}{2}$ and
$0<h<-2\frac{L_{v}-\alpha}{L_{b}+L_{\hat{b}}}$ then $\beta<0$ and thus the SSM
is Mean-square contractive.

The proof of this result is postponed to Section 4.2. We illustrate
immediately with several numerical results the scope of our findings. The
reasoning behind the specification of the Lipschitz constants of $b$ and
$\sigma$ will become apparent in the examples Section 3.4. The main motivation
is to account for the contribution of the spatial and measure components
separately as a way to explore the flexibility allowed by the scheme in
choosing the drift coefficients $v$ and $b$.

The equivalent definition to Mean-square contractivity (Definition 2.7) for
the initial MV-SDE (1.1) is the so-called Exponential mean-square stability
inequality defined next. We will make use of this definition in Section 3.4
below.

######  Definition 2.9 (Exponential mean-square contractive solutions).

Let $X,Y$ be two solutions to (1.1) with initial conditions $X_{0},Y_{0}\in
L_{0}^{m}(\mathbb{R}^{d})$ respectively. If $X,Y$ satisfy

$\displaystyle\mathbb{E}\big{[}\,|X_{t}-Y_{t}|^{2}\big{]}\leq e^{\eta
t}\mathbb{E}\big{[}|X_{0}-Y_{0}|^{2}\big{]}$

for some real number $\eta<0$, then the MV-SDE (1.1) is said to generate
exponential mean-square contractive solutions.

##  3 Examples of interest

For the convenience of the reader, we provide in A a brief description of the
taming scheme [18] and time-adaptive method [45] we use next. In the following
examples, we apply our settings of $v,b,\sigma$ in (1.1) and $i\in\llbracket
1,N\rrbracket$. As the “true” solution is unknown, and to compare the
convergence rates for these examples, the result of smaller timestep $h$
simulation will be used as a proxy to the true solution. Within each example,
all the methods will use same Brownian motion paths. We consider the weak
error ($k=1$) and the strong error ($k=2$) between the true solution $X_{T}$
and the approximation $\hat{X}_{T}$ as follow

$\displaystyle\epsilon_{k}=\Bigg{(}\mathbb{E}\Big{[}\,(X_{T}-\hat{X}_{T})^{k}\Big{]}\Bigg{)}^{1/k}\approx\Big{(}\frac{1}{N}\sum_{j=1}^{N}(X_{T}^{j}-\hat{X}_{T}^{j})^{k}\Big{)}^{\frac{1}{k}},\qquad
k=1,2.$

Our main theorem covers only the strong convergence result, nonetheless, we
also present the weak convergence rate estimation. We leave the weak error
analysis for future research.

We study several examples. The stochastic Ginzburg Landau example is a well-
studied one and it provides a comparison example to other methods. The second
example is a multi-dimensional FitzHugh-Nagumo model of McKean-Vlasov type
from neuroscience which shows that the split-step method can properly deal
with the super-linear term in a complex system. For these two examples we
discuss as well the simulation runtime in serial and parallel implementation.
The third example lies outside the scope of our assumptions by featuring a
non-Lipschitz measure dependency, however, the split-step method still works.
Proving the convergence of the method under this setting is left for future
research.

In the last part, we discuss the stability of the SSM as understood by Theorem
2.8. We first look at a linear case to compare the conditions for mean-square
contractivity between MV-SDEs and the numerical scheme. Then, the non-linear
Ginzburg Landau type equation gives some results to show the mean-square
contractivity for the split-step method. At last, the Cucker-Smale flocking
model (a degenerate MV-SDE) shows the split-step method has better properties
compared to the other methods under larger choices of timestep $h$.

######  Remark 3.1 (Parallel implementation).

To implement the split-step method (2.3)-(2.4) in parallel, at each timestep,
we first solve step (2.3) of the SSM distributed by the cores, then calculate
the empirical measure of the particle system in the first core, and finally
the second step (2.4) of the SSM is executed also in parallel. To implement
the taming method (A.1) and the time adaptive method(A.2) in parallel, in each
timestep $h$, one only need to first communicate the empirical measure of the
particle system between different cores, then each core can calculate its
particles dynamics independently.

A priori one should expect the taming to be the fastest of the three
algorithms. We highlight that due to the nature of the empirical measure,
after each time-step the processors need to communicate to update the
calculation of the empirical measure, this leads to a well-known loss of
parallelization power [9].

###  3.1 Example: the stochastic Ginzburg Landau equation

The first example we consider is the mean-field perturbation of the classic
one-dimensional stochastic Ginzburg Landau equation, namely, for all
$t\in[0,T]$,

$\displaystyle dX_{t}$
$\displaystyle=\Big{(}\frac{(\sigma^{\prime})^{2}}{2}X_{t}-X_{t}^{3}+c\mathbb{E}[X_{t}]\Big{)}dt+\sigma^{\prime}X_{t}dW_{t},\quad
X_{0}=x_{0}\in\mathbb{R}$ (3.1) $\displaystyle\textrm{where we set:}\quad
v(t,x)$ $\displaystyle=-x^{3},\quad
b(t,x,\mu)=\frac{(\sigma^{\prime})^{2}}{2}x+c\int_{\mathbb{R}}x\mu(dx),\quad\sigma(t,x,\mu)=\sigma^{\prime}x,$
(3.2)

where $\sigma^{\prime},c$ are constants.

This equation is a toy one-dimensional MV-SDE that we use for a methodological
comparison analysis. It features in [18, Section 4.1] and [45, Section 4.2],
thus a comparison is insightful under the choice of the same coefficients.
Namely, take $\sigma^{\prime}=1.5$, $x_{0}=1$, $c=0.5$ , $T=1$ and $N=1000$
particles. The result is shown in Figure 3.1. The timesteps are
$h\in\\{10^{-4},2\times 10^{-4},5\times 10^{-4},10^{-3},\dots,10^{-1}\\}$. The
true solution is calculated under $h=10^{-6}$.

######  Remark 3.2.

The map $v$ is a cubic polynomial and by Lemma 4.1 the implicit equation of
(2.3) has a unique solution. We determine it via an explicit formula for cubic
polynomials yielding a unique real root.

(a) Weak Error

(b) Strong Error (RMSE)

(c) Algorithm runtime

Figure 3.1:  Simulations of the stochastic Ginzburg Landau equation with
$N=1000$ particles. (a) Weak error of different methods. (b) Strong error of
different methods. (c) Runtime of different methods in serial and parallel
implementations.

Taming method is applied with $\alpha=0.5$ while the Time Adaptive method
under the choice $\mathbf{h}^{\delta}(x)=h\min(1,|x|^{-2})$. Fig 3.1(a) shows
the Weak error of the Taming method to roughly be $1/2$ with other methods
being 1.0. Fig 3.1(b) shows the strong error of all the methods to be around
$1/2$ but the taming method’s error is about one-order of magnitude higher
than the other two methods (this was also observed in [45, Section 4]). Fig
3.1(c) depicts running times of the 3 methods, the top 3 lines are the
standard implementations (non-parallel) and the bottom 3 lines are the
parallel implementation with 4 cores. In both parallel and non-parallel
implementation the taming scheme is the fastest while the SSM takes slightly
longer than the other methods but with a performance comparable to the time-
adaptive scheme. In the parallel implementation with 4 cores we reach a
reduction to nearly $25\%$ of the non-parallel implementation’s computation
time.

###  3.2 Example: the mean field FitzHugh-Nagumo model from Neuroscience

This is a three-dimensional ($d=3$) MV-SDE (1.1) defined with
$v:[0,T]\times\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$,
$b:[0,T]\times\mathbb{R}^{3}\times\mathcal{P}_{2}(\mathbb{R}^{3})\rightarrow\mathbb{R}^{3}$,
$\sigma:[0,T]\times\mathbb{R}^{3}\times\mathcal{P}_{2}(\mathbb{R}^{3})\rightarrow\mathbb{R}^{3\times
3}$ for $x=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$, $z\in\mathbb{R}$ and
$\mu^{x_{3}}$ is the marginal measure of $x_{3}$

$\displaystyle v(t,x)=\begin{pmatrix}-(x_{1})^{3}/3\\\ 0\\\
0\end{pmatrix},\quad
b(t,x,\mu)=\begin{pmatrix}x_{1}-x_{2}+I-\int_{\mathbb{R}}J(x_{1}-V_{rev})zd\mu^{x_{3}}(z)\\\
c(x_{1}+a-bx_{2})\\\
a_{r}\frac{T_{max}(1-x_{3})}{1+\exp(-\lambda(x_{1}-V_{T}))}-a_{d}x_{3}\end{pmatrix}$
$\displaystyle\sigma(t,x,\mu)=\begin{pmatrix}\sigma_{ext}&0&-\int_{\mathbb{R}}\sigma_{J}(x_{1}-V_{rev})zd\mu^{x_{3}}(z)\\\
0&0&0\\\ 0&\sigma_{32}(x)&0\end{pmatrix},\quad
x_{0}\sim\mathcal{N}\begin{pmatrix}\begin{pmatrix}V_{0}\\\ w_{0}\\\
y_{0}\end{pmatrix},\begin{pmatrix}&\sigma_{V_{0}}&0&0\\\
&0&\sigma_{w_{0}}&0\\\ &0&0&\sigma_{y_{0}}\end{pmatrix}\end{pmatrix}$

where $\mathcal{N}(\cdot,\cdot)$ denote the usual multivariate normal-
distribution,

$\displaystyle\sigma_{32}(x):=\mathds{1}_{\\{x_{3}\in(0,1)\\}}\sqrt{a_{r}\frac{T_{max}(1-x_{3})}{1+\exp(-\lambda(x_{1}-V_{T}))}+a_{d}x_{3}}\times\Gamma\exp\Bigg{[}-\frac{\Lambda}{1-(2x_{3}-1)^{2}}\Bigg{]},$

and
$I,J,V_{rev},V_{T},T_{\max},\Gamma,\Lambda,\sigma_{ext},\sigma_{J},\sigma_{32},a,b,c,a_{r},a_{d},\lambda$
are constants.

All the parameters are the same as in [18, Section 4.3]
($V_{0}=0,w_{0}=$1/2$,\dots$), also [45, Section 4.4]. The split structure of
the SSM enable us to flexibly deal with the super-linear term and the
Lipschitz terms separately (we make judicious choices regarding $v$ and $b$),
thus we can optimize the solver of the implicit part. Had we included the
terms of $b$ into $v$, then we would not be able to appeal to an explicit
cubic root formulae to solve (2.3), see Remark 3.2. In other words, we would
have been forced to use a general purpose solver which is computationally
costlier.

(a) Weak Error

(b) Strong Error (RMSE)

(c) Algorithm runtime

Figure 3.2:  Simulations of the FitzHugh-Nagumo model with $N=1000$ particles.
(a) Weak error of different methods. (b) Strong error of different methods.
(c) Runtime of different methods in serial and parallel.

For the simulation, we take $N=1000$, $T=2$, the time step is taken from
$h\in\\{10^{-4},2\times 10^{-4},5\times 10^{-4},10^{-3},\dots,10^{-1}\\}$ and
the true solution is calculated with $h=10^{-6}$. The Taming method is
implemented with $\alpha=1/2$ and the Adaptive method with
$\mathbf{h}^{\delta}(x)=h\min(1,|x|^{2}/|b(t,x,\mu)|^{2})$. Fig 3.2(a) shows
the weak error of the Taming method to roughly be $1/2$ with other methods
being 1.0. Fig 3.2(b) shows the strong error of all the methods to roughly be
$1/2$ and again the taming method’s error is one order of magnitude higher
than the others. Fig 3.2(c) shows the split-step method taking similar
computation time as the taming method and half the time of the adaptive
method. The implementation in parallel decreased the runtime of all schemes to
nearly $25\%$ of their non-parallel implementation runtimes when using 4
cores.

###  3.3 Example: Polynomial drift (non-Lipschitz measure dependency)

We present an example from [8, Section 3.3] which falls outside the theoretic
framework of this work. Take the one-dimensional MV-SDE for $t\in[0,T]$ and
$\gamma\in\mathbb{R}$,

$\displaystyle dX_{t}$ $\displaystyle=\Big{(}\gamma
X_{t}+\mathbb{E}\left[X_{t}\right]-X_{t}\mathbb{E}\big{[}\,|X_{t}|^{2}\big{]}\Big{)}dt+X_{t}dW_{t},\quad\text{with
}X(0)=x_{0}\in\mathbb{R},$ (3.3) $\displaystyle\textrm{where we set:}\quad
v(t,x)$ $\displaystyle=\gamma x-x\int_{\mathbb{R}}|x|^{2}\mu(dx),\quad
b(t,x,\mu)=\int_{\mathbb{R}}x\mu(dx),\quad\sigma(t,x,\mu)=x.$

The dynamics of the interacting particle system (1.2), for ${i\in\llbracket
1,N\rrbracket}$, $X^{i,N}\in\mathbb{R}$, is

$dX_{t}^{i,N}=\Big{[}\gamma
X_{t}^{i,N}+\frac{1}{N}\sum_{j=1}^{N}X_{t}^{j,N}-X_{t}^{i,N}\frac{1}{N}\sum_{j=1}^{N}|X_{t}^{j,N}|^{2}\Big{]}dt+X_{t}^{i,N}dW_{t}^{i}$

where $(W^{i})_{i}$ are independent Brownian motions. Take
$\gamma=-1,\leavevmode\nobreak\ T=1$. We have a super-linear measure component
in (3.3) and none of the schemes is applicable (insofar as existing
theoretical results allow). If the measure component was fixed (with finite
2nd moment), then, the drift would satisfy a one-sided Lipschitz condition –
the convergence of this type of scheme for this equation is left for future
work.

(a) Weak Error

(b) Strong Error (RMSE)

Figure 3.3: Simulations of the polynomial drift model with $N=1000$ particles.
(a) Weak error of different methods. (b) Strong error of different methods.

In terms of implementation of this example, at each timestep $t_{n}$ the
empirical measure $\hat{\mu}^{X,N}_{n}$ is fixed, then we proceed to solve the
split-step part for $Y^{\star}$. Specifically, we substitute (2.3) by

$\displaystyle Y_{n}^{i,\star,N}$
$\displaystyle=\hat{X}_{n}^{i,N}+hv\left(t_{n},Y_{n}^{i,\star,N},\hat{\mu}^{X,N}_{n}\right),\quad\quad\hat{\mu}^{X,N}_{n}(dx):=\frac{1}{N}\sum_{j=1}^{N}\delta_{\hat{X}_{n}^{j,N}}(dx),$

and afterwards apply step (2.4) as before.

The results are shown in Figure 3.3. We take $N=1000,\leavevmode\nobreak\ T=2$
and the timestep is taken from $h\in\\{10^{-4},2\times 10^{-4},5\times
10^{-4},10^{-3},\dots,10^{-1}\\}$. The true solution is calculated under
$h=10^{-6}$. Taming method is applied with $\alpha=0.5$ while the Time
Adaptive method under the choice $\mathbf{h}^{\delta}(x)=h\min(1,|x|^{-2})$.
Figure 3.3(a) shows Weak error of the Taming method with $\alpha=0.5$ to
roughly be $1/2$ with other methods being 1.0. Figure 3.3(b) shows strong
error of all the methods to roughly be $1/2$.

###  3.4 Stability of the SSM: linear, non-linear and the Cucker-Smale model
case

In the following cases, let $t\in[0,T]$, $X_{0},Z_{0}\in
L_{0}^{m}(\mathbb{R}^{d})$ with
$\hat{X}_{0}^{i},\hat{Z}_{0}^{i},\leavevmode\nobreak\ i\in\llbracket
1,N\rrbracket$ i.i.d. samples from $X_{0}$ and $Z_{0}$. The collections of
random variables $X_{n}^{i,N}$ and $Z_{n}^{i,N}$ are defined in Theorem 2.8
and follow from the SSM.

####  3.4.1 Linear case: Ornstein–Uhlenbeck McKean-Vlasov SDE

For the MV-SDE (see e.g. [8, Section 2.1]), for all $t\in[0,T]$ and
$x_{0}\in\mathbb{R}$

$\displaystyle dX_{t}$ $\displaystyle=\Big{(}\rho
X_{t}+\lambda\mathbb{E}[X_{t}]\Big{)}dt+\nu dW_{t},\quad X_{0}=X_{0},$ (3.4)
$\displaystyle\textrm{where we set:}\quad v(t,x)$ $\displaystyle=\rho x,\quad
b(t,x,\mu)=\lambda\int_{\mathbb{R}}x\mu(dx),\quad\sigma(t,x,\mu)=\nu,$ (3.5)

where $\rho,\lambda,\nu$ are constants. The first and second moments of $X$
are respectively given by
$\mathbb{E}[X_{t}]=\mathbb{E}[X_{0}]\exp((\rho+\lambda)t)$ and
$\mathbb{E}[X_{t}^{2}]=\mathbb{E}[X_{0}^{2}]\exp(2(\rho+\lambda)t)+\frac{\nu^{2}}{2\rho}(\exp(2\rho
t)-1)$.

Let $X,Z$ be two solution of (3.4) with $X_{0}$ and $Z_{0}$ as initial
condition respectively, then by direct calculation

$\displaystyle\mathbb{E}\Big{[}|X_{t}-Z_{t}|^{2}\Big{]}$
$\displaystyle=\frac{1}{2\lambda}e^{2(\rho+\lambda)t}+\mathbb{E}\Big{[}|X_{0}-Z_{0}|^{2}\Big{]}e^{2\rho
t}.$

Under $\rho\leq 0$, if $\rho+\lambda<0$ then from Definition 2.9, MV-SDE (3.4)
generates exponential mean-square contractive solutions. Also, the parameters
of this example are $L_{v}=\rho,\leavevmode\nobreak\ L_{\hat{b}}=\lambda^{2}$,
$L_{b}=L_{{g}}=L_{\hat{\sigma}}=0$. Plug these into (2.6) and in order to make
$\beta<0$, we need to choose $h$ such that

$\displaystyle
h<-\frac{2\rho+1+\lambda^{2}}{\lambda^{2}}=-1-\frac{2\rho+1}{\lambda^{2}}.$

Then, from Definition 2.7 the SSM (2.3)–(2.4) is mean-square contractive. So,
for the SSM, $h$ exists when $1+\lambda^{2}<-2\rho$. Both this condition and
the condition for the SDE to generate exponential mean-square contractive
solutions need $\rho<0$, thus when $\lambda\neq 1$, the condition for a mean-
square stable numerical solution is stronger than the condition for the MV-SDE
to generate exponential mean-square contractive solutions. When $\lambda=1$
then the two conditions coincide.

####  3.4.2 Nonlinear case I: a stochastic Ginzburg Landau type equation

To illustrate the stability of SSM scheme with a stochastic Ginzburg Landau
type equation (see Section 3.1), we consider the following one-dimensional MV-
SDE for all $t\in[0,T]$,

$\displaystyle dX_{t}$
$\displaystyle=\Big{(}-\frac{5}{2}X_{t}-\frac{1}{4}X_{t}^{3}+\mathbb{E}[X_{t}]\Big{)}dt+X_{t}dW_{t},\qquad
X_{0}=1,$ (3.6) $\displaystyle\textrm{where we set:}\quad v(t,x)$
$\displaystyle=-\frac{5}{2}x-\frac{1}{4}x^{3}+\gamma x,\quad
b(t,x,\mu)=\int_{\mathbb{R}}x\mu(dx)-\gamma x\ \textrm{ for }\
\gamma\in\mathbb{R},\quad\sigma(t,x,\mu)=x.$

The parameters of this example are $L_{v}=-{5}/{2}+\gamma$,
$L_{\hat{b}}=L_{\sigma}=1$, and $L_{b}=-\gamma,L_{\hat{\sigma}}=0$. Plug these
into (2.6) and when $\gamma=0$, in order to make $\beta<0$ one must have
$0<h<2$. We now employ the SSM to simulate $N=1000$ particles, under different
choices of $h$ and different initial values and provide the results in Figure
3.4. The SSM’s step (2.3) is solved via Remark 3.2. Figure 3.4(a) shows, when
$\gamma=0$, that $t\mapsto\sup_{i\in\llbracket
1,N\rrbracket}|X_{t}^{i,N}|^{2}$ decreases to zero under the different choices
of $h$, which is as expected. In Figure 3.4(b), under $\gamma=0$ and
$X_{0}=1$, the Mean square differences between $X$ and $Z$ also decrease
exponentially as time goes by.

Figure 3.4(c) shows that modifying (3.6) with $\gamma\neq 0$, and thus
$L_{v}=\gamma-5/2$, does not significantly change the result even when
$1-2hL_{v}>0$ is not satisfied. Nonetheless, when $1-2hL_{v}$ is near $0$ the
scheme becomes more unstable deviation. This indicates the SSM can be flexibly
applied and allows for modified Lipschitz terms in $v,b$ but one needs to
consider the restriction carefully.

(a) $\sup_{i\in\llbracket 1,N\rrbracket}|X_{t}^{i,N}|^{2}$

(b) $\sup_{i\in\llbracket 1,N\rrbracket}|X_{t}^{i,N}-Z_{t}^{i,N}|^{2}$

(c) $\gamma\mapsto\sup_{i\in\llbracket 1,N\rrbracket}|X_{T}^{i,N}|^{2}$

Figure 3.4: Simulations of the stochastic Ginzburg Landau type equation with
$N=1000$ particles. (a) shows $\sup_{i\in\llbracket
1,N\rrbracket}|X_{t}^{i,N}|^{2}$ under three different timesteps
$h\in\\{10^{-2},10^{-3},10^{-4}\\}$. (b) shows the mean square differences
between $X,Z$ when $h=0.01$ for fixed $X_{0}=1$ and $Z_{0}\in\\{5,10,100\\}$.
(c) shows $[-100,100]\ni\gamma\mapsto\sup_{i\in\llbracket
1,N\rrbracket}|X_{T}^{i,N}|^{2}$ under $h=0.01,\leavevmode\nobreak\ T=3$.

####  3.4.3 Nonlinear case II: the two-dimensional Cucker-Smale flocking model

This example, from [22, Section 2], highlights the stability of the SSM under
larger timestep $h$. The SSM is stable under larger timesteps $h$ by using the
implicit step to control for the super-linear part. The explicit methods (the
taming and the adaptive) fail to have acceptable results at this level of $h$.

Applied our settings, this is a two-dimensional MV-SDE defined under
$v:[0,T]\times\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$,
$b,\sigma:[0,T]\times\mathbb{R}^{2}\times\mathcal{P}_{2}(\mathbb{R}^{2})\rightarrow\mathbb{R}^{2}$
for $x=(V,X)\in\mathbb{R}^{2}\leavevmode\nobreak\ ,z\in\mathbb{R}$, $\mu^{V}$
is the measure of $V$ as:

$\displaystyle v(t,x)=\begin{pmatrix}-(V)^{3}\\\ 0\end{pmatrix},\ \
b(t,x,\mu)=\begin{pmatrix}\lambda\int_{\mathbb{R}}(V-z)d\mu^{V}(z)\\\
V\end{pmatrix},\ \textrm{ and }\
\sigma(t,x,\mu)=\begin{pmatrix}\sigma^{\prime}\int_{\mathbb{R}}(V-z)d\mu^{V}(z)\\\
0\end{pmatrix},$

where $\lambda,\sigma^{\prime}$ are constants.

The dynamics of the particle system follows easily

$\displaystyle dV_{t}^{i,N}$
$\displaystyle=\Big{[}-(V_{t}^{i,N})^{3}+\frac{\lambda}{N}\sum_{j=1}^{N}(V_{t}^{j,N}-V_{t}^{i,N})\Big{]}dt+\frac{\sigma^{\prime}}{N}\sum_{j=1}^{N}(V_{t}^{j,N}-V_{t}^{i,N})dW_{t}^{i},\quad
dX_{t}^{i,N}=V_{t}^{i,N}dt,$

where $i\in\llbracket 1,N\rrbracket,V^{i,N},X^{i,N}\in\mathbb{R}$,
$(W^{i})_{i}$ are independent Brownian motions. Take
$\lambda=2,\leavevmode\nobreak\ \sigma^{\prime}=4,\leavevmode\nobreak\ T=5$,
$V_{0}^{i,N}\sim\mathcal{N}(0,1)$ are i.i.d. samples from standard normal
distribution, the timestep takes values from $h\in\\{10^{-3},2\times
10^{-3},5\times 10^{-3},10^{-3},\dots,10^{-1}\\}$, and we simulate $N=1000$
particles. The SSM’s step (2.3) is solved via Remark 3.2. The true solution is
calculated under $h=10^{-4}$. The Taming method with $\alpha=0.5$ and $1$, the
Adaptive method with $\mathbf{h}^{\delta}(x)=h\min(1,|x|^{-2})$.

Figure 3.5(a) and (b) show that for larger timesteps $h=0.1$ the SSM is stable
whereas both the Adaptive and the Taming method blow up (error large than
$10^{4}$). The strong error of the SSM is consistently small across the
timesteps but only for smaller timesteps, $h<0.05$, do the time-adaptive and
the Taming method with $\alpha=1$ have low errors. The taming method with
$\alpha=0.5$ failed to have a good result. This seems to indicate that the SSM
has better properties on long time simulation under larger timestep than the
other two methods.

(a) Strong Error (RMSE) of $V$

(b) Strong Error (RMSE) of $X$

Figure 3.5:  Simulations of the Cucker-Smale flocking model with $N=1000$
particles. (a) Strong error of $V$ under different methods. (b) Strong error
of $X$ under different methods.

##  4 Proof of the convergence result for the split-step method (SSM)

Throughout this section Assumption 2.1 is assumed to hold for all results.

###  4.1 Proof of the main convergence result, Theorem 2.6

For all auxiliary results next, we assume the conditions of Theorem 2.6 are in
force and we thus do not state them.

#### Preliminary results

As a first step, we state a result that allows us to re-write (2.3) and (2.4)
as a map of $\hat{X}^{i,N}$ without the presence of the $Y^{i,\star,N}$. We
present first a new general version of [28, Lemma 3.4] where the
differentiability assumption is lifted and the maps are allowed to depend on
time.

######  Lemma 4.1.

Choose $h>0$ satisfying $1-2hL_{v}>0$. Given $t\in[0,T],\leavevmode\nobreak\
c,d\in\mathbb{R}^{d}$ the implicit equation:

$\displaystyle c=d+hv(t,c),$ (4.1)

has a unique solution $c$. Define the functions

$\displaystyle[0,T]\times\mathbb{R}^{d}\ni(t,d)\mapsto
F_{h}(t,d)=c\in\mathbb{R}^{d},\quad v_{h}(t,d)=v(t,F_{h}(t,d)).$ (4.2)

We then have for all $t\in[0,T]$, $x,y\in\mathbb{R}^{d}$ the following three
inequalities

$\displaystyle|v_{h}(t,x)|\leq$ $\displaystyle\frac{|v(t,x)|}{1-hL_{v}},$
(4.3) $\displaystyle|F_{h}(t,x)-F_{h}(t,y)|^{2}\leq$
$\displaystyle\frac{|x-y|^{2}}{1-2hL_{v}},$ (4.4) $\displaystyle\langle
x-y,v_{h}(t,x)-v_{h}(t,y)\rangle\leq$
$\displaystyle\frac{L_{v}}{1-2hL_{v}}|x-y|^{2}.$ (4.5)

For all $t\in[0,T]$, $x_{i},y_{i}\in\mathbb{R}^{d},\leavevmode\nobreak\
i\in\llbracket 1,N\rrbracket$ and
$\mu^{x},\mu^{y}\in\mathcal{P}_{2}(\mathbb{R}^{d})$ are empirical measure
associated with $\\{x_{i}\\}_{i},\\{y_{i}\\}_{i}$. Define the maps

$\displaystyle
b_{h}(t,x_{i},\mu^{F_{h,x}})=b(t,F_{h}(t,x_{i}),\mu^{F_{h,x}}),\qquad\sigma_{h}(t,x_{i},\mu^{F_{h,x}})=\sigma(t,F_{h}(t,x_{i}),\mu^{F_{h,x}}),$
$\displaystyle\mu^{F_{h,x}}(dx)=\frac{1}{N}\
\sum_{j=1}^{N}\delta_{F_{h}(t,x_{j})}(dx),\quad\mu^{F_{h,y}}(dx)=\frac{1}{N}\
\sum_{j=1}^{N}\delta_{F_{h}(t,y_{j})}(dx).$

Then, $b_{h}$ and $\sigma_{h}$ are satisfy

$\displaystyle\Big{|}b_{h}(t,x_{i},\mu^{F_{h,x}})-b_{h}(t,y_{i},\mu^{F_{h,y}})\Big{|}^{2}\leq\frac{L_{b}}{1-2hL_{v}}\Bigg{[}|x_{i}-y_{i}|^{2}+\Big{(}\frac{1}{N}\sum_{j=1}^{N}|x_{j}-y_{j}|^{2}\Big{)}\Bigg{]},$
(4.6)
$\displaystyle\Big{|}\sigma_{h}(t,x_{i},\mu^{F_{h,x}})-\sigma_{h}(t,y_{i},\mu^{F_{h,y}})\Big{|}^{2}\leq\frac{L_{\sigma}}{1-2hL_{v}}\Bigg{[}|x_{i}-y_{i}|^{2}+\Big{(}\frac{1}{N}\sum_{j=1}^{N}|x_{j}-y_{j}|^{2}\Big{)}\Bigg{]}.$
(4.7)

Lastly, $v_{h}\rightarrow v$, $b_{h}\rightarrow b$, and
$\sigma_{h}\rightarrow\sigma$ uniformly over the compacts of
$[0,T]\times\mathbb{R}^{d}$ as $h\rightarrow 0^{+}$.

We observe that to establish (4.3) one only needs $1-hL_{v}>0$, in other words
the condition $1-2hL_{v}>0$ is not sharp. Nonetheless, inequalities (4.4) and
(4.5) are critical for our work and hence we right one single constraint.

###### Proof.

Existence and uniqueness for (4.1) can be proved from a via a strict
monotonicity contraction argument. Namely, fix some $t\in[0,T]$, from (4.1)
one defines the operator $A:\mathbb{R}^{d}\to\mathbb{R}^{d}$ as
$A(u)=u-hv(t,u)$ for $u\in\mathbb{R}^{d}$. Following [48, Definition 25.2
(p500)], the operator $A$ is continuous and strongly monotone (uniformly in
$t$) under the restrictions $h>0$ and $1-hL_{v}>0$. This follows by directly
injecting the one-sided condition of $v$ (from Assumption 2.1) in the
definition of strongly monotone operator. Finally, from [48, Theorem 26.A
(p557)] we conclude that the operator $A$ is invertible and the inverse map is
Lipschitz continuous. Thus, (4.1) has a unique measurable inverse given by
$F_{h}$ from (4.2). See also [38, (p2596)].

We now determine the Lipschitz constant of $v_{h}$ and $F_{h}$ of (4.2). For
(4.3), suppose $c,d\in\mathbb{R}^{d}$ satisfy $c=d+hv(t,c)$ then:

$\displaystyle c$ $\displaystyle=d+h\Big{(}v(t,c)-v(t,d)\Big{)}+hv(t,d),$
$\displaystyle|c-d|^{2}$
$\displaystyle=\big{\langle}c-d,v(t,c)-v(t,d)\big{\rangle}h+\langle
c-d,v(t,d)\rangle h\leq hL_{v}|c-d|^{2}+\langle c-d,v(t,d)\rangle h.$

Since $c=d+hv_{h}(t,d)$ we have

$\displaystyle|v_{h}(t,d)|=\frac{1}{h}|c-d|\leq\frac{|v(t,d)|}{1-hL_{v}}.$

For (4.4), suppose $c_{1},c_{2},d_{1},d_{2}\in\mathbb{R}^{d}$ satisfy
$c_{1}=d_{1}+hv(t,c_{1})$ and $c_{2}=d_{2}+hv(t,c_{2})$, then

$\displaystyle\Big{|}F_{h}(t,d_{1})-F_{h}(t,d_{2})\Big{|}^{2}=|c_{1}-c_{2}|^{2}$
$\displaystyle=\langle
c_{1}-c_{2},d_{1}-d_{2}\rangle+\big{\langle}c_{1}-c_{2},v(t,c_{1})-v(t,c_{2})\big{\rangle}h$
$\displaystyle\leq\frac{1}{2}|c_{1}-c_{2}|^{2}+\frac{1}{2}|d_{1}-d_{2}|^{2}+hL_{v}|c_{1}-c_{2}|^{2}.$

For (4.5), we employ the results above and find that

$\displaystyle\Big{\langle}d_{1}-d_{2},(d_{1}-d_{2})+h\big{(}v_{h}(t,d_{1})-v_{h}(t,d_{2})\big{)}\Big{\rangle}$
$\displaystyle=\langle
d_{1}-d_{2},c_{1}-c_{1}\rangle\leq\frac{1}{2}|d_{1}-d_{2}|^{2}+\frac{1}{2}\frac{|d_{1}-d_{2}|^{2}}{1-2hL_{v}}$

and thus

$\displaystyle\big{\langle}d_{1}-d_{2},v_{h}(t,d_{1})-v_{h}(t,d_{2})\big{\rangle}$
$\displaystyle=\frac{L_{v}}{1-2hL_{v}}|d_{1}-d_{2}|^{2}.$

Since $b_{h}$ and $\sigma_{h}$ are of the same nature, we show the proof for
$b_{h}$ as the proof for $\sigma_{h}$ is identical. Then, let $i\in\llbracket
1,N\rrbracket$, using the definition of $b_{h}$, the Lipschitz property of $b$
(Assumption 2.1), and (4.4) we have

$\displaystyle\Big{|}b_{h}(t,x_{i},\mu^{F_{h,x}})-b_{h}(t,y_{i},\mu^{F_{h,y}})\Big{|}^{2}$
$\displaystyle\leq
L_{b}\Big{(}|F_{h}(t,x_{i})-F_{h}(t,y_{i})|^{2}+W^{(2)}(\mu^{F_{h,x}},\mu^{F_{h,y}})^{2}\Big{)}$
$\displaystyle\leq
L_{b}\Big{(}|F_{h}(t,x_{i})-F_{h}(t,y_{i})|^{2}+\frac{1}{N}\sum_{j=1}^{N}|F_{h}(t,x_{j})-F_{h}(t,y_{j})|^{2}\Big{)}$
$\displaystyle\leq\frac{L_{b}}{1-2hL_{v}}\Big{(}|x_{i}-y_{i}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|x_{j}-y_{j}|^{2}\Big{)}.$

The convergence result is classical.

∎

After having introduced Lemma 4.1, we can finally address the continuous-time
extension of the SSM (2.3)-(2.4) has referenced in Theorem 2.6. The SSM can be
written as a continuous time SDE via linear interpolation of the iterates,
namely, for $t\in[t_{n},{t_{n+1}}]$, $i\in\llbracket
1,N\rrbracket,\leavevmode\nobreak\ \hat{X}_{0}^{i}\in
L_{0}^{m}(\mathbb{R}^{d})$:

$\displaystyle d\hat{X}_{t}^{i,N}$
$\displaystyle=\Big{(}v_{h}(\kappa(t),\hat{X}_{\kappa(t)}^{i,N})+b_{h}(\kappa(t),\hat{X}_{\kappa(t)}^{i,N},\widetilde{\mu}_{\kappa(t)}^{N})\Big{)}dt+\sigma_{h}(\kappa(t),\hat{X}_{\kappa(t)}^{i,N},\widetilde{\mu}_{\kappa(t)}^{N})dW_{t}^{i},$
(4.8) $\displaystyle\quad\widetilde{\mu}^{N}_{\kappa(t)}(dx)$
$\displaystyle:=\frac{1}{N}\
\sum_{j=1}^{N}\delta_{F_{h}(\kappa(t),\hat{X}_{\kappa(t)}^{j,N})}(dx),\leavevmode\nobreak\
\kappa(t)=\inf\big{\\{}t_{n}:t_{n}\leq t,\ n\in\llbracket
0,M\rrbracket\big{\\}},\leavevmode\nobreak\
\widetilde{\mu}^{N}_{t_{n}}=\widetilde{\mu}^{N}_{n}.$

#### Moment bounds

We now employ the results of Lemma 4.1 to establishing a domination of
$|Y_{n}^{i,\star,N}|$ by $|\hat{X}_{n}^{i,N}|$.

######  Lemma 4.2.

Choose $h$ as in Theorem 2.6 and recall $C_{T},\widehat{L}_{v}$ as defined in
Remark 2.2.

Then, $|Y_{n}^{i,\star,N}|$ in (2.3) satisfies

$\displaystyle|Y_{n}^{i,\star,N}|^{2}$
$\displaystyle\leq|\hat{X}_{n}^{i,N}|^{2}(1+\frac{2\widehat{L}_{v}}{1-2\widehat{L}_{v}h}h)+\frac{2C_{T}}{1-2\widehat{L}_{v}h}h,\quad\textrm{for
any }i\in\llbracket 1,N\rrbracket,\ n\in\llbracket 0,M-1\rrbracket.$ (4.9)

######  Remark 4.3.

Recall the stepsize constraint on $h$ in Theorem 2.6. By inspection of the
proof of this Lemma and the definition of $\widehat{L}_{v}$ in Remark 2.2, the
reader will find that the constraint $L_{v}<-1/2$ needed to ensure
$1-h\widehat{L}_{v}>0$ is not sharp and can be replaced by some number
$\varepsilon\in(0,1)$, i.e., $L_{v}<-\varepsilon$.

The lack of sharpness arises in the choice of $\widehat{L}_{v}$ in Remark 2.2.
There we used the Cauchy-Schwarz inequality where we could have used a Young
type inequality from which the parameter $\varepsilon$ would have arisen. We
choose $1/2$ for ease of presentation.

###### Proof.

From (2.3) and Remark 2.2, for any $i,n$ and any $t_{n}\in\pi$ we have using
Cauchy-Schwarz and the properties of $v$

$\displaystyle|Y_{n}^{i,\star,N}|^{2}$ $\displaystyle=\langle
Y_{n}^{i,\star,N},\hat{X}_{n}^{i,N}\rangle+\langle
Y_{n}^{i,\star,N},v(t_{n},Y_{n}^{i,\star,N})\rangle h$
$\displaystyle\leq\frac{1}{2}|Y_{n}^{i,\star,N}|^{2}+\frac{1}{2}|\hat{X}_{n}^{i,N}|^{2}+h(C_{T}+\widehat{L}_{v}|Y_{n}^{i,\star,N}|^{2})\leq|\hat{X}_{n}^{i,N}|^{2}(1+\frac{2\widehat{L}_{v}}{1-2\widehat{L}_{v}h}h)+\frac{2C_{T}}{1-2\widehat{L}_{v}h}h.$

∎

######  Remark 4.4.

Under Assumption 2.1 and $h$ as in Theorem 2.6, there exist constants
$C>0,\widetilde{C}\in\mathbb{R}$ depending on $L_{v},C_{T},\widehat{L}_{v}$
but independent of $h$ such that in Lemma 4.1 and Lemma 4.2 we have

$\displaystyle|v_{h}(t,x)|\leq C{|v(t,x)|},\quad$
$\displaystyle|F_{h}(t,x)-F_{h}(t,y)|^{2}\leq C|x-y|^{2},$
$\displaystyle\big{\langle}x-y,v_{h}(t,x)-v_{h}(t,y)\big{\rangle}\leq\widetilde{C}|x-y|^{2},\quad$
$\displaystyle|Y_{k}^{i,\star,N}|^{2}\leq|\hat{X}_{k}^{i,N}|^{2}(1+Ch)+Ch.$

After this remark emphasizing the independence of the constants in $h$, we are
in a position to finally prove the moment bounds for the output of the SSM.

######  Proposition 4.5 (Moment bounds of SSM).

Choose $h$ as in Theorem 2.6. Then for any $i\in\llbracket 1,N\rrbracket,\
n\in\llbracket 0,M\rrbracket$, $\hat{X}_{n}^{i,N}$ of (2.3)-(2.4) satisfies,
for $p\geq 1$ (for an $m$ sufficiently large depending on $p$), there exist a
constant $C>0$ independent of $h,N,M$, such that

$\displaystyle\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\sup_{n\in\llbracket
0,M\rrbracket}|\hat{X}_{n}^{i,N}|^{2p}\Big{]}<C\Big{(}1+\mathbb{E}[\,|\hat{X}_{0}|^{2p}]\Big{)}<\infty.$

###### Proof.

Let $i\in\llbracket 1,N\rrbracket$, ${n\in\llbracket 0,M-1\rrbracket}$ and
recall (2.4). Using Assumption 2.1, Remark 4.4, letting the constant $C>0$ be
independent of $h,N,M$ and possibly changing from line by line, we have by
taking squares

$\displaystyle|\hat{X}_{n+1}^{i,N}|^{2}$
$\displaystyle=|Y_{n}^{i,\star,N}|^{2}+\Big{|}b(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})h+\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})\Delta
W_{n}^{i}\Big{|}^{2}$
$\displaystyle+2\Big{\langle}Y_{n}^{i,\star,N},b(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})h+\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})\Delta
W_{n}^{i}\Big{\rangle}$
$\displaystyle\leq|\hat{X}_{n}^{i,N}|^{2}(1+Ch)+Ch+2|b(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})|^{2}h^{2}+2|\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})|^{2}(\Delta
W_{n}^{i})^{2}$
$\displaystyle+2\big{\langle}Y_{n}^{i,\star,N},\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})\Delta
W_{n}^{i}\big{\rangle}+2C\Big{(}1+|Y_{n}^{i,\star,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|Y_{n}^{j,\star,N}|^{2}\Big{)}h$
$\displaystyle\leq|\hat{X}_{n}^{i,N}|^{2}+C\Big{(}1+|\hat{X}_{n}^{i,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{n}^{j,N}|^{2}\Big{)}\Big{(}(\Delta
W_{k}^{i})^{2}+h\Big{)}+2\big{\langle}Y_{n}^{i,\star,N},\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}_{n}^{N})\Delta
W_{n}^{i}\big{\rangle}$

By backward induction from $n+1$ to zero, we have (after some simplification)

$\displaystyle|\hat{X}_{n+1}^{i,N}|^{2}\leq$
$\displaystyle|\hat{X}_{0}^{i,N}|^{2}+2\sum_{k=0}^{n}\Big{\\{}\big{\langle}Y_{k}^{i,\star,N},\sigma(t_{k},Y_{k}^{i,\star,N},\hat{\mu}_{k}^{N})\Delta
W_{k}^{i}\big{\rangle}\Big{\\}}$
$\displaystyle+C\sum_{k=0}^{M-1}\Big{\\{}(\Delta
W_{k}^{i})^{2}\Big{\\}}+C\sum_{k=0}^{M-1}\Big{\\{}|\hat{X}_{k}^{i,N}|^{2}(\Delta
W_{k}^{i})^{2}\Big{\\}}+C\sum_{k=0}^{M-1}\Big{\\{}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}(\Delta
W_{k}^{i})^{2}\Big{\\}}$
$\displaystyle+C\sum_{k=0}^{M-1}\Big{\\{}h\Big{\\}}+C\sum_{k=0}^{M-1}\Big{\\{}|\hat{X}_{k}^{i,N}|^{2}h\Big{\\}}+C\sum_{k=0}^{M-1}\Big{\\{}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}h\Big{\\}}$

Taking power $p\geq 2$ on both sides, expectations and re-organizing the
terms, we have (with $C_{p}$ independent of $h,N$)

$\displaystyle\mathbb{E}\Big{[}|\hat{X}_{n+1}^{i,N}|^{2p}\Big{]}$
$\displaystyle\leq$ $\displaystyle
C_{p}\Bigg{[}\mathbb{E}\Big{[}|\hat{X}_{0}^{i,N}|^{2p}\Big{]}+\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{n}\Big{\\{}\langle
Y_{k}^{i,\star,N},\sigma(t_{k},Y_{k}^{i,\star,N},\hat{\mu}_{k}^{N})\Delta
W_{k}^{i}\rangle\Big{\\}}\Big{)}^{p}\Big{]}+\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\Big{\\{}(\Delta
W_{k}^{i})^{2}\Big{\\}}\Big{)}^{p}\Big{]}$
$\displaystyle+\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\Big{\\{}|\hat{X}_{k}^{i,N}|^{2}(\Delta
W_{k}^{i})^{2}\Big{\\}}\Big{)}^{p}\Big{]}+\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\Big{\\{}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}(\Delta
W_{k}^{i})^{2}\Big{\\}}\Big{)}^{p}\Big{]}$
$\displaystyle+\mathbb{E}\Big{[}\Big{(}p\sum_{k=0}^{M-1}\Big{\\{}h\Big{\\}}\Big{)}^{p}\Big{]}+\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\Big{\\{}|\hat{X}_{k}^{i,N}|^{2}h\Big{\\}}\Big{)}^{p}\Big{]}+\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\Big{\\{}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}h\Big{\\}}\Big{)}^{p}\Big{]}\Bigg{]}.$

There are 8 terms to be estimated, but in essence only 3 arguments are needed.
We present them only for the most complex terms since for the remaining ones
it is just a simplification of the arguments presented. We start with the last
term of the 2nd line: apply Jensen’s inequality twice after scaling the outer
summation and then tower property to take advantage of the independence
between $\hat{X}_{k}^{j,N}$ and $\Delta W_{k}^{i}$, namely, (recall $h=T/M$)

$\displaystyle\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\Big{\\{}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}(\Delta
W_{k}^{i})^{2}\Big{\\}}\Big{)}^{p}\Big{]}$
$\displaystyle\leq\mathbb{E}\Big{[}\frac{1}{M}\sum_{k=0}^{M-1}\Big{\\{}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2p}(\Delta
W_{k}^{i})^{2p}\Big{\\}}M^{p}\Big{]}$
$\displaystyle=h\mathbb{E}\Big{[}\sum_{k=0}^{M-1}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2p}\Big{]}CT^{p}.$

We now address the second term in the 1st line. Using Burkholder–Davis–Gundy
(BDG) inequality, the Lipschitz property of $\sigma$, Jensen’s inequality and
the domination of $Y_{k}^{i,\star,N}$ by $\hat{X}_{k}^{i,N}$ reported on Lemma
4.2 (and Remark 4.4) yields

$\displaystyle\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{n}\Big{\\{}$
$\displaystyle\langle
Y_{k}^{i,\star,N},\sigma(t_{k},Y_{k}^{i,\star,N},\hat{\mu}_{k}^{N})\Delta
W_{k}^{i}\rangle\Big{\\}}\Big{)}^{p}\Big{]}$ $\displaystyle\leq$
$\displaystyle C\mathbb{E}\Big{[}\sum_{k=0}^{M-1}\Big{\\{}|\langle
Y_{k}^{i,\star,N},\sigma(t_{k},Y_{k}^{i,\star,N},\hat{\mu}_{k}^{N})\rangle|^{2}h\Big{\\}}\Big{]}^{\frac{p}{2}}$
$\displaystyle\leq$ $\displaystyle
C\mathbb{E}\Big{[}\frac{1}{M}\sum_{k=0}^{M-1}\Big{\\{}1+|Y_{k}^{i,\star,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|Y_{k}^{j,\star,N}|^{2}\Big{\\}}h\Big{]}^{\frac{p}{2}}M^{\frac{p}{2}}$
$\displaystyle\leq$ $\displaystyle
C\mathbb{E}\Bigg{[}1+\Big{[}\sum_{k=0}^{M-1}|\hat{X}_{k}^{i,N}|^{p}\Big{]}h+\frac{1}{M}\Big{[}\sum_{k=0}^{M-1}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{p}\Big{]}h\Bigg{]}$

Finally, we address the last term in the 3rd line of the inequality. The
result follows by applying Jensen’s inequality

$\displaystyle\mathbb{E}\Big{[}\Big{(}\sum_{k=0}^{M-1}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}h\Big{)}^{p}\Big{]}$
$\displaystyle=\mathbb{E}\Big{[}\Big{(}\frac{1}{M}\sum_{k=0}^{M-1}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2}h\Big{)}^{p}M^{p}\Big{]}\leq
h\mathbb{E}\Big{[}\sum_{k=0}^{M-1}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2p}\Big{]}T^{p}.$

Collecting the several inequalities and injecting them back in the initial
inequality, we conclude that

$\displaystyle\mathbb{E}\Big{[}\,|\hat{X}_{n+1}^{i,N}|^{2p}\Big{]}\leq$
$\displaystyle
C\Bigg{[}1+\mathbb{E}\Big{[}\,|\hat{X}_{0}^{i,N}|^{2p}\Big{]}+\mathbb{E}\Big{[}\sum_{k=0}^{M-1}|\hat{X}_{k}^{i,N}|^{2p}h\Big{]}+\mathbb{E}\Big{[}\sum_{k=0}^{M-1}\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{k}^{j,N}|^{2p}h\Big{]}\Bigg{]}.$

Taking supremum and using that the particles are conditional i.d.d. (for fixed
$k$)

$\displaystyle\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\sup_{n\in\llbracket
0,M\rrbracket}|\hat{X}_{n}^{i,N}|^{2p}\Big{]}$ $\displaystyle\leq
C\Bigg{[}1+\mathbb{E}\Big{[}\,|\hat{X}_{0}^{i,N}|^{2p}\Big{]}+\sum_{k=0}^{M-1}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq n\leq
k}|\hat{X}_{n}^{i,N}|^{2p}\Big{]}h\Bigg{]}.$

The proof finishes after applying the discrete Grönwall’s inequality to the
inequality. ∎

We now provide (moment) estimates for the continuous-time extension of the
SSM.

######  Proposition 4.6.

Choose $h$ as in Theorem 2.6. Take $(\hat{X}_{t}^{i,N})_{t\in[0,T]}$ as the
map satisfying (4.8), i.e., the continuous time extension of the SSM.,

Then, for any $p\geq 1$ (with $m$ sufficiently large depending on $p$), there
exist a constant $C>0$ independent of $h,N,M$, such that:

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq
T}|\hat{X}_{t}^{i,N}|^{2p}\Big{]}<C\Big{(}1+\mathbb{E}\Big{[}\,|\hat{X}_{0}|^{2p}\Big{]}\Big{)}<\infty.$
(4.10)

###### Proof.

Let $i\in\llbracket 1,N\rrbracket$ and ${n\in\llbracket 0,M-1\rrbracket}$.
From (4.8), set $t_{n}+s=t$, for all $t\in[0,T],n\in\llbracket
0,M-1\rrbracket$, then

$\displaystyle\hat{X}_{t}^{i,N}=\hat{X}_{n}^{i,N}+v_{h}(t_{n},\hat{X}_{n}^{i,N})s+b_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})s+\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})(W_{t_{n}+s}^{i}-W_{t_{n}}^{i})$
(4.11)

Plugging (2.3)

in (4.11) gives

$\displaystyle\hat{X}_{t}^{i,N}=\hat{X}_{n}^{i,N}(1-\frac{s}{h})+\frac{s}{h}Y_{n}^{i,\star,N}+b_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})s+\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})(W_{t_{n}+s}^{i}-W_{t_{n}}^{i}).$

Since $s\leq h$, by the Lipschitz conditions on $b_{h}$ and $\sigma_{h}$, and
Lemma 4.2, we have

$\displaystyle|\hat{X}_{t}^{i,N}|^{2}$ $\displaystyle\leq
C\Big{(}1+|\hat{X}_{n}^{i,N}|^{2}+|b_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})|^{2}+\Big{|}\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})(W_{t_{n}+s}^{i}-W_{t_{n}}^{i})\Big{|}^{2}\Big{)}$
$\displaystyle\leq
C\Big{(}1+|\hat{X}_{n}^{i,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{n}^{j,N}|^{2}+\Big{|}\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})(W_{t_{n}+s}^{i}-W_{t_{n}}^{i})\Big{|}^{2}\Big{)}.$

Taking supremum over time and expectations on both sides, yields

$\displaystyle\mathbb{E}\Big{[}\sup_{0\leq t\leq
T}|\hat{X}_{t}^{i,N}|^{2p}\Big{]}$ $\displaystyle=\mathbb{E}\Big{[}\sup_{0\leq
n\leq M}\sup_{0\leq s\leq h}|\hat{X}_{t_{n}+s}^{i,N}|^{2p}\Big{]}$
$\displaystyle\leq C\mathbb{E}\Bigg{[}1+\sup_{0\leq n\leq
M}\Big{\\{}\Big{(}|\hat{X}_{n}^{i,N}|^{2p}+\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{n}^{j,N}|^{2}\Big{)}+I_{h}\Big{\\}}\Bigg{]},$

where $I_{h}$ is given by $I_{h}:=\sup_{0\leq s\leq
h}\big{|}\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})(W_{t_{n}+s}^{i}-W_{t_{n}}^{i})\big{|}^{2p}$.
Using the BDG inequality, the Lipschitz condition on $\sigma_{h}$ and
Proposition 4.5, gives

$\displaystyle\mathbb{E}[I_{h}]$ $\displaystyle\leq
C\mathbb{E}\Big{[}\Big{(}|\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})|^{2}h\Big{)}^{p}\Big{]}\leq
Ch^{p}\mathbb{E}\Big{[}1+|\hat{X}_{n}^{i,N}|^{2p}+\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{n}^{j,N}|^{2p}\Big{]}\leq
Ch^{p}.$

Take supremum over $i\in\llbracket 1,N\rrbracket$, by Proposition 4.5 once
again there is a constant $C$ such that

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq T}|\hat{X}_{t}^{i,N}|^{2p}\Big{]}$
$\displaystyle\leq\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}C(1+\sup_{0\leq n\leq
M}|\hat{X}_{n}^{i,N}|^{2p}+Ch^{p})\Big{]}\leq
C\Big{(}1+\mathbb{E}\Big{[}\,|\hat{X}_{0}|^{2p}\Big{]}\Big{)}<\infty.$

∎

The last result in this block concerns incremental (in time) moment bounds of
$\hat{X}^{i,N}$.

######  Proposition 4.7.

Let the assumptions of Proposition 4.6 hold. Then, there exists a constant
$C>0$ such that for any $p\geq 2$ (with $m$ sufficiently large depending on
$p$)

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq T}|\hat{X}_{t}^{i,N}-\hat{X}_{\kappa(t)}^{i,N}|^{p}\Big{]}\leq
Ch^{\frac{p}{2}}.$ (4.12)

###### Proof.

From (4.8), set for all $t\in[0,T],\leavevmode\nobreak\ t_{n}+s=t$

$\displaystyle\hat{X}_{t}^{i,N}=\hat{X}_{n}^{i,N}+v_{h}(t_{n},\hat{X}_{n}^{i,N})s+b_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})s+\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})(W_{t}^{i}-W_{t_{n}}^{i})$

Thus, there is a constant $C_{p}$ independent of $h,N$ such that:

$\displaystyle|\hat{X}_{t}^{i,N}-\hat{X}_{n}^{i,N}|^{p}\leq$ $\displaystyle
C_{p}\Big{[}|v_{h}(t_{n},\hat{X}_{n}^{i,N})|^{p}h^{p}+|b_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})|^{p}h^{p}+|\sigma_{h}(t_{n},\hat{X}_{n}^{i,N},\widetilde{\mu}_{n}^{N})|^{p}|(W_{t}^{i}-W_{t_{n}}^{i})|^{p}\Big{]}.$

Take take sup over time and expectation on both sides. Using Assumption 2.1
one deals with the last term with the BDG inequality (using conditional
expectation $\mathbb{E}[\cdot|\mathcal{F}_{t_{n}}]$). From Lemma 4.1 and
Propositions 4.5 and 4.6, there exists $C$ independent of $h,N,M$

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq T}|\hat{X}_{t}^{i,N}-\hat{X}_{\kappa(t)}^{i,N}|^{p}\Big{]}\leq
C(h^{p}+h^{\frac{p}{2}})\Big{(}1+\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq n\leq
M}|\hat{X}_{n}^{i,N}|^{(q+1)p}\Big{]}\Big{)}\leq Ch^{\frac{p}{2}}.$

where $q$ follows from the polynomial growth property of $v$. ∎

#### Local errors

After having discussed moment bounds, we now discuss the local error.

######  Proposition 4.8.

Choose $h$ as in Theorem 2.6. Take the functions $v,b,\sigma$ and the
corresponding functions $v_{h},b_{h},\sigma_{h}$ as defined as Lemma 4.1.

Then, there are positive constants $C_{1},C_{2},C_{3}$ and
$q^{\prime}\in\mathbb{N}$, for all $t\in[0,T]$,
$z_{i}\in\mathbb{R}^{d},\leavevmode\nobreak\ i\in\llbracket 1,N\rrbracket$ and

$\displaystyle|v_{h}(t,z_{i})-v(t,z_{i})|^{2}$ $\displaystyle\leq
C_{1}\Big{(}1+|z_{i}|^{q^{\prime}}\Big{)}^{2}h^{2},$ (4.13)
$\displaystyle|b_{h}(t,z_{i},\mu_{z}^{h})-b(t,z_{i},\mu_{z})|^{2}$
$\displaystyle\leq
C_{2}\Big{(}1+|z_{i}|^{q^{\prime}}+\frac{1}{N}\sum_{j=1}^{N}|z_{j}|^{q^{\prime}}\Big{)}^{2}h^{2},$
(4.14)
$\displaystyle|\sigma_{h}(t,z_{i},\mu_{z}^{h})-\sigma(t,z_{i},\mu_{z})|^{2}$
$\displaystyle\leq
C_{3}\Big{(}1+|z_{i}|^{q^{\prime}}+\frac{1}{N}\sum_{j=1}^{N}|z_{j}|^{q^{\prime}}\Big{)}^{2}h^{2}.$
(4.15)

where $\mu_{z}$ and $\mu_{z}^{h}$ the two empirical measures associated with
$\\{z_{i}\\}_{i},\\{F_{h}(t,Z_{i})\\}_{i}$ respectively, namely,

$\displaystyle\mu_{z}(dz)=\frac{1}{N}\
\sum_{j=1}^{N}\delta_{z_{j}}(dz),\qquad\mu_{z}^{h}(dz)=\frac{1}{N}\
\sum_{j=1}^{N}\delta_{F_{h}(t,z_{j})}(dz).$

###### Proof.

Recall the estimates given in Lemma 4.1. Using
$F_{h}(t,z_{i})=z_{i}+hv_{h}(t,z_{i})$ and (4.3)

$\displaystyle\big{|}v_{h}(t,z_{i})-v(t,z_{i})\big{|}^{2}$
$\displaystyle=\big{|}v\big{(}t,z_{i}+hv_{h}(t,z_{i})\big{)}-v(t,z_{i})\big{|}^{2}$
$\displaystyle\leq
C\Big{(}1+|z_{i}+hv_{h}(t,z_{i})|^{q}+|z_{i}|^{q}\Big{)}^{2}h^{2}|v_{h}(t,z_{i})|^{2}$
$\displaystyle\leq
C\Big{(}1+|z_{i}|^{2q}+|z_{i}|^{2q(q+1)}h^{2q}\Big{)}\frac{h^{2}}{(1-hL_{v})^{2}}\Big{(}1+|z_{i}|^{2q+2}\Big{)}$
$\displaystyle\leq
C\Big{(}1+|z_{i}|^{4q+2}+|z_{i}|^{(2q+2)(q+1)}h^{2q}\Big{)}h^{2}$

As in Lemma 4.1 we show only the result for $b_{h}$ as the computation is the
same for $\sigma_{h}$ (and overall very close to that for $v_{h}$). Using the
definition of $b_{h}$, the Lipschitz property of $b$ and the definition of
$\mu_{z}^{h},\mu_{z}$ we have

$\displaystyle\big{|}b_{h}(t,z_{i},\mu_{z}^{h})-b(t,z_{i},\mu_{z})\big{|}^{2}$
$\displaystyle\leq
L_{b}\Big{(}h^{2}|v_{h}(t,z_{i})|^{2}+W^{(2)}(\mu_{z}^{h},\mu_{z})\Big{)}$
$\displaystyle\leq
L_{b}\Big{(}h^{2}|v_{h}(t,z_{i})|^{2}+\frac{1}{N}\sum_{j=1}^{N}h^{2}|v_{h}(t,z_{j})|^{2}\Big{)}\leq
C\Big{(}1+|z_{i}|^{2q+2}\Big{)}h^{2}$

Applying Inequality (4.3) and the growth in $v$ from assumption 2.1 yields the
claim.

∎

######  Proposition 4.9.

Choose $h$ as in Theorem 2.6. Take $h$ as in Theorem Let $i\in\llbracket
1,N\rrbracket$ and take $X^{i,N}$ to be the solution to the interacting
particle system (1.2) and let $\hat{X}^{i,N}$ be the continuous-time extension
of the SSM given by (4.8). We then have

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq T}|X_{t}^{i,N}-\hat{X}_{t}^{i,N}|^{2}\Big{]}=\mathcal{O}(h).$

###### Proof.

Let $i\in\llbracket 1,N\rrbracket$, $t\in[0,T]$ and $p\geq 2$. From (2.2) and
(4.10), both process $X^{i,N}$ and $\hat{X}^{i,N}$ have bounded $2p$-moments.
Define the auxiliary quantity $\Delta X^{i}:=X^{i,N}-\hat{X}^{i,N}$. Itô’s
formula applied to $|X_{t}^{i,N}-\hat{X}_{t}^{i,N}|^{2}=|\Delta
X^{i}_{t}|^{2}$ yields

$\displaystyle|\Delta X^{i}_{t}|^{2}=$ $\displaystyle
2\int_{0}^{t}\Big{\langle}v(s,X_{s}^{i,N})-v_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}ds$ (4.16)
$\displaystyle+2\int_{0}^{t}\Big{\langle}b(s,X_{s}^{i,N},\mu_{s}^{N})-b_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N}),\Delta
X^{i}_{s}\Big{\rangle}ds$ (4.17)
$\displaystyle+\int_{0}^{t}\Big{|}\sigma(s,X_{s}^{i,N},\mu_{s}^{N})-\sigma_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}ds$
(4.18) $\displaystyle+2\int_{0}^{t}\Big{\langle}\Delta
X^{i}_{s},\Big{(}\sigma(s,X_{s}^{i,N},\mu_{s}^{N})-\sigma_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{)}dW^{i}_{s}\Big{\rangle}.$
(4.19)

We analyse the components term by term. Namely, for (4.16)

$\displaystyle\Big{\langle}v(s,X_{s}^{i,N})-v_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}=$
$\displaystyle\Big{\langle}v(s,X_{s}^{i,N})-v(s,\hat{X}_{s}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}$
$\displaystyle+\Big{\langle}v(s,\hat{X}_{s}^{i,N})-v(\kappa(s),\hat{X}_{\kappa(s)}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}$
$\displaystyle+\Big{\langle}v(\kappa(s),\hat{X}_{\kappa(s)}^{i,N})-v_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}.$

We estimate the expectation of each term. From the one-sided Lipschitz
condition on $v$ we have

$\displaystyle\mathbb{E}\Big{[}\Big{\langle}v(s,X_{s}^{i,N})-v(s,\hat{X}_{s}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}\Big{]}$ $\displaystyle\leq
L_{v}\mathbb{E}\Big{[}\,|\Delta X^{i}_{s}|^{2}\Big{]}.$

Using the Cauchy-Schwarz inequality, the $1/2$-Hölder regularity in time and
the one-sided Lipschitz condition on $v$, yields

$\displaystyle\mathbb{E}\Big{[}\Big{\langle}v(s,\hat{X}_{s}^{i,N})-v(\kappa(s),\hat{X}_{\kappa(s)}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}\Big{]}$
$\displaystyle\leq\frac{1}{2}\mathbb{E}\Big{[}\Big{|}v(s,\hat{X}_{s}^{i,N})-v(\kappa(s),\hat{X}_{\kappa(s)}^{i,N})\Big{|}^{2}\Big{]}+\frac{1}{2}\mathbb{E}\Big{[}|\Delta
X^{i}_{s}|^{2}\Big{]}$ $\displaystyle\leq
Ch+C\mathbb{E}\Big{[}\Big{(}1+|\hat{X}_{s}^{i,N}|^{2q}+|\hat{X}_{\kappa(s)}^{i,N}|^{2q}\Big{)}|\hat{X}_{s}^{i,N}-\hat{X}_{\kappa(s)}^{i,N}|^{2}\Big{]}+\frac{1}{2}\mathbb{E}\Big{[}|\Delta
X^{i}_{s}|^{2}\Big{]}\leq Ch+C\mathbb{E}\Big{[}|\Delta X^{i}_{s}|^{2}\Big{]}.$

For the last inequality we used Hölder’s inequality on the product term in
combination with Proposition 4.7. We now make use of Proposition 4.8 and
arguments similar to those above to deal with the last term of the initial
inequality

$\displaystyle\mathbb{E}\Big{[}\Big{\langle}v(\kappa(s),\hat{X}_{\kappa(s)}^{i,N})-v_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N}),\Delta
X^{i}_{s}\Big{\rangle}\Big{]}\leq\frac{1}{2}\mathbb{E}\Big{[}\Big{|}v(\kappa(s),\hat{X}_{\kappa(s)}^{i,N})-v_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N})\Big{|}^{2}+|\Delta
X^{i}_{s}|^{2}\Big{]}$ $\displaystyle\leq
C\mathbb{E}\Bigg{[}1+|\hat{X}_{\kappa(s)}^{i,N}|^{q^{\prime}}\Bigg{]}h^{2}+\frac{1}{2}\mathbb{E}\Big{[}|\Delta
X^{i}_{s}|^{2}\Big{]}\leq Ch^{2}+\frac{1}{2}\mathbb{E}\Big{[}|\Delta
X^{i}_{s}|^{2}\Big{]}.$

We now proceed to estimate the $b$ components. Using Cauchy-Schwarz
inequality, for (4.17)

$\displaystyle\mathbb{E}\Big{[}\Big{\langle}b(s,X_{s}^{i,N},\mu_{s}^{N})-b_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N}),\Delta
X^{i}_{s}\Big{\rangle}\Big{]}$
$\displaystyle\leq\mathbb{E}\Big{[}\frac{1}{2}\Big{|}b(s,X_{s}^{i,N},\mu_{s}^{N})-b(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}+\frac{1}{2}|\Delta
X^{i}_{s}|^{2}\Big{]}$
$\displaystyle\leq\frac{1}{2}\mathbb{E}\Big{[}\leavevmode\nobreak\
\Big{|}b(s,X_{s}^{i,N},\mu_{s}^{N})-b(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}+\Big{|}b(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})-b_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}+|\Delta
X^{i}_{s}|^{2}\leavevmode\nobreak\ \Big{]}$

From the Lipschitz condition on $b$, by (4.4) and Proposition 4.7,

$\displaystyle\mathbb{E}\Big{[}\leavevmode\nobreak\
\Big{|}b(s,X_{s}^{i,N},\mu_{s}^{N})-b(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}\leavevmode\nobreak\
\Big{]}$ $\displaystyle\leq C\mathbb{E}\Big{[}\leavevmode\nobreak\
h+|X_{s}^{i,N}-\hat{X}_{s}^{i,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|X_{s}^{j,N}-\hat{X}_{s}^{j,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}\Big{|}F_{h}(s,\hat{X}_{s}^{j,N})-F_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{j,N})\Big{|}^{2}\leavevmode\nobreak\
\Big{]}$ $\displaystyle\leq C\mathbb{E}\Big{[}\leavevmode\nobreak\
h+|X_{s}^{i,N}-\hat{X}_{s}^{i,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|X_{s}^{j,N}-\hat{X}_{s}^{j,N}|^{2}\leavevmode\nobreak\
\Big{]}.$

From the Lipschitz condition on $b$

$\displaystyle\mathbb{E}\Big{[}\leavevmode\nobreak\
\Big{|}b(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})-b_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}\leavevmode\nobreak\
\Big{]}\leq C\mathbb{E}\Big{[}\leavevmode\nobreak\
\Big{(}1+|\hat{X}_{\kappa(s)}^{i,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|\hat{X}_{\kappa(s)}^{j,N}|^{2}\Big{)}h^{2}\leavevmode\nobreak\
\Big{]}\leq Ch^{2},$

and similarly we obtain
$\mathbb{E}\Big{[}\,\Big{|}\sigma(s,X_{s}^{i,N},\mu_{s}^{N})-\sigma_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}\,\Big{]}\leq
Ch^{2}$.

Consider the last term (4.19) , take expectation, using the BDG inequality and
above result, by Cauchy-Schwarz inequality:

$\displaystyle\mathbb{E}\Big{[}\sup_{0\leq u\leq
t}\int_{0}^{t}\Big{\langle}\Delta
X^{i}_{s},\Big{(}\sigma(s,X_{s}^{i,N},\mu_{s}^{N})-\sigma_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{)}dW^{i}_{s}\Big{\rangle}\leavevmode\nobreak\
\Big{]}$ $\displaystyle\leq\mathbb{E}\Big{[}\leavevmode\nobreak\
C\Big{(}\int_{0}^{T}|\Delta X^{i}_{s}|^{2}\leavevmode\nobreak\
\Big{|}\sigma(s,X_{s}^{i,N},\mu_{s}^{N})-\sigma_{h}(\kappa(s),\hat{X}_{\kappa(s)}^{i,N},\widetilde{\mu}_{\kappa(s)}^{N})\Big{|}^{2}ds\Big{)}^{\frac{1}{2}}\leavevmode\nobreak\
\Big{]}$ $\displaystyle\leq\frac{1}{2}\mathbb{E}\Big{[}\sup_{0\leq t\leq
T}|\Delta X^{i}_{t}|^{2}\leavevmode\nobreak\ \Big{]}+Ch^{2}.$

Gathering all inequalities together and taking supremum on $i$, we conclude
(where $h$ is the leading term)

$\displaystyle\sup_{i\in\llbracket 1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq
t\leq T}|X_{t}^{i,N}-\hat{X}_{t}^{i,N}|^{2}\leavevmode\nobreak\ \Big{]}$
$\displaystyle\leq\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\leavevmode\nobreak\
C\int_{0}^{T}\Big{(}h+|X_{s}^{i,N}-\hat{X}_{s}^{i,N}|^{2}+\frac{1}{N}\sum_{j=1}^{N}|X_{s}^{j,N}-\hat{X}_{s}^{j,N}|^{2}\Big{)}ds\leavevmode\nobreak\
\Big{]}$ $\displaystyle\leq Ch+C\int_{0}^{T}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}\sup_{0\leq u\leq
s}|X_{u}^{i,N}-\hat{X}_{u}^{i,N}|^{2}\leavevmode\nobreak\ \Big{]}ds.$

Grönwall’s lemma delivers the final result. ∎

We can now conclude with the proof of the main result.

###### Proof of Theorem 2.6.

In relation to the Points 1, 2 and 3 in the theorem’s statement: Point 1
follows from Proposition 4.6; Point 2 follows from Proposition 4.9; the last
point follows by a straightforward combination of Proposition 2.4 and
Proposition 4.9. ∎

###  4.2 Proof of the stability theorem, Theorem 2.8

###### Proof of Theorem 2.8.

Let $i\in\llbracket 1,N\rrbracket$ and $n\in\mathbb{N}$. From Lemma 4.1, we
have immediately that

$\displaystyle|Y_{n}^{i,\star,N}-G_{n}^{i,\star,N}|^{2}\leq\frac{|\hat{X}_{n}^{i,N}-\hat{Z}_{n}^{i,N}|^{2}}{1-2hL_{v}}.$

By definition of the SSM one also has

$\displaystyle(Y_{n}^{i,\star,N}-\hat{X}_{n}^{i,N})=hv\left(t_{n},Y_{n}^{i,\star,N}\right)\quad\textrm{and}\quad(G_{n}^{i,\star,N}-\hat{Z}_{n}^{i,N})=hv\left(t_{n},G_{n}^{i,\star,N}\right).$

Thus, by Assumption 2.1 and the Cauchy-Schwarz inequality, for some $n$,

$\displaystyle\mathbb{E}\Big{[}|\hat{X}_{n+1}^{i,N}-\hat{Z}_{n+1}^{i,N}|^{2}\Big{]}$
$\displaystyle\leq$ $\displaystyle\mathbb{E}\bigg{[}\leavevmode\nobreak\
|Y_{n}^{i,\star,N}-G_{n}^{i,\star,N}|^{2}+2\Big{\langle}Y_{n}^{i,\star,N}-G_{n}^{i,\star,N},b(t_{n},Y_{n}^{i,\star,N},\hat{\mu}^{Y,N}_{n})-b(t_{n},G_{n}^{i,\star,N},\hat{\mu}^{G,N}_{n})\Big{\rangle}h$
$\displaystyle+\Big{|}b(t_{n},Y_{n}^{i,\star,N},\hat{\mu}^{Y,N}_{n})-b(t_{n},G_{n}^{i,\star,N},\hat{\mu}^{G,N}_{n})\Big{|}^{2}h^{2}+\Big{|}\sigma(t_{n},Y_{n}^{i,\star,N},\hat{\mu}^{Y,N}_{n})-\sigma(t_{n},G_{n}^{i,\star,N},\hat{\mu}^{G,N}_{n})\Big{|}^{2}(\Delta
W_{n}^{i})^{2}\bigg{]}$ $\displaystyle\leq$
$\displaystyle\mathbb{E}\bigg{[}(1+h+L_{b}h+L_{b}h^{2}+L_{\sigma}h)|Y_{n}^{i,\star,N}-G_{n}^{i,\star,N}|^{2}+(L_{\hat{b}}h^{2}+L_{\hat{\sigma}}h+L_{\hat{b}}h)\Big{[}\frac{1}{N}\sum_{j=1}^{N}|Y_{n}^{j,\star,N}-G_{n}^{j,\star,N}|^{2}\Big{]}\bigg{]}.$

where we used the tower property to deal with the Brownian increment. Taking
supremum over $i$ and using (4.4) yields

$\displaystyle\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}_{n+1}^{i,N}-\hat{Z}_{n+1}^{i,N}|^{2}\Big{]}\leq(1+\beta
h)\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}_{n}^{i,N}-\hat{Z}_{n}^{i,N}|^{2}\Big{]}.$

where $\beta$ and $\alpha$ are exactly given by (2.6). A straightforward
induction arguments leads to

$\displaystyle\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}_{n}^{i,N}-\hat{Z}_{n}^{i,N}|^{2}\Big{]}\leq$
$\displaystyle(1+\beta h)^{n}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}^{i,N}_{0}-\hat{Z}^{i,N}_{0}|^{2}\Big{]}.$

Recall (2.6) for $\beta$ and $\alpha$. If $L_{v}\geq\alpha$ then $1+\beta h>1$
and hence $\lim_{n\rightarrow\infty}\big{(}1+\beta)^{n}\neq 0$ implying that
the SSM is not Mean-square contractive. On the other hand, since $(1+\beta h)$
is always positive then when $\beta<0,(1+\beta h)\in(0,1)$ $\Leftrightarrow
L_{v}<\alpha<-\frac{1}{2}$ with the restriction
$h<-2\frac{L_{v}-\alpha}{L_{b}+L_{\hat{b}}}$ and consequently the SSM is Mean-
square contractive

$\displaystyle\lim_{n\rightarrow\infty}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}_{n}^{i,N}-\hat{Z}_{n}^{i,N}|^{2}\Big{]}\leq$
$\displaystyle\lim_{n\rightarrow\infty}\big{(}1+\beta
h)^{n}\sup_{i\in\llbracket
1,N\rrbracket}\mathbb{E}\Big{[}|\hat{X}^{i,N}_{0}-\hat{Z}^{i,N}_{0}|^{2}\Big{]}=0.$

∎

##  Appendix A Short description of the Taming and time-adaptive method

We provide a brief review of the Taming method [18] and the time adaptive
method [45] for super-linear growth MV-SDEs in the context of the notation set
in Section 2 for $N,M\in\mathbb{N}$:
$[0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\ni(t,x,\mu)\mapsto\widehat{b}(t,x,\mu)=v(t,x)+b(t,x,\mu)$,
$Mh=T,\leavevmode\nobreak\ i\in\llbracket 1,N\rrbracket$,
$\leavevmode\nobreak\ n\in\llbracket 0,M-1\rrbracket$. Each method
approximates (1.1) through the interacting particle system (1.2) as described
next.

###  A.1 Explicit Euler particle system for the super-linear case: taming

The taming method [18] approximates (1.2) as follows (see also [45, Section
4]):

$\displaystyle\bar{X}_{n+1}^{i,N,M}=\bar{X}_{n}^{i,N,M}$
$\displaystyle+\dfrac{\widehat{b}\Big{(}t_{n},\bar{X}_{n}^{i,N,M},\bar{\mu}^{X,N}_{n}\Big{)}}{1+M^{-\alpha}\Big{|}\widehat{b}\Big{(}t_{n},\bar{X}_{n}^{i,N,M},\bar{\mu}^{X,N}_{n}\Big{)}\Big{|}}h+\sigma\Big{(}t_{n},\bar{X}_{n}^{i,N,M},\bar{\mu}^{X,N}_{n}\Big{)}\Delta
W_{n}^{i},$ (A.1)

where
$\bar{\mu}^{X,N}_{n}(\mathrm{d}x)=\frac{1}{N}\sum_{j=1}^{N}\delta_{\bar{X}_{n}^{j,N,M}}(\mathrm{d}x)$,$\leavevmode\nobreak\
\Delta W_{n}^{i}=W_{t_{n+1}}^{i}-W_{t_{n}}^{i}$ with
$\bar{X}_{0}^{i,N,M}=X_{0}^{i}$. The parameter $\alpha\in(0,1]$ where setting
$\alpha=1/2$ delivers a rMSE convergence rate of order $1/2$ while setting
$\alpha=1$ delivers a rMSE convergence rate of order $1$ (for a constant
diffusion $\sigma$).

###  A.2 Time adaptive schemes

The adaptive method from [45] approximates (1.2) as follows for
$t_{n}\in[k_{n}h,(k_{n}+1)h),k_{n}\in\mathbb{N}$ and

$\displaystyle\bar{X}_{t_{n+1}}^{i,N}=\bar{X}_{t_{n}}^{i,N}$
$\displaystyle+\widehat{b}\Big{(}t_{n},\bar{X}_{t_{n}}^{i,N},\bar{\mu}^{X,N}_{k_{n}h}\Big{)}h_{n}^{i}+\sigma\Big{(}t_{n},\bar{X}_{t_{n}}^{i,N},\bar{\mu}^{X,N}_{k_{n}h}\Big{)}\Delta
W_{t_{n}}^{i},$ (A.2)

where
$\bar{\mu}^{X,N}_{k_{n}h}(\mathrm{d}x)=\frac{1}{N}\sum_{j=1}^{N}\delta_{\bar{X}_{k_{n}h}^{j,N,M}}(\mathrm{d}x)$,
$t_{n+1}=t_{n}+h_{n}^{i}$, $\leavevmode\nobreak\ \Delta
W_{t_{n}}^{i}=W_{t_{n+1}}^{i}-W_{t_{n}}^{i}$ with
$\bar{X}_{0}^{i,N,M}=X_{0}^{i}$ and for a map
$\mathbf{h}^{\delta}(x):\mathbb{R}^{d}\rightarrow[0,h]$

$h_{n}^{i}=\min\\{\mathbf{h}^{\delta}(\bar{X}_{t_{n}}^{i,N}),(k_{n}+1)h-t_{n}\\}.$

The function $\mathbf{h}^{\delta}$ is specified at each example and is to be
understood as similar technique to the taming method. There the drift
$\widehat{b}$ is modified to control the growth across the application of the
scheme, for the adaptive scheme, one modifies instead the time-step $h$ in a
dynamic fashion to control the growth of $\widehat{b}$, see [23] & [45].

## References

  * [1] D. Adams, G. dos Reis, R. Ravaille, W. Salkeld, and J. Tugaut.  Large Deviations and Exit-times for reflected McKean-Vlasov equations with self-stabilizing terms and superlinear drifts.  arXiv preprint arXiv:2005.10057, 2020.
  * [2] A. Agarwal and S. Pagliarani.  A Fourier-based Picard-iteration approach for a class of McKean-Vlasov SDEs with Lévy jumps.  Stochastics, 93(4):592–624, 2021.
  * [3] J. Baladron, D. Fasoli, O. Faugeras, and J. Touboul.  Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons.  J. Math. Neurosci., 2:Art. 10, 50, 2012.
  * [4] J. Bao, C. Reisinger, P. Ren, and W. Stockinger.  Milstein schemes for delay McKean equations and interacting particle systems.  arXiv preprint arXiv:2005.01165, 2020.
  * [5] J. Bao, C. Reisinger, P. Ren, and W. Stockinger.  First-order convergence of Milstein schemes for McKean-Vlasov equations and interacting particle systems.  Proceedings of the Royal Society A, 477(2245):20200258, 2021.
  * [6] M. Bauer, T. Meyer-Brandis, and F. Proske.  Strong solutions of mean-field stochastic differential equations with irregular drift.  Electron. J. Probab., 23:Paper No. 132, 35, 2018.
  * [7] D. Belomestny and J. Schoenmakers.  Projected particle methods for solving McKean-Vlasov stochastic differential equations.  SIAM Journal on Numerical Analysis, 56(6):3169–3195, 2018.
  * [8] O. Bencheikh and B. Jourdain.  Bias behaviour and antithetic sampling in mean-field particle approximations of SDEs nonlinear in the sense of McKean.  ESAIM: Proceedings and Surveys, 65:219–235, 2019.
  * [9] F. Bernal, G. dos Reis, and G. Smith.  Hybrid PDE solver for data-driven problems and modern branching.  European J. Appl. Math., 28(6):949–972, 2017.
  * [10] S. Biswas, C. Kumar, G. d. Reis, C. Reisinger, et al.  Well-posedness and tamed Euler schemes for McKean-Vlasov equations driven by Lévy noise.  arXiv preprint arXiv:2010.08585, 2020.
  * [11] F. Bolley, J. A. Cañizo, and J. A. Carrillo.  Stochastic mean-field limit: non-Lipschitz forces and swarming.  Math Models Methods Appl Sci, 21, 2011.
  * [12] M. Bossy, O. Faugeras, and D. Talay.  Clarification and complement to “Mean-Field description and propagation of chaos in networks of Hodgkin–Huxley and FitzHugh–Nagumo Neurons”.  Journal of mathematical neuroscience, 5(1):1–23, 2015.
  * [13] M. Bossy and D. Talay.  A stochastic particle method for the McKean-Vlasov and the Burgers equation.  Math. Comp., 66(217):157–192, 1997.
  * [14] E. Buckwar, A. Samson, M. Tamborrino, and I. Tubikanec.  Splitting methods for SDEs with locally Lipschitz drift. an illustration on the Fitzhugh-Nagumo model.  arXiv:2101.01027, 2021.
  * [15] R. Carmona.  Lectures on BSDEs, stochastic control, and stochastic differential games with financial applications, volume 1 of Financial Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2016.
  * [16] F. Delarue, D. Lacker, and K. Ramanan.  From the master equation to mean field game limit theory: a central limit theorem.  Electron. J. Probab., 24:Paper No. 51, 54, 2019.
  * [17] X. Ding and H. Qiao.  Euler–Maruyama Approximations for Stochastic McKean-Vlasov Equations with Non-Lipschitz Coefficients.  Journal of Theoretical Probability, pages 1–18, 2020.
  * [18] G. dos Reis, S. Engelhardt, and G. Smith.  Simulation of McKean–Vlasov SDEs with super-linear growth.  IMA Journal of Numerical Analysis, 01 2021.  draa099.
  * [19] G. dos Reis, W. Salkeld, and J. Tugaut.  Freidlin-Wentzell LDP in path space for McKean-Vlasov equations and the functional iterated logarithm law.  Ann. Appl. Probab., 29(3):1487–1540, 2019.
  * [20] G. dos Reis, G. Smith, and P. Tankov.  Importance sampling for McKean-Vlasov SDEs.  arXiv preprint arXiv:1803.09320, 2018.
  * [21] W. Dreyer, M. Gaberšček, C. Guhlke, R. Huth, and J. Jamnik.  Phase transition in a rechargeable lithium battery.  European J. Appl. Math., 22(3):267–290, 2011.
  * [22] R. Erban, J. Haškovec, and Y. Sun.  A Cucker–Smale Model with Noise and Delay.  SIAM Journal on Applied Mathematics, 76(4):1535–1557, 2016.
  * [23] W. Fang and M. B. Giles.  Adaptive Euler–Maruyama method for SDEs with nonglobally Lipschitz drift.  Annals of Applied Probability, 30(2):526–560, 2020.
  * [24] J. Garnier, G. Papanicolaou, and T.-W. Yang.  Mean field model for collective motion bistability.  Discrete Contin. Dyn. Syst. Ser. B, 24(2):851–879, 2019.
  * [25] S. N. Gomes, G. A. Pavliotis, and U. Vaes.  Mean field limits for interacting diffusions with colored noise: phase transitions and spectral numerical methods.  Multiscale Model. Simul., 18(3):1343–1370, 2020.
  * [26] C. Guhlke, P. Gajewski, M. Maurelli, P. K. Friz, and W. Dreyer.  Stochastic many-particle model for LFP electrodes.  Contin. Mech. Thermodyn., 30(3):593–628, 2018.
  * [27] S. Herrmann, P. Imkeller, I. Pavlyukevich, and D. Peithmann.  Stochastic resonance, volume 194 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2014.  A mathematical approach in the small noise limit.
  * [28] D. J. Higham, X. Mao, and A. M. Stuart.  Strong convergence of Euler-type methods for nonlinear stochastic differential equations.  SIAM Journal on Numerical Analysis, 40(3):1041–1063, 2002.
  * [29] X. Huang, P. Ren, and F.-Y. Wang.  Distribution dependent stochastic differential equations.  arXiv preprint arXiv:2012.13656, 2020.
  * [30] M. Hutzenthaler, A. Jentzen, and P. E. Kloeden.  Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with non-globally Lipschitz continuous coefficients.  Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 467(2130):1563–1576, 2011.
  * [31] M. Hutzenthaler, T. Kruse, and T. A. Nguyen.  Multilevel Picard approximations for McKean-Vlasov stochastic differential equations.  arXiv preprint arXiv:2103.10870, 2021.
  * [32] M. Kac.  Foundations of kinetic theory.  In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1954–1955, vol. III, pages 171–197. University of California Press, Berkeley and Los Angeles, 1956.
  * [33] C. Kumar et al.  On Explicit Milstein-type scheme for McKean-Vlasov Stochastic Differential Equations with Super-linear Drift Coefficient.  arXiv preprint arXiv:2004.01266, 2020.
  * [34] C. Kumar, C. Reisinger, W. Stockinger, et al.  Well-posedness and tamed schemes for McKean-Vlasov Equations with Common Noise.  arXiv preprint arXiv:2006.00463, 2020.
  * [35] D. Lacker.  On a strong form of propagation of chaos for McKean-Vlasov equations.  Electronic Communications in Probability, 23, 2018.
  * [36] D. Lacker.  Hierarchies, entropy, and quantitative propagation of chaos for mean field diffusions.  arXiv e-prints, pages arXiv–2105, 2021.
  * [37] G. Leobacher, C. Reisinger, and W. Stockinger.  Well-posedness and numerical schemes for McKean-Vlasov equations and interacting particle systems with discontinuous drift.  arXiv preprint arXiv:2006.14892, 2020.
  * [38] A. Lionnet, G. dos Reis, and L. Szpruch.  Time discretization of FBSDE with polynomial growth drivers and reaction-diffusion PDEs.  Ann. Appl. Probab., 25(5):2563–2625, 2015.
  * [39] F. Malrieu.  Convergence to equilibrium for granular media equations and their Euler schemes.  Ann. Appl. Probab., 13(2):540–560, 2003.
  * [40] F. Malrieu and D. Talay.  Concentration inequalities for Euler schemes.  In Monte Carlo and quasi-Monte Carlo methods 2004, pages 355–371. Springer, 2006.
  * [41] V. Malyshkin, editor.  Parallel Computing Technologies 15th International Conference, PaCT 2019, Almaty, Kazakhstan, August 19–23, 2019, Proceedings.  Theoretical Computer Science and General Issues ; 11657. Springer International Publishing : Imprint: Springer, Cham, 1st ed. 2019.. edition, 2019\.
  * [42] H. P. McKean, Jr.  A class of Markov processes associated with nonlinear parabolic equations.  Proc. Nat. Acad. Sci. U.S.A., 56:1907–1911, 1966.
  * [43] S. Méléard.  Asymptotic behaviour of some interacting particle systems; McKean-Vlasov and Boltzmann models.  In Probabilistic models for nonlinear partial differential equations, pages 42–95. Springer, 1996.
  * [44] Y. S. Mishura and A. Y. Veretennikov.  Existence and uniqueness theorems for solutions of McKean–Vlasov stochastic equations.  arXiv preprint arXiv:1603.02212, 2020.
  * [45] C. Reisinger and W. Stockinger.  An adaptive Euler-maruyama scheme for McKean SDEs with super-linear growth and application to the mean-field Fitzhugh-Nagumo model.  arXiv preprint arXiv:2005.06034, 2020.
  * [46] M. Röckner and X. Zhang.  Well-posedness of distribution dependent SDEs with singular drifts.  arXiv preprint arXiv:1809.02216, 2018.
  * [47] A.-S. Sznitman.  Topics in propagation of chaos.  Ecole d’Eté de Probabilités de Saint-Flour XIX — 1989, pages 165–251, 1991.
  * [48] E. Zeidler.  Nonlinear functional analysis and its applications. II/A.  Springer-Verlag, New York, 1990.  Linear monotone operators, Translated from the German by the author and Leo F. Boron.
  * [49] X. Zong, F. Wu, and C. Huang.  Convergence and stability of the semi-tamed Euler scheme for stochastic differential equations with non-Lipschitz continuous coefficients.  Appl. Math. Comput., 228:240–250, 2014.