# An explicit Milstein-type scheme for interacting particle systems and
McKean–Vlasov SDEs with common noise and non-differentiable drift coefficients

Sani<EMAIL_ADDRESS>[ Chaman
<EMAIL_ADDRESS>[ Neelima
<EMAIL_ADDRESS>[ Gonçalo dos
<EMAIL_ADDRESS>[ Christoph
<EMAIL_ADDRESS>[ India Institute of
Technology Roorkee, India presep=, ]e1  India Institute of Technology Roorkee,
Indiapresep=, ]e2  Delhi University, Indiapresep=, ]e3  Edinburgh University,
United Kingdom presep=, ]e4  Oxford University, United Kingdom presep=, ]e5

###### Abstract

We propose an explicit drift-randomised Milstein scheme for both McKean–Vlasov
stochastic differential equations and associated high dimensional interacting
particle systems with common noise. By using a drift randomisation step in
space and measure, we establish the scheme’s strong convergence rate of $1$
under reduced regularity assumptions on the drift coefficient: no classical
(Euclidean) derivatives in space or measure derivatives (e.g., Lions/Fréchet)
are required. The main result is established by enriching the concepts of
bistability and consistency of numerical schemes used previously for standard
SDE. We introduce certain Spijker-type norms (and associated Banach spaces) to
deal with the interaction of particles present in the stochastic systems being
analysed. A discussion of the scheme’s complexity is provided.

65C30, 60H35,

65C05, 65C35,

stochastic interacting particle systems,

McKean–Vlasov equations,

common noise,

Milstein scheme,

non-differentiable drift,

drift randomisation,

bistability,

###### keywords:

[class=MSC]

###### keywords:

, and

##  1 Introduction

For a given $T>0$, consider the following stochastic differential equation
(SDE) of McKean–Vlasov type and with common noise,

$\displaystyle X_{t}=X_{0}+$
$\displaystyle\int_{0}^{t}b(s,X_{s},\mathcal{L}^{1}(X_{s}))ds+\sum_{\ell=1}^{m_{1}}\int_{0}^{t}\sigma_{1}^{\ell}(s,X_{s},\mathcal{L}^{1}(X_{s}))dW_{s}^{\ell}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{0}^{t}\sigma_{0}^{\ell}(s,X_{s},\mathcal{L}^{1}(X_{s}))dW_{s}^{0,\ell},$
(1)

almost surely for all $t\in[0,T]$, where $W:=\\{W_{t}\\}_{t\in[0,T]}$ and
$W^{0}:=\\{W_{t}^{0}\\}_{t\in[0,T]}$ are, respectively, $m_{1}$ and $m_{0}$
dimensional independent Wiener processes and
$\\{\mathcal{L}^{1}(X_{t})\\}_{t\in[0,T]}$ denotes the stochastic flow of
conditional marginal laws of $X:=\\{X_{t}\\}_{t\in[0,T]}$ given $W^{0}$. The
initial value $X_{0}$ is an $\mathscr{F}^{0}$\- measurable random variable,
independent of $W$ and $W^{0}$. The McKean–Vlasov SDE (1) can be viewed as an
infinite-dimensional system of particles with $W$ representing the randomness
inherent in the individual particle and $W^{0}$ the randomness common to all
the particles. When $\sigma_{0}\equiv 0$, the particles are governed by only
one source of randomness, $W$, and the stochastic flow
$\\{\mathcal{L}^{1}(X_{t})\\}_{t\in[0,T]}$ becomes a deterministic one. Notice
that McKean–Vlasov SDEs are different from standard SDEs due to the dependence
of the coefficients on the (conditional) marginal law $\mathcal{L}^{1}(X_{t})$
of $X_{t}$ given $W^{0}$, which brings additional difficulties.

Due to their wide applications in areas such as Finance, mathematical
neuroscience and biology, machine learning and physics — animal swarming, cell
movement induced by chemotaxis, opinion dynamics, particle movement in porous
media and electrical battery modelling, self-assembly of particles and
dynamical density functional theory (see for example [31, 16, 8, 25, 2, 27,
33, 9, 13, 14, 28, 32])— McKean–Vlasov equations and associated interacting
particle systems, with or without common noise, addressed via stochastic
systems or associated Fokker Plank equations ([38, 22, 21]) have gained
immense popularity.

As in the case of SDEs, explicit solutions of McKean–Vlasov SDEs are typically
not available, which necessitates the development of numerical schemes to
approximate them. The numerical approximation of McKean–Vlasov SDEs can be
carried out in two steps, as explained below.

  * •

As a first step one builds the so-called interacting particle system,
$\\{X^{i,N}\\}_{i\in\\{1,\ldots,N\\}}$, where one replaces the (conditional)
marginal law appearing in the coefficients of (1) by the empirical law
obtained from the particles. Concretely, taking $N$ i.i.d. copies
$\\{W^{i}\\}_{i\in\\{1,\ldots,N\\}}$ of $W$ and
$\\{X^{i}_{0}\\}_{i\in\\{1,\ldots,N\\}}$ of $X_{0}$,

one defines the interacting particle system associated with the above
McKean–Vlasov SDE by

$\displaystyle X_{t}^{i,N}=X_{0}^{i}$
$\displaystyle+\int_{0}^{t}b(s,X_{s}^{i,N},\mu_{s}^{X,N})\,ds+\sum_{\ell=1}^{m_{1}}\int_{0}^{t}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})\,dW_{s}^{i,\ell}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{0}^{t}\sigma_{0}^{\ell}(s,X^{i,N}_{s},\mu_{s}^{X,N})\,dW_{s}^{0,\ell},$
(2)

almost surely for any $t\in[0,T]$ and $i\in\\{1,\ldots,N\\}$, where

$\displaystyle\mu_{t}^{X,N}(\cdot):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i,N}}(\cdot)$

is the empirical measure of $N$ particles. Subsequently, one needs to show,
roughly put, that $\textrm{Law}(X^{i,N}_{\cdot})$ for some $i$ (fixed)
convergences to $\textrm{Law}(X_{\cdot})$ of (1) as $N\to\infty$; see the
seminal work by Sznitman [53] (and Proposition 2.5 below).

  * •

In the second step of approximation, the temporal discretization of the
interacting particle system is performed to obtain fully implementable
numerical schemes for the McKean–Vlasov SDEs such as Euler-type schemes and
Milstein-type schemes. The main difficultly is to show that all estimates are
independent of the number of particles $N$ (see Theorem 3.3 and Corollary 3.4
below). The direct application of results from classic SDE theory do not
deliver this independence.

_It is critical to note that the results of this manuscript double for either
the numerical approximation of McKean–Vlasov SDEs if one’s starting point is (
1), or for stand-alone systems of interacting particle SDE systems if one’s
starting point is (2)._

Results on the strong well-posedness and propagation of chaos for
McKean–Vlasov SDEs are being extensively researched and we cannot possibly do
justice to that growing body of literature. Nonetheless, we mention some of
the milestones and more recent results fitting thematically with our
manuscript: for PoC results, one starts from Sznitman’s seminal work [53] to
the monographs [13, 14] and the recent developments by [23, 39, 40] to mention
a few – overall, there are still gaps in the existing PoC results, especially
across dimensions. Concretely, in [20, Section 3.5] the PoC rates across the
dimension $d$ is estimated numerically for a diffusion $\sigma$ of polynomial
growth and a drift $b$ of polynomial growth (a setting outside the scope of
this work) and the rates estimated are better than the rates found in any
theoretical results at present. For a focus on well-posedness, one starts from
the McKean’s seminal work [43], to again [13, 14] and the recent developments
[6, 17, 37, 44, 45, 52, 20] (and references therein).

The second step of numerical approximation of McKean–Vlasov SDEs saw
relatively little development after the first Euler-type scheme was proposed
and analysed (in a weak sense) in [10], but experience rapid advances in a
string of recent papers, [3, 4, 5, 18, 19, 36, 37, 41, 24, 49]. In particular,
[24] proposed an Euler-type numerical scheme for the interacting particle
system associated with the McKean–Vlasov SDE, and its strong convergence is
investigated when the drift coefficient grows super-linearly in the state
variable. More precisely, the drift is assumed to be one-sided Lipschitz
continuous in the state variable and Lipschitz continuous in the measure
variable, while the diffusion coefficient is assumed to be Lipschitz
continuous in both the state and measure variables. The rate of strong
convergence of the scheme is shown to be equal to $1/2$. In [4, 36], a
Milstein-type scheme for the interacting particle system associated with the
McKean–Vlasov SDEs is proposed using the notion of Lions derivatives,
introduced by P.-L. Lions in his lectures at the Collège de France and
presented in [12], and its strong convergence is shown with rate $1$. The
drift coefficient is assumed to satisfy a one-sided Lipschitz condition in the
state variable and a polynomial Lipschitz condition, and the diffusion
coefficient to satisfy Lipschitz condition; both are assumed to be Lipschitz
continuous in the measure variable. Furthermore, only the diffusion
coefficient is required to be once differentiable (in both state and measure
variables). In [4], the authors additionally require second order
differentiability of the coefficients. In [37], an Euler-type scheme and a
Milstein-type scheme are developed for the interacting particle system
connected with McKean–Vlasov SDEs with common noise, where all the
coefficients are allowed to grow super-linearly in the state variable.

In this article, we develop a Milstein-type scheme for the interacting
particle system corresponding to the McKean–Vlasov SDE without assuming the
differentiability of the drift coefficient (in space or measure component),
which is therefore more relaxed than the corresponding results in [4, 36, 37].
It is out of the scope of this work to lift the differentiability conditions
on the diffusion coefficient and hence those assumptions match those already
existing in the most recent literature. The relaxation of the regularity
requirement of the drift coefficient is achieved by a certain randomisation
strategy that needs to be applied to _both state and measure components_ : the
technical developments necessary to deal with this difficulty are the second
contribution of this manuscript. In the case of SDEs (when the coefficients do
not depend on the law of the solution process), the technique of randomisation
has been studied in [35] (also [7, 34] and and more recently in [46, 48]) to
construct a Milstein-type scheme without assuming the first order
differentiability of the drift coefficient. As the coefficients in our
settings depend on the law of the solution process as well, we require a two-
fold randomisation – one with respect to the state variable and the other with
respect to the measure variable. For this, we use a uniform random variable to
generate a random point in each sub-interval of the time mesh and the Euler
scheme is used to obtain values of the particles’ states at these random
points, which are then used in the drift coefficient of the Milstein scheme,
both for the state and empirical measure. The precise details of this
randomisation can be found in Section 3. Critically, the technique developed
in [35] for the analysis cannot be used directly in our settings and a novel
approach is required. We observe the following.

  * •

The interacting particle system associated with the McKean–Vlasov SDE can be
treated as an $(\mathbb{R}^{d})^{N}$-dimensional SDE and thus the results of
[35] could be applied directly. However, all estimates would depend on $N$ and
hence “explode” as $N$ tends to infinity. This implies that to establish our
results a new tool must be developed in order to show the independence on $N$.

  * •

We propose a new notion of bistability and consistency of the numerical scheme
that is appropriate for the context of high-dimensional interacting particle
systems. Inspired by [35] we propose suitable stochastic Spijker norms capable
of dealing with the interaction component of the particles. Further details
can be found in Section 5.

  * •

A discussion on the practicalities of implementing our scheme is given in
Section 3.1.1 which includes a critical view on complexity and the consequence
of having common-noise.

  * •

In the simplest version possible of (1), three drift functions are well within
the scope of our work are (linear interaction, convolution-functionals and
linear interaction kernels)

$\displaystyle b(s,X_{s},\mathcal{L}^{1}(X_{s}))$
$\displaystyle=f_{1}(X_{s})+f_{2}\big{(}\tilde{\mathbb{E}}^{1}(g(X_{s})\big{)},$
$\displaystyle\widehat{b}(s,X_{s},\mathcal{L}^{1}(X_{s}))$
$\displaystyle=f_{1}(X_{s})+\int_{\mathbb{R^{d}}}g(X_{s}-y)\mathcal{L}^{1}(X_{s})(dy),$
$\displaystyle\tilde{b}(s,X_{s},\mathcal{L}^{1}(X_{s}))$
$\displaystyle=f_{1}(X_{s})+\int_{\mathbb{R^{d}}}\tilde{g}(X_{s},y)\mathcal{L}^{1}(X_{s})(dy),$

for any real-valued functions $f_{1},f_{2},g,\tilde{g}$ that satisfy a
standard Lipschitz condition in space (but are not differentiable), e.g.,
$g(x)=-|x|$ and more complex examples for $g,\tilde{g}$ can be found in [16,
13, 14, 30, 32] and $\tilde{\mathbb{E}}^{1}$ represents the conditional
expectation given the common noise $W^{0}$. The 2nd drift example, $\hat{b}$,
corresponds to the usual convolution operator fairly common in modelling with
McKean–Vlasov SDE and associated interacting particle systems (with or without
common noise) [15, 30, 32, 1]. Lastly, we point the reader to Example 3.15 in
[1] that intuitively highlights why $\hat{b}$ is Lipschitz in the Wasserstein
metric but is not Lions differentiable.

### Organization

The main framework of the McKean–Vlasov equation and the interacting particle
system including well-posedness and propagation of chaos is given in Section
2. The numerical scheme focusing on the approximation of the interacting
particle system is found in Section 3 as are the main convergence results. All
proofs are given in Section 5.

###  1.1 Notations

Both the Euclidean norm on $\mathbb{R}^{d}$ and the standard matrix norm on
$\mathbb{R}^{d\times m}$ are denoted by $|\cdot|$. The notation $\delta_{x}$
stands for the Dirac measure centred at $x\in\mathbb{R}^{d}$. We use the same
notation $a^{\ell}$ to denote the $\ell$-th column of a matrix
$a\in\mathbb{R}^{d\times m}$ and the $\ell$-th element of a vector
$a\in\mathbb{R}^{d}$. $\mathscr{B}(\chi)$ stands for the Borel
$\sigma$-algebra on a topological space $\chi$. Further,
$\mathscr{P}_{2}(\mathbb{R}^{d})$ denotes the space of all probability
measures on $(\mathbb{R}^{d},\mathscr{B}(\mathbb{R}^{d}))$ having finite
second moment and equipped with the $\mathscr{L}^{2}$-Wasserstein metric given
by

$\displaystyle\mathcal{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{]}^{1/2},$

for any $\mu$, $\nu\in\mathscr{P}_{2}(\mathbb{R}^{d})$, where $\Pi(\mu,\nu)$
denotes the set of couplings of $\mu$ and $\nu$. Clearly,
$\mathscr{P}_{2}(\mathbb{R}^{d})$ is a Polish space under this metric. We use
$\mathscr{L}^{p}(\Omega)$ to denote the Banach space of all
$\mathbb{R}^{d}$-valued random variables $Y$ defined on a probability space
$(\Omega,\mathscr{F},\mathbb{P})$ and satisfies
$\|Y\|_{\mathscr{L}^{p}(\Omega)}:=\big{[}\mathbb{E}|Y|^{p}\big{]}^{1/p}<\infty$,
where $\mathbb{E}$ stands for expectation with respect to $\mathbb{P}$.
Similarly, we use $\mathscr{L}^{p}([0,T]\times\Omega)$ to denote the Banach
space of processes $Y:[0,T]\times\Omega\mapsto\mathbb{R}^{d}$ having
$\|Y\|_{\mathscr{L}^{p}([0,T]\times\Omega)}=\big{[}\int_{0}^{T}\|Y(s)\|_{\mathscr{L}^{p}(\Omega)}^{p}ds\big{]}^{1/p}<\infty$.
Also, $\mathscr{C}^{\alpha}([0,T],\,\mathscr{L}^{p}(\Omega))$ stands for the
space of all $\alpha$-Hölder continuous functions
$Y:[0,T]\mapsto\mathscr{L}^{p}(\Omega)$ with the following norm,

$\displaystyle\|Y\|_{\mathscr{C}^{\alpha}([0,T],\,\mathscr{L}^{p}(\Omega))}=\sup_{t\in[0,T]}\|Y(t)\|_{\mathscr{L}^{p}(\Omega)}+\sup_{\underset{t\neq
t^{\prime}}{t,t^{\prime}\in[0,T]}}\frac{\|Y(t)-Y(t^{\prime})\|_{\mathscr{L}^{p}(\Omega)}}{|t-t^{\prime}|^{\alpha}}.$
(3)

For a function
$f:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\rightarrow\mathbb{R}$,
$\partial_{x}f:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\rightarrow\mathbb{R}^{d}$
is the derivative of $f$ with respect to the space variable and
$\partial_{\mu}f:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
is the Lions derivative with respect to the measure variable. $I_{A}$ stands
for the indicator function of a set $A$ and
$\bar{\mathbb{N}}=\mathbb{N}\cup\\{0\\}$. The constants that appear in the
paper vary from line to line, will depend on the problems data, for instance
$T$, $m_{0},m_{1}$, etc., but critically are independent of the number of
particle $N$ and the schemes timestep $h$ (given below).

##  2 McKean–Vlasov Stochastic Differential Equations and the interacting
particle systems

Let $T>0$ be a fixed constant. Consider probability spaces
$(\Omega^{1},\mathscr{F}^{1},\mathbb{P}^{1})$ and
$(\Omega^{0},\mathscr{F}^{0},\mathbb{P}^{0})$ equipped with filtrations
$\mathbb{F}^{1}:=\\{\mathscr{F}_{t}^{1}\\}_{t\in[0,T]}$ and
$\mathbb{F}^{0}:=\\{\mathscr{F}_{t}^{0}\\}_{t\in[0,T]}$, respectively. The
filtrations $\mathbb{F}^{1}$ and $\mathbb{F}^{0}$ satisfy the usual
conditions, i.e., they are complete and right continuous. Assume that
$W:=\\{W_{t}{\in\mathbb{R}^{m_{1}}}\\}_{t\in[0,T]}$ and
$W^{0}:=\\{W_{t}^{0}{\in\mathbb{R}^{m_{0}}}\\}_{t\in[0,T]}$ are independent
Brownian motions defined on
$(\Omega^{1},\mathscr{F}^{1},\mathbb{F}^{1},\mathbb{P}^{1})$ and
$(\Omega^{0},\mathscr{F}^{0},\mathbb{F}^{0},\mathbb{P}^{0})$, respectively. In
what follows, the interacting particles are governed by i.i.d. copies of $W$
and $W^{0}$ represents the noise common to all the particles. Let us define a
product probability space
$(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathbb{P}})$, where
$\tilde{\Omega}:=\Omega^{1}\times\Omega^{0}$,
$(\tilde{\mathscr{F}},\tilde{\mathbb{P}})$ is the completion of
$(\mathscr{F}^{1}\otimes\mathscr{F}^{0},\mathbb{P}^{1}\otimes\mathbb{P}^{0})$
and $\tilde{\mathbb{F}}:=\\{\tilde{\mathscr{F}}_{t}\\}_{t\in[0,T]}$ is the
completion and right-continuous augmentation of
$\\{\mathscr{F}_{t}^{1}\otimes\mathscr{F}_{t}^{0}\\}_{t\in[0,T]}$. The
expectation with respect to $\tilde{\mathbb{P}}$ is denoted by
$\tilde{\mathbb{E}}$.

Let
$b:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\mapsto\mathbb{R}^{d}$,
$\sigma_{1}:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\mapsto\mathbb{R}^{d\times
m_{1}}$ and
$\sigma_{0}:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\mapsto\mathbb{R}^{d\times
m_{0}}$ be
$\mathscr{B}([0,T])\otimes\mathscr{B}(\mathbb{R}^{d})\otimes\mathscr{B}(\mathscr{P}_{2}(\mathbb{R}^{d}))$-measurable
functions.

In this article, we consider the following $\mathbb{R}^{d}$-valued
McKean–Vlasov stochastic differential equations (SDEs) with common noise
defined on
$(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathbb{F}},\tilde{\mathbb{P}})$,

$\displaystyle X_{t}=X_{0}+$
$\displaystyle\int_{0}^{t}b(s,X_{s},\mathcal{L}^{1}(X_{s}))ds+\sum_{\ell=1}^{m_{1}}\int_{0}^{t}\sigma_{1}^{\ell}(s,X_{s},\mathcal{L}^{1}(X_{s}))dW_{s}^{\ell}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{0}^{t}\sigma_{0}^{\ell}(s,X_{s},\mathcal{L}^{1}(X_{s}))dW_{s}^{0,\ell},$
(4)

almost surely for all $t\in[0,T]$, where
$\\{\mathcal{L}^{1}(X_{t})\\}_{t\in[0,T]}$ is the flow of conditional laws of
$X_{t}$ given $W^{0}$ and $X_{0}$. A priori, it is not certain that the flow
of conditional marginals $\\{\mathcal{L}^{1}(X_{t})\\}_{t\in[0,T]}$ is an
$\mathbb{F}^{0}$-adapted continuous process. However, due to Lemma 2.5 in
[14], when the McKean–Vlasov SDE (2) has a unique $\tilde{\mathbb{F}}$-adapted
continuous solution with uniformly bounded second moment, then
$\\{\mathcal{L}^{1}(X_{t})\\}_{t\in[0,T]}$ is an $\mathbb{F}^{0}$-adapted
continuous process.

We make the following assumptions.

######  Assumption H– 1.

$X_{0}\in\mathscr{L}^{\bar{p}}(\tilde{\Omega})$ for some $\bar{p}\geq 2$.

######  Assumption H– 2.

There exists a constant $L>0$ such that

$\displaystyle|b(t,x,\mu)-b(t,x^{\prime},\mu^{\prime})|+\sum_{u=0}^{1}|\sigma_{u}(t,x,\mu)-\sigma_{u}(t,x^{\prime},\mu^{\prime})|\leq
L\\{|x-x^{\prime}|+\mathcal{W}_{2}(\mu,\mu^{\prime})\\},$

for all $t\in[0,T]$, $x,x^{\prime}\in\mathbb{R}^{d}$ and
$\mu,\mu^{\prime}\in\mathscr{P}_{2}(\mathbb{R}^{d})$.

######  Assumption H– 3.

There exists a constant $L>0$ such that

$\displaystyle|b(t,x,\mu)-b(t^{\prime},x,\mu)|$ $\displaystyle\leq
L\\{1+|x|+\mathcal{W}_{2}(\mu,\delta_{0})\\}|t-t^{\prime}|^{1/2},$
$\displaystyle\sum_{u=0}^{1}|\sigma_{u}(t,x,\mu)-\sigma_{u}(t^{\prime},x,\mu)|$
$\displaystyle\leq
L\\{1+|x|+\mathcal{W}_{2}(\mu,\delta_{0})\\}|t-t^{\prime}|,$

for all $t,t^{\prime}\in[0,T]$, $x\in\mathbb{R}^{d}$ and
$\mu\in\mathscr{P}_{2}(\mathbb{R}^{d})$.

######  Remark 2.1.

From Assumptions H–2 and H–3, for all $t\in[0,T]$, $x\in\mathbb{R}^{d}$ and
$\mu\in\mathscr{P}_{2}(\mathbb{R}^{d})$,

$\displaystyle|b(t,x,\mu)|+\sum_{u=0}^{1}|\sigma_{u}(t,x,\mu)|$
$\displaystyle\leq\bar{L}\\{1+|x|+\mathcal{W}_{2}(\mu,\delta_{0})\\},$

where
$\bar{L}=\max\big{\\{}L,LT,L\sqrt{T},b(0,0,\delta_{0}),\sigma_{0}(0,0,\delta_{0}),\sigma_{1}(0,0,\delta_{0})\big{\\}}$.

The proof of the following proposition can be found in [14, 37] and Appendix
A.1.

######  Proposition 2.2 (Well-posedness and Moment Bounds).

If Assumptions H–1 with $\bar{p}\geq 2$, H–2 and H–3 are satisfied, then the
McKean–Vlasov SDE (2) has a unique $\tilde{\mathbb{F}}$-adapted solution
$\\{X_{t}\\}_{t\in[0,T]}$ and

$\displaystyle\big{\|}\sup_{t\in[0,T]}|X_{t}|\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\leq
C_{1}(1+\|X_{0}\|_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}^{\bar{p}}),$

for
$C_{1}:=4^{\bar{p}-1}\max\big{\\{}1,3^{\bar{p}-1}\bar{L}^{\bar{p}}\big{(}T^{\bar{p}}+2\big{(}\frac{\bar{p}^{3}T}{2(\bar{p}-1)}\big{)}^{\bar{p}/2}\big{)}\big{\\}}\exp\big{(}12^{\bar{p}-1}\bar{L}^{\bar{p}}\big{(}2T^{%̄
̉\bar{p}}+4\big{(}\frac{\bar{p}^{3}T}{2(\bar{p}-1)}\big{)}^{\bar{p}/2}\big{)}\big{)}$.

To introduce the interacting particle system connected with the McKean–Vlasov
SDEs (2), let us consider $N\in\mathbb{N}$ i.i.d. copies of $W$ and $X_{0}$,
denoted by $W^{i}$ and $X_{0}^{i}$ for $i\in\\{1,\ldots,N\\}$, respectively.
Define the following system of equations,

$\displaystyle X^{i}_{t}=X^{i}_{0}$
$\displaystyle+\int_{0}^{t}b(s,X^{i}_{s},\mathcal{L}^{1}(X_{s}^{i}))ds+\sum_{\ell=1}^{m_{1}}\int_{0}^{t}\sigma_{1}^{\ell}(s,X^{i}_{s},\mathcal{L}^{1}(X_{s}^{i}))dW^{i,\ell}_{s}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{0}^{t}\sigma_{0}^{\ell}(s,X^{i}_{s},\mathcal{L}^{1}(X_{s}^{i}))dW^{0,\ell}_{s},$
(5)

almost surely for any $t\in[0,T]$ and $i\in\\{1,\ldots,N\\}$. Notice that due
to Proposition 2.11 in [14],
$\mathbb{P}^{0}\big{(}\mathcal{L}^{1}(X_{t}^{1})=\mathcal{L}^{1}(X_{t}^{i})\mbox{
for all }t\in[0,T]\big{)}=1$. We remark that the proof of Proposition 2.11 in
[14], which establishes the result in a more restrictive setting, uses (only)
the well-posedness of the system (2) and Theorem 1.33 from [14] (Yamada
Watanabe Theorem) and thus particles have the same law under our settings
also. The system (2) is popularly known as the conditional non-interacting
particle system. On approximating $\mathcal{L}^{1}(X_{t}^{1})$ by the
empirical measure of the states $\\{X_{t}^{i}\\}_{i\in\\{1,\ldots,N\\}}$ of
$N$ particles at time $t$, and denoted $\mu_{t}^{X,N}$, one obtains the
following system of interacting particles,

$\displaystyle X_{t}^{i,N}=X_{0}^{i}$
$\displaystyle+\int_{0}^{t}b(s,X_{s}^{i,N},\mu_{s}^{X,N})ds+\sum_{\ell=1}^{m_{1}}\int_{0}^{t}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})dW_{s}^{i,\ell}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{0}^{t}\sigma_{0}^{\ell}(s,X^{i,N}_{s},\mu_{s}^{X,N})dW_{s}^{0,\ell},$
(6)

almost surely for any $t\in[0,T]$ and $i\in\\{1,\ldots,N\\}$, where

$\displaystyle\mu_{t}^{X,N}(\cdot):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i,N}}(\cdot),$

is the empirical measure of $N$ particles.

######  Remark 2.3.

The system (2) can be understood as an $\mathbb{R}^{d\times N}$-dimensional
SDE and thus its well-posedness and moment stability up to order $\bar{p}$
follow from [29] under Assumptions H–1 (with $\bar{p}\geq 2$), H–2 and H–3. In
other words, the interacting particle system (2) connected with the
McKean–Vlasov SDE (2) has a unique strong solution
$\\{X_{t}^{i,N}\\}_{t\in[0,T]}$ adapted to the filtration
$\\{\tilde{\mathscr{F}}_{t}\\}_{t\in[0,T]}$ and one can show that

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}\sup_{t\in[0,T]}|X_{t}^{i,N}|\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\leq
C_{1}\Big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}^{\bar{p}}\Big{)},$

where $C_{1}$ is the same constant as in Proposition 2.2. Note that the RHS of
the moment estimate above is _critically_ independent of $N$ – this follows
from [37] but does not from [29].

The following lemma gives the time-regularity of the interacting particle
system (2) and its proof is given in Appendix A.2.

######  Lemma 2.4.

Let Assumptions H–1 with $\bar{p}\geq 2$, H–2 and H–3 hold. Then,

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{t}^{i,N}-X_{t^{\prime}}^{i,N}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}$
$\displaystyle\leq
C_{2}|t-t^{\prime}|^{\bar{p}/2}\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i,N}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\big{)},$

for all $t>t^{\prime}\in[0,T]$ and $N\in\mathbb{N}$, where
$C_{2}:=(1+2C_{1})9^{\bar{p}-1}L^{\bar{p}}\big{(}T^{\bar{p}/2}+2\big{(}\frac{\bar{p}(\bar{p}-1)}{2}\big{)}^{\bar{p}/2}\big{)}$
and the constant $C_{1}$ is defined in Proposition 2.2.

The convergence of the interacting particle system (2) to the non-interacting
particle system (2) is popularly known in the literature as the propagation of
chaos and is stated in the following proposition (see Theorem 2.12 in [14] for
details or Appendix A.3). For this, let us define the empirical measure of the
non-interacting particle system (2) as

$\displaystyle\mu_{t}^{X}(\cdot):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i}}(\cdot)$

almost surely for any $t\in[0,T]$ and $N\in\mathbb{N}$. Further, use Theorem
5.8 in [13] and Proposition 2.2 to obtain the following estimate,

$\displaystyle\big{\|}\mathcal{W}_{2}(\mathcal{L}^{1}(X_{t}^{1}),\mu_{t}^{X})\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}$
$\displaystyle\leq C_{3}\begin{cases}N^{-1/2},&\mbox{ if }d<4,\\\
N^{-1/2}\log_{2}N,&\mbox{ if }d=4,\\\ N^{-2/d},&\mbox{ if }d>4,\end{cases}$
(7)

for any $t\in[0,T]$ and $N\in\mathbb{N}$, where
$C_{3}:=C_{3}(d,\bar{p},\|X_{0}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}})$
is a positive constant.

######  Proposition 2.5 (Propagation of Chaos).

Let Assumptions H–1 with $\bar{p}>4$, H–2 and H–3 hold. Then,

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}\sup_{t\in[0,T]}$
$\displaystyle|X_{t}^{i}-X_{t}^{i,N}|\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}\leq
C_{4}\begin{cases}N^{-1/2},&\mbox{ if }d<4,\\\ N^{-1/2}\ \log_{2}N,&\mbox{ if
}d=4,\\\ N^{-2/d},&\mbox{ if }d>4,\end{cases}$

where
$C_{4}:=4(3T+24)TL^{2}e^{6(3T+24)TL^{2}}C_{3}(d,\bar{p},\|X_{0}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}})$,
where the constant $C_{3}$ appears in (7).

The proof of this result can be found in Appendix A.3.

All in all, the PoC rates here are in line with the PoC results found in the
numerics for McKean Vlasov SDE literature [3, 4, 5, 18, 19, 52, 36, 37, 41,
24, 49]. There are improved results obtaining rate of $1/N$ and $1/N^{2}$
(e.g., [23, 39]) that hold under stricter conditions that do not fit the scope
of the work presented in the manuscript. We have added text to this effect in
the main body of the paper.

##  3 Main Result

The drift-randomised Milstein scheme was originally proposed and analysed for
standard SDEs with non-differentiable drift coefficient in [35]. In the case
of McKean–Vlasov SDEs, for which the coefficients depend on the law of the
solution process as well, one needs additional randomisation of the drift
coefficient with respect to the measure component and a way to deal with the
implications from interacting particles. The associated difficulties are
tackled in this paper.

###  3.1 The Scheme

In order to propose the randomised Milstein scheme for the interacting
particle system (2) connected with the McKean–Vlasov SDEs (2), the map
$(x,\mu)\mapsto(\sigma_{1}(t,x,\mu),\sigma_{0}(t,x,\mu))$ is assumed to be
continuously differentiable for every $t\in[0,T]$. The notion of Lions
derivative (see Appendix A.4 for details) is used to differentiate these
functions with respect to the measure component. For $u\in\\{0,1\\}$, let us
define $d\times d$ matrices $\Lambda_{\sigma_{u}\sigma_{1}}$ and
$\Lambda_{\sigma_{u}\sigma_{0}}$ whose $\ell$-th column are given by,

$\displaystyle\Lambda_{\sigma_{u}\sigma_{1}}^{\ell}(t,s,x^{i},\mu):=$
$\displaystyle\sum_{\ell_{1}=1}^{m_{1}}\int_{s}^{t}\partial_{x}\sigma_{u}^{\ell}(s,x^{i},\mu)\sigma_{1}^{\ell_{1}}(s,x^{i},\mu)dW^{i,\ell_{1}}_{r}$
$\displaystyle+\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell_{1}=1}^{m_{1}}\int_{s}^{t}\partial_{\mu}\sigma_{u}^{\ell}(s,x^{i},\mu,x^{k})\sigma_{1}^{\ell_{1}}(s,x^{k},\mu)dW^{k,\ell_{1}}_{r},$
(8) $\displaystyle\Lambda^{\ell}_{\sigma_{u}\sigma_{0}}(t,s,x^{i},\mu):=$
$\displaystyle\sum_{\ell_{1}=1}^{m_{0}}\int_{s}^{t}\partial_{x}\sigma_{u}^{\ell}(s,x^{i},\mu)\sigma_{0}^{\ell_{1}}(s,x^{i},\mu)dW_{r}^{0,\ell_{1}}$
$\displaystyle+\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell_{1}=1}^{m_{0}}\int_{s}^{t}\partial_{\mu}\sigma_{u}^{\ell}(s,x^{i},\mu,x^{k})\sigma_{0}^{\ell_{1}}(s,x^{k},\mu)dW_{r}^{0,\ell_{1}},$
(9)

for all $s,t\in[0,T]$, $\ell\in\\{1,\ldots,m_{u}\\}$,
$x^{i}\in\mathbb{R}^{d}$, $i\in\\{1,\ldots,N\\}$,
$\mu\in\mathscr{P}_{2}(\mathbb{R}^{d})$. Also, define the $d$-dimensional
vectors $\tilde{\sigma}_{1}^{\ell}$ and $\tilde{\sigma}_{0}^{\ell_{1}}$ for
any $\ell\in\\{1,\ldots,m_{1}\\}$ and $\ell_{1}\in\\{1,\ldots,m_{0}\\}$ as

$\displaystyle\left.\begin{aligned}
\tilde{\sigma}_{1}^{\ell}(t,s,x^{i},\mu):=&\sigma_{1}^{\ell}(s,x^{i},\mu)+\Lambda^{\ell}_{\sigma_{1}\sigma_{1}}(t,s,x^{i},\mu)+\Lambda^{\ell}_{\sigma_{1}\sigma_{0}}(t,s,x^{i},\mu)\\\
\tilde{\sigma}_{0}^{\ell_{1}}(t,s,x^{i},\mu):=&\sigma^{\ell_{1}}_{0}(s,x^{i},\mu)+\Lambda^{\ell_{1}}_{\sigma_{0}\sigma_{1}}(t,s,x^{i},\mu)+\Lambda^{\ell_{1}}_{\sigma_{0}\sigma_{0}}(t,s,x^{i},\mu)\end{aligned}\right\\}$
(10)

for all $i\in\\{1,\ldots,n_{h}\\}$.

Now, consider a sequence $\eta:=\\{\eta_{j}\\}_{j\in\mathbb{N}}$ of i.i.d.
standard uniformly distributed random variables defined on a probability space
$(\Omega^{\eta},\mathscr{F}^{\eta},\mathbb{P}^{\eta})$, equipped with the
natural filtration
$\mathbb{F}^{\eta}:=\\{\mathscr{F}_{j}^{\eta}\\}_{j\in\mathbb{N}}$ of
$\\{\eta_{j}\\}_{j\in\mathbb{N}}$. Let $\mathbb{E}^{\eta}$ stand for the
expectation with respect to $\mathbb{P}^{\eta}$. The random variables
$\\{\eta_{j}\\}_{j\in\mathbb{N}}$ are assumed to be independent of $W$,
$W^{0}$, $W^{i}$ and $X_{0}^{i}$ for all $i\in\\{1,\ldots,N\\}$. Now consider
the general non-equidistant temporal grid $\varrho_{h}$ of $[0,T]$ with
$n_{h}$ subintervals,

$\displaystyle\varrho_{h}=\\{(t_{0},t_{1},\ldots,t_{n_{h}}):0=t_{0}<t_{1}<\cdots<t_{n_{h}}=T\\},$
(11)

with $h_{j}:=(t_{j}-t_{j-1})>0$ for $j\in\\{1,\ldots,n_{h}\\}$,
$n_{h}\in\mathbb{N}$ and
$h:=\max_{j\in\\{1,\ldots,n_{h}\\}}h_{j}\leq\min(1,T)$. Now, consider a new
probability space
$(\Omega,\mathscr{F},\mathbb{P})=(\tilde{\Omega}\times\Omega^{\eta},\tilde{\mathscr{F}}\otimes\mathscr{F}^{\eta},\tilde{\mathbb{P}}\otimes\mathbb{P}^{\eta})$
equipped with a filtration
$\mathbb{F}:=\\{\mathscr{F}_{j}^{h}\\}_{j\in\\{0,\ldots,n_{h}\\}}$; where
$\mathscr{F}_{j}^{h}=\tilde{\mathscr{F}}_{t_{j}}\otimes\mathscr{F}_{j}^{\eta}$.
We will denote the expectation with respect to $\mathbb{P}$ by $\mathbb{E}$.

The drift-randomised Milstein scheme for the interacting particle system (2)
of the McKean–Vlasov SDE (2) is given by

$\displaystyle X_{j,\eta}^{i,N,h}=$ $\displaystyle
X_{j-1}^{i,N,h}+\eta_{j}h_{j}b(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\int_{t_{j-1}}^{t_{j-1}+\eta_{j}h_{j}}dW_{s}^{i,\ell}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\sigma_{0}^{\ell}(t_{j-1},X^{i,N,h}_{j-1},\mu_{j-1}^{X,N,h})\int_{t_{j-1}}^{t_{j-1}+\eta_{j}h_{j}}dW_{s}^{0,\ell},$
(12) $\displaystyle X_{j}^{i,N,h}=$ $\displaystyle
X_{j-1}^{i,N,h}+h_{j}b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h})$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell}_{s}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW_{s}^{0,\ell},$
(13)

almost surely for all $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$
with the initial value $X_{0}^{i,N,h}=X_{0}^{i}$, where the empirical measures
$\mu_{j-1}^{X,N,h}$ and $\mu_{j,\eta}^{X,N,h}$ are defined as

$\displaystyle\mu_{j-1}^{X,N,h}(\cdot)$
$\displaystyle:=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{j-1}^{i,N,h}}(\cdot)\quad\mbox{
and
}\quad\mu_{j,\eta}^{X,N,h}(\cdot):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{j,\eta}^{i,N,h}}(\cdot).$
(14)

######  Remark 3.1 (Comparison to prior art: classical case).

When $\sigma_{0}\equiv 0$, the conditional law $\mathcal{L}^{1}(\cdot)$
becomes the unconditional law of the solution process. If either
$\mathcal{L}^{1}(\cdot)$ is known or $b$ and $\sigma_{1}$ do not depend on
$\mathcal{L}^{1}(\cdot)$, then the McKean–Vlasov SDE (2) becomes a standard
SDE. In such a case, the randomised Milstein scheme (13) considered here
reduces to the one considered in [35]. Indeed, the terms involving measure
derivatives in (8) and hence in (10), (12) and (13) vanish in this case.
Additional terms that appear in (12) and (13) are due to the dependence of the
coefficients $b$, $\sigma_{1}$ and $\sigma_{0}$ on the law of the solution
process (measure variable).

####  3.1.1 Practical Implementation

In this section, we comment on implementation issues, specifically the
sampling of the random time mesh and simulation of the Lévy area. The latter
results explicitly from the presence of common noise and is not specific to
our randomisation method, the former is an essential aspect of the scheme.

We first note that in Equations (12) and (13) the same uniform random
variables $\eta_{1}$, $\ldots$, $\eta_{n_{h}}$ are used for each particle in
the system to identify random points in each sub-interval of the temporal grid
$\rho_{h}$. This approach is similar to the adaptive time-stepping Euler
scheme of [49] where the same random points (arising due to adaptive step-
sizes) are used for each particle of the system. A refined version uses
different time meshes for individual particles, but common, uniform timesteps
for the definition of the empirical measure. In our setting, for each
realisation of the common noise, a single path of random timesteps is
simulated, so that the computational effort of the randomisation is comparable
to the simulation of a discrete Brownian path, but negligible compared to the
simulation of iterated stochastic integrals as required for the Milstein
scheme, which we discuss further now.

If the following, commutative conditions are imposed on the diffusion
coefficients,

$\displaystyle\partial_{x}\sigma_{1}^{\ell}(t,x,\mu)\sigma_{1}^{\ell_{1}}(t,x,\mu)$
$\displaystyle=\partial_{x}\sigma_{1}^{\ell_{1}}(t,x,\mu)\sigma_{1}^{\ell}(t,x,\mu),$

for any $\ell,\ell_{1}=1,\ldots,m_{1}$, $t\in[0,T]$, $x\in\mathbb{R}^{d}$ and
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$, then

$\displaystyle\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell_{1}}_{r}dW^{i,\ell}_{s}$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell}_{r}dW^{i,\ell}_{s}$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}>\ell}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\int_{t_{j-1}}^{s}\Big{(}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell_{1}}_{r}dW^{i,\ell}_{s}$
$\displaystyle\qquad+\partial_{x}\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell}_{r}dW^{i,\ell_{1}}_{s}\Big{)}$
$\displaystyle=\sum_{\ell=1}^{m_{1}}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\frac{1}{2}\big{(}(\Delta
W^{i,\ell})^{2}-h_{j}\big{)}$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}>\ell}^{m_{1}}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\Delta
W^{i,\ell}\Delta W^{i,\ell_{1}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}=1}^{m_{1}}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})(\Delta
W^{i,\ell}\Delta W^{i,\ell_{1}}-h_{j}I_{\\{\ell=\ell_{1}}\\}),$

almost surely for any $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$,
and thus one can write the third term on the right-hand side of (13) as

$\displaystyle\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell}_{s}=\sum_{\ell=1}^{m_{1}}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\Delta
W^{i,\ell}$
$\displaystyle+\frac{1}{2}\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}=1}^{m_{1}}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})(\Delta
W^{i,\ell}\Delta W^{i,\ell_{1}}-h_{j}I_{\\{\ell=\ell_{1}}\\})$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell_{1}=1}^{m_{1}}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j-1}^{k,N,h})$
$\displaystyle\qquad\qquad\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{k,N,h},\mu_{j-1}^{X,N,h})dW^{k,\ell_{1}}_{r}dW^{i,\ell}_{s}$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{0}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{0,\ell_{1}}_{r}dW^{i,\ell}_{s}$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell_{1}=1}^{m_{0}}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j-1}^{k,N,h})$
$\displaystyle\qquad\qquad\sigma_{0}^{\ell_{1}}(t_{j-1},X_{j-1}^{k,N,h},\mu_{j-1}^{X,N,h})dW^{0,\ell_{1}}_{r}dW^{i,\ell}_{s}.$

In the above, the last three terms require the approximation of the Lévy area,
which can be done with the help of the techniques developed in [54]. A similar
conclusion holds for the fourth term on the right-hand side of (13).
Furthermore, the terms involving Lions derivatives are of order
$\mathcal{O}(1/N)$, as shown in [4, Proof of Proposition 2.3] for the regular
case, and hence can be ignored when $N$ is large. In addition, if
$\sigma_{0}\equiv 0$, i.e., the common noise term is not present, the fourth
term on the right-hand side of the above equation can also be dropped and thus
we have

$\displaystyle\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i,\ell}_{s}=\sum_{\ell=1}^{m_{1}}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\Delta
W^{i,\ell}$
$\displaystyle+\frac{1}{2}\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}=1}^{m_{1}}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})(\Delta
W^{i,\ell}\Delta W^{i,\ell_{1}}-h_{j}I_{\\{\ell=\ell_{1}}\\}),$

which leads to a fully implementable randomised Milstein scheme. Summing up,
when $N$ is large and $\sigma_{0}\equiv 0$, the randomised Milstein scheme can
be reduced to

$\displaystyle
X_{j,\eta}^{i,N,h}=X_{j-1}^{i,N,h}+\eta_{j}h_{j}b(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})+\sum_{\ell=1}^{m_{1}}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})(W^{i,\ell}_{t_{j-1}+\eta_{j}h_{j}}-W^{i,\ell}_{t_{j-1}}),$
$\displaystyle
X_{j}^{i,N,h}=X_{j-1}^{i,N,h}+h_{j}b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h})+\sum_{\ell=1}^{m_{1}}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\Delta
W^{i,\ell}$
$\displaystyle+\frac{1}{2}\sum_{\ell=1}^{m_{1}}\sum_{\ell_{1}=1}^{m_{1}}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})(\Delta
W^{i,\ell}\Delta W^{i,\ell_{1}}-h_{j}I_{\\{\ell=\ell_{1}}\\}),$

almost surely for any $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$.

###  3.2 The Main Convergence Result and Its Assumptions

In order to investigate the rate of convergence of the randomised Milstein
scheme (13), we make the following additional assumptions.

######  Assumption H– 4.

There exists a constant $L>0$ such that

$\displaystyle|\partial_{x}\sigma_{u}^{\ell}(t,x,\mu)-\partial_{x}\sigma_{u}^{\ell}(t,x^{\prime},\mu^{\prime})|$
$\displaystyle\leq
L\big{\\{}|x-x^{\prime}|+\mathcal{W}_{2}(\mu,\mu^{\prime})\big{\\}},$
$\displaystyle|\partial_{\mu}\sigma_{u}^{\ell}(t,x,\mu,y)-\partial_{\mu}\sigma_{u}^{\ell}(t,x^{\prime},\mu^{\prime},y^{\prime})|$
$\displaystyle\leq
L\big{\\{}|x-x^{\prime}|+|y-y^{\prime}|+\mathcal{W}_{2}(\mu,\mu^{\prime})\big{\\}},$

for all $u\in\\{0,1\\}$, $\ell\in\\{1,\ldots,m_{u}\\}$, $t\in[0,T]$,
$x,x^{\prime},y,y^{\prime}\in\mathbb{R}^{d}$ and
$\mu,\mu^{\prime}\in\mathcal{P}_{2}(\mathbb{R}^{d})$.

######  Assumption H– 5.

There exists a constant $L>0$ such that

$\displaystyle|\partial_{x}\sigma_{u}^{\ell}(t,x,\mu)\sigma_{v}^{\ell_{1}}(t,x,\mu)-\partial_{x}\sigma_{u}^{\ell}(t,x^{\prime},\mu^{\prime})\sigma_{v}^{\ell_{1}}(t,x^{\prime},\mu^{\prime})|\leq$
$\displaystyle L\\{|x-x^{\prime}|+\mathcal{W}_{2}(\mu,\mu^{\prime})\\},$
$\displaystyle|\partial_{\mu}\sigma_{u}^{\ell}(t,x,\mu,y)\sigma_{v}^{\ell_{1}}(t,y,\mu)-\partial_{\mu}\sigma_{u}^{\ell}(t,x^{\prime},\mu^{\prime},y^{\prime})\sigma_{v}^{\ell_{1}}(t,y^{\prime},\mu^{\prime})|\leq$
$\displaystyle
L\\{|x-x^{\prime}|+|y-y^{\prime}|+\mathcal{W}_{2}(\mu,\mu^{\prime})\\},$

for all $u,v\in\\{0,1\\}$, $\ell\in\\{1,\ldots,m_{u}\\}$,
$\ell_{1}\in\\{1,\ldots,m_{v}\\}$, $t\in[0,T]$,
$x,x^{\prime},y,y^{\prime}\in\mathbb{R}^{d}$ and
$\mu,\mu^{\prime}\in\mathcal{P}_{2}(\mathbb{R}^{d})$.

######  Remark 3.2.

Due to Assumption H–2, we have

$\displaystyle|\partial_{x}\sigma_{u}^{\ell}(t,x,\mu)|+|\partial_{\mu}\sigma_{u}^{\ell}(t,x,\mu,y)|$
$\displaystyle\leq L,$

for all $u\in\\{0,1\\}$, $\ell\in\\{1,\ldots,m_{u}\\}$, $t\in[0,T]$,
$x,y\in\mathbb{R}^{d}$ and $\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$.

Below, we state the main result of this paper containing the error rate for
the approximation of the scheme to the interacting particle system (2). For
notational simplicity, we use $X_{j}^{i,N}$ to represent $X_{t_{j}}^{i,N}$ for
any $i\in\\{1,\ldots,N\\}$ and $j\in\\{1,\ldots,n_{h}\\}$. This result is
proved in Section 5.4.

######  Theorem 3.3.

Let Assumptions H–1 with $\bar{p}\geq 4$, H–2 to H–5 hold. Then, the drift-
randomised Milstein scheme (13) converges in the strong sense to the true
solution of the interacting particle system (2) with order 1. Concretely, for
$q=\bar{p}/2$ and $h\leq\min(1,T)$ we have

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{j\in\\{0,\ldots,n_{h}\\}}|X_{j}^{i,N}-X_{j}^{i,N,h}|\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq
C_{7}C_{10}h,$

where the positive constants $C_{7}$ and $C_{10}$ appear in Proposition 5.5
and Proposition 5.9, respectively.

A direct combination of Proposition 2.5 and Theorem 3.3 delivers the control
on the error for the numerical approximation of the McKean–Vlasov Equation SDE
(2) (and (2)). For convenience of notation, $X_{j}^{i}$ is used to denote
$X_{t_{j}}^{i}$ for any $i\in\\{1,\ldots,N\\}$ and $j\in\\{1,\ldots,n_{h}\\}$.

######  Corollary 3.4.

Let assumptions of Theorem 3.3 hold with $\bar{p}>4$. Then, the drift
randomized Milstein scheme (13) converges in the strong sense to the true
solution of the McKean–Vlasov SDE (2),

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{j\in\\{0,\ldots,n_{h}\\}}|X_{j}^{i}-X_{j}^{i,N,h}|\big{\|}_{\mathscr{L}^{2}(\Omega)}\leq\begin{cases}\sqrt{C_{4}}N^{-1/4}+C_{7}C_{10}h,&\mbox{
if }d<4\\\ \sqrt{C_{4}}N^{-1/4}\sqrt{\log_{2}N}+C_{7}C_{10}h,&\mbox{ if
}d=4\\\ \sqrt{C_{4}}N^{-1/d}+C_{7}C_{10}h,&\mbox{ if }d>4\end{cases}$

with $h\leq\min\\{1,T\\}$ and
$C_{4}:=2(3T+24)TL^{2}e^{4(3T+24)TL^{2}}C_{3}(d,\bar{p},\|X_{0}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}})$
where the constant $C_{3}$ appears in (7), and $C_{7}$ and $C_{10}$ come from
Proposition 5.5 and Proposition 5.9 respectively.

######  Remark 3.5 (Application to optimal control and machine learning).

In an $n$-player stochastic differential game with an
$\left(\mathbb{R}^{d}\right)^{n}$-valued state process
$\boldsymbol{X}=\left(X^{1},\ldots,X^{n}\right)$, agent $i$ chooses a control
process $(\alpha^{i}_{t})$ with values in an action space $A$ so as to
minimize some target functional (see, e.g., the framework considered in [23,
(1.1)]). For a _Markovian control_
$\alpha^{i}_{t}=\alpha^{i}\left(t,\boldsymbol{X}_{t}\right)$, the resulting
dynamics can be described by a (controlled) McKean–Vlasov SDE,

$\displaystyle
dX_{t}^{i}=\widehat{b}\big{(}X_{t}^{i},m_{\boldsymbol{X}_{t}}^{n},\alpha^{i}\left(t,\boldsymbol{X}_{t}\right)\big{)}\,dt+\sigma\,dB_{t}^{i}+\sigma_{0}\,dW_{t}.$
(15)

Throughout, $W$ and $B^{1},\ldots,B^{n}$ are independent Wiener processes, and
we write

$\displaystyle m_{\boldsymbol{x}}^{n}=\frac{1}{n}\sum_{k=1}^{n}\delta_{x_{k}}$

to denote the empirical measure of a vector
$\boldsymbol{x}=\left(x_{1},\ldots,x_{n}\right)$ in
$\left(\mathbb{R}^{d}\right)^{n}$. A similar situation arises in mean-field
control, where a central agent chooses the same feedback control
$\alpha\left(t,\boldsymbol{X}_{t}\right)$ for each agent so as to minimise
their (the central agent’s) objective.

In both these situations, for a non-differentiable control of (Markovian)
type, $\alpha(t,X_{t})$, appearing in $\widehat{b}$ and not appearing in
$\sigma$ or $\sigma_{0}$ (as Assumption H-4 and H-5 require
differentiability), then our approximation scheme will be applicable to the
simulation of the controlled (mean-field) SDE and still produce an
approximation of strong order $1$ as long as one can establish sufficient
regularity of $\alpha^{\cdot}$ such that Assumption H-2 and H-3 holds for the
modified drift $b$

$[0,T]\times\mathbb{R}^{d}\times\mathcal{P}_{2}(\mathbb{R}^{d})\ni(t,x,\mu)\mapsto
b(t,x,\mu):=\widehat{b}\big{(}t,x,\mu,\alpha(t,x)\big{)}.$

This exact same argument would work for controlled SDE in classical settings
where the control is non-differentiable, e.g., when (adapting from (15))

$\displaystyle
dX_{t}=\widehat{b}\big{(}t,X_{t},\alpha\left(t,{X}_{t}\right)\big{)}\,dt+\sigma(t,X_{t})\,dW_{t}.$

In situations when $\alpha$ is a random field our theory would not apply
directly (Kruse et al.’s [7, 34, 35] or [46, 48] as well). It might be
possible to address space-measure mean-field controls
$\alpha(t,X_{t},\mu_{t})$ as in [47] but [47] also shows that proving
regularity properties for $\alpha$ in its measure component is involved.

Lastly, our method also fits into a situation where machine learning is
applied, via tools like reinforcement learning or policy iteration, to solve
the optimal control problem. The requirement is a suitable choice of
control/policy iteration class that would ensure Assumption H-2 and H-3 holds.
A popular choice, especially in the moderate- to high-dimensional context, are
deep neural networks, and a commonly used activation function therein is a
ReLU, which makes the resulting parametric ansatz function Lipschitz but not
everywhere differentiable; see [50] for applications of such a policy gradient
method to non-smooth mean-field control, and to [51] for a proof that the
resulting feedback control remains uniformly Lipschitz over the iterations, in
a setting with controlled drift but without mean-field interaction.

##  4 Moment Bound

In this section, we assume throughout that the conditions of Theorem 3.3 are
in force. Here, we establish moment bounds for the scheme (13), but before
proving it (Lemma 4.3), we state and prove the following auxiliary result.

######  Lemma 4.1.

Let Assumptions  H–2 and H–3 be satisfied. For some $p\geq 2$, if
$X_{j-1}^{i,N,h}\in\mathscr{L}^{p}(\Omega)$ for any $i\in\\{1,\ldots,N\\}$ and
$j\in\\{1,\ldots,n_{h}\\}$, then

$X_{j,\eta}^{i,N,h}\in\mathscr{L}^{p}(\Omega)\mbox{ and
}X_{j}^{i,N,h}\in\mathscr{L}^{p}(\Omega),$

for all $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$.

###### Proof.

Notice that for every $i\in\\{1,\ldots,N\\}$, $X_{j-1}^{i,N,h}$ is
$\mathscr{F}_{j-1}^{h}$-measurable and Assumption H–2 gives the continuity of
$b$ and $\sigma_{u}$ for $u\in\\{0,1\\}$, which in turn implies that
$X_{j,\eta}^{i,N,h}$ defined in (12) is $\mathscr{F}_{j}^{h}$-measurable for
any $j\in\\{1,\ldots,n_{h}\\}$. Also, continuity of
$\partial_{x}\sigma_{u}^{\ell}$ and $\partial_{\mu}\sigma_{u}^{\ell}$ for
$u\in\\{0,1\\}$ and $\ell\in\\{1,2,\ldots,m_{u}\\}$ implies $X_{j}^{i,N,h}$ in
(13) is $\mathscr{F}_{j}^{h}$-measurable for any $j\in\\{1,\ldots,n_{h}\\}$
and $i\in\\{1,\ldots,N\\}$.

As $X_{j-1}^{i,N,h}\in\mathscr{L}^{p}(\Omega)$, we have from Remark 2.1 and
Minkowski’s inequality that

$\displaystyle\big{\|}b(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq\bar{L}\big{\\{}1+\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j-1}^{X,N,h},\delta_{0})\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}$
$\displaystyle\leq\bar{L}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}<\infty,$
(16)

where the last inequality is obtained by using

$\displaystyle\big{\|}\mathcal{W}_{2}(\mu_{j-1}^{X,N,h},\delta_{0})\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq\frac{1}{N}\sum_{i=1}^{N}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}.$

Similarly, we get for $u\in\\{0,1\\}$,

$\displaystyle\big{\|}\sigma_{u}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq\bar{L}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}<\infty,$
(17)

which along with (12) further implies

$\displaystyle\big{\|}X_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq$
$\displaystyle\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}+h_{j}\big{\|}b(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle+h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}\sum_{u=0}^{1}\big{\|}\sigma_{u}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}<\infty,$

for any $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$. Thus, on using
Remark 2.1, one also obtains

$\displaystyle\big{\|}b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq\bar{L}\big{\\{}1+\big{\|}X_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j,\eta}^{X,N,h},\delta_{0})\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}$
$\displaystyle\leq\bar{L}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}<\infty,$
(18)

for all $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$. Moreover,
recall (8) and use Remarks 2.1 and 3.2 along with (4) and (17) to obtain the
following,

$\displaystyle\sum_{\ell=1}^{m_{u}}\big{\|}\Lambda_{\sigma_{u}\sigma_{1}}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq
h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}\sum_{\ell=1}^{m_{u}}\sum_{\ell_{1}=1}^{m_{1}}\big{\|}\partial_{x}\sigma_{u}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle+h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}\sum_{\ell=1}^{m_{u}}\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell_{1}=1}^{m_{1}}\big{\|}\partial_{\mu}\sigma_{u}^{\ell}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j-1}^{k,N,h})$
$\displaystyle\qquad\qquad\sigma_{1}^{\ell_{1}}(t_{j-1},X_{j-1}^{k,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq
2h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}m_{u}m_{1}\bar{L}^{2}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}<\infty,$
(19)

for all $u\in\\{0,1\\}$, $s\in[t_{j-1},t_{j}]$, $i\in\\{1,\ldots,N\\}$ and
$j\in\\{1,\ldots,n_{h}\\}$. Similarly,

$\displaystyle\sum_{\ell=1}^{m_{u}}\big{\|}\Lambda_{\sigma_{u}\sigma_{0}}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq
2h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}m_{u}m_{0}\bar{L}^{2}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}<\infty,$
(20)

for all $u\in\\{0,1\\}$, $s\in[t_{j-1},t_{j}]$, $i\in\\{1,\ldots,N\\}$ and
$j\in\\{1,\ldots,n_{h}\\}$.

Recalling the expressions of $\tilde{\sigma}_{1}$ and $\tilde{\sigma}_{0}$
from (10) and then applying Remark 2.1 along with Equations (20) and (19), we
have

$\displaystyle\sum_{u=0}^{1}\sum_{\ell=1}^{m_{u}}\big{\|}\tilde{\sigma}_{u}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq\sum_{u=0}^{1}\sum_{\ell=1}^{m_{u}}\big{\|}\sigma^{\ell}_{u}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}+\sum_{u=0}^{1}\sum_{\ell=1}^{m_{1}}\big{\|}\Lambda^{\ell}_{\sigma_{1}\sigma_{u}}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\hskip
28.45274pt+\sum_{u=0}^{1}\sum_{\ell=1}^{m_{0}}\big{\|}\Lambda^{\ell}_{\sigma_{0}\sigma_{u}}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle\leq(\bar{L}(m_{0}+m_{1})+2h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}(m_{0}+m_{1})^{2}\bar{L}^{2})\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\big{\\}}<\infty,$
(21)

for all $i\in\\{1,\ldots,N\\}$, $s\in[t_{j-1},t_{j}]$, $i\in\\{1,\ldots,N\\}$
and $j\in\\{1,\ldots,n_{h}\\}$. Thus, by (13) and Theorem 7.1 in [42],

$\displaystyle\big{\|}X_{j}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq$
$\displaystyle\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{p}(\Omega)}+h_{j}\big{\|}b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}$
$\displaystyle+h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}\sum_{u=0}^{1}\sum_{\ell=1}^{m_{u}}\Big{(}\int_{t_{j-1}}^{t_{j}}\big{\|}\tilde{\sigma}_{u}^{\ell}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{p}(\Omega)}^{p}ds\Big{)}^{1/p},$

for all $i\in\\{1,\ldots,N\\}$ and $j\in\\{1,\ldots,n_{h}\\}$, which on using
(4) and (21) completes the proof. ∎

As a consequence of the above lemma, we obtain the following corollary.

######  Corollary 4.2.

Let Assumptions  H–1 with $\bar{p}\geq 2$, H–2 and H–3 hold. Then,
$X_{j}^{i,N,h}\in\mathscr{L}^{\bar{p}}(\Omega)$ for all $i\in\\{1,\ldots,N\\}$
and $j\in\\{1,\ldots,n_{h}\\}$.

###### Proof.

First, let us recall the randomised Milstein scheme given in (13). For $j=1$,
the result holds due to Assumption H–1 and Lemma 4.1. Assume now that the
result is true for $j=k$ for some $k\in\\{1,\ldots,n_{h}\\}$, i.e.,
$X_{k}^{i,N,h}\in\mathscr{L}^{\bar{p}}(\Omega)$ for all
$i\in\\{1,\ldots,N\\}$. Then, Lemma 4.1 yields
$X_{k+1}^{i,N,h}\in\mathscr{L}^{\bar{p}}(\Omega)$. An inductive argument
complete the proof. ∎

The following lemma gives the moment bound of the randomised Milstein scheme
(13).

######  Lemma 4.3.

Let Assumptions  H–1 with $\bar{p}\geq 2$,  H–2 and H–3 hold. Then, for any
time grid $\varrho_{h}$ (11) with $h\leq\min(1,T)$,

$\displaystyle\sup_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{j\in\\{1,\ldots,n_{h}\\}}\big{|}X_{j}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq
C_{5}(1+\sup_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}),$

where

$C_{5}:=\max\\{2,2\bar{L}T+4\bar{L}^{2}T+8\bar{L}^{2}T\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}+2C_{0}\sqrt{T}\\}e^{4T(\bar{L}+2\bar{L}^{2}+4\bar{L}^{2}(\frac{\bar{p}(\bar{p}-1)}{2})^{1/2}+C_{0}^{2})}$

with
$C_{0}:=\bar{C}_{\bar{p}}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}(\bar{L}(m_{0}+m_{1})+2h_{j}^{1/2}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}(m_{0}+m_{1})^{2}\bar{L}^{2})$
and $\bar{C}_{\bar{p}}$ defined in Lemma A.1.

###### Proof.

Recall (12) and use Minkowski’s inequality, Remark 2.1 and $h\leq 1$ to obtain
the following,

$\displaystyle\big{\|}$ $\displaystyle
X_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+h_{j}\big{\|}\eta_{j}b(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+h_{j}^{1/2}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}\sum_{u=0}^{1}\big{\|}\sigma_{u}(t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+(\bar{L}+2\bar{L}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}\big{(}1+\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j-1}^{X,N,h},\delta_{0})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{)}$
$\displaystyle\leq$
$\displaystyle\bar{L}+2\bar{L}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2}+(1+2\bar{L}+4\bar{L}\Big{(}\frac{p(p-1)}{2}\Big{)}^{1/2})\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)},$
(22)

for all $i\in\\{1,\ldots,N\\}$, $N\in\mathbb{N}$ and
$j\in\\{1,\ldots,n_{h}\\}$. Due to (13) and Minkowski’s inequality, for any
$\bar{k}\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$,

$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}$
$\displaystyle\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}h_{j}b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h})\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i}_{s}|\|_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{0}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW_{s}^{0}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)},$

which on application of Lemma A.1 and Remark 2.1 yields

$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}$
$\displaystyle\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+\sum_{j=1}^{\bar{k}}h_{j}\big{\|}b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\bar{C}_{\bar{p}}\big{\|}\big{(}\sum_{j=1}^{\bar{k}}\big{|}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i}_{s}\big{|}^{2}\big{)}^{1/2}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\bar{C}_{\bar{p}}\big{\|}\big{(}\sum_{j=1}^{\bar{k}}\big{|}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{0}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW_{s}^{0}\big{|}^{2}\big{)}^{1/2}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+\sum_{j=1}^{\bar{k}}h_{j}\bar{L}\big{\\{}1+\big{\|}X_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j,\eta}^{X,N,h},\delta_{0})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{\\}}$
$\displaystyle+\bar{C}_{\bar{p}}\big{\|}\sum_{j=1}^{\bar{k}}|\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i}_{s}|^{2}\big{\|}_{\mathscr{L}^{\bar{p}/2}(\Omega)}^{1/2}$
$\displaystyle+\bar{C}_{\bar{p}}\big{\|}\sum_{j=1}^{\bar{k}}|\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{0}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW_{s}^{0}|^{2}\big{\|}_{\mathscr{L}^{\bar{p}/2}(\Omega)}^{1/2}$
$\displaystyle\leq$
$\displaystyle\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+\sum_{j=1}^{\bar{k}}h_{j}\bar{L}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{\\}}$
$\displaystyle+\bar{C}_{\bar{p}}\Big{(}\sum_{j=1}^{\bar{k}}\big{\|}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW^{i}_{s}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}^{2}\Big{)}^{1/2}$
$\displaystyle+\bar{C}_{\bar{p}}\Big{(}\sum_{j=1}^{\bar{k}}\big{\|}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{0}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})dW_{s}^{0}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}^{2}\Big{)}^{1/2},$

for any $i\in\\{1,\ldots,N\\}$ and $\bar{k}\in\\{1,\ldots,n_{h}\\}$.
Furthermore, the application of (22) and Theorem 7.1 in [42] give

$\displaystyle\big{\|}$
$\displaystyle\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\sum_{j=1}^{\bar{k}}h_{j}\bar{L}\big{\\{}1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}$
$\displaystyle+2(1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2})\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{\\}}$
$\displaystyle+\bar{C}_{\bar{p}}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}\sum_{u=0}^{1}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}^{(\bar{p}-2)/\bar{p}}\Big{(}\int_{t_{j-1}}^{t_{j}}\big{\|}\tilde{\sigma}_{u}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}^{\bar{p}}ds\Big{)}^{2/\bar{p}}\Big{)}^{1/2}$
(23)

for any $\bar{k}\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$.

Also, by adapting an argument similar to the one used in (21), we have

$\displaystyle\sum_{u=0}^{1}\big{\|}\tilde{\sigma}_{u}(s,t_{j-1},X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h})\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle\leq(\bar{L}(m_{0}+m_{1})+2h_{j}^{1/2}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}(m_{0}+m_{1})^{2}\bar{L}^{2})\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{\\}},$

for any $s\in[t_{j-1},t_{j}]$, $i\in\\{1,\ldots,N\\}$ and
$j\in\\{1,\ldots,n_{h}\\}$. On substituting the above in Equation (23), one
obtains

$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}$
$\displaystyle\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\bar{L}(1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2})\sum_{j=1}^{\bar{k}}h_{j}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{\\}}$
$\displaystyle+C_{0}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\big{\\}}^{2}\Big{)}^{1/2},$

where
$C_{0}:=\bar{C}_{\bar{p}}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}(\bar{L}(m_{0}+m_{1})+2h_{j}^{1/2}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}(m_{0}+m_{1})^{2}\bar{L}^{2})$,
which on using Young’s inequality yields

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\bar{L}T(1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2})+C_{0}\sqrt{T}$
$\displaystyle+2\bar{L}(1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2})\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+2C_{0}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}|X_{\bar{n}}^{i,N,h}|\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}^{1/2}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\Big{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+\bar{L}T(1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2})+C_{0}\sqrt{T}$
$\displaystyle+2\bar{L}(1+2\bar{L}+4\bar{L}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2})\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+\frac{1}{2}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+2C_{0}^{2}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}.$

Thus, we have

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq
2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}$
$\displaystyle+2\big{(}\bar{L}T+2\bar{L}^{2}T+4\bar{L}^{2}T\Big{(}\frac{\bar{p}(\bar{p}-1}{2}\Big{)}^{1/2})+C_{0}\sqrt{T}\big{)}$
$\displaystyle+4\big{(}\bar{L}+2\bar{L}^{2}++4\bar{L}^{2}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}+C_{0}^{2}\big{)}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}<\infty,$

for any $\bar{k}\in\\{1,\ldots,n_{h}\\}$. The application of the discrete
Grönwall inequality (see Lemma A.2) yields

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,n_{h}\\}}\big{|}X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}\leq\big{(}2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i}\big{\|}_{\mathscr{L}^{\bar{p}}(\Omega)}+2\bar{L}T+4\bar{L}^{2}T$
$\displaystyle+8\bar{L}^{2}T\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{1/2}+2C_{0}\sqrt{T}\big{)}e^{4T(\bar{L}+2\bar{L}^{2}+4\bar{L}^{2}(\frac{\bar{p}(\bar{p}-1)}{2})^{1/2}+C_{0}^{2})}.$

This completes the proof. ∎

##  5 Proof of Main Result

The proof mechanism we use builds from the notions of “bistability” and
“consistency” introduced in [7] (and further explored in [34, 35]). We
introduce a notion of bistability and consistency for the numerical scheme
(13) of the interacting particle system (2) (and associated with the
McKean–Vlasov SDE (2)) by choosing suitable norms and spaces. This choice of
Banach spaces and norms is designed to capture the underlying key feature of
the systems being analysed, namely, that we deal with interacting particle
systems – and to the best of our knowledge, are new.

Throughout this section and in line with the statement of Theorem 3.3, we take
$q=\bar{p}/2$, where $\bar{p}$ comes from Assumption H–1 and assumed to
satisfy $\bar{p}\geq 4$. We next introduce the required notation and
definitions to prove our main result Theorem 3.3.

###  5.1 Quantities of Interest, Norms, Banach Spaces and Residuals

For the time grid $\varrho_{h}$ given in (11), define the Banach spaces
$(\mathscr{G}^{h}_{q},\|\cdot~{}\|_{\mathscr{G}^{h}_{q}})$,
$(\mathscr{G}^{h}_{S,q},\|\cdot\|_{\mathscr{G}^{h}_{S,q}})$ of stochastic grid
processes
$Y^{h}:=\big{\\{}Y_{j}^{i,N,h}\in\mathscr{L}^{q}(\Omega,\mathscr{F}_{j}^{h},\mathbb{P};\,\mathbb{R}^{d});j\in\\{0,1,\ldots,n_{h}\\}\mbox{
and }i\in\\{1,\ldots,N\\}\big{\\}}$ as

$\displaystyle\|Y^{h}\|_{\mathscr{G}^{h}_{q}}:=\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\>j\in\\{0,\ldots,n_{h}\\}}|Y_{j}^{i,N,h}|\big{\|}_{\mathscr{L}^{q}(\Omega)}<\infty,$

and

$\displaystyle\|Y^{h}\|_{\mathscr{G}^{h}_{S,q}}:=$
$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}Y_{0}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{j\in\\{1,\ldots,n_{h}\\}}|\sum_{k=1}^{j}Y_{k}^{i,N,h}|\big{\|}_{\mathscr{L}^{q}(\Omega)}<\infty,$
(24)

respectively. Also, define

$\displaystyle\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j}):=h_{j}b(t_{j-1}+\eta_{j}h_{j},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h})$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})dW^{i,\ell}_{s}+\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})dW_{s}^{0,\ell},$
(25)

almost surely for any $j\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$,
where $Y_{j,\eta}^{i,N,h}$ is defined using (12), and the empirical measures
$\mu_{j-1}^{Y,N,h}$ and $\mu_{j,\eta}^{Y,N,h}$ are defined using (14).

Define $\mathcal{R}[Y^{h}]\in\mathscr{G}^{h}_{S,q}$ for $q\geq 2$ as the
collection of the pointwise residuals $\mathcal{R}^{i,N}_{j}[Y^{h}]$
(associated to executing the scheme with $Y^{h}$) by

$\displaystyle\mathcal{R}^{i,N}_{0}[Y^{h}]$
$\displaystyle=Y^{i,N,h}_{0}-X_{0}^{i,N,h},$
$\displaystyle\mathcal{R}^{i,N}_{j}[Y^{h}]$
$\displaystyle=Y^{i,N,h}_{j}-Y^{i,N,h}_{j-1}-\Gamma^{h}_{j}(Y^{i,N,h}_{j-1},\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j}),\,j\in\\{1,\ldots,n_{h}\\}$
(26)

almost surely for any $i\in\\{1,\ldots,N\\}$.

### Randomised Quadrature Rule for Stochastic Processes.

In this small section we discuss the randomised quadrature rule for stochastic
processes developed in [35]. For this, let us recall the sequence of i.i.d.
uniform random variables $\eta:=\\{\eta_{j}\\}_{j\in\mathbb{N}}$ and the
temporal grid $\varrho_{h}$ from Section 3.1. Consider a stochastic process
$V:[0,T]\times\tilde{\Omega}\mapsto\mathbb{R}^{d}$ on
$(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathbb{P}})$ satisfying
$V\in\mathscr{L}^{p}([0,T]\times\tilde{\Omega})$ for $p\geq 2$. For each
$\bar{n}\in\\{1,\ldots,n_{h}\\}$, $\displaystyle\int_{0}^{t_{\bar{n}}}V(s)ds$
is approximated by the randomised Riemann sum

$\displaystyle\Theta_{\bar{n},\eta}^{h}\big{[}V\big{]}=\displaystyle\sum_{j=1}^{\bar{n}}h_{j}V(t_{j-1}+\eta_{j}h_{j}),$
(27)

which is a random variable on $(\Omega,\mathscr{F},\mathbb{P})$ and an
unbiased estimator of $\displaystyle\int_{0}^{t_{\bar{n}}}V(s)ds$, i.e.,
$\mathbb{E}^{\eta}\Theta_{\bar{n},\eta}^{h}\big{[}V\big{]}=\displaystyle\int_{0}^{t_{\bar{n}}}V(s)ds\in\mathscr{L}^{p}(\tilde{\Omega})$.
Moreover, due to Theorem 4.1 in [35], we have

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}|\Theta_{\bar{n},\eta}^{h}\big{[}V\big{]}-\int_{0}^{t_{\bar{n}}}V(s)ds|\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq
2{\color[rgb]{0,0,0}\bar{C}_{p}}T^{(p-2)/(2p)}\|V\|_{\mathscr{L}^{p}([0,T]\times\tilde{\Omega})}h^{1/2}.$

Additionally, if
$V\in\mathscr{C}^{\alpha}([0,T],\,\mathscr{L}^{p}(\tilde{\Omega}))$ for some
$\alpha\in(0,1]$, then

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\Theta_{\bar{n},\eta}^{h}[V]-\int_{0}^{t_{\bar{n}}}V(s)ds\big{|}\big{\|}_{\mathscr{L}^{p}(\Omega)}\leq{\color[rgb]{0,0,0}\bar{C}_{p}}\sqrt{T}\|V\|_{\mathscr{C}^{\alpha}([0,T],\,\mathscr{L}^{p}(\tilde{\Omega}))}h^{\alpha+1/2}$
(28)

where $\bar{C}_{p}$ is defined in Lemma A.1.

###  5.2 Bistability of the Scheme

We now specify the notion of a scheme’s bistability and show that the proposed
randomised Milstein scheme (13) of the interacting particle system (2) is
bistable (see Proposition 5.5 below). For this, recall the scheme from (13)
and define
$X^{h}:=\\{X_{j}^{i,N,h}:j\in\\{0,\ldots,n_{h}\\},i\in\\{1,\ldots,N\\}\\}$.
Clearly, $X^{h}\in\mathscr{G}_{q}^{h}$ due to Lemma 4.3. Also, recall the
definition of residuals $\mathcal{R}[Y^{h}]$ from (5.1) for some set
$Y^{h}\in\mathscr{G}_{q}^{h}$.

######  Definition 5.1 (Stochastic Bistability).

The randomised Milstein scheme (13) associated to the interacting particle
system (2) is called stochastically bistable if there exist constants
$C_{6},C_{7}>0$, independent of $h$ and $N\in\mathbb{N}$, such that for any
arbitrary time grid $\varrho_{h}$ as given in (11) and for any
$Y^{h}\in\mathscr{G}_{q}^{h}$, the following holds,

$\displaystyle
C_{6}\|\mathcal{R}[Y^{h}]\|_{\mathscr{G}_{S,q}^{h}}\leq\|Y^{h}-X^{h}\|_{\mathscr{G}_{q}^{h}}\leq
C_{7}\|\mathcal{R}[Y^{h}]\|_{\mathscr{G}_{S,q}^{h}},$

for $q\geq 2$.

Let us first establish some useful lemmas.

######  Lemma 5.2.

Let Assumptions H–2 and H–3 hold. Then, for any $q\geq 2$ and
$Y^{h},Z^{h}\in\mathscr{G}_{q}^{h}$,

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}$
$\displaystyle\big{|}\sum_{j=1}^{\bar{n}}h_{j}\big{[}b(t_{j-1}+\eta_{j}h_{j},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h})-b(t_{j-1}+\eta_{j}h_{j},Z_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Z,N,h})\big{]}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq
C_{8}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y^{i,N,h}_{\bar{n}}-Z^{i,N,h}_{\bar{n}}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)},$

for any $\bar{k}\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$ where
$C_{8}:=2L+4L^{2}T+8L^{2}\sqrt{Tq(q-1)/2}$.

###### Proof.

Notice that all the terms in the statement of this lemma are well-defined due
to Remark 2.1. By using (12) with $Y^{h},Z^{h}\in\mathscr{G}_{q}^{h}$ along
with Minkowski’s inequality and Assumption H–2, we have

$\displaystyle\big{\|}Y_{j,\eta}^{i,N,h}-Z_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq\big{\|}Y_{j-1}^{i,N,h}-Z_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+h_{j}\big{\|}\eta_{j}(b(t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-b(t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h}))\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\sqrt{h_{j}}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\sum_{u=0}^{1}\big{\|}\sqrt{\eta_{j}}(\sigma_{u}(t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\hskip
170.71652pt-\sigma_{u}(t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h}))\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\big{\|}Y_{j-1}^{i,N,h}-Z_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\Big{(}LT+2L\sqrt{T}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\Big{)}$
$\displaystyle\qquad\times\big{\\{}\big{\|}Y_{j-1}^{i,N,h}-Z_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j-1}^{Y,N,h},\mu_{j-1}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)}\big{\\}}$
$\displaystyle\leq\Big{(}1+2LT+4L\sqrt{T}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\Big{)}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-Z_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)},$
(29)

where the last inequality is obtained by using

$\displaystyle\mathcal{W}_{2}(\mu_{j-1}^{Y,N,h},\mu_{j-1}^{Z,N,h})\leq\Big{(}\frac{1}{N}\sum_{i=1}^{N}\big{|}Y_{j-1}^{i,N,h}-Z_{j-1}^{i,N,h}\big{|}^{2}\Big{)}^{1/2},$

for all $i\in\\{1,\ldots,N\\}$, $N\in\mathbb{N}$ and
$j\in\\{1,\ldots,n_{h}\\}$.

Now, the application Assumption H–2 yields

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}$
$\displaystyle
h_{j}\big{[}b(t_{j-1}+\eta_{j}h_{j},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h})-b(t_{j-1}+\eta_{j}h_{j},Z_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Z,N,h})\big{]}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq
L\sum_{j=1}^{\bar{k}}h_{j}\big{\\{}\big{\|}Y_{j,\eta}^{i,N,h}-Z_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j,\eta}^{Y,N,h},\mu_{j,\eta}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)}\big{\\}}$
$\displaystyle\leq
2L\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}Y_{j,\eta}^{i,N,h}-Z_{j,\eta}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)},$

and then one uses (29) to complete the proof. ∎

######  Lemma 5.3.

Let Assumptions H–2, H–3, H–5 hold. Then, for any $q\geq 2$ and
$Y^{h},Z^{h}\in\mathscr{G}_{q}^{h}$,

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\qquad\qquad\qquad\qquad-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW_{s}^{i,\ell}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\qquad\qquad\qquad\qquad-\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW_{s}^{0,\ell}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\qquad\leq
C_{9}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y^{i,N,h}_{\bar{n}}-Z^{i,N,h}_{\bar{n}}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}^{2}\Big{)}^{1/2},$

for all $i\in\\{1,\ldots,N\\}$ and $\bar{k}\in\\{1,\ldots,n_{h}\\}$, where
$C_{9}:=(m_{0}+m_{1})\bar{C}_{q}\sqrt{q(q-1)/2}\big{(}2L+5(m_{0}+m_{1})L\sqrt{Tq(q-1)/2}\big{)}$
and $\bar{C}_{q}$ is given in Lemma A.1.

###### Proof.

Notice that

$\displaystyle\Big{\\{}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW_{s}^{i,\ell}\Big{\\}}_{\bar{n}},$

with ${\bar{n}\in\\{0,\ldots,n_{h}\\}}$ is an
$\\{\mathscr{F}_{\bar{n}}^{h}\\}_{\bar{n}\in\\{0,\ldots,n_{h}\\}}$-adapted
standard martingale. Thus, on using Lemma A.1 and Minkowski’s inequality, one
obtains, for all $i\in\\{1,\ldots,N\\}$ and $\bar{k}\in\\{1,\ldots,n_{h}\\}$,

$\displaystyle\big{\|}$
$\displaystyle\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\hskip
142.26378pt-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW_{s}^{i,\ell}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq\bar{C}_{q}\big{\|}\sum_{j=1}^{\bar{k}}\big{|}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\hskip
142.26378pt-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW_{s}^{i,\ell}\big{|}^{2}\big{\|}_{\mathscr{L}^{q/2}(\Omega)}^{1/2}$
$\displaystyle\leq\bar{C}_{q}\Big{(}\sum_{j=1}^{\bar{k}}\big{\|}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\hskip
142.26378pt-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})dW_{s}^{i,\ell}\big{\|}^{2}_{\mathscr{L}^{q}(\Omega)}\Big{)}^{1/2}$
$\displaystyle\leq\sqrt{m}_{1}\bar{C}_{q}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}^{(q-2)/q}\sum_{\ell=1}^{m_{1}}\big{[}\int_{t_{j-1}}^{t_{j}}\big{\|}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle\hskip
142.26378pt-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{\|}^{q}_{\mathscr{L}^{q}(\Omega)}ds\big{]}^{2/q}\Big{)}^{1/2},$
(30)

where the last inequality is obtained due to Theorem 7.1 in [42].

Now, recall (8), (9) and (10) and use Assumption H–2 to get the following,

$\displaystyle\big{\|}\tilde{\sigma}_{1}^{\ell}(s$
$\displaystyle,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\|_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\|\sigma_{1}^{\ell}(t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\sigma_{1}^{\ell}(t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\|_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\sum_{u=0}^{1}\big{\|}\Lambda_{\sigma_{1}\sigma_{u}}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\Lambda_{\sigma_{1}\sigma_{u}}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
L\big{\\{}\big{\|}Y_{j-1}^{i,N,h}-Z_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{j-1}^{Y,N,h},\mu_{j-1}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)}\big{\\}}$
$\displaystyle+h_{j}^{1/2}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\sum_{u=0}^{1}\sum_{\ell_{1}=1}^{m_{u}}\big{\|}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})\sigma_{u}^{\ell_{1}}(t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})$
$\displaystyle-\partial_{x}\sigma_{1}^{\ell}(t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\sigma_{u}^{\ell_{1}}(t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+h_{j}^{1/2}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\frac{1}{N}\sum_{u=0}^{1}\sum_{\ell_{1}=1}^{m_{u}}\sum_{k=1}^{N}\big{\|}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},Y_{j-1}^{k,N,h})$
$\displaystyle\times\sigma_{u}^{\ell_{1}}(t_{j-1},Y_{j-1}^{k,N,h},\mu_{j-1}^{Y,N,h})-\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h},Z_{j-1}^{k,N,h})$
$\displaystyle\times\sigma_{u}^{\ell_{1}}(t_{j-1},Z_{j-1}^{k,N,h},\mu_{j-1}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)},$

which due to Assumption H–5 yields

$\displaystyle\big{\|}\tilde{\sigma}_{1}^{\ell}(s$
$\displaystyle,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\|_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\Big{(}L+2(m_{0}+m_{1})L\Big{(}\frac{Tq(q-1)}{2}\Big{)}^{1/2}\Big{)}$
$\displaystyle\times\big{\\{}\big{\|}Y_{j-1}^{i,N,h}-Z_{j-1}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\big{\|}W_{2}(\mu_{j-1}^{Y,N,h},\mu_{j-1}^{Z,N,h})\big{\|}_{\mathscr{L}^{q}(\Omega)}\big{\\}}$
$\displaystyle+(m_{0}+m_{1})L\Big{(}\frac{Tq(q-1)}{2}\Big{)}^{1/2}\frac{1}{N}\sum_{k=1}^{N}\big{\|}Y_{j-1}^{k,N,h}-Z_{j-1}^{k,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\Big{(}2L+5(m_{0}+m_{1})L\Big{(}\frac{Tq(q-1)}{2}\Big{)}^{1/2}\Big{)}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-Z_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)},$

for any $s\in[t_{j-1},t_{j}]$, $j\in\\{1,\ldots,n_{h}\\}$ and
$i\in\\{1,\ldots,N\\}$. The proof of the first part of the lemma is completed
by substituting the above in (5.2). A bound for the terms involving
$\tilde{\sigma}_{0}$ follows by similar arguments. ∎

As a consequence of Lemmas 5.2 and 5.3, we obtain the following corollary.

######  Corollary 5.4.

If Assumptions H–2, H–3 and H–5 are satisfied. Then, for any $q\geq 2$ and
$Y^{h},Z^{h}\in\mathscr{G}_{q}^{h}$,

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\big{[}\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j})$
$\displaystyle\hskip
142.26378pt-\Gamma^{h}_{j}(Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h},Z_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Z,N,h},\eta_{j})\big{]}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
C_{8}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y^{i,N,h}_{\bar{n}}-Z^{i,N,h}_{\bar{n}}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+C_{9}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y^{i,N,h}_{\bar{n}}-Z^{i,N,h}_{\bar{n}}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}^{2}\Big{)}^{1/2}$
$\displaystyle\leq$
$\displaystyle(C_{8}T+C_{9}\sqrt{T})\big{\|}Y^{h}-Z^{h}\big{\|}_{\mathscr{G}_{q}^{h}},$

for all $\bar{k}\in\\{1,\ldots,n_{h}\\}$ and $i\in\\{1,\ldots,N\\}$, where the
constants $C_{8}$ and $C_{9}$ appear in Lemmas 5.2 and 5.3, respectively.

###### Proof.

Recall $\Gamma^{h}_{j}$ from (5.1) and write

$\displaystyle\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},$
$\displaystyle\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j})-\Gamma^{h}_{j}(Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h},Z_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Z,N,h},\eta_{j})$
$\displaystyle=h_{j}\big{[}b(t_{j-1}+\eta_{j}h_{j},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h})-b(t_{j-1}+\eta_{j}h_{j},Z_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Z,N,h})\big{]}$
$\displaystyle\quad+\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW^{i,\ell}_{s}$
$\displaystyle\quad+\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\big{[}\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h})-\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},Z_{j-1}^{i,N,h},\mu_{j-1}^{Z,N,h})\big{]}dW_{s}^{0,\ell}.$

Then, on using Lemmas 5.2 and 5.3, we get the required result. ∎

We now prove the bistability of the scheme (13) in the following proposition.

######  Proposition 5.5.

Let Assumptions H–1 with $\bar{p}\geq 4$ and set $q=\bar{p}/2\geq 2$, H–2, H–3
and H–5 hold. Then, for any $Y^{h},X^{h}\in\mathscr{G}_{q}^{h}$, we have

$\displaystyle
C_{6}\|\mathcal{R}[Y^{h}]\|_{\mathscr{G}_{S,q}^{h}}\leq\|Y^{h}-X^{h}\|_{\mathscr{G}_{q}^{h}}\leq
C_{7}\|\mathcal{R}[Y^{h}]\|_{\mathscr{G}_{S,q}^{h}},$

in other words, the randomised Milstein scheme given in (13) is stochastically
bistable in the sense of Definition 5.1 with
$C_{6}:=\frac{1}{3+C_{8}T+C_{9}\sqrt{T}}$ and
$C_{7}:=2e^{(2C_{8}+C_{9}^{2})T}$, where the constants $C_{8}$ and $C_{9}$
appear in Lemmas 5.2 and 5.3, respectively.

###### Proof.

Recall Equations (13), (5.1) and (5.1) to write

$\displaystyle X_{\bar{n}}^{i,N,h}$ $\displaystyle-
X_{0}^{i,N,h}-\sum_{j=1}^{\bar{n}}\Gamma^{h}_{j}(X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h},\eta_{j})=0,$
$\displaystyle\sum_{j=1}^{\bar{n}}\mathcal{R}_{j}^{i,N}[Y^{h}]$
$\displaystyle=Y_{\bar{n}}^{i,N,h}-Y_{0}^{i,N,h}-\sum_{j=1}^{\bar{n}}\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j})$
$\displaystyle=\big{(}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{)}-\big{(}Y_{0}^{i,N,h}-X_{0}^{i,N,h}\big{)}$
$\displaystyle-\sum_{j=1}^{\bar{n}}\big{(}\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},$
$\displaystyle
Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j})-\Gamma^{h}_{j}(X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h},\eta_{j})\big{)},$
(31)

for any $\bar{n}\in\\{1,\ldots,n_{h}\\}$, $i\in\\{1,\ldots,N\\}$ and then use
the Spijker norm (24) to get the following,

$\displaystyle\big{\|}\mathcal{R}[Y$
${}^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}}=\max_{i\in\\{1\ldots,N\\}}\big{\|}\big{|}\mathcal{R}_{0}^{i,N}[Y^{h}]\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\mathcal{R}_{j}^{i,N}[Y^{h}]\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
2\max_{i\in\\{1,\ldots,N\\}}\big{\|}Y_{0}^{i,N,h}-X_{0}^{i,N,h}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\big{(}\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j})$
$\displaystyle\qquad-\Gamma^{h}_{j}(X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h},\eta_{j})\big{)}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)},$

which on using Corollary 5.4 yields

$\displaystyle\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}}\leq(3+C_{8}T+C_{9}\sqrt{T})\|Y^{h}-X^{h}\|_{\mathscr{G}_{q}^{h}}.$
(32)

Further, by rearranging the terms of (31) and using Minkowski’s inequality,
one obtains

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$
$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\big{(}\Gamma^{h}_{j}(Y_{j-1}^{i,N,h},\mu_{j-1}^{Y,N,h},Y_{j,\eta}^{i,N,h},\mu_{j,\eta}^{Y,N,h},\eta_{j})$
$\displaystyle\qquad-\Gamma^{h}_{j}(X_{j-1}^{i,N,h},\mu_{j-1}^{X,N,h},X_{j,\eta}^{i,N,h},\mu_{j,\eta}^{X,N,h},\eta_{j})\big{)}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\max_{i\in\\{1,\ldots,N\\}}\big{\|}|\mathcal{R}_{0}^{i,N}[Y^{h}]|\big{\|}_{\mathscr{L}^{q}(\Omega)}+\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{1,\ldots,\bar{k}\\}}\big{|}\sum_{j=1}^{\bar{n}}\mathcal{R}_{j}^{i,N}[Y^{h}]\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)},$

and then application of Corollary 5.4 gives the following,

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle\leq$ $\displaystyle
C_{8}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+C_{9}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}^{2}\Big{)}^{1/2}+\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}}$
$\displaystyle\leq$ $\displaystyle
C_{8}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}^{1/2}$
$\displaystyle\qquad\times
C_{9}\Big{(}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\Big{)}^{1/2}+\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}},$

which due to Young’s inequality yields

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq\frac{1}{2}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\big{(}C_{8}+\frac{C_{9}^{2}}{2}\big{)}\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}}.$

This further implies

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq
2\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}}$
$\displaystyle+(2C_{8}+C_{9}^{2})\sum_{j=1}^{\bar{k}}h_{j}\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{\bar{n}\in\\{0,\ldots,j-1\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}.$

Due to Lemma A.2, we get

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}$
$\displaystyle\big{\|}\max_{\bar{n}\in\\{0,\ldots,\bar{k}\\}}\big{|}Y_{\bar{n}}^{i,N,h}-X_{\bar{n}}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq
2e^{(2C_{8}+C_{9}^{2})T}\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}},$

for any $\bar{k}\in\\{0,\ldots,n_{h}\\}$, which further implies

$\displaystyle\|Y^{h}-X^{h}\|_{\mathscr{G}_{q}^{h}}\leq
2e^{(2C_{8}+C_{9}^{2})T}\big{\|}\mathcal{R}[Y^{h}]\big{\|}_{\mathscr{G}_{S,q}^{h}},$

and the proof is completed by combining the above with (32). ∎

###  5.3 Consistency of the Scheme

Recall the temporal grid $\varrho_{h}$ from (11) and set the values of the
interacting particle system (2) over the grid points of $\varrho_{h}$ as the
set $X^{\varrho_{h}}:=\\{X^{i,N}_{j}\\}_{j\in\\{0,\ldots,n_{h}\\}}$ where
$X_{j}^{i,N}=X_{t_{j}}^{i,N}$ for any $j\in\\{0,\ldots,n_{h}\\}$ and
$i\in\\{1,\ldots,N\\}$. Notice that $X^{\varrho_{h}}$ is different from
$X^{h}$, where the later stands for the randomised Milstein scheme (13).
Further, $\mathcal{R}[X^{\varrho_{h}}]$ is defined using $X^{\varrho_{h}}$ in
(5.1). Next, we introduce the notion of a scheme’s consistency and establish
that the randomised Milstein scheme (13) is consistent (see below Proposition
5.9).

######  Definition 5.6 (Consistency).

The randomised Milstein scheme (13) for the interacting particle system (2) is
called consistent of order $\gamma>0$ if there exists a constant $C_{10}>0$,
independent of $h$ and $N\in\mathbb{N}$, such that

$\displaystyle\|\mathcal{R}[X^{\varrho_{h}}]\|_{\mathscr{G}_{S,q}^{h}}\leq
C_{10}h^{\gamma},$

for $q\geq 2$.

In the context of the measure dependent drift coefficient, we obtain the
following randomised quadrature rule.

######  Corollary 5.7.

Let Assumptions H–1 with $\bar{p}\geq 4$ and set $q=\bar{p}/2\geq 2$, H–2, H–
3 hold. Then, for all $i\in\\{1,\ldots,N\\}$, the operator
$\Theta_{\bar{n},\eta}^{h}[\cdot]$ defined in (27) and applied to $b$
satisfies

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\Theta_{\bar{n},\eta}^{h}[b]-\int_{0}^{t_{\bar{n}}}b(s,X_{s}^{i,N},\mu_{s}^{X,N})ds\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq
C_{11}h,$

where
$C_{11}:=\sqrt{T}{\color[rgb]{0,0,0}\bar{C}_{q}}\big{(}2\bar{L}+4\bar{L}C_{1}^{1/\bar{p}}+2{\color[rgb]{0,0,0}\bar{L}}C_{2}^{1/\bar{p}}\big{)}\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}\big{)}^{1/\bar{p}}$
and the positive constants $C_{1}$, $C_{2}$, $\bar{C}_{q}$ appear in
Proposition 2.2, Lemma 2.4 and Lemma A.1, respectively.

###### Proof.

We take $V(s)=b\big{(}s,X_{s}^{i,N},\mu_{s}^{X,N}\big{)}$ for any $s\in[0,T]$
and $i\in\\{1,\ldots,N\\}$, which is defined on
$(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathbb{P}})$. Notice that Remarks
2.1 and 2.3 along with Hölder’s inequality (as $q=\bar{p}/2$) imply

$\displaystyle\|V\|_{\mathscr{L}^{q}([0,T]\times\tilde{\Omega})}$
$\displaystyle=\Big{(}\int_{0}^{T}\|b\big{(}s,X_{s}^{i,N},\mu_{s}^{X,N}\big{)}\|_{\mathscr{L}^{q}(\tilde{\Omega})}^{q}ds\Big{)}^{1/q}$
$\displaystyle\leq\bar{L}\Big{(}\int_{0}^{T}\big{(}1+\|X_{s}^{i,N}\|_{\mathscr{L}^{q}(\tilde{\Omega})}+\|\mathcal{W}_{2}(\mu_{s}^{X,N},\delta_{0})\|_{\mathscr{L}^{q}(\tilde{\Omega})}\big{)}^{q}ds\Big{)}^{1/q}$
$\displaystyle\leq\bar{L}\Big{(}\int_{0}^{T}(1+2\max_{i\in\\{1,\ldots,N\\}}\|X_{s}^{i,N}\|_{\mathscr{L}^{q}(\tilde{\Omega})})^{q}ds\Big{)}^{1/q}$
$\displaystyle\leq\bar{L}\Big{(}\int_{0}^{T}(1+2\max_{i\in\\{1,\ldots,N\\}}\|X_{s}^{i,N}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})})^{q}ds\Big{)}^{1/q}$
$\displaystyle\leq\bar{L}\Big{(}\int_{0}^{T}(1+2C_{1}^{1/\bar{p}}(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})})^{1/\bar{p}})^{q}ds\Big{)}^{1/q}$
$\displaystyle\leq\bar{L}(T)^{1/q}(1+2C_{1}^{1/\bar{p}}(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})})^{1/\bar{p}})<\infty.$

Also, due to Assumptions H–2, H–3 and Lemma 2.4, for $t,t^{\prime}\in[0,T]$,

$\displaystyle\big{\|}V(t)-V(t^{\prime})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}=\big{\|}b(t,X_{t}^{i,N},\mu_{t}^{X,N})-b(t^{\prime},X_{t^{\prime}}^{i,N},\mu_{t^{\prime}}^{X,N})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}\leq
L\big{\\{}\big{\|}X_{t}^{i,N}-X_{t^{\prime}}^{i,N}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle\quad+\big{\|}\mathcal{W}_{2}(\mu_{t}^{X,N},\mu_{t^{\prime}}^{X,N})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}+|t-t^{\prime}|^{1/2}\big{(}1+\big{\|}X^{i,N}_{t}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}+\big{\|}\mathcal{W}_{2}(\mu_{t}^{X,N},\delta_{0})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}\big{)}\\}$
$\displaystyle\leq
2L\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{t}^{i,N}-X_{t^{\prime}}^{i,N}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}+|t-t^{\prime}|^{1/2}L\\{1+2C_{1}^{1/\bar{p}}(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}})^{1/\bar{p}}\\}$
$\displaystyle\leq\\{2L(C_{1}^{1/\bar{p}}+C_{2}^{1/\bar{p}})(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}})^{1/\bar{p}}+L\\}|t-t^{\prime}|^{1/2},$
(33)

for all $i\in\\{1,\ldots,N\\}$. Due to (3),
$\|V\|_{\mathscr{C}^{1/2}([0,T],\,\mathscr{L}^{q}(\tilde{\Omega}))}\leq
2\bar{L}+\bar{L}(4C_{1}^{1/\bar{p}}+2C_{2}^{1/\bar{p}}\big{)}\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}\big{)}^{1/\bar{p}}$,
as $L\leq\bar{L}$. Then, (28) completes the proof. ∎

######  Lemma 5.8.

Let Assumptions H–1 with $\bar{p}\geq 4$ and set $q=\bar{p}/2\geq 2$, H–2 to
H–4 hold. Let $h\leq\min(1,T)$. Then,

$\displaystyle\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]}dW^{i,\ell}_{s}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle+\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\big{[}\sigma_{0}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]}dW^{0,\ell}_{s}\big{|}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle\hskip 113.81102pt\leq C_{12}h,$

for all $i\in\\{1,\ldots,N\\}$, where

$\displaystyle C_{12}$
$\displaystyle:=(m_{0}+m_{1})\bar{C}_{q}\sqrt{{Tq(q-1)}/{2}}\Big{\\{}L+2\bar{L}^{2}+4L^{2}\sqrt{{q(q-1)}/{2}}+\Big{(}2(L+2\bar{L}^{2})C_{1}^{1/\bar{p}}$
$\displaystyle\qquad+8L^{2}\sqrt{{q(q-1)}{2}}\big{(}C_{1}^{1/\bar{p}}+C_{2}^{1/\bar{p}}\big{)}+4LC_{2}^{2/\bar{p}}\Big{)}\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}^{\bar{p}}\big{)}^{2/\bar{p}}\Big{\\}}$

and the positive constants $C_{1}$, $C_{2}$ and $\bar{C}_{q}$ appear in
Proposition 2.2, Lemma 2.4 and Lemma A.1, respectively.

###### Proof.

For notational simplicity, define

$\zeta^{\ell}(s):=\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})$

for all $\ell\in\\{1,\ldots,m_{1}\\}$, $i\in\\{1,\ldots,N\\}$,
$j\in\\{1,\ldots,n_{h}\\}$ and $s\in[t_{j-1},t_{j}]$. Notice that

$\Big{\\{}\displaystyle\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j}-1}^{t_{j}}\zeta^{\ell}(s)dW^{i,\ell}_{s}\Big{\\}}_{\bar{n}\in\\{0,\ldots,n_{h}\\}}$

is an
$\\{\tilde{\mathscr{F}}_{t_{\bar{n}}}\\}_{\bar{n}\in\\{0,\ldots,n_{h}\\}}$-adapted
standard martingale. Thus, by using Lemma A.1 and Theorem 7.1 in [42], we get

$\displaystyle\big{\|}$
$\displaystyle\max_{\bar{n}\in\\{0,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\zeta^{\ell}(s)dW^{i,\ell}_{s}\big{|}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}\leq\bar{C}_{q}\big{\|}\big{(}\sum_{j=1}^{n_{h}}\big{|}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\zeta^{\ell}(s)dW^{i,\ell}_{s}\big{|}^{2}\big{)}^{1/2}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle=\bar{C}_{q}\big{\|}\sum_{j=1}^{n_{h}}\big{|}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\zeta^{\ell}(s)dW^{i,\ell}_{s}\big{|}^{2}\big{\|}_{\mathscr{L}^{q/2}(\tilde{\Omega})}^{1/2}\leq\sqrt{m}_{1}\bar{C}_{q}\Big{(}\sum_{j=1}^{n_{h}}\sum_{\ell=1}^{m_{1}}\big{\|}\int_{t_{j-1}}^{t_{j}}\zeta^{\ell}(s)dW^{i,\ell}_{s}\big{\|}^{2}_{\mathscr{L}^{q}(\tilde{\Omega})}\Big{)}^{1/2}$
$\displaystyle\leq\sqrt{m}_{1}\bar{C}_{q}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\Big{(}\sum_{j=1}^{n_{h}}\sum_{\ell=1}^{m_{1}}h_{j}^{(q-2)/{q}}\big{(}\int_{t_{j-1}}^{t_{j}}\big{\|}\zeta^{\ell}(s)\big{\|}^{q}_{\mathscr{L}^{q}(\tilde{\Omega})}ds\big{)}^{2/q}\Big{)}^{1/2},$
(34)

for all $i\in\\{1,\ldots,N\\}$. Now, recall $\tilde{\sigma}_{1}^{\ell}$ from
(10) to write the following for any $s\in[t_{j-1},t_{j}]$,

$\displaystyle\zeta^{\ell}(s)=$
$\displaystyle\big{[}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\sigma_{1}^{\ell}(t_{j-1},X_{s}^{i,N},\mu_{s}^{X,N})\big{]}+\big{[}\sigma_{1}^{\ell}(t_{j-1},X_{s}^{i,N},\mu_{s}^{X,N})$
$\displaystyle-\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})-\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})(X_{s}^{i,N}-X_{j-1}^{i,N})$
$\displaystyle-\frac{1}{N}\sum_{k=1}^{N}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})(X_{s}^{k,N}-X_{j-1}^{k,N})\big{]}+\big{[}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})$
$\displaystyle\times(X_{s}^{i,N}-X_{j-1}^{i,N})+\frac{1}{N}\sum_{k=1}^{N}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})(X_{s}^{k,N}-X_{j-1}^{k,N})$
$\displaystyle\hskip
28.45274pt-\Lambda_{\sigma_{1}\sigma_{1}}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})-\Lambda_{\sigma_{1}\sigma_{0}}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]},$

which on the application of Assumption H–3, Lemma A.3 along with (2) and the
values of $\Lambda_{\sigma_{1}\sigma_{1}}^{\ell}$,
$\Lambda_{\sigma_{1}\sigma_{0}}^{\ell}$ from (8), (9) yields

$\displaystyle\zeta^{\ell}$ $\displaystyle(s)\leq
L(1+|X_{s}^{i,N}|+\mathcal{W}_{2}(\mu_{s}^{X,N},\delta_{0}))|s-t_{j-1}|+\frac{3L}{2}|X_{s}^{i,N}-X_{j-1}^{i,N}|^{2}$
$\displaystyle+\frac{5L}{2}\frac{1}{N}\sum_{k=1}^{N}|X_{s}^{k,N}-X_{j-1}^{k,N}|^{2}+\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})b(r,X_{r}^{i,N},\mu_{r}^{X,N})dr$
$\displaystyle+\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})(\sigma_{1}(r,X_{r}^{i,N},\mu_{r}^{X,N})-\sigma_{1}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N}))dW^{i}_{r}$
$\displaystyle+\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})(\sigma_{0}(r,X_{r}^{i,N},\mu_{r}^{X,N})-\sigma_{0}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N}))dW^{0}_{r}$
$\displaystyle+\frac{1}{N}\sum_{k=1}^{N}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})b(r,X_{r}^{k,N},\mu_{r}^{X,N})dr$
$\displaystyle+\frac{1}{N}\sum_{k=1}^{N}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})$
$\displaystyle\hskip
85.35826pt\times(\sigma_{1}(r,X_{r}^{k,N},\mu_{r}^{X,N})-\sigma_{1}(t_{j-1},X_{j-1}^{k,N},\mu_{j-1}^{X,N}))dW^{k}_{r}$
$\displaystyle+\frac{1}{N}\sum_{k=1}^{N}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})$
$\displaystyle\hskip
85.35826pt\times(\sigma_{0}(r,X_{r}^{k,N},\mu_{r}^{X,N})-\sigma_{0}(t_{j-1},X_{j-1}^{k,N},\mu_{j-1}^{X,N}))dW^{0}_{r},$

for all $i\in\\{1,\ldots,N\\}$, $j\in\\{1,\ldots,n_{h}\\}$ and
$\ell\in\\{1,\ldots,m_{1}\\}$. This further implies, due to Remarks 2.1 and
3.2,

$\displaystyle\|$
$\displaystyle\zeta^{\ell}(s)\|_{\mathscr{L}^{q}(\tilde{\Omega})}\leq
L(1+2\max_{i\in\\{1,\ldots,N\\}}\|X_{s}^{i,N}\|_{\mathscr{L}^{q}(\tilde{\Omega})})h_{j}+4L\max_{i\in\\{1,\ldots,N\\}}\|X_{s}^{i,N}-X_{j-1}^{i,N}\|_{\mathscr{L}^{2q}(\tilde{\Omega})}^{2}$
$\displaystyle+h_{j}^{(q-1)/{q}}2\bar{L}^{2}\big{(}\int_{t_{j-1}}^{s}(1+2\max_{i\in\\{1,\ldots,N\\}}\|X_{r}^{i,N}\|_{\mathscr{L}^{q}(\tilde{\Omega})})^{q}dr\big{)}^{1/q}$
$\displaystyle+\big{\|}\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})(\sigma_{1}(r,X_{r}^{i,N},\mu_{r}^{X,N})-\sigma_{1}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N}))dW^{i}_{r}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle+\big{\|}\int_{t_{j-1}}^{s}\partial_{x}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})(\sigma_{0}(r,X_{r}^{i,N},\mu_{r}^{X,N})-\sigma_{0}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N}))dW^{0}_{r}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle+\big{\|}\frac{1}{N}\sum_{k=1}^{N}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})$
$\displaystyle\hskip
85.35826pt\times(\sigma_{1}(r,X_{r}^{k,N},\mu_{r}^{X,N})-\sigma_{1}(t_{j-1},X_{j-1}^{k,N},\mu_{j-1}^{X,N}))dW^{k}_{r}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle+\big{\|}\frac{1}{N}\sum_{k=1}^{N}\int_{t_{j-1}}^{s}\partial_{\mu}\sigma_{1}^{\ell}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j-1}^{k,N})$
$\displaystyle\hskip
85.35826pt\times(\sigma_{0}(r,X_{r}^{k,N},\mu_{r}^{X,N})-\sigma_{0}(t_{j-1},X_{j-1}^{k,N},\mu_{j-1}^{X,N}))dW^{0}_{r}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}.$

Moreover, the application of Hölder’s inequality (as $q=\bar{p}/2)$, Remarks
2.3 and 3.2, Lemma 2.4 and Theorem 7.1 in [42] yield

$\displaystyle\|$
$\displaystyle\zeta^{\ell}(s)\|_{\mathscr{L}^{q}(\tilde{\Omega})}\leq\big{\\{}(L+2\bar{L}^{2})\big{(}1+2C_{1}^{1/\bar{p}}(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i,N}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}})^{1/\bar{p}}\big{)}$
$\displaystyle+4LC_{2}^{2/\bar{p}}\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i,N}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\big{)}^{2/\bar{p}}\big{\\}}h_{j}$
$\displaystyle+\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}Lh_{j}^{(q-2)/{2q}}\sum_{u=0}^{1}\big{(}\int_{t_{j-1}}^{s}\big{\|}\sigma_{u}(r,X_{r}^{i,N},\mu_{r}^{X,N})-\sigma_{u}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}^{q}dr\big{)}^{1/{q}}$
$\displaystyle+\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}Lh_{j}^{(q-2)/{2q}}\sum_{u=0}^{1}\frac{1}{N}\sum_{k=1}^{N}\big{(}\int_{t_{j-1}}^{s}\big{\|}\sigma_{u}(r,X_{r}^{k,N},\mu_{r}^{X,N})-\sigma_{u}(t_{j-1},X_{j-1}^{k,N},\mu_{j-1}^{X,N})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}^{q}dr\big{)}^{1/{q}},$
(35)

for any $s\in[t_{j-1},t_{j}]$, $i\in\\{1,\ldots,N\\}$,
$j\in\\{1,\ldots,n_{h}\\}$ and $l\in\\{1,\ldots,m_{1}\\}$. Further, by using
Assumptions H–2 and H–3, Remark 2.3, Hölder’s inequality, Lemma 2.4 and $h\leq
1$, one obtains

$\displaystyle\big{\|}\sigma_{u}(r,X_{r}^{i,N},\mu_{r}^{X,N})-\sigma_{u}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}^{q}\leq
L\big{\\{}\big{\|}X_{r}^{i,N}-X_{j-1}^{i,N}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}$
$\displaystyle\quad+\big{\|}\mathcal{W}_{2}(\mu_{r}^{X,N},\mu_{j-1}^{X,N})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}+|r-t_{j-1}|\big{(}1+\big{\|}X^{i,N}_{j-1}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}+\big{\|}\mathcal{W}_{2}(\mu_{j-1}^{X,N},\delta_{0})\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}\big{)}\\}$
$\displaystyle\leq
2L\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{r}^{i,N}-X_{j-1}^{i,N}\big{\|}_{\mathscr{L}^{q}(\tilde{\Omega})}+|r-t_{j-1}|L\\{1+2C_{1}^{1/\bar{p}}(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}^{\bar{p}})^{1/\bar{p}}\\}$
$\displaystyle\leq\\{2L(C_{1}^{1/\bar{p}}+C_{2}^{1/\bar{p}})(1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}^{\bar{p}})^{1/\bar{p}}+L\\}h_{j}^{1/2},$
(36)

for all $i\in\\{1,\ldots,N\\}$, $j\in\\{1,\ldots,n_{h}\\}$, $u\in\\{0,1\\}$
and $r\in[t_{j-1},t_{j}]$. On substituting (36) in (35), one obtains

$\displaystyle\|\zeta^{\ell}(s)\|_{\mathscr{L}^{q}(\tilde{\Omega})}\leq\Big{\\{}L+2\bar{L}^{2}+4L^{2}\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}+\Big{(}2(L+2\bar{L}^{2})C_{1}^{1/\bar{p}}+8L^{2}\big{(}C_{1}^{1/\bar{p}}+C_{2}^{1/\bar{p}}\big{)}$
$\displaystyle\times\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}+4LC_{2}^{2/\bar{p}}\Big{)}\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|_{{}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}^{\bar{p}}\big{)}^{2/\bar{p}}\Big{\\}}h_{j},$

which, substituting in (34), gives the first estimate. The second estimate can
be proved by similar arguments. ∎

######  Proposition 5.9.

Let Assumptions H–1 with $\bar{p}\geq 4$ and set $q=\bar{p}/2\geq 2$, H–2 to
H–4 hold. Then, the drift-randomised Milstein scheme (13) is consistent of
order 1, i.e.,

$\displaystyle\big{\|}\mathcal{R}[X^{\varrho_{h}}]\big{\|}_{\mathscr{G}_{S,q}^{h}}\leq
C_{10}h,$

where

$\displaystyle C_{10}:=2L\big{(}1+2\sqrt{{q(q-1)}/{2}}\big{)}\big{\\{}L+$
$\displaystyle 2L(C_{1}^{1/\bar{p}}+C_{2}^{1/\bar{p}})$
$\displaystyle\times\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}\big{)}^{1/\bar{p}}\big{\\}}+C_{11}+C_{12}$

and the positive constants $C_{1}$, $C_{2}$, $C_{11}$ and $C_{12}$ appear in
Proposition 2.2, Lemma 2.4, Corollary 5.7 and Lemma 5.8, respectively.

###### Proof.

Let us recall the residual $\mathcal{R}_{j}^{i,N}[X^{\varrho_{h}}]$ from
(5.1), the interacting particle system from (2) and $\Gamma^{h}_{j}$ from
(5.1) to write $\mathcal{R}_{0}^{i,N}[X^{\varrho_{h}}]=0$ and

$\displaystyle\mathcal{R}_{j}^{i,N}[X^{\varrho_{h}}]$
$\displaystyle=X_{j}^{i,N}-X_{j-1}^{i,N}-\Gamma^{h}_{j}(X_{j-1}^{i,N},\mu_{j-1}^{X,N},X_{j,\eta}^{i,N},\mu_{j,\eta}^{X,N},\eta_{j})$
$\displaystyle=$
$\displaystyle\int_{t_{j-1}}^{t_{j}}\big{[}b(s,X_{s}^{i,N},\mu_{s}^{X,N})-b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N},\mu_{j,\eta}^{X,N})\big{]}ds$
$\displaystyle+\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]}dW^{i,\ell}_{s}$
$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\big{[}\sigma_{0}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]}dW_{s}^{0,\ell},$
(37)

for all $i\in\\{1,\ldots,N\\}$ and $j\in\\{1,\ldots,n_{h}\\}$. For any
$\bar{n}\in\\{1,\ldots,n_{h}\\}$, from (27),

$\Theta_{\bar{n},\eta}^{h}\big{[}b\big{]}=\displaystyle\sum_{j=1}^{\bar{n}}h_{j}b(t_{j-1}+\eta_{j}h_{j},X^{i,N}_{t_{j-1}+\eta_{j}h_{j}},\mu_{t_{j-1}+\eta_{j}h_{j}}^{X,N}),$

which along with Assumption H–2 yields

$\displaystyle\sum_{j=1}^{\bar{n}}$
$\displaystyle\int_{t_{j-1}}^{t_{j}}\big{[}b(s,X_{s}^{i,N},\mu_{s}^{X,N})-b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N},\mu_{j,\eta}^{X,N})\big{]}ds$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{\bar{n}}\int_{t_{j-1}}^{t_{j}}\big{[}b(s,X_{s}^{i,N},\mu_{s}^{X,N})-b(t_{j-1}+\eta_{j}h_{j},X^{i,N}_{t_{j-1}+\eta_{j}h_{j}},\mu_{t_{j-1}+\eta_{j}h_{j}}^{X,N})\big{]}ds$
$\displaystyle+\sum_{j=1}^{\bar{n}}\int_{t_{j-1}}^{t_{j}}\big{[}b(t_{j-1}+\eta_{j}h_{j},X^{i,N}_{t_{j-1}+\eta_{j}h_{j}},\mu_{t_{j-1}+\eta_{j}h_{j}}^{X,N})-b(t_{j-1}+\eta_{j}h_{j},X_{j,\eta}^{i,N},\mu_{j,\eta}^{X,N})\big{]}ds$
$\displaystyle\leq$
$\displaystyle\int_{0}^{t_{\bar{n}}}b(s,X_{s}^{i,N},\mu_{s}^{X,N})ds-\Theta_{\bar{n},\eta}^{h}\big{[}b\big{]}$
$\displaystyle+L\sum_{j=1}^{\bar{n}}h_{j}\big{[}\big{|}X^{i,N}_{t_{j-1}+\eta_{j}h_{j}}-X_{j,\eta}^{i,N}\big{|}+\big{|}\mathcal{W}_{2}(\mu_{t_{j-1}+\eta_{j}h_{j}}^{X,N},\mu_{j,\eta}^{X,N})\big{|}\big{]},$

and thus from (37) one obtains

$\displaystyle\big{\|}$
$\displaystyle\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\mathcal{R}^{i,N}_{j}[X^{\varrho_{h}}]\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\Theta^{h}_{\bar{n},\eta}[b]-\int_{0}^{t_{\bar{n}}}b(s,X_{s}^{i,N},\mu_{s}^{X,N})ds\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+L\sum_{j=1}^{n_{h}}h_{j}\big{[}\big{\|}X^{i,N}_{t_{j-1}+\eta_{j}h_{j}}-X_{j,\eta}^{i,N}\big{\|}_{\mathscr{L}^{q}(\Omega)}+\big{\|}\mathcal{W}_{2}(\mu_{t_{j-1}+\eta_{j}h_{j}}^{X,N},\mu_{j,\eta}^{X,N})\big{\|}_{\mathscr{L}^{q}(\Omega)}\big{]}$
$\displaystyle+\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{1}}\int_{t_{j-1}}^{t_{j}}\big{[}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{1}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]}dW^{i,\ell}_{s}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}$
$\displaystyle+\big{\|}\max_{\bar{n}\in\\{1,\ldots,n_{h}\\}}\big{|}\sum_{j=1}^{\bar{n}}\sum_{\ell=1}^{m_{0}}\int_{t_{j-1}}^{t_{j}}\big{[}\sigma_{0}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\tilde{\sigma}_{0}^{\ell}(s,t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{]}dW_{s}^{0,\ell}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)},$
(38)

for all $i\in\\{1,\ldots,N\\}$. Furthermore, Theorem 7.1 in [42] yields

$\displaystyle\big{\|}X^{i,N}_{t_{j-1}+\eta_{j}h_{j}}-X_{j,\eta}^{i,N}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq
h_{j}^{(q-1)/q}\Big{(}\int_{t_{j-1}}^{t_{j}}\big{\|}b(s,X_{s}^{i,N},\mu_{s}^{X,N})-b(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{\|}_{\mathscr{L}^{q}(\Omega)}^{q}ds\Big{)}^{1/q}$
$\displaystyle\quad+\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}h_{j}^{(q-2)/2q}\sum_{u=0}^{1}\Big{(}\int_{t_{j-1}}^{t_{j}}\big{\|}\sigma_{u}(s,X_{s}^{i,N},\mu_{s}^{X,N})-\sigma_{u}(t_{j-1},X_{j-1}^{i,N},\mu_{j-1}^{X,N})\big{\|}_{\mathscr{L}^{q}(\Omega)}^{q}ds\Big{)}^{1/q},$

which on using Equations (33) and (36),

$\displaystyle\big{\|}X^{i,N}_{t_{j-1}+\eta_{j}h_{j}}-X_{j,\eta}^{i,N}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq$
$\displaystyle\Big{(}1+2\Big{(}\frac{q(q-1)}{2}\Big{)}^{1/2}\Big{)}\Big{\\{}L+2L(C_{1}^{1/\bar{p}}+C_{2}^{1/\bar{p}})$
$\displaystyle\times\big{(}1+\max_{i\in\\{1,\ldots,N\\}}\|X_{0}^{i}\|^{\bar{p}}_{{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}}\big{)}^{1/\bar{p}}\Big{\\}}h_{j},$
(39)

for all $i\in\\{1,\ldots,N\\}$ and $j\in\\{1,\ldots,n_{h}\\}$. Also,

$\displaystyle\big{\|}\mathcal{W}_{2}(\mu_{t_{j-1}+\eta_{j}h_{j}}^{X,N},\mu_{j,\eta}^{X,N})\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq\max_{i\in\\{1,\ldots,N\\}}\big{\|}X^{i,N}_{t_{j-1}+\eta_{j}h_{j}}-X_{j,\eta}^{i,N}\big{\|}_{\mathscr{L}^{q}(\Omega)},$
(40)

for all $j\in\\{1,\ldots,n_{h}\\}$. The proof is completed by substituting
(5.3) and (40) in (5.3) and using Corollary 5.7 and Lemma 5.8. ∎

###  5.4 Rate of Convergence of the Scheme

After proving bistability and consistency of the scheme (13), the main result,
Theorem 3.3, follows immediately.

###### Proof of Theorem 3.3.

Set $Y^{h}$ in Proposition 5.5 as
$Y^{h}=X^{\varrho_{h}}=\\{X^{i,N}_{j}\\}_{j\in\\{0,\ldots,n_{h}\\}}\in{\mathscr{G}_{q}^{h}}$
(as introduced in Section 5.3) to obtain,

$\displaystyle\|X^{\varrho_{h}}-X^{h}\|_{\mathscr{G}_{q}^{h}}$
$\displaystyle=\max_{i\in\\{1,\ldots,N\\}}\big{\|}\max_{j\in\\{0,1,\ldots,n_{h}\\}}\big{|}X_{j}^{i,N}-X_{j}^{i,N,h}\big{|}\big{\|}_{\mathscr{L}^{q}(\Omega)}\leq
C_{7}\big{\|}\mathcal{R}[X^{\varrho_{h}}]\big{\|}_{\mathscr{G}_{S,q}^{h}},$

which on using Proposition 5.9 completes the proof. ∎

##  Appendix A Auxiliary results

The following lemma is the discrete Burkholder–Davis–Gundy inequality, see
[11].

######  Lemma A.1 (Discrete Burkholder–Davis–Gundy Inequality).

Let $\\{M_{n}\\}_{n\in\bar{\mathbb{N}}}$ be a discrete-time martingale on a
probability space $(\Omega,\mathscr{F},\mathbb{P})$ with respect to the
filtration $\\{\mathscr{F}_{n}\\}_{n\in\bar{\mathbb{N}}}$. Then, there exist
constants $\underline{C}_{p},\bar{C}_{p}>0$ such that,

$\displaystyle\underline{C}_{p}\|[M]_{n}^{1/2}\|_{\mathscr{L}^{p}(\Omega)}\leq\|\max_{j\in\\{0,1,\cdots,n\\}}|M_{j}|\|_{\mathscr{L}^{p}(\Omega)}\leq\bar{C}_{p}\|[M]_{n}^{1/2}\|_{\mathscr{L}^{p}(\Omega)},$

for any $p\geq 1$ where $\\{[M]_{n}\\}_{n\in\bar{\mathbb{N}}}$ is the
quadratic variation process of $\\{M_{n}\\}_{n\in\bar{\mathbb{N}}}$.

The discrete version of the Grönwall’s inequality is given below, see
Proposition 4.1 in [26].

######  Lemma A.2 (Discrete Grönwall’s Inequality).

Let $\\{y_{n}\\}_{n\in\mathbb{N}}$ and $\\{z_{n}\\}_{n\in\mathbb{N}}$ be
sequences of non-negative real numbers satisfying $y_{n}\leq
a+\displaystyle\sum_{j=1}^{n-1}z_{j}y_{j}$ for all $n\in\mathbb{N}$ where
$a>0$ is a constant. Then, $y_{n}\leq
a\exp\big{(}\displaystyle\sum_{j=1}^{n-1}z_{j}\big{)}$ for all
$n\in\mathbb{N}$.

###  A.1 Proof of Proposition 2.2

Let $t^{\prime}\in[0,T]$. Using (2) and Hölder’s inequality, we get

$\displaystyle\sup_{t\in[0,t^{\prime}]}|X_{t}|^{\bar{p}}$ $\displaystyle\leq
4^{\bar{p}-1}\Big{\\{}|X_{0}|^{\bar{p}}+T^{\bar{p}-1}\int_{0}^{t^{\prime}}|b(s,X_{s},\mathcal{L}^{1}(X_{s}))|^{\bar{p}}ds$
$\displaystyle+\sup_{t\in[0,t^{\prime}]}\big{|}\int_{0}^{t}\sigma_{1}(s,X_{s},\mathcal{L}^{1}(X_{s}))dW_{s}\big{|}^{\bar{p}}+\sup_{t\in[0,t^{\prime}]}\big{|}\int_{0}^{t}\sigma_{0}(s,X_{s},\mathcal{L}^{1}(X_{s}))dW_{s}^{0}\big{|}^{\bar{p}}\Big{\\}},$

which on using Burkholder–Gundy–Davis inequality (Theorem 7.2 in [42]) and
Remark 2.1 yields

$\displaystyle\big{\|}\sup_{t\in[0,t^{\prime}]}|X_{t}|\big{\|}^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}\leq
4^{\bar{p}-1}\Big{\\{}\big{\|}X_{0}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}+T^{\bar{p}-1}\int_{0}^{t^{\prime}}\big{\|}b(s,X_{s},\mathcal{L}^{1}(X_{s}))\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}ds$
$\displaystyle\quad+\Big{(}\frac{\bar{p}^{3}}{2(\bar{p}-1)}\Big{)}^{\bar{p}/2}T^{(\bar{p}-2)/2}\sum_{u=0}^{1}\int_{0}^{t^{\prime}}\big{\|}\sigma_{u}(s,X_{s},\mathcal{L}^{1}(X_{s}))\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}ds\Big{\\}}$
$\displaystyle\leq
4^{\bar{p}-1}\big{\|}X_{0}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}+4^{\bar{p}-1}3^{\bar{p}-1}\bar{L}^{\bar{p}}\Big{(}T^{%̄
̉\bar{p}-1}+2\Big{(}\frac{\bar{p}^{3}}{2(\bar{p}-1)}\Big{)}^{\bar{p}/2}T^{(\bar{p}-2)/2}\Big{)}$
$\displaystyle\qquad\times\int_{0}^{t^{\prime}}\big{\\{}1+\big{\|}X_{s}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}+\big{\|}\mathcal{W}_{2}(\mathcal{L}^{1}(X_{s}),\delta_{0})\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\big{\\}}ds$
$\displaystyle\leq
4^{\bar{p}-1}\big{\|}X_{0}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}+12^{\bar{p}-1}\bar{L}^{\bar{p}}\Big{(}T^{\bar{p}}+2\Big{(}\frac{\bar{p}^{3}}{2(\bar{p}-1)}\Big{)}^{\bar{p}/2}T^{\bar{p}/2}\Big{)}$
$\displaystyle\qquad+12^{\bar{p}-1}\bar{L}^{\bar{p}}\Big{(}T^{%̄
̉\bar{p}-1}+2\Big{(}\frac{\bar{p}^{3}}{2(\bar{p}-1)}\Big{)}^{\bar{p}/2}T^{(\bar{p}-2)/2}\Big{)}2\int_{0}^{t^{\prime}}\big{\|}\sup_{r\in[0,s]}|X_{r}|\big{\|}^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}ds$

for all $t^{\prime}\in[0,T]$. Then, Grönwall’s inequality gives

$\displaystyle\big{\|}\sup_{t\in[0,T]}|X_{t}|\big{\|}^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}\leq$
$\displaystyle\Big{(}4^{\bar{p}-1}\big{\|}X_{0}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}+12^{\bar{p}-1}\bar{L}^{\bar{p}}\Big{(}T^{\bar{p}}%̄+2\Big{(}\frac{\bar{p}^{3}}{2(\bar{p}-1)}\Big{)}^{\bar{p}/2}T^{\bar{p}/2}\Big{)}\Big{)}$
$\displaystyle\times\exp\Big{(}12^{\bar{p}-1}\bar{L}^{\bar{p}}\Big{(}T^{%̄̉\bar{p}}+2\Big{(}\frac{\bar{p}^{3}}{2(\bar{p}-1)}\Big{)}^{\bar{p}/2}T^{\bar{p}/2}\Big{)}2\Big{)},$

which completes the proof. $\square$

###  A.2 Proof of Lemma 2.4

Recall Equation (2) and use Theorem 7.1 in [42] along with Remark 2.1 to get
the following,

$\displaystyle\big{\|}X_{t}^{i,N}-X_{t^{\prime}}^{i,N}\big{\|}^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}\leq
3^{\bar{p}-1}(t-t^{\prime})^{\bar{p}-1}\int_{t^{\prime}}^{t}\big{\|}b(s,X_{s}^{i,N},\mu_{s}^{X,N})\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}ds$
$\displaystyle+3^{\bar{p}-1}\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{\bar{p}/2}(t-t^{\prime})^{(\bar{p}-2)/2}\sum_{u=0}^{1}\int_{t^{\prime}}^{t}\big{\|}\sigma_{u}(s,X_{s}^{i,N},\mu_{s}^{X,N})\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}ds$
$\displaystyle\leq$ $\displaystyle
9^{\bar{p}-1}L^{\bar{p}}\Big{(}(t-t^{\prime})^{\bar{p}-1}+2\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{\bar{p}/2}(t-t^{\prime})^{(\bar{p}-2)/2}\Big{)}$
$\displaystyle\qquad\times\int_{t^{\prime}}^{t}\big{(}1+\big{\|}X_{s}^{i,N}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}+\big{\|}\mathcal{W}_{2}(\mu_{s}^{X,N},\delta_{0})\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\big{)}ds,$

for all $t>t^{\prime}\in[0,T]$ and $i\in\\{1,\ldots,N\\}$, which due to Remark
2.3 yields,

$\displaystyle\big{\|}$ $\displaystyle
X_{t}^{i,N}-X_{t^{\prime}}^{i,N}\big{\|}^{\bar{p}}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}\leq
9^{\bar{p}-1}L^{\bar{p}}\Big{(}(t-t^{\prime})^{\bar{p}-1}+2\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{\bar{p}/2}(t-t^{\prime})^{(\bar{p}-2)/2}\Big{)}$
$\displaystyle\times\int_{t^{\prime}}^{t}\big{\\{}1+2\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{s}^{i,N}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\big{\\}}ds$
$\displaystyle\leq$ $\displaystyle
9^{\bar{p}-1}L^{\bar{p}}\Big{(}(t-t^{\prime})^{\bar{p}/2}+2\Big{(}\frac{\bar{p}(\bar{p}-1)}{2}\Big{)}^{\bar{p}/2}\Big{)}(t-t^{\prime})^{\bar{p}/2}\Big{(}1+2C_{1}\Big{(}1+\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{0}^{i,N}\big{\|}_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}}\Big{)}\Big{)}$

and thus the proof is completed. $\square$

###  A.3 Proof of Proposition 2.5

By using (2) and (2) along with Hölder’s inequality, one obtains

$\displaystyle\sup_{t\in[0,T]}$ $\displaystyle|X_{t}^{i}-X_{t}^{i,N}|^{2}\leq
3\Big{\\{}T\int_{0}^{T}|b(s,X_{s}^{i},\mathcal{L}^{1}(X_{s}^{1}))-b(s,X_{s}^{i,N},\mu_{s}^{X,N})|^{2}ds$
$\displaystyle+\sup_{t\in[0,T]}\big{|}\int_{0}^{t}\big{(}\sigma_{1}(s,X_{s}^{i},\mathcal{L}^{1}(X_{s}^{1}))-\sigma_{1}(s,X_{s}^{i,N},\mu_{s}^{X,N})\big{)}dW_{s}\big{|}^{2}$
$\displaystyle+\sup_{t\in[0,T]}\big{|}\int_{0}^{t}\big{(}\sigma_{0}(s,X_{s}^{i},\mathcal{L}^{1}(X_{s}))-\sigma_{0}(s,X_{s}^{i,N},\mu_{s}^{X,N})\big{)}dW_{s}^{0}\big{|}^{2}\Big{\\}},$

which on using the martingale inequality (Theorem 7.2 in [42]) and Assumption
H–2 yields

$\displaystyle\big{\|}$
$\displaystyle\sup_{t\in[0,t^{\prime}]}|X_{t}^{i}-X_{t}^{i,N}|\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}\leq
3T\int_{0}^{t^{\prime}}\big{\|}b(s,X_{s}^{i},\mathcal{L}^{1}(X_{s}^{1}))-b(s,X_{s}^{i,N},\mu_{s}^{X,N})\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}ds$
$\displaystyle\quad+12\sum_{u=0}^{1}\int_{0}^{t^{\prime}}\big{\|}\sigma_{u}(s,X_{s}^{i},\mathcal{L}^{1}(X_{s}^{1}))-\sigma_{u}(s,X_{s}^{i,N},\mu_{s}^{X,N})\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}ds$
$\displaystyle\leq
2(3T+24)L^{2}\int_{0}^{t^{\prime}}\big{\\{}\big{\|}X_{s}^{i}-X_{s}^{i,N}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}+\big{\|}\mathcal{W}_{2}(\mathcal{L}^{1}(X_{s}^{1}),\mu_{s}^{X,N}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}\big{\\}}ds$
$\displaystyle\leq
2(3T+24)L^{2}\int_{0}^{t^{\prime}}\big{\\{}\big{\|}X_{s}^{i}-X_{s}^{i,N}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}+2\big{\|}\mathcal{W}_{2}(\mathcal{L}^{1}(X_{s}^{1}),\mu_{s}^{X}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}+2\big{\|}\mathcal{W}_{2}(\mathcal{\mu}_{s}^{X},\mu_{s}^{X,N}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}\big{\\}}ds$
$\displaystyle\leq
2(3T+24)L^{2}\int_{0}^{t^{\prime}}\big{\\{}3\max_{i\in\\{1,\ldots,N\\}}\big{\|}X_{s}^{i}-X_{s}^{i,N}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}+2\big{\|}\mathcal{W}_{2}(\mathcal{L}^{1}(X_{s}^{1}),\mu_{s}^{X}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}\big{\\}}ds$

for any $t^{\prime}\in[0,T]$ and $N\in\mathbb{N}$. An application of
Grönwall’s inequality yields

$\displaystyle\max_{i\in\\{1,\ldots,N\\}}\big{\|}\sup_{t\in[0,T]}|X_{t}^{i}-X_{t}^{i,N}|\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}\leq
4(3T+24)L^{2}e^{6(3T+24)TL^{2}}\int_{0}^{T}\big{\|}\mathcal{W}_{2}(\mathcal{L}^{1}(X_{s}^{1}),\mu_{s}^{X}\big{\|}_{\mathscr{L}^{2}(\tilde{\Omega})}^{2}ds,$

for any $N\in\mathbb{N}$. The proof is completed by using (7). $\square$

###  A.4 Lions Derivative and Particle Projection Function Inequalities

Given a function $f:\mathscr{P}_{2}(\mathbb{R}^{d})\to\mathbb{R}$ and
$\nu_{0}\in\mathscr{P}_{2}(\mathbb{R}^{d})$, we say that $f$ is Lions
differentiable at $\nu_{0}$ if we can find an atomless, Polish probability
space $(\tilde{\Omega},\tilde{\mathcal{F}},\tilde{P})$ and a random variable
$Y_{0}\in\mathscr{L}^{2}(\tilde{\Omega})$ with law
$\mathcal{L}^{1}({Y_{0}}):=\tilde{P}\circ Y_{0}^{-1}=\nu_{0}$ such that the
“lift” function $F:\mathscr{L}^{2}(\tilde{\Omega})\to\mathbb{R}$ defined by
$F(Z):=f(\mathcal{L}^{1}(Z))$ has Fréchet derivative $F^{\prime}[Y_{0}]$ at
$Y_{0}\in\mathscr{L}^{2}(\tilde{\Omega})$. By Riesz representation theorem, we
can find $DF(Y_{0})\in\mathscr{L}^{2}(\tilde{\Omega})$ such that
$F^{\prime}[Y_{0}](Z)=\tilde{E}\langle DF[Y_{0}],Z\rangle$ for all
$Z\in\mathscr{L}^{2}(\tilde{\Omega};\mathbb{R}^{d})$. Further, Theorem 6.5
(structure of the gradient) in [12] guarantees the existence of a function
$\partial_{\mu}f(\nu_{0}):\mathbb{R}^{d}\to\mathbb{R}^{d}$, independent of the
choice of the probability space
$(\tilde{\Omega},\tilde{\mathcal{F}},\tilde{P})$ and the random variable
$Y_{0}$, satisfying
$\displaystyle\int_{\mathbb{R}^{d}}|\partial_{\mu}f(\nu_{0})(x)|^{2}\nu_{0}(dx)<\infty$
such that $DF(Y_{0})=\partial_{\mu}f(\nu_{0})(Y_{0})$. Then, we call
$\partial_{\mu}f(\nu_{0})$ as Lions’ derivative of $f$ at
$\nu_{0}=\mathcal{L}^{1}({Y_{0}})$.

By using similar arguments used in Lemma 4.2 in [36], one can prove the
following lemma.

######  Lemma A.3.

Let $g$ be a real valued function defined on
$[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})$ such that
derivative with respect to the state variable
$\partial_{x}g(\cdot,\cdot,\cdot)$ and with respect to the measure variable
$\partial_{\mu}g(\cdot,\cdot,\cdot,\cdot)$ satisfy Lipschitz condition
uniformly in time, i.e., there exists a constant $L>0$ such that,

$\displaystyle|\partial_{x}g(t,x,\mu)-\partial_{x}g(t,\bar{x},\bar{\mu})|$
$\displaystyle\leq L\\{|x-\bar{x}|+\mathcal{W}_{2}(\mu,\bar{\mu})\\},$
$\displaystyle|\partial_{\mu}g(t,x,\mu,y)-\partial_{\mu}g(t,\bar{x},\bar{\mu},\bar{y})|$
$\displaystyle\leq
L\\{|x-\bar{x}|+\mathcal{W}_{2}(\mu,\bar{\mu})+|y-\bar{y}|\\},$

for all $t\in[0,T]$, $x,\bar{x},y,\bar{y}\in\mathbb{R}^{d}$ and
$\mu,\bar{\mu}\in\mathscr{P}_{2}(\mathbb{R}^{d})$. Then, for all $t\in[0,T]$,
$x^{i}$, $\bar{x}^{i}\in\mathbb{R}^{d}$ and $i\in\\{1,\ldots,N\\}$, the
following holds,

$\displaystyle\big{|}g\big{(}t,x^{i},\frac{1}{N}\sum_{k=1}^{N}\delta_{x^{k}}\big{)}-g\big{(}t,\bar{x}^{i},\frac{1}{N}\sum_{k=1}^{N}\delta_{\bar{x}^{k}}\big{)}-\partial_{x}g\big{(}t,\bar{x}^{i},\frac{1}{N}\sum_{k=1}^{N}\delta_{\bar{x}^{k}}\big{)}(x^{i}-\bar{x}^{i})$
$\displaystyle-\frac{1}{N}\sum_{k=1}^{N}\partial_{\mu}g\big{(}t,\bar{x}^{i},\frac{1}{N}\sum_{k=1}^{N}\delta_{\bar{x}^{k}},\bar{x}^{k}\big{)}(x^{k}-\bar{x}^{k})\big{|}\leq\frac{3L}{2}|x^{i}-\bar{x}^{i}|^{2}+\frac{5L}{2}\frac{1}{N}\sum_{k=1}^{N}|x^{k}-\bar{x}^{k}|^{2}.$

[Acknowledgments] The authors would like to thank the anonymous referees whose
comments improved the quality of this paper.

G. dos Reis acknowledges support from the _Fundação para a Ci $\hat{e}$ncia e
a Tecnologia_ (Portuguese Foundation for Science and Technology) through the
project UIDB/00297/2020 and UIDP/00297/2020 (Center for Mathematics and
Applications).

## References

  * [1] Adams, D., Dos Reis, G., Ravaille, R., Salkeld, W. and Tugaut, J. (2022). Large Deviations and Exit-times for reflected McKean–Vlasov equations with self-stabilising terms and superlinear drifts. _Stochastic Process. Appl._ 146 264–310.
  * [2] Baladron, J., Fasoli, D., Faugeras, O. and Touboul, J. (2012). Mean-field description and propagation of chaos in networks of Hodgkin–Huxley and FitzHugh–Nagumo neurons. _J. Math. Neurosci._ 2(10)
  * [3] Bao, J., and Huang, X. (2022). Approximations of McKean–Vlasov stochastic differential equations with irregular coefficients. _J. Theoret. Probab._ 35 1187–1215.