∎

11institutetext: Tan Zhang 22institutetext: Department of Mathematics, The
University of Hong Kong, Pokfulam, Hong Kong.
22email<EMAIL_ADDRESS>33institutetext: Zhongjian Wang 44institutetext:
Division of Mathematical Sciences, School of Physical and Mathematical
Sciences, Nanyang Technological University, Singapore 637371.
44email<EMAIL_ADDRESS>55institutetext: Jack Xin 66institutetext:
Department of Mathematics, University of California at Irvine, Irvine, CA
92697, USA.
66email<EMAIL_ADDRESS>77institutetext: Zhiwen Zhang 88institutetext:
Department of Mathematics, The University of Hong Kong, Pokfulam, Hong Kong.
88email<EMAIL_ADDRESS>
# A structure-preserving scheme for computing effective diffusivity and
anomalous diffusion phenomena of random flows

Tan Zhang  Zhongjian Wang  Jack Xin  Zhiwen Zhang

(Received: date / Accepted: date)

###### Abstract

This paper aims to investigate the diffusion behavior of particles moving in
stochastic flows under a structure-preserving scheme. We compute the effective
diffusivity for normal diffusive random flows and establish the power law
between spatial and temporal variables for cases with anomalous diffusion
phenomena. From a Lagrangian approach, we separate the corresponding
stochastic differential equations (SDEs) into sub-problems and construct a
one-step structure-preserving method to solve them. Then by modified equation
systems, the convergence analysis in calculating the effective diffusivity is
provided and compared between the structure-preserving scheme and the Euler-
Maruyama scheme. Also, we provide the error estimate for the structure-
preserving scheme in calculating the power law for a series of super-diffusive
random flows. Finally, we calculate the effective diffusivity and anomalous
diffusion phenomena for a series of 2D and 3D random fields.

###### Keywords:

Convection-enhanced diffusionrandom flowsstructure-preserving schemecorrector
problemergodic theoryMarkov processbackward error analysis.

###### MSC:

37M2560J6060H3565P10, 65M7576M50.

##  1 Introduction

The study of diffusion enhancement phenomena in fluid advection has been
ongoing since the groundbreaking research conducted by Sir G. TaylorTaylor1922
. A key challenge lies in accurately characterizing and quantifying the
effective diffusion at large scales in complex fluid flows, as this holds
significant theoretical and practical significance. The applications of this
research span various scientific and technological domains, including
atmospheric science, oceanography, chemical engineering, and combustion.

In this work, we aim to compute the effective diffusivity of the passive
tracer models and investigate anomalous diffusion phenomena in chaotic and
random flows. For the passive tracer models in random flows, it can be
described by the following stochastic differential equations (SDEs).

$d\bm{x}(t)=\bm{v}(\bm{x},t,\omega)dt+\sigma d\bm{w}(t),\
\bm{x}\in\mathbb{R}^{d},\ t\geq 0,$ (1)

where $\bm{x}(t)=(x_{1}(t),...,x_{d}(t))^{T}\in\mathbb{R}^{d}$ is the spatial
position of the passive tracer particle, $\bm{v}(\bm{x},t,\omega)$ is the
incompressible velocity field, $\sigma>0$ is the molecular diffusivity, and
$\bm{w}(t)$ is the standard $d$-dimensional Brownian motion. The velocity
field $\bm{v}(\bm{x},t,\omega)$ is a mean zero, stationary, ergodic random
vector field over a certain probability space. Here $\omega$ is an element of
the probability space describing all possible environments which is exclusive
of the Brownian motion, $\bm{w}(t)$. We assume that the realizations of
$\bm{v}(\bm{x},t,\omega)$ are almost surely locally Lipschitz in $\bm{x}$.
This is to guarantee the existence of the solution of equation (1). In
practice, for a given energy spectrum density, we follow the generation method
in Kraichnan1970 Majda1999  to model the random field
$\bm{v}(\bm{x},t,\omega)$.

One of the primary challenges associated with the SDEs (1) involves
investigating the behavior of the particles $\bm{x}(t)$ in (1) over an
extended time and on a large scale. For the rescaled process,
$\epsilon\bm{x}(t/\epsilon^{2})$, it is crucial to verify whether it will
converge in distribution to a new Brownian motion $\tilde{\bm{w}}(t)$ with a
covariance matrix $\mathbf{D}^{E}$. If so, then $\mathbf{D}^{E}$ is denoted as
the effective diffusivity matrix. This problem can be viewed as the
homogenization of the passive tracer model and is highly nontrivial. Extensive
research has been conducted to compute the effective diffusivity matrix
$\mathbf{D}^{E}$ under various conditions on the flows. In situations where
there are spatial-temporal periodic fields or random fields with short-range
correlations, the homogenization theory can be applied to calculate the
effective diffusivity matrix $\mathbf{D}^{E}$ Panpanicolaou2011 ,Garnier1997
,Jikov1994 ,Pavliotis2008 . This enables us to determine the behavior of the
diffusion process in such systems.

In the case of random flows characterized by long-range correlations, studying
the behavior of particle motion over an extended time and on a large scale
poses significant challenges due to its intricate and general complexity. The
presence of various forms of anomalous diffusion, including super-diffusion
and sub-diffusion, further adds to the difficulty of the analysis.

Multiple theoretical studies have been conducted on the homogenization of
time-dependent random flows Carmona1997 ,Fannjiang1999 ,Fannjiang1996
,Komorowski2001 ,Landim1998 . One notable example of these results is the work
by the authors in Carmona1997 , where they successfully demonstrated the
existence of the effective diffusivity for the 2D incompressible time-
dependent Gaussian field. In Fannjiang1999 ,Fannjiang1996 , the authors
demonstrated certain necessary conditions in turbulent flows that should be
satisfied for the behavior of advection-diffusion over an extended time to
exhibit normal diffusion behavior. For random flows that exhibit normal
diffusion behavior, we have proposed a structure-preserving method to compute
the effective diffusivity in the Lagrangian approach Junlong2021 . In our
previous studies, the focus has primarily been on random flows that exhibit
diffusive behavior, rather than on random fields that display abnormal
diffusion. Moreover, the numerical analysis for calculating the effective
diffusivity is unexplored when the diffusion term’s coefficient $\sigma$
vanishes. This motivates us to provide the following contributions in this
paper:

  * $\bullet$

Provide the convergence analysis of our structure-preserving method in
computing effective diffusivity when the diffusion coefficient $\sigma$
vanishes;

  * $\bullet$

Compare the error estimate of our structure-preserving method with the
traditional Euler-Maruyama scheme;

  * $\bullet$

Investigate the anomalous diffusion phenomena and verify the accuracy of our
structure-preserving method in a series of random flows.

The whole procedure is as follows. We first apply an operator splitting method
to separate the SDE system (1) into two sub-problems. In our structure-
preserving method, we adopt an implicit structure-preserving scheme to solve
the advection term and adopt the Euler-Maruyama scheme to solve the diffusion
term. From the perspective of the modified equation, we obtain the continuous
modified equation systems (54) which correspond to our structure-preserving
method and Euler-Maruyama scheme respectively. Then by applying Ito’s formula,
we can analyze and obtain the error estimate of two modified equation systems
in computing effective diffusivity (Theorem 4.1). We have also analyzed that
the structure-preserving method’s modified system can maintain the diffusion
behavior of the original system better when $\sigma$ is relatively small,
regardless of whether it is sub-diffusive (Theorem 4.2) or normal diffusive
(Corollary 1). While the Euler-Maruyama scheme does not hold this property.
Meanwhile, for a group of random velocity fields that exhibit super-diffusion
as mentioned in Fannjiang1998 , we prove that our structure-preserving method
can maintain its original super-diffusion behavior (Theorem 4.3). Finally, we
present numerical results to verify the theorems, demonstrate the convection-
enhanced diffusion phenomenon, and compute the effective diffusivity and
anomalous diffusion phenomena in 2D and 3D random flows.

The rest of the paper is organized as follows. In Section 2, we introduce the
definition of effective diffusivity and briefly review some existing results
for diffusion in random flows. We also provide the sharp condition for random
fields to exhibit normal diffusive phenomena. In Section 3, we develop the
implicit structure-preserving method for computing the effective diffusivity
of passive tracer models in random flows. We also make a comparison of the
volume-preserving property between the structure-preserving method and the
Euler-Maruyama method. Moreover, we introduce the modified equation systems of
two numerical schemes. In Section 4, we provide the convergence analysis for
our structure-preserving method and Euler-Maruyama method based on the
modified equation system. Error analysis of the performance of the structure-
preserving scheme under random flows with super-diffusive behavior is also
carried out. In Section 5, we conduct numerical experiments and present
results to verify the error analyses we proved, demonstrate the convection-
enhanced diffusion phenomenon, and compute the effective diffusivity and
anomalous diffusion phenomena in 2D and 3D random flows.

##  2 Preliminaries

To ensure the self-contained nature of this paper, we provide a concise
overview of the current research in the area of advection-enhanced diffusion
in stochastic turbulent flows and the definition of effective diffusivity of
the associate passive tracer models (1). Some results for the general cases we
put in Appendix A, and in the main body we use these results in a more
specific system (Subsection 2.1). We will also go through the criterion to
determine whether a stochastic turbulent flow will perform normal diffusive
behavior or not (Subsection 2.2).

###  2.1 Continuous corrector problem and the definition of effective
diffusivity

To study the effective diffusivity of the stochastic flows connected to (1),
we first present a concise overview of the current findings regarding
convection-enhanced diffusion in random flows.

Let $(\mathcal{X},\mathcal{H},P_{0})$ be a probability space. Here
$\mathcal{H}$ is a $\sigma$-algebra over $\mathcal{X}$, $P_{0}$ is a
probability measure. Let $\tau_{\bm{x}}$, $\bm{x}\in\mathbb{R}^{d}$, be an
almost surely continuous, jointly measurable group of measure-preserving
transformation on $\mathcal{X}$. The properties of this system are illustrated
in Appendix A and here we follow its defined notations. Then let $\Omega$ be
the space of $\mathcal{X}$-valued continuous function
$C([0,\infty);\mathcal{X})$ and let $\mathcal{B}$ be its Borel
$\sigma$-algebra. Then we denote $\omega(\cdot+t),t\geq 0$ as the standard
shift operator on the path space $(\Omega,\mathcal{B})$.

From a Lagrangian particle approach, we can define the following environment
process $\eta:[0,\infty)\times\Omega\to\mathcal{X}$ at time $t$ as

$\eta(t)=\tau_{\bm{x}_{t}^{\omega}}\omega(t),\ \text{where}\
\eta(0)=\omega(0).$ (2)

Here $\bm{x}_{t}^{\omega}\in\mathbb{R}^{d}$ is the position of the particle at
time $t$ with the realization of the environment $\omega$. The random field
$\bm{v}(\bm{x},t,\omega)$ is stationary provided $\tau_{\bm{x}}$ preserves the
probability measure $P_{0}$, for every $\bm{x}\in\mathbb{R}^{d}$. And the
ergodicity of $\bm{v}(\bm{x},t,\omega)$ is also provided by the transformation
$\tau_{\bm{x}},\bm{x}\in\mathbb{R}^{d}$ is ergodic. Furthermore, inconvenience
to the subsequent analysis of the effective diffusivity, we require that the
incompressible random velocity field $\bm{v}(\bm{x},t,\omega)$ satisfies the
following assumptions. The incompressible random field
$\bm{v}(\bm{x},t,\omega)$ is jointly continuous in $(\bm{x},t)$, bounded and
with finite second moments.

$\bm{v}(\bm{x},t,\omega)$ is locally Lipschitz in the spatial variable and has
certain regularity, i.e. $\bm{v}\in(C_{b}^{m}(\mathcal{X}))^{d}$ for some
$m\geq 1$.

By denoting M as the expectation corresponding to the Brownian motion term in
the SDEs, the environment process will generate a semigroup of transformation

$\bm{S}^{t}f(\tilde{\bm{x}})=\textbf{ME}_{\tilde{\bm{x}}}f(\eta(t)),t\geq 0,\
f\in L^{\infty}(\mathbb{R}^{d}),$ (3)

For clarity, we denote $\mathbb{E}$ as the total expectation with respect to
all randomness, i.e., $\mathbb{E}=\mathbf{ME}$, in the remaining part of this
paper. From Proposition 3 and Proposition 6 in Fannjiang1999 , we know that
$\bm{S}^{t}$ is a strongly continuous, measure-preserving, Markov semigroup of
contraction on $L^{2}(\mathbb{R}^{d})$. Then we let $\mathcal{L}_{0}$ denote
the generator of the semigroup $\bm{S}^{t},t\geq 0$, i.e.

$\mathcal{L}_{0}f=\frac{\sigma^{2}}{2}\Delta f+\bm{v}\cdot\nabla f.$ (4)

$\mathcal{L}_{0}$ is a concretization of generator $L$ which is mentioned in
Appendix A. So $\bm{S}^{t}$ also holds the spectral gap condition (128), i.e.

$||\bm{S}^{t}f||_{L^{2}(\mathbb{R}^{d})}\leq\exp(-2c_{1}t)||f||_{L^{2}(\mathbb{R}^{d})},\
c_{1}>0,\ f\in L_{0}^{2}(\mathbb{R}^{d}).$ (5)

The proof of (5) can be seen from Proposition 6 in Fannjiang1999 . Then we can
define $\bm{G}=\int_{0}^{\infty}\bm{S}^{t}\bm{v}dt$ and it satisfies the
following continuous corrector problem,

$\mathcal{L}_{0}\bm{G}=-\bm{v}.$ (6)

By solving the continuous corrector problem (6), we can define the effective
diffusivity.

######  Proposition 1 (Fannjiang1999 , Lemma 1)

Let $\bm{x}(t)$ be the solution to the passive tracer model (1) and
$\bm{x}_{\epsilon}(t)\equiv\epsilon\bm{x}(t/\epsilon^{2})$. For any unit
vector $\bm{e}\in\mathbb{R}^{d}$, let $\bm{G}_{\bm{e}}$ denote the projection
of the vector solution $\bm{G}$ in (6) along the direction $\bm{e}$. Then the
law of the process of the projection of $\bm{x}_{\epsilon}(t)$ along direction
$\bm{e}$, i.e. $\bm{x}_{\epsilon}(t)\cdot\bm{e}$ will weakly converge to a
Brownian motion with a diffusion coefficient which is given by

$\bm{e}^{T}\mathbf{D}^{E}\bm{e}=D_{0}+(-\mathcal{L}_{0}\bm{G}_{\bm{e}},\bm{G}_{\bm{e}})_{L^{2}(\mathbb{R}^{d})},$
(7)

where $\mathbf{D}^{E}$ is the effective diffusivity associated with (1) and
$D_{0}=\sigma^{2}/2$.

We can see that the correction to $D_{0}$ in (7) is non-negative definite and
for any unit vector $\bm{e}\in\mathbb{R}^{d}$,
$\bm{e}^{T}\mathbf{D}^{E}\bm{e}\geq D_{0}$. This is called convection-enhanced
diffusion. The effective diffusivity fundamentally captures the long-time,
large-scale behaviors of passive tracer models. In practice, as $D_{0}$
approaches extremely small values, the solutions of (6) exhibit sharp
gradients, requiring a significant number of Fourier modes to resolve. Also
here we consider the flow with unbounded domain setting. These characteristics
render the Eulerian approach computationally demanding. Alternatively, the
effective diffusivity matrix can be computed using the Lagrangian approach,
defined by,

$D^{E}_{ij}=\lim_{t\to\infty}\frac{\big{\langle}\big{(}x_{i}(t)-x_{i}(0)\big{)}\big{(}x_{j}(t)-x_{j}(0)\big{)}\big{\rangle}}{2t},\
1\leq i,j\leq d,$ (8)

where $\langle\cdot\rangle$ is the average taken over a collection of
particles. If $D_{ij}^{E}$ exists and is a nonzero value, then it illustrates
that the transport of the particle under this random field exhibits a
diffusive behavior on a long-time scale.

If $\bm{v}(\bm{x},t,\omega)$ satisfies Assumption 2 and 2.1, then the
dispersion of the particles
$\big{\langle}\big{(}x_{i}(t)-x_{i}(0)\big{)}\big{(}x_{j}(t)-x_{j}(0)\big{)}\big{\rangle}$
will grow linearly with respect to the time $t$. In our previous study
Junlong2021 , we have given the convergence analysis of the effective
diffusivity calculated by the structure-preserving scheme for random fields
with normal diffusive performance.

There are also instances where the dispersion of the transport of the particle
does not exhibit linear growth with respect to the time $t$, but instead
follows a power law of $\gamma$, i.e. $\mathbb{E}[x_{i}^{2}(t)]\sim
t^{\gamma}$, where $1<\gamma<2$ and $0<\gamma<1$ are corresponding to super-
diffusion and sub-diffusion respectively. Such random fields with abnormal
diffusion behavior are what we want to focus on in this paper, and we also
want to compare the structure-preserving scheme with other methods of
computing effective diffusivity. Due to our method being implemented by the
Lagrangian particle approach, the computational cost is relatively small in
the high-dimensional cases and when $D_{0}$ is extremely small. The major
difficulty for the Lagrangian approach is to guarantee the accuracy during the
evolution which is long enough to approach the diffusion time scale. This
motivates us to develop a structure-preserving scheme that is accurate for
long-time integration and compare it with the explicit Euler-Maruyama scheme
in this problem.

###  2.2 Sharp conditions for random flows to exhibit normal diffusive
behavior

In this subsection, we discuss the criterion of random fields to exhibit
diffusion behavior from their energy spectrum density. Notice that we have
suppressed the dependence of the velocity on $\omega$ in $\bm{v}(\bm{x},t)$
for notation simplicity here. In our following analysis and experiments, we
will consider the particles moving in an isotropic velocity field
$\bm{v}(\bm{x},t)$ satisfies Assumption 2. Thus, the random field
$\bm{v}(\bm{x},t)$ is fully specified by its covariance function
$C(r,t)=Cov(\bm{v}(\bm{x},t),\bm{v}(\bm{0},0))$ with $r=||\bm{x}||_{2}$. We
assume the covariance functions have the following forms

$C(r,t)=2D(t)\int_{0}^{\infty}E(k)\frac{\sin(kr)}{kr}dk,\ \text{for 3D
flows},$ (9)

and

$C(r,t)=2D(t)\int_{0}^{\infty}E(k)J_{0}(kr)dk,\ \text{for 2D flows},$ (10)

where $D(t)$ is the time-correlation function, $E(k)$ is the energy spectrum
function, and $J_{0}(\cdot)$ is the Bessel functions of the first kind. Given
the covariance function (9) or (10), we use the randomization method in
Majda1999  to generate realizations of the random velocity field.
Specifically, we generate the velocity field $\bm{v}(\bm{x},t)$ with the form

$\bm{v}(\bm{x},t)=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\big{(}\bm{u}_{n}\cos(\bm{k}_{n}\cdot\bm{x}+\theta_{n}t)+\bm{w}_{n}\sin(\bm{k}_{n}\cdot\bm{x}+\theta_{n}t)\big{)}.$
(11)

For two-dimensional random flows, we have

$\bm{u}_{n}=\xi_{n}\bm{k}_{n}^{\perp},\
\bm{w}_{n}=\zeta_{n}\bm{k}_{n}^{\perp},\
\bm{k}_{n}^{\perp}=(-k_{n}^{2},k_{n}^{1}),\ n=1,2,...,N,$

where $\xi_{n}$ and $\zeta_{n}$ are independent 1D Gaussian variables. For
three-dimensional random flows, we have

$\bm{u}_{n}=\bm{\xi}_{n}\times\bm{k}_{n},\
\bm{w}_{n}=\bm{\zeta}_{n}\times\bm{k}_{n},\
\bm{k}_{n}=(k_{n}^{1},k_{n}^{2},k_{n}^{3}),\ n=1,2,...,N,$

where the samples $\bm{\xi}_{n}$ and $\bm{\zeta}_{n}$ are independent 3D
Gaussian variables. Sample points for each component of $\bm{k}_{n}$ are
independently chosen according to the energy spectrum function $E(k)$.

Here we assume that the velocity field is independent of the time $t$ in the
following analysis, i.e. $D(t)=1$. For a 2D random velocity field with an
isotropic energy spectrum $E(k)$, a sharp condition, which has been proved in
Fannjiang1996 , to determine whether it will perform normal diffusion behavior
is that

$\int\frac{1}{k^{2}}\hat{C}(k_{1},k_{2})dk_{1}dk_{2}<\infty.$ (12)

Here $k=\sqrt{k_{1}^{2}+k_{2}^{2}}$ and $\hat{C}(k_{1},k_{2})$ is the Fourier
transform of the covariance function $C(x_{1},x_{2})$ defined in (10). Here
$r=\sqrt{x_{1}^{2}+x_{2}^{2}}$. Then together with the correlation function
for 2D random flows (10), we can further obtain that

$\displaystyle\hat{C}(k_{1},k_{2})$ $\displaystyle=2\int
E(k^{\prime})J_{0}(k^{\prime}r)e^{-i(x_{1}k_{1}+x_{2}k_{2})}dx_{1}dx_{2}dk^{\prime}$
$\displaystyle=2\int E(k^{\prime})J_{0}(k^{\prime}r)e^{-ikr\cos\theta}rd\theta
drdk^{\prime}$ $\displaystyle=4\pi\int
E(k^{\prime})J_{0}(k^{\prime}r)J_{0}(kr)rdrdk^{\prime}$
$\displaystyle=4\pi\int
E(k^{\prime})\frac{1}{k^{\prime}}\delta(k^{\prime}-k)dk^{\prime}=4\pi
E(k)\frac{1}{k}.$ (13)

Then the sharp condition (12) for 2D cases can be transformed into

$\int\frac{1}{k^{2}}\hat{C}(k_{1},k_{2})dk_{1}dk_{2}=2\pi\int\frac{1}{k^{2}}4\pi
E(k)\frac{1}{k}kdk=8\pi^{2}\int\frac{E(k)}{k^{2}}dk<\infty.$ (14)

A similar analysis can be done for 3D random flows. For a 3D random velocity
field with energy spectrum $E(k)$, the sharp condition (14) to determine
whether it will perform normal diffusion behavior will now be transformed into

$\displaystyle\hat{C}(k_{1},$ $\displaystyle
k_{2},k_{3})=2\int_{0}^{\infty}\int_{\mathbb{R}^{3}}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}r}e^{-i(k_{1}x_{1}+k_{2}x_{2}+k_{3}x_{3})}dx_{1}dx_{2}dx_{3}dk^{\prime}$
$\displaystyle=2\int_{0}^{\infty}\int_{0}^{\infty}\int_{0}^{2\pi}\int_{0}^{\pi}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}r}e^{-ikr\sin\theta\cos\phi}r^{2}\sin\theta
d\theta d\phi drdk^{\prime}$
$\displaystyle=4\pi\int_{0}^{\infty}\int_{0}^{\infty}\int_{0}^{\pi}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}r}J_{0}(kr\sin\theta)r^{2}\sin\theta
d\theta drdk^{\prime}$
$\displaystyle=8\pi\int_{0}^{\infty}\int_{0}^{\infty}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}r}\frac{\sin(kr)}{kr}r^{2}drdk^{\prime}$
$\displaystyle=8\pi\int_{0}^{\infty}\int_{0}^{\infty}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}}\frac{\sin(kr)}{k}drdk^{\prime}.$
(15)

$\displaystyle\int\frac{1}{k^{2}}\hat{C}(k_{1},$ $\displaystyle
k_{2},k_{3})dk_{1}dk_{2}dk_{3}$
$\displaystyle=4\pi\int_{0}^{\infty}\frac{1}{k^{2}}(8\pi\int_{0}^{\infty}\int_{0}^{\infty}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}}\frac{\sin(kr)}{k}drdk^{\prime})k^{2}dk$
$\displaystyle=32\pi^{2}\int_{0}^{\infty}\int_{0}^{\infty}\int_{0}^{\infty}E(k^{\prime})\frac{\sin(k^{\prime}r)}{k^{\prime}}\frac{\sin(kr)}{k}drdk^{\prime}dk$
$\displaystyle=16\pi^{3}\int_{0}^{\infty}\int_{0}^{\infty}\frac{E(k^{\prime})\sin(k^{\prime}r)}{k^{\prime}}drdk^{\prime}<\infty.$
(16)

With the criterion of $2$-dimensional and $3$-dimensional normal diffusion
obtained, we can make the following assumption about the random velocity
fields for the subsequent study of diffusive and sub-diffusive cases. The
random velocity field $\bm{v}(\bm{x},t)$ possesses an isotropic energy
spectral density $E(k)$ which satisfies (14) in 2D cases or (16) in 3D cases.

##  3 Structure-preserving scheme for computing effective diffusivity

In this section, we would like to introduce our structure-preserving scheme
for $2$-dimensional and $d$-dimensional ($d\geq 3$) random flows (Subsection
3.1 and 3.2). Due to the difference in nature between $2$-dimensional and
$d$-dimensional random flows, the conditions required for constructing the
corresponding numerical schemes are not identical, so we discuss the
construction of the numerical schemes for the $2$-dimensional and
$d$-dimensional cases separately. Meanwhile, we verify the volume-preserving
property of our method and make a comparison with the Euler-Maruyama scheme
(Subsection 3.3). We also introduce the modified equation systems of two
numerical schemes which will be used in the convergence analysis section
(Subsection 3.4).

###  3.1 Construction of the scheme for $2$-dimensional random flows

Here we first introduce the numerical scheme for computing the effective
diffusivity of $2$-dimensional random flows. For a given random realization of
the environment, $\omega$, we first split the passive tracer model (1) into
two sub-problems based on the operator splitting method.

$d\bm{x}(t)=\bm{v}(\bm{x}(t),\omega)dt,\ \text{and}\ \ d\bm{x}(t)=\sigma
d\bm{w}(t).$ (17)

By applying the Lie-Trotter splitting method and discretizing the sub-
problems, we can develop our numerical schemes. Specifically, from time
$t_{n}$ to $t_{n+1}=t_{n}+\Delta t$, where $\Delta t>0,t_{0}=0$, assuming the
numerical solution $\bm{x}_{n}$ for time $t_{n}$ is given, we discretize the
subproblems and obtain

$\displaystyle\bm{x}_{n}^{*}=\bm{x}_{n}+\Delta
t\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega),$ (18)
$\displaystyle\bm{x}_{n+1}=\bm{x}_{n}^{*}+\sigma\sqrt{\Delta t}\bm{w}_{n},$
(19)

here $\bm{w}_{n}$ is a $d$-dimensional standard Gaussian random variable. By
solving the numerical scheme (18) (19), we obtain the numerical solution
$\bm{x}_{n+1}$, which approximates the exact solution $\bm{x}(t)$ to the SDE
system (1) at time $t=t_{n+1}$.

Here we can see that we apply the Euler midpoint rule to discrete the first
sub-problem and obtain numerical scheme (18). One natural question arises: Is
the numerical scheme (18) volume-preserving? Then we can derive that

$\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}}=\mathbf{I}+\Delta
t\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)(\frac{1}{2}\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}}+\frac{1}{2}\mathbf{I}),$
(20)

where $\mathbf{I}$ is the identity matrix of size $d$, and
$\mathbf{D}\bm{v}(\cdot)$ is the Jacobin matrix. Equation (20) implies that

$\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}}=\big{(}\mathbf{I}-\frac{\Delta
t}{2}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)\big{)}^{-1}\big{(}\mathbf{I}+\frac{\Delta
t}{2}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)\big{)}.$
(21)

To ensure that the numerical scheme (18) is volume-preserving, it is necessary
to require the condition that the determinant of
$\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}}$ should equal to one, i.e.
$\det(\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}})=1$. Inspired by the
approach in KangFeng1989 , here we denote
$P(\lambda)=\det\big{(}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)-\lambda\mathbf{I}\big{)}$
as the characteristic determinant of
$\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)$. Then we can
derive that

$\det(\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}})=\frac{(\frac{\Delta
t}{2})^{d}\det\big{(}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)+\frac{2}{\Delta
t}\mathbf{I}\big{)}}{(\frac{-\Delta
t}{2})^{d}\det\big{(}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)-\frac{2}{\Delta
t}\mathbf{I}\big{)}}=(-1)^{d}\frac{P(-\frac{2}{\Delta t})}{P(\frac{2}{\Delta
t})}.$ (22)

######  Proposition 2 (KangFeng1989 , Theorem 24)

When $d=2$, the numerical scheme (18) is volume-preserving if the velocity
field $\bm{v}$ is divergence-free.

###### Proof

When $d=2$, equation (22) will become into
$\det(\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}})=P(-\frac{2}{\Delta
t})/P(\frac{2}{\Delta t})$. Notice that

$P(\lambda)=\lambda^{2}-\text{Tr}\big{(}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)\big{)}\lambda+\det\big{(}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)\big{)}=\lambda^{2}.$
(23)

When $\bm{v}$ is divergence-free, we can further derive that
$P(\lambda)=P(-\lambda)$. Then we can obtain that
$\det(\frac{\partial\bm{x}_{n}^{*}}{\partial\bm{x}_{n}})=1$ which guarantees
the scheme (18) is volume-preserving in this case. $\square$

###  3.2 Construction of the scheme for $d$-dimensional separable random flows

For the commonly cases with $d\geq 3$, we now require the condition that
$P(\lambda)=(-1)^{d}P(-\lambda)$. To satisfy this, it requires the velocity
field $\bm{v}$ is divergence-free and $\det\big{(}\mathbf{D}\bm{v}$
$(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)\big{)}=0$ when $d=3$. It appears
that the condition
$\det\big{(}\mathbf{D}\bm{v}(\frac{\bm{x}_{n}^{*}+\bm{x}_{n}}{2},\omega)\big{)}=0$
is highly restrictive and may not be satisfied in general cases. Conditions
that are difficult to satisfy also appear when $d>3$. To overcome this issue,
we employ the splitting method proposed in Kang1995 . We split the divergence-
free vector field $\bm{v}(\bm{x},\omega):\mathbb{R}^{d}\to\mathbb{R}^{d}$ into
$d-1$ vector fields, i.e.
$\bm{v}(\bm{x},\omega)=\sum_{k=1}^{d-1}\bm{v}_{k}(\bm{x},\omega)$, where each
$\bm{v}_{k}(\bm{x},\omega)=\bm{v}_{k}(x_{k},x_{k+1},\omega)$ is Hamiltonian in
the variables $(x_{k},x_{k+1})$. Then we can split the equation (1) into $d$
sub-problems as following:

$\left\\{\begin{aligned}
&\bm{x}^{1}(t_{n})=\bm{x}(t_{n-1})+\bm{v}_{1}(\frac{\bm{x}(t_{n-1})+\bm{x}^{1}(t_{n})}{2},\omega)\Delta
t,\\\
&\bm{x}^{2}(t_{n})=\bm{x}^{1}(t_{n})+\bm{v}_{2}(\frac{\bm{x}^{1}(t_{n})+\bm{x}^{2}(t_{n})}{2},\omega)\Delta
t,\\\ &\cdots\cdots\\\
&\bm{x}^{d-1}(t_{n})=\bm{x}^{d-2}(t_{n})+\bm{v}_{d-1}(\frac{\bm{x}^{d-2}(t_{n})+\bm{x}^{d-1}(t_{n})}{2},\omega)\Delta
t,\\\ &\bm{x}(t_{n})=\bm{x}^{d-1}(t_{n})+\sigma\sqrt{\Delta
t}\bm{w}_{n-1}.\end{aligned}\right.$

For example, following the Randomization Method discussed in Majda1999 , a 3D
random velocity field can be generated by (11) with $\theta_{n}=0$. Then we
can separate this velocity field into

$\displaystyle\bm{v}_{1}(\bm{x})=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\big{(}\bm{u}_{n}^{1,2}(\bm{k}_{n})\cos(\bm{k}_{n}\cdot\bm{x})+\bm{w}_{n}^{1,2}(\bm{k}_{n})\sin(\bm{k}_{n}\cdot\bm{x})\big{)},$
$\displaystyle\bm{v}_{2}(\bm{x})=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\big{(}\bm{u}_{n}^{2,3}(\bm{k}_{n})\cos(\bm{k}_{n}\cdot\bm{x})+\bm{w}_{n}^{2,3}(\bm{k}_{n})\sin(\bm{k}_{n}\cdot\bm{x})\big{)},$

here

$\displaystyle\bm{u}_{n}^{1,2}(\bm{k}_{n})=(\xi_{n,2}k_{n,3}-\xi_{n,3}k_{n,2},\
\xi_{n,3}k_{n,1}-\xi_{n,2}\frac{k_{n,1}k_{n,3}}{k_{n,2}},\ 0),$ (24)
$\displaystyle\bm{u}_{n}^{2,3}(\bm{k}_{n})=(0,\
\xi_{n,2}\frac{k_{n,1}k_{n,3}}{k_{n,2}}-\xi_{n,1}k_{n,3},\
\xi_{n,1}k_{n,2}-\xi_{n,2}k_{n,1}),$ (25)

and by replacing $\xi_{n,i}$ with $\zeta_{n,i}$, where $i=1,2,3$ in (24),
(25), we get $\bm{w}_{n}^{1,2}(\bm{k}_{n})$ and
$\bm{w}_{n}^{2,3}(\bm{k}_{n})$. Then we adopt the first-order Lie-Trotter
splitting method to discrete these three sub-problems where the subproblems
are discretized by the Euler midpoint rule, separately. By applying the
splitting method in McLachlan2002 , we can split a $d$-dimensional velocity
field into a summation of $d-1$ velocity fields and construct a structure-
preserving scheme. More details can be found in Kang1995 , Hairer2006 . This
allows us to construct the corresponding numerical scheme for each concrete
realization of the $d$-dimensional random field.

###  3.3 Volume-preserving property of the structure-preserving scheme

In this subsection, we would like to show the volume-preserving property of
our implicit structure-preserving scheme in the stochastic sense and make a
comparison with the Euler-Maruyama scheme. Here we use 2D random flows for
illustration since, as mentioned in the previous subsection, for high
dimensional random flows we are also separating them into multiple stochastic
flows which are all Hamiltonian in two spatial variables. The analysis of high
dimensional cases can be proved by following the idea of Theorem 4.1 and
Theorem 4.2 in the next section.

By following the numerical scheme (18) (19), we can generate a solution series
$\\{\bm{x}_{n}^{\omega}\\},n=1,2,3,\cdots$. We can see that the solution
$\bm{x}_{n+1}^{\omega}$ is only determined by the solution at the previous
time step, $\bm{x}_{n}^{\omega}$. So it can be viewed as a Markov process.
Here we demonstrate some important properties of the random flows are
preserved even after numerical discretization. As an analogy to the
continuous-time case (2), we can define the environment process as viewed from
the numerical solution $\bm{x}_{n}^{\omega}$ at different discrete time steps

$\eta_{n}=\tau_{\bm{x}_{n}^{\omega}}\omega(n\Delta
t),\quad\eta_{0}=\omega(0).$ (26)

The above environment process is defined on the space of trajectories
$(\tilde{\Omega},\mathcal{B})$ with $\tilde{\Omega}=C([0,\infty)\cap\Delta
t\mathbb{Z};\mathbb{R}^{d})$. It is a subspace of $\Omega$ with time parameter
lies only on $\Delta t\mathbb{Z}$ and we still use
$\mathbf{E}_{\tilde{\bm{x}}}$ to denote the corresponding expectation
operator. Then, we can also define

$\bm{S}_{n}f(\tilde{\bm{x}})=\mathbf{ME}_{\tilde{\bm{x}}}f(\eta_{n}).$ (27)

###### Proposition 3

Under the structure-preserving numerical scheme setting, $P_{0}$ is an
invariant measure of $\eta_{n}$, i.e., $P_{0}$ is an invariant measure of the
Markov semigroup $\\{\bm{S}_{n}\\}$.

###### Proposition 4

Under the structure-preserving numerical scheme setting, $\bm{S}_{n}$ has the
property that

$||\bm{S}_{n}f||_{L^{2}(\mathbb{R}^{d})}\leq\exp\big{(}-2c_{1}(\sigma)n\Delta
t\big{)}||f||_{L^{2}(\mathbb{R}^{d})}$ (28)

The proof of Proposition 3 and Proposition 4 can be found in Junlong2021 . We
can obtain from Proposition 3 that, in the structure-preserving scheme,
$\mathbf{E}\bm{S}_{n}f=\mathbf{E}\bm{S}_{n-1}f$ for all $n$. Then we use
$\bm{\Theta}_{\Delta t}^{\omega,S}$ to denote the numerical integrator
associated with the structure-preserving scheme during $\Delta t$ time and
define the transition probability density function (PDF) of the corresponding
solution process,

$p^{S}_{\omega}(\bm{x},\bm{y})=p_{0}(\bm{\Theta}_{\Delta
t}^{\omega,S}(\bm{x}),\bm{y}),$ (29)

where $p_{0}$ is the kernel of the sub-process corresponding $d\bm{w}_{t}$,
i.e.

$p_{0}(\bm{x},\bm{y})=\frac{1}{2\pi\sigma^{2}\Delta
t}\exp\Big{(}-\frac{||\bm{y}-\bm{x}||^{2}}{2\sigma^{2}\Delta t}\Big{)}.$ (30)

###### Definition 1

A numerical scheme is volume-preserving if the integral of its transition PDF
is always equal to 1, i.e.

$\int p^{S}_{\omega}(\bm{x},\bm{y})d\bm{x}\equiv 1.$ (31)

Then we can find that

$\int p^{S}_{\omega}(\bm{x},\bm{y})d\bm{x}=\int p_{0}(\bm{\Theta}_{\Delta
t}^{\omega,S}(\bm{x}),\bm{y})d\bm{x}=\int
p_{0}(\bm{z},\bm{y})\det(\mathbf{D}\bm{\Theta}_{\Delta
t}^{\omega,S})^{-1}d\bm{z}.$ (32)

Here $\bm{z}=\bm{\Theta}_{\Delta t}^{\omega,S}(\bm{x})$. It can be verified
that in the structure-preserving scheme,

$\det(\mathbf{D}\bm{\Theta}_{\Delta
t}^{\omega,S})=\det(\frac{\partial\bm{x}^{*}}{\partial\bm{x}})\equiv 1,$ (33)

always holds from (22). Together with the fact that $\int
p_{0}(\bm{z},\bm{y})d\bm{z}=1$, we can prove that

$\int p^{S}_{\omega}(\bm{x},\bm{y})d\bm{x}\equiv 1.$ (34)

This implies that $\bm{S}_{n}$ is volume-preserving for the structure-
preserving scheme. While for the Euler-Maruyama scheme, we use
$\bm{\Theta}_{\Delta t}^{\omega,E}$ to denote the corresponding numerical
integrator and $p^{E}_{\omega}(\bm{x},\bm{y})=p_{0}(\bm{\Theta}_{\Delta
t}^{\omega,E},\bm{y})$ to denote the corresponding transition PDF. To be
specific

$\bm{\Theta}_{\Delta t}^{\omega,E}(\bm{x})=\bm{x}+\Delta t\bm{v}(\bm{x}),$
(35)

and it can be computed that

$\displaystyle\det(\mathbf{D}\bm{\Theta}_{\Delta
t}^{\omega,E})=\begin{vmatrix}1-\Delta t\Psi_{x_{1}x_{2}}&-\Delta
t\Psi_{x_{2}x_{2}}\\\ \Delta t\Psi_{x_{1}x_{1}}&1+\Delta
t\Psi_{x_{1}x_{2}}\end{vmatrix}=1+(\Delta t)^{2}\det(\mathbf{H}(\bm{x})).$
(36)

Here $\mathbf{H}(\bm{x})$ is the Hessian matrix of the stream function $\Psi$.

Following the randomization method in Kraichnan1970 , we can generate 2D
random velocity fields in the form (11) with $\theta_{n}=0$. With Assumption 2
and 2.1, here we consider $\bm{v}\in(C_{b}^{1}(\mathcal{X}))^{2}$, then we can
find the upper bound for $|\frac{\partial v_{i}}{\partial x_{j}}|,\
i,j\in\\{1,2\\},$ and denote it as $M_{1}$. Then we can obtain that

$|det(\mathbf{H}(\bm{x}))|=|\Psi_{x_{1}x_{1}}\Psi_{x_{2}x_{2}}-\Psi_{x_{1}x_{2}}^{2}|\leq
2M_{1}^{2}.$ (37)

Here $\Psi$ is the stream function (Hamiltion) of the velocity field $\bm{v}$
and we use $\Psi_{x_{i}x_{j}}$ to denote $\frac{\partial^{2}\Psi}{\partial
x_{i}\partial x_{j}}$. It can be represented by the following formula,

$\Psi(\bm{x})=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\xi_{n}\sin(\bm{k}_{n}\cdot\bm{x})-\zeta_{n}\cos(\bm{k}_{n}\cdot\bm{x})=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}\Psi_{n}(\bm{x}),$
(38)

and here we define
$\Psi_{n}(\bm{x})=\xi_{n}\sin(\bm{k}_{n}\cdot\bm{x})-\zeta_{n}\cos(\bm{k}_{n}\cdot\bm{x})$,
$n=1,\cdots,N$.

It can be obtained from (37) that $det(\mathbf{H}(\bm{x}))$ is also bounded.
Then we require $\Delta t<\frac{1}{2M_{1}}$ and then can obtain that
$\frac{1}{2}\leq 1+(\Delta t)^{2}\det(\mathbf{H}(\bm{x}))\leq\frac{3}{2}$.
This guarantees the fraction in (42) always meaningful. Then we can further
compute that

$\displaystyle det(\mathbf{H}(\bm{x}))$
$\displaystyle=\Psi_{x_{1}x_{1}}\Psi_{x_{2}x_{2}}-\Psi_{x_{1}x_{2}}^{2}$
$\displaystyle=\frac{1}{N}\big{(}(\sum_{i=1}^{N}k_{i,1}^{2}\Psi_{i}(\bm{x}))(\sum_{j=1}^{N}k_{j,2}^{2}\Psi_{j}(\bm{x}))-(\sum_{m=1}^{N}k_{m,1}k_{m,2}\Psi_{m}(\bm{x}))^{2}\big{)}$
$\displaystyle=\frac{1}{N}\sum_{i>j}(k_{i,1}k_{j,2}-k_{i,2}k_{j,1})^{2}\Psi_{i}(\bm{x})\Psi_{j}(\bm{x}).$
(39) $\displaystyle\Psi_{i}(\bm{x})\Psi_{j}(\bm{x})=$
$\displaystyle\sqrt{(\xi_{i}^{2}+\zeta_{i}^{2})(\xi_{j}^{2}+\zeta_{j}^{2})}\sin(\bm{k}_{i}\cdot\bm{x}+\alpha_{i})\sin(\bm{k}_{j}\cdot\bm{x}+\alpha_{j})$
$\displaystyle=$
$\displaystyle\frac{\sqrt{(\xi_{i}^{2}+\zeta_{i}^{2})\big{(}\xi_{j}^{2}+\zeta_{j}^{2})}}{2}(\cos((\bm{k}_{i}-\bm{k}_{j})\cdot\bm{x}+\alpha_{i}-\alpha_{j})$
$\displaystyle-\cos((\bm{k}_{i}+\bm{k}_{j})\cdot\bm{x}+\alpha_{i}+\alpha_{j})\big{)}.$
(40)

Here $\alpha_{1},\cdots,\alpha_{N}$ are obtained by the auxiliary angle
formula. Due to $k_{n,1},k_{n,2}$, $\xi_{n},\zeta_{n},n=1,\cdots N$ are
independent random variables, w.l.o.g we can assume that
$k_{i,1}k_{j,2}-k_{i,2}k_{j,1}\neq 0$ and $\bm{k}_{i}\pm\bm{k}_{j}$ are
different from each other for any $i,j$. Then the whole
$det(\mathbf{H}(\bm{x}))$ can be treated as the summation of trigonometric
functions and no terms can be eliminated from each other. So
$det(\mathbf{H}(\bm{x}))$ is not a constant function and it oscillates up and
down around 0 as $\bm{x}$ changes.

###### Lemma 1

For 2D random fields generated by (11), with $\theta_{n}=0,$ and $\Delta
t<\frac{1}{2M_{1}}$ small enough, the Euler-Maruyama scheme is not a volume-
preserving scheme, i.e.

$\int p^{E}_{\omega}(\bm{x},\bm{y})d\bm{x}\not\equiv 1.$ (41)

###### Proof

For the Euler-Maruyama scheme, the equation (32) will become

$\int p^{E}_{\omega}(\bm{x},\bm{y})d\bm{x}=\int
p_{0}(\bm{z},\bm{y})\frac{1}{1+(\Delta
t)^{2}det(\mathbf{H}(\bm{x}))}d\bm{z}=\int
p_{0}(\bm{z},\bm{y})h(\bm{z})d\bm{z}.$ (42)

Here $\bm{z}=\bm{\Theta}_{\Delta t}^{\omega,E}(\bm{x})$ and we use $h(\bm{z})$
to denote $\frac{1}{1+(\Delta t)^{2}det(\mathbf{H}(\bm{x}))}$. Then we can
find two points $\bm{z}_{0}=\bm{\Theta}_{\Delta
t}^{\omega,E}(\bm{x}_{0}),\bm{z}_{1}=\bm{\Theta}_{\Delta
t}^{\omega,E}(\bm{x}_{1})$, s.t. $det(\mathbf{H}(\bm{x}_{0}))>0$,
$det(\mathbf{H}(\bm{x}_{1}))<0$, i.e. we can find

$h(\bm{z}_{0})<1<h(\bm{z}_{1}),\ \bm{z}_{0},\bm{z}_{1}\in\mathbb{R}^{2}.$ (43)

Due to the form of $p_{0}(\bm{x},\bm{y})$ defined in (30) and $h(\bm{z})$ is
continuous, the integral $\int p_{0}(\bm{z},\bm{y})h(\bm{z})d\bm{z}$ will be
dominated by the neighborhood of $\bm{y}$ when $\Delta t$ is sufficient small.
In this situation, by setting $\bm{y}_{0}=\bm{x}_{0}$ and
$\bm{y}_{1}=\bm{x}_{1}$, we can obtain that

$\int p^{E}_{\omega}(\bm{x},\bm{y}_{0})d\bm{x}\neq\int
p^{E}_{\omega}(\bm{x},\bm{y}_{1})d\bm{x}$ (44)

So that $\int p^{E}_{\omega}(\bm{x},\bm{y})d\bm{x}\not\equiv 1$ in the
explicit Euler-Maruyama scheme. This shows that the explicit Euler-Maruyama
scheme is not volume-preserving. $\square$

###  3.4 Modified equation systems of two numerical schemes

In this subsection, we would like to first construct the modified equation
systems of two numerical schemes that will be used in the convergence
analysis. In the subsequent analysis, all variables depend on the
implementation of random velocity fields, i.e. are dependent on $\omega$, and
for writing convenience we have all ignored the $\omega$ in the writing. With
operator $\mathcal{L}_{0}$ defined in (4), we consider the backward Kolmogorov
equation defined as follows,

$u_{t}=\mathcal{L}_{0}u,\ \text{and}\ u(\bm{x},0)=\phi(\bm{x}).$ (45)

Then integrate (45) from $t=0$ to $t=\Delta t$, we can obtain that

$u(\bm{x},\Delta t)=\phi(\bm{x})+\int_{0}^{\Delta
t}\mathcal{L}_{0}u(\bm{x},s)ds.$ (46)

Then by applying weak Taylor expansion on $u(\bm{x},s)$ and denoting
$\frac{\partial^{k}}{\partial
s^{k}}u(\bm{x},0)=\mathcal{L}_{0}^{k}\phi(\bm{x})$, we can obtain that

$u(\bm{x},\Delta t)=\phi(\bm{x})+\sum_{k=0}^{N}\frac{\Delta
t^{k+1}}{(k+1)!}\mathcal{L}_{0}^{k+1}\phi(\bm{x})+O(\Delta t^{N+2}).$ (47)

For $d$-dimensional cases ($d\geq 2$), we study the flow generated by the
method mentioned in subsection 3.1 and 3.2. We define
$\mathcal{G}_{1}=\bm{v}_{1}\cdot\nabla$,
$\mathcal{G}_{2}=\bm{v}_{2}\cdot\nabla$, …,
$\mathcal{G}_{d-1}=\bm{v}_{d-2}\cdot\nabla$,
$\mathcal{G}_{d}=\frac{\sigma^{2}}{2}\Delta$ and compute

$\left\\{\begin{aligned} &\partial_{t}u^{1}=\mathcal{G}_{1}u^{1},\
u^{1}(\bm{x},0)=u(\bm{x},0),\\\ &\partial_{t}u^{2}=\mathcal{G}_{2}u^{2},\
u^{2}(\bm{x},0)=u^{1}(\bm{x},\Delta t),\\\ &\cdots\\\
&\partial_{t}u^{d-1}=\mathcal{G}_{d-1}u^{d-1},\
u^{d-1}(\bm{x},0)=u^{d-2}(\bm{x},\Delta t),\\\
&\partial_{t}u^{d}=\mathcal{G}_{d}u^{d},\
u^{d}(\bm{x},0)=u^{d-1}(\bm{x},\Delta t),\\\
&u(\bm{x},0)=\phi(\bm{x}).\end{aligned}\right.$ (48)

Then $u^{d}(\cdot,\Delta t)$ will be the flow at time $t=\Delta t$ generated
by our scheme and be an approximation to $u(\cdot,\Delta t)$. To analyze the
difference between (45) and (48), we shall resort to the Baker-Campbell-
Hausdorff (BCH) formula and replace the matrices with differential operators.
Let $\mathcal{I}_{\Delta t}=\exp(\Delta t\mathcal{L}^{num})$ denote the
composite flow operator associated with (48) and $\mathcal{L}^{num}$ be the
corresponding generator, i.e.,

$\displaystyle\mathcal{I}_{\Delta t}$ $\displaystyle u(\bm{x},0):=\exp(\Delta
t\mathcal{G}_{d})\exp(\Delta t\mathcal{G}_{d-1})\cdots\exp(\Delta
t\mathcal{G}_{1})\phi(\bm{x}),$ $\displaystyle=\exp(\Delta
t(\mathcal{G}_{d}+\cdots+\mathcal{G}_{1}))\phi(\bm{x})+\frac{\Delta
t^{2}}{2}\sum_{i>j}^{d}[\mathcal{G}_{i},\mathcal{G}_{j}]\phi(\bm{x})+O(\Delta
t^{3}),$ $\displaystyle=\exp(\Delta t\mathcal{L}_{0})\phi(\bm{x})+\frac{\Delta
t^{2}}{2}\sum_{i>j}^{d}[\mathcal{G}_{i},\mathcal{G}_{j}]\phi(\bm{x})+O(\Delta
t^{3}).$ (49)

Here $[,]$ is the Lie-Bracket. Then we apply weak Taylor expansion and obtain
that

$u^{d}(\bm{x},\Delta t)=\phi(\bm{x})+\Delta
t\mathcal{L}_{0}\phi(\bm{x})+\Delta t^{2}\mathcal{A}_{1}\phi(\bm{x})+O(\Delta
t^{3}).$ (50)

Here
$\mathcal{A}_{1}=\frac{1}{2}\sum_{i>j}^{d}[\mathcal{G}_{i},\mathcal{G}_{j}]+\frac{1}{2}\mathcal{L}_{0}^{2}$.
We can also obtain an asymptotic form of the generator of this modified
equation in terms of $\Delta t$

$\mathcal{L}^{num}=\mathcal{L}_{0}+\Delta t\mathcal{L}_{1}+\Delta
t^{2}\mathcal{L}_{2}+\cdots.$ (51)

We denote the truncated generator as

$\mathcal{L}^{\Delta t,k}=\mathcal{L}_{0}+\Delta
t\mathcal{L}_{1}+\cdots+\Delta t^{k}\mathcal{L}_{k}.$ (52)

Substitute (52) into (47). Then by comparing the corresponding term’s
coefficient with (50), we can obtain that

$\mathcal{L}_{1}=\mathcal{A}_{1}-\frac{1}{2}\mathcal{L}_{0}^{2}=\frac{1}{2}\sum_{i>j}^{d}[\mathcal{G}_{i},\mathcal{G}_{j}].$
(53)

Here we give a specific form for the 2-dimensional cases, while for the high-
dimensional cases, there are concrete discussions in the subsequent proofs.
Hence, the modified flow of $\bm{x}^{\Delta t,k}$ under the structure-
preserving scheme can be written as

$\left\\{\begin{aligned} dx_{1}=&\big{[}-\Psi_{x_{2}}+\Delta
t\big{(}\frac{1}{2}(\Psi_{x_{2}x_{2}}\Psi_{x_{1}}+\Psi_{x_{1}x_{2}}\Psi_{x_{2}})+\frac{\sigma^{2}}{4}(\Psi_{x_{1}x_{1}x_{2}}+\Psi_{x_{2}x_{2}x_{2}})\big{)}\big{]}dt\\\
&+\sigma dW_{1}+\frac{\sigma}{2}\Delta
t\Psi_{x_{1}x_{2}}dW_{1}+\frac{\sigma}{2}\Delta t\Psi_{x_{1}x_{1}}dW_{2},\\\
dx_{2}=&\big{[}\Psi_{x_{1}}-\Delta
t\big{(}\frac{1}{2}(\Psi_{x_{1}x_{2}}\Psi_{x_{1}}+\Psi_{x_{1}x_{1}}\Psi_{x_{2}})+\frac{\sigma^{2}}{4}(\Psi_{x_{1}x_{2}x_{2}}+\Psi_{x_{1}x_{1}x_{1}})\big{)}\big{]}dt\\\
&+\sigma dW_{2}-\frac{\sigma}{2}\Delta
t\Psi_{x_{1}x_{2}}dW_{2}-\frac{\sigma}{2}\Delta
t\Psi_{x_{2}x_{2}}dW_{1}.\end{aligned}\right.$ (54)

Here $\Psi$ is the stream function and we use $\Psi_{x_{1}},\Psi_{x_{2}}$ to
denote $\frac{\partial\Psi}{\partial x_{1}}$ and $\frac{\partial\Psi}{\partial
x_{2}}$ respectively. Similar notations are also used in the remaining of this
paper. Then we can denote

$\tilde{\bm{v}}_{S}:=\begin{pmatrix}\frac{1}{2}(\Psi_{x_{2}x_{2}}\Psi_{x_{1}}+\Psi_{x_{1}x_{2}}\Psi_{x_{2}})+\frac{\sigma^{2}}{4}(\Psi_{x_{1}x_{1}x_{2}}+\Psi_{x_{2}x_{2}x_{2}})\\\
-\frac{1}{2}(\Psi_{x_{1}x_{2}}\Psi_{x_{1}}+\Psi_{x_{1}x_{1}}\Psi_{x_{2}})-\frac{\sigma^{2}}{4}(\Psi_{x_{1}x_{2}x_{2}}+\Psi_{x_{1}x_{1}x_{1}})\end{pmatrix}.$
$\tilde{\mathbf{D}}_{S}:=\big{(}(\mathbf{I}_{2}+\frac{\Delta
t}{2}\begin{pmatrix}\Psi_{x_{1}x_{2}}&\Psi_{x_{1}x_{1}}\\\
-\Psi_{x_{2}x_{2}}&-\Psi_{x_{1}x_{2}}\end{pmatrix})(\mathbf{I}_{2}+\frac{\Delta
t}{2}\begin{pmatrix}\Psi_{x_{1}x_{2}}&\Psi_{x_{1}x_{1}}\\\
-\Psi_{x_{2}x_{2}}&-\Psi_{x_{1}x_{2}}\end{pmatrix})^{T}-\mathbf{I}_{2}\big{)}/\Delta
t.$ (55)

Then the density function of particles $u(\bm{x},t)$ obtained from the
structure-preserving scheme satisfies a modified Fokker-Planck equation given
by

$u_{t}=-\nabla((\bm{v}+\Delta
t\tilde{\bm{v}}_{S})u)+D_{0}\nabla\nabla:(\mathbf{I}_{2}+\Delta
t\tilde{\mathbf{D}}_{S})u.$ (56)

The inner product between matrices is denoted by
$A:B=tr(A^{T}B)=\sum_{i,j}a_{(i,j)}$ $b_{(i,j)}$. It follows that
$\Delta=\nabla\nabla:\mathbf{I}_{2}$ and $\nabla\nabla:\tilde{\mathbf{D}}_{S}$
are defined accordingly. Similar analysis can be done for the Euler-Maruyama
scheme and we can obtain and denote $\tilde{\bm{v}}_{E}$ and
$\tilde{\mathbf{D}}_{E}$ correspondingly, where

$\tilde{\bm{v}}_{E}:=\begin{pmatrix}\frac{1}{2}(\Psi_{x_{2}x_{2}}\Psi_{x_{1}}-\Psi_{x_{1}x_{2}}\Psi_{x_{2}})+\frac{\sigma^{2}}{4}(\Psi_{x_{1}x_{1}x_{2}}+\Psi_{x_{2}x_{2}x_{2}})\\\
-\frac{1}{2}(\Psi_{x_{1}x_{2}}\Psi_{x_{1}}-\Psi_{x_{1}x_{1}}\Psi_{x_{2}})-\frac{\sigma^{2}}{4}(\Psi_{x_{1}x_{2}x_{2}}+\Psi_{x_{1}x_{1}x_{1}})\end{pmatrix}.$
$\tilde{\mathbf{D}}_{E}:=\big{(}(\mathbf{I}_{2}+\frac{\Delta
t}{2}\begin{pmatrix}\Psi_{x_{1}x_{2}}&-\Psi_{x_{1}x_{1}}\\\
\Psi_{x_{2}x_{2}}&-\Psi_{x_{1}x_{2}}\end{pmatrix})(\mathbf{I}_{2}+\frac{\Delta
t}{2}\begin{pmatrix}\Psi_{x_{1}x_{2}}&-\Psi_{x_{1}x_{1}}\\\
\Psi_{x_{2}x_{2}}&-\Psi_{x_{1}x_{2}}\end{pmatrix})^{T}-\mathbf{I}_{2}\big{)}/\Delta
t.$ (57)

##  4 Convergence analysis

In this section, we shall provide the convergence analysis of our stochastic
structure-preserving scheme in computing effective diffusivity and anomalous
diffusion phenomena. It has been mentioned in Fannjiang1998  that the
diffusivity is always enhanced in incompressible flows when $\sigma>0$. So the
analysis of the diffusive and super-diffusive cases is straightforward, while
the sub-diffusive cases are investigated with $\sigma\to 0$.

With the modified flow system, we can observe and analyze the difference
between our structure-preserving scheme and the Euler-Maruyama scheme when
computing random flows with a specific energy spectrum density (Subsection
4.1). Meanwhile, the error caused by the structure-preserving scheme in
computing the power law of random fields with super-diffusive behavior is also
estimated (Subsection 4.2).

###  4.1 Convergence analysis of two schemes in computing effective
diffusivity

In this subsection, we would like to provide the convergence analysis of two
numerical schemes in computing the effective diffusivity and make a comparison
between them. To make the following analysis more concise, here we consider
the 2-dimensional random flows with stream functions defined in (38), where
$\xi_{n},\zeta_{n}$ are i.i.d. standard Gaussian variables and
$\bm{k}_{n}=(k_{n,1},k_{n,2})$ are random variables with mean zero and finite
variance. Here we require $k^{2}_{n,1}+k_{n,2}^{2}>0$. For 3-dimensional
cases, as mentioned in subsection 3.2, we can separate the problem into sub-
problems and the velocity field in each sub-problem is Hamiltonian in two
spatial variables. This is equivalent to each sub-problem being a
2-dimensional case, and then we combine these sub-problems together.

From the perspective of the stream function (Hamiltonian) $\Psi$, we can
easily see the numerical diffusion error caused by the Euler-Maruyama scheme,
but this is not relevant to our subsequent convergence analysis in computing
effective diffusivity, so we put the analysis of this part in Appendix B,
which is available to interested readers.

Then we denote and consider $\Phi(\bm{x}(t))=\frac{||\bm{x}(t)||^{2}}{d},\
\bm{x}(t)=(x_{1}(t),\cdots,x_{d}(t))^{T}$ and the effective diffusivity of the
modified flow system can be computed by
$\frac{\mathbb{E}[\Phi(\bm{x}(t))]}{2t}$.

###### Remark 1

Since the random velocity fields we study are all isotropic, and thus in
calculating the effective diffusivity, the value of
$\frac{\mathbb{E}[\Phi(t)]}{2t}$ is equivalent to the value of $D^{11}_{E}$.

Due to the spectral gap property and the operator for these two numerical
schemes being uniform elliptic when $\Delta t$ is small enough, we can
consider for the corresponding corrector problem of these two schemes and
define

$\text{Euler-Maruyama:}\
\hat{\bm{\mu}}_{E}(t):=\int_{0}^{\infty}(\bm{v}(\bm{x}_{E}(t+s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(t+s)))ds.$ (58) $\text{Structure-Preserving:}\
\hat{\bm{\mu}}_{S}(t):=\int_{0}^{\infty}(\bm{v}(\bm{x}_{S}(t+s))+\Delta
t\tilde{\bm{v}}_{S}(\bm{x}_{S}(t+s)))ds.$ (59)

Here we use $\bm{x}_{E}$ and $\bm{x}_{S}$ to denote the different paths of
these two numerical schemes’ modified equation system. And
$\hat{\bm{\mu}}_{E}$ satisfies the following continuous-type corrector problem
which corresponds to the modified equation system.

$\mathcal{L}^{E}\hat{\bm{\mu}}_{E}=-(\bm{v}+\Delta t\tilde{\bm{v}}_{E}).$ (60)

Here $\mathcal{L}^{E}$ is the generator of the semigroup corresponding to the
Euler-Maruyama scheme. Then we can deduce that

$\frac{d\hat{\bm{\mu}}_{E}}{dt}=-(\bm{v}(\bm{x}_{E}(t))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(t))).$ (61)

Similar results can also be obtained for $\hat{\bm{\mu}}_{S}$. With
$\bm{x}_{E}$ and $\bm{x}_{S}$, we can define the corresponding functions
$\Phi^{E}$ and $\Phi^{S}$ for two numerical schemes, i.e.

$\Phi^{E}(t)=\Phi(\bm{x}_{E}(t)),\ \text{and}\
\Phi^{S}(t)=\Phi(\bm{x}_{S}(t)).$ (62)

Then we also denote $\mathbf{G}_{E}$ and $\mathbf{G}_{S}$ as the following
formulas.

$\mathbf{G}_{E}:=\sigma(\mathbf{I}_{2}+\frac{\Delta
t}{2}\begin{pmatrix}\Psi_{x_{1}x_{2}}&-\Psi_{x_{1}x_{1}}\\\
\Psi_{x_{2}x_{2}}&-\Psi_{x_{1}x_{2}}\end{pmatrix}),\quad\mathbf{G}_{S}:=\sigma(\mathbf{I}_{2}+\frac{\Delta
t}{2}\begin{pmatrix}\Psi_{x_{1}x_{2}}&\Psi_{x_{1}x_{1}}\\\
-\Psi_{x_{2}x_{2}}&-\Psi_{x_{1}x_{2}}\end{pmatrix}).$ (63)

###### Lemma 2

Let $\bm{x}_{E}(t)$ and $\bm{x}_{S}(t)$ be the solution of the modified flow
system of two schemes with $\bm{v}$ satisfies Assumption 2.2. Then together
with $\hat{\bm{\mu}}_{E}$, we have the following formula for $\Phi^{E}(t)$ in
the Euler-Maruyama scheme.

$\frac{d\Phi^{E}}{dt}=-\frac{d(\bm{x}_{E}\cdot\hat{\bm{\mu}}_{E})}{dt}+(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})\cdot\hat{\bm{\mu}}_{E}+\sigma^{2}+O(\sigma^{2}\Delta
t^{2})+\bm{x}_{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt},$ (64)

where the integral of $\bm{x}_{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt}$ is
a mean zero martingale, i.e.
$\mathbf{E}[\int_{0}^{t}\bm{x}_{E}\cdot\mathbf{G}_{E}d\bm{w}_{t}]$ $=0$.
Similar results can also be obtained for the modified flow system of the
structure-preserving scheme, i.e.

$\frac{d\Phi^{S}}{dt}=-\frac{d(\bm{x}_{S}\cdot\hat{\bm{\mu}}_{S})}{dt}+(\bm{v}+\Delta
t\tilde{\bm{v}}_{S})\cdot\hat{\bm{\mu}}_{S}+\sigma^{2}+O(\sigma^{2}\Delta
t^{2})+\bm{x}_{S}\cdot\mathbf{G}_{S}\frac{d\bm{w}_{t}}{dt},$ (65)

###### Theorem 4.1

Let $\bm{x}_{t},\ t>0$ be the solutions of the passive tracer model (1) with
the random field satisfies Assumption 2.2. Let $\hat{\bm{\mu}}$ be the
solution of the corrector problem corresponding to (1). Then we have both the
error of two schemes’ modified equation systems in computing effective
diffusivity converge to $0$ as $\Delta t\to 0$. To be specific, we have

$\frac{\mathbb{E}[\Phi^{*}(t)]}{2t}=\frac{\sigma^{2}}{2}+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}+\rho(\Delta
t)+O(\sigma^{2}\Delta t^{2})+O(\frac{1}{\sqrt{t\Delta t}}).$ (66)

Here $\rho(\Delta t)=O(\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}})$ is a function
which satisfies $\lim_{\Delta t\to 0}\rho(\Delta t)=0$ and is independent of
$t$. Equation (66) holds for $\Phi^{*}(t)$ taken as $\Phi^{E}(t)$ or
$\Phi^{S}(t)$. Moreover, for $d$-dimensional ($d\geq 3$) cases, equation (66)
also holds.

From the Assumption 2 and spectral gap condition, $\bm{v}$ and
$\hat{\bm{\mu}}$ are both bounded. Then for a given time $t$,
$\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]$ will be a
finite value. Then for the structure-preserving scheme, we can split
$\rho(\Delta t)$ into two parts, i.e. $\rho(\Delta t)=\rho_{1}^{S}(\Delta
t)+\rho_{2}^{S}(\Delta t)$ and denote

$\rho_{1}^{S}(\Delta
t):=\big{|}\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\bm{v}(\bm{x}^{S}_{s+\tau})\rangle
d\tau
ds-\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}_{s})\cdot\bm{v}(\bm{x}_{s+\tau})\rangle
d\tau ds\big{|}\Big{/}2t.$ (67) $\displaystyle\rho_{2}^{S}(\Delta t):=$
$\displaystyle\Delta
t\big{|}\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds\
+\int_{0}^{t}\int_{0}^{\infty}\langle\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s})\cdot\bm{v}(\bm{x}^{S}_{s+\tau})\rangle$
$\displaystyle d\tau ds+\Delta
t\int_{0}^{t}\int_{0}^{\infty}\langle\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds\big{|}\Big{/}2t.$ (68)

Similar definitions can be done for the Euler-Maruyama scheme and we can
obtain $\rho_{1}^{E}(\Delta t),\ \rho_{2}^{E}(\Delta t)$. Here we use
$\bm{x}_{t},\bm{x}^{S}_{t}$ and $\bm{x}^{E}_{t}$ to denote
$\bm{x}(t),\bm{x}_{S}(t)$ and $\bm{x}_{E}(t)$. $\rho_{1}^{S}(\Delta t)$ and
$\rho_{1}^{E}(\Delta t)$ will keep the original diffusion behavior, which is
common to both numerical schemes. Then we consider for the difference of two
numerical schemes between $\rho_{2}^{S}(\Delta t)$ and $\rho_{2}^{E}(\Delta
t)$.

###### Theorem 4.2

For the structure-preserving scheme, $\rho_{2}^{S}(\Delta t)$ can be bounded
by the summation of one higher-order term of $\Delta t$ and one term
correlated with the original convection enhanced behavior, i.e.

$\rho_{2}^{S}(\Delta t)=O(\Delta
t\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t})+O(\Delta
t^{\varrho}).$ (69)

Here $\varrho>1$. Moreover, when $\sigma^{2}$ has the same scale as a fixed
$\Delta t$, $\rho_{2}^{S}(\Delta t)$ in the structure-preserving method can
keep the original sub-diffusion behavior of the system when $\sigma$ vanishes,
i.e.

$\lim_{t\to\infty}\frac{\rho_{2}^{S}(\Delta
t)}{\frac{\sigma^{2}}{2}+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}}=O(\Delta
t^{\varrho-1})<\infty$ (70)

While the Euler-Maruyama method does not hold this property.

The proof of Lemma 2, Theorem 4.1 and Theorem 4.2 can be seen in the
subsection 4.3. Theorem 4.2 and the following Corollary 1 also holds for
3-dimensional cases, the detailed explanation we put together in the proof of
Theorem 4.2.

###### Corollary 1

When diffusion coefficient $\sigma^{2}$ has the same scale as a fixed $\Delta
t$, the error terms $\rho_{2}^{S}(\Delta t)$ caused by the structure-
preserving will keep the behavior of the original random velocity field with
energy spectral $E(k)$ which will perform the diffusion behavior, i.e.

$\lim_{t\to\infty}\frac{\rho_{2}^{S}(\Delta
t)}{\frac{\sigma^{2}}{2}+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}}<O(\Delta
t)+O(\Delta t^{\varrho-1})<\infty.$ (71)

###### Proof

for the field with diffusion behavior, the term $\lim_{t\to\infty}$
$\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}$
will be a non-zero constant and is much larger than $\sigma^{2}$ when
$\sigma^{2}$ has the same scale as a fixed $\Delta t$. Then we can obtain that

$\displaystyle\lim_{t\to\infty}\frac{\rho_{2}^{S}(\Delta
t)}{\frac{\sigma^{2}}{2}+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}}$
$\displaystyle=\lim_{t\to\infty}\frac{h_{1}\Delta
t\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}+h_{2}\Delta
t^{\varrho}}{h_{3}\Delta
t+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}}$
$\displaystyle<O(\Delta t)+O(\Delta t^{\varrho-1})<\infty.$ (72)

$\square$

###  4.2 Error estimate for structure-preserving scheme in super-diffusive
cases

For random fields with energy spectral $E(k)$ which satisfies (14), we can
define the corresponding corrector problem well due to the spectral gap
condition (28). Instead, for those random fields with energy spectral that
does not satisfy condition (14) and exhibit super diffusion performance,
$\hat{\bm{\mu}}_{E}$ and $\hat{\bm{\mu}}_{S}$ may diverge and can not be well
defined.

###### Proposition 5

For 2D isotropic random flows with velocity spectrum

$\hat{R}_{i,j}(\bm{k})\sim\frac{1}{|\bm{k}|^{2\alpha}}(\delta_{i,j}-\frac{k_{i}k_{j}}{|\bm{k}|^{2}}),$
(73)

here $\ 0<\alpha<1,\ |\bm{k}|\leq L$. Then it can be obtained that the
dispersion of this system will have the following scaling law after $t$ is
large enough,

$\mathbb{E}[(\bm{x}_{t})^{2}]\sim t^{\frac{2}{2-\alpha}}.$ (74)

The proof of Proposition 5 can be found in Fannjiang1998  , then the two-point
correlation function $R_{i,j}(\bm{x})$ is asymptotically
$R_{i,j}(\bm{x})\sim\frac{1}{|\bm{x}|^{2(1-\alpha)}},\ |\bm{x}|\gg 1$. The
random velocity field $\bm{v}$ which has the velocity spectrum mentioned in
Proposition 5 and its corresponding stream function can be realized by

$\Psi(\bm{x})=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\frac{1}{|\bm{k}_{i}|}(\xi_{i}\sin(\bm{k}_{i}\cdot\bm{x})-\zeta_{i}\cos(\bm{k}_{i}\cdot\bm{x})),\
and\ \bm{v}(\bm{x})=\nabla^{\perp}\Psi(\bm{x}).$ (75)

Here $\forall\ 1\leq i\leq N$, $\xi_{i},\zeta_{i}$ are i.i.d standard Gaussian
random variables, the direction of $\bm{k}_{i}$ is uniformly distributed on a
unit circle and the length of $\bm{k}_{i}$ satisfies the density function
$\rho(|\bm{k}_{i}|)\sim 1/|\bm{k}_{i}|^{2\alpha-1},0\leq|\bm{k}_{i}|\leq L$.
Here the form of $\Psi(\bm{x})$ is similar but different from (38). Then we
have the following Theorem 4.3 for the structure-preserving scheme and its
proof can be seen from subsection 4.3.

###### Theorem 4.3

Let $\bm{x}_{t}^{S}$ be the solution of the modified flow system of the
structure-preserving scheme at time $t$ with the original random field
possessing an isotropic velocity spectrum satisfying (73). Then for
$\frac{1}{4}\leq\alpha<1$ and $t$ large enough, the modified flow system will
exhibit super-diffusion performance, and the scaling law of its dispersion is
consistent with the original system, i.e.

$\mathbf{E}[(\bm{x}_{t}^{S})^{2}]\sim t^{\frac{2}{2-\alpha}},\
\frac{1}{4}\leq\alpha<1.$ (76)

###  4.3 Proof of Lemma and Theorems

#### Proof of Lemma 2

###### Proof

By applying Ito’s formula on $\Phi^{E}(\bm{x}(t))$ of the modified system of
the Euler-Maruyama scheme, we can obtain that

$d\Phi^{E}=(\bm{x}_{E}\cdot(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})+\frac{\sigma^{2}}{2}\text{Tr}(\mathbf{G}_{E}^{T}\mathbf{G}_{E}))dt+\bm{x}_{E}\cdot\mathbf{G}_{E}d\bm{w}_{t}.$
(77)

Then we can compute that

$\frac{\sigma^{2}}{2}\text{Tr}(\mathbf{G}_{E}^{T}\mathbf{G}_{E})=\sigma^{2}+\frac{\sigma^{2}\Delta
t^{2}}{8}(\Psi_{x_{1}x_{1}}^{2}+2\Psi_{x_{1}x_{2}}^{2}+\Psi_{x_{2}x_{2}}^{2}),$
(78)

and denote the following integral as

$\displaystyle\rho_{E}$
$\displaystyle:=\int_{0}^{t}\bm{x}_{E}(s)\cdot(\bm{v}(\bm{x}_{E}(s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s)))ds$
$\displaystyle=\int_{0}^{t}(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})\cdot\hat{\bm{\mu}}_{E}ds-[\bm{x}_{E}\cdot\hat{\bm{\mu}}_{E}]_{0}^{t}.$
(79)

A similar definition can also be obtained for $\rho_{S}$.

$\displaystyle\rho_{S}$
$\displaystyle:=\int_{0}^{t}\bm{x}_{S}(s)\cdot(\bm{v}(\bm{x}_{S}(s))+\Delta
t\tilde{\bm{v}}_{S}(\bm{x}_{S}(s)))ds$
$\displaystyle=\int_{0}^{t}(\bm{v}+\Delta
t\tilde{\bm{v}}_{S})\cdot\hat{\bm{\mu}}_{S}ds-[\bm{x}_{S}\cdot\hat{\bm{\mu}}_{S}]_{0}^{t}.$
(80)

So consider the derivative to time $t$ and move the Brownian motion related
term inside the martingale related term, we can get

$\frac{d\Phi^{E}}{dt}=-\frac{d(\bm{x}_{E}\cdot\hat{\bm{\mu}}_{E})}{dt}+(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})\cdot\hat{\bm{\mu}}_{E}+\sigma^{2}+O(\sigma^{2}\Delta
t^{2})+\bm{x}_{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt}$ (81)

where $\mathbb{E}[\int_{0}^{t}\bm{x}_{E}\cdot\mathbf{G}_{E}d\bm{w}_{t}]=0$.
For the structure-preserving scheme, we can similarly obtain that

$\frac{d\Phi^{S}}{dt}=-\frac{d(\bm{x}_{S}\cdot\hat{\bm{\mu}}_{S})}{dt}+(\bm{v}+\Delta
t\tilde{\bm{v}}_{S})\cdot\hat{\bm{\mu}}_{S}+\sigma^{2}+O(\sigma^{2}\Delta
t^{2})+\bm{x}_{S}\cdot\mathbf{G}_{S}\frac{d\bm{w}_{t}}{dt},$ (82)

here $\mathbb{E}[\int_{0}^{t}\bm{x}_{S}\cdot\mathbf{G}_{S}d\bm{w}_{t}]=0$.
$\square$

#### Proof of Theorem 4.1

###### Proof

For the Euler-Maruyama scheme, we can obtain from Lemma 2 that

$\frac{d\Phi^{E}}{dt}=-\frac{d(\bm{x}_{E}\cdot\hat{\bm{\mu}}_{E})}{dt}+(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})\cdot\hat{\bm{\mu}}_{E}+\sigma^{2}+O(\sigma^{2}\Delta
t^{2})+\bm{x}_{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt},$ (83)

It has been proved in Junlong2021  that the term
$-[\bm{x}_{E}\cdot\hat{\bm{\mu}}_{E}]_{0}^{t}/2t$ is $O(\frac{1}{t\sqrt{\Delta
t}})$. It is negligible when $t$ is large enough. And then we consider the
term $(\bm{v}+\Delta t\tilde{\bm{v}}_{E})\cdot\hat{\bm{\mu}}_{E}$ at time $t$.
We want to first prove that the solution $\hat{\bm{\mu}}_{E}$ converges to the
solution of the original system’s corrector problem $\hat{\bm{\mu}}$ as
$\Delta t\to 0$.

Considering the spectral gap condition, we can find a truncation time
$T_{0}=N\Delta t$ and obtain the following inequalities

$||\int_{T_{0}}^{\infty}(\bm{v}(\bm{x}_{E}(s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s)))ds||_{L^{2}}\leq\frac{1}{2c_{1}}\exp(-2c_{1}T_{0}),$
(84)
$||\int_{T_{0}}^{\infty}\bm{v}(\bm{x}(s))ds||_{L^{2}}\leq\frac{1}{2c_{1}}\exp(-2c_{1}T_{0}).$
(85)

Then for any $\epsilon>0$, we can set the truncation time $T_{0}$ big enough
to let the following inequality hold

$||\int_{T_{0}}^{\infty}(\bm{v}(\bm{x}_{E}(s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s)))ds-\int_{T_{0}}^{\infty}\bm{v}(\bm{x}(s))ds||_{L^{2}}\leq\frac{1}{c_{1}}\exp(-c_{1}T_{0})<\epsilon.$
(86)

Then we consider the error before the truncation time $T_{0}$

$\displaystyle||\int_{0}^{T_{0}}$ $\displaystyle(\bm{v}(\bm{x}_{E}(s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s)))ds-\int_{0}^{T_{0}}\bm{v}(\bm{x}(s))ds||_{L^{2}}$
$\displaystyle\leq||\int_{0}^{T_{0}}\bm{v}(\bm{x}_{E}(s))ds-\int_{0}^{T_{0}}\bm{v}(\bm{x}(s))ds||_{L^{2}}+||\int_{0}^{T_{0}}\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s))ds||_{L^{2}}$
$\displaystyle=||\int_{0}^{T_{0}}(\bm{v}(\bm{x}_{E}(s))-\bm{v}(\bm{x}(s)))ds||_{L^{2}}+\Delta
t||\int_{0}^{T_{0}}\tilde{\bm{v}}_{E}(\bm{x}_{E}(s))ds||_{L^{2}}.$ (87)

From Assumption 2 that $\bm{v}$ is bounded, we can obtain that $\Delta
t\tilde{\bm{v}}_{E}$ is also bounded. Then we can obtain that $\Delta
t||\int_{0}^{T_{0}}\tilde{\bm{v}}_{E}(\bm{x}_{E}(s))ds||_{L^{2}}$ is
$O(T_{0}\Delta t)=O(N\Delta t^{2})$. Also from the error estimate for Euler-
Maruyama scheme in Kloeden1992 , we know that
$||\int_{0}^{T_{0}}(\bm{v}(\bm{x}_{E}(s))-\bm{v}(\bm{x}(s)))ds||_{L^{2}}$ is
$O(T_{0}\Delta t)=O(N\Delta t^{2})$. Then we can combine these two terms
together and obtain that

$||\int_{0}^{T_{0}}(\bm{v}(\bm{x}_{E}(s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s)))ds-\int_{0}^{T_{0}}\bm{v}(\bm{x}(s))ds||_{L^{2}}\leq
c_{2}T_{0}\Delta t.$ (88)

here $c_{2}$ is a constant. Then we can choose $\Delta
t<\epsilon/(c_{2}T_{0})$. And we can choose
$\epsilon\approx\frac{1}{c_{1}}\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}}$ to
balance each value of $\epsilon$ where $c_{2}=\exp(c_{3}T_{0})$.

Now we finish the proof that the solution $\hat{\bm{\mu}}_{E}$ converges to
the solution of the original system’s corrector problem $\hat{\bm{\mu}}$ as
$\Delta t\to 0$. A similar analysis can be done for the solution
$\hat{\bm{\mu}}_{S}$. Then we go back to the equation (81). We can obtain that
the long-time behavior of $(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})\cdot\hat{\bm{\mu}}_{E}$ will converge to the original
system’s behavior as $\Delta t\to 0$. Also, we can further get

$\displaystyle\frac{\mathbb{E}[\int_{0}^{t}(\bm{v}(\bm{x}_{E}(s))+\Delta
t\tilde{\bm{v}}_{E}(\bm{x}_{E}(s)))\cdot\hat{\bm{\mu}}_{E}(s)ds]}{2t}=\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}(s))\cdot\hat{\bm{\mu}}(s)ds]}{2t}+\rho(\Delta
t).$ (89)

Here $\rho(\Delta t)=O(\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}})$ is a function
which satisfies $\lim_{\Delta t\to 0}\rho(\Delta t)=0$ and is independent of
$t$.

Then we consider d-dimensional cases ($d\geq 3$) processed under the method
mentioned in subsection 3.2. If we separate the velocity field
$\bm{v}(x_{1},x_{2},\cdots,x_{d})=\sum_{i=1}^{d-1}\bm{v}_{i}(x_{i},x_{i+1})$
and use the Euler-Maruyama scheme in the advection sub-problems, we can
compute and obtain the following modified equation system

$d\bm{x}_{E}=\big{(}\sum_{i=1}^{d-1}\bm{v}_{i}+\Delta
t\sum_{j=1}^{d-1}\tilde{\bm{v}}_{j}^{E}\big{)}dt+(\sigma\mathbf{I}_{d}+\Delta
t\tilde{\mathbf{D}}_{E})d\bm{w}_{t}$ (90)

Here

$\tilde{\bm{v}}_{j}^{E}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\
\frac{1}{2}(\Psi^{(j)}_{x_{j+1}x_{j+1}}\Psi^{(j)}_{x_{j}}-\Psi^{(j)}_{x_{j}x_{j+1}}\Psi^{(j)}_{x_{j+1}})+\frac{\sigma^{2}}{4}(\Psi^{(j)}_{x_{j}x_{j}x_{j+1}}+\Psi^{(j)}_{x_{j+1}x_{j+1}x_{j+1}})\\\
-\frac{1}{2}(\Psi^{(j)}_{x_{j}x_{j+1}}\Psi^{(j)}_{x_{j}}-\Psi^{(j)}_{x_{j}x_{j}}\Psi^{(j)}_{x_{j+1}})-\frac{\sigma^{2}}{4}(\Psi^{(j)}_{x_{j}x_{j+1}x_{j+1}}+\Psi^{(j)}_{x_{j}x_{j}x_{j}})\\\
0\\\ \vdots\\\ 0\\\ \end{pmatrix}.$

Here $\Psi^{(j)}$ is the Hamiltonian of $\bm{v}_{j}$ in its corresponding two
spatial variables $x_{j},x_{j+1}$. The first $j-1$ components and the last
$d-j-1$ components of $\tilde{\bm{v}}_{j}^{E}$ are all $0$. From the previous
analysis, the additional coefficient of the Brownian motion term, $\Delta
t\tilde{\mathbf{D}}_{E}$, will result in an error term that can be included by
$O(\sigma^{2}\Delta t)$ in computing the effective diffusivity. Then we just
treat $\sum_{i=1}^{d-1}\bm{v}_{i}=\bm{v}$ and $\Delta
t\sum_{j=1}^{d-1}\tilde{\bm{v}}_{j}^{E}=\Delta t\tilde{\bm{v}}_{E}$ and then
applying Ito’s formula on $\Phi^{E}(t)$. The subsequent proof is independent
of dimensionality and, following the idea of the previous proof, we obtain
that (66) also holds for the Euler-Maruyama scheme in $d$-dimensional cases.

If we use the structure-preserving scheme in the advection sub-problems, we
can compute and obtain the following modified equation system

$d\bm{x}_{S}=\big{(}\sum_{i=1}^{d-1}\bm{v}_{i}+\Delta
t\sum_{j=1}^{d-1}\tilde{\bm{v}}_{j}^{S}\big{)}dt+(\sigma\mathbf{I}_{d}+\Delta
t\tilde{\mathbf{D}}_{S})d\bm{w}_{t}$ (91)

Here

$\tilde{\bm{v}}_{j}^{S}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\
\frac{1}{2}(\Psi^{(j)}_{x_{j+1}x_{j+1}}\Psi^{(j)}_{x_{j}}+\Psi^{(j)}_{x_{j}x_{j+1}}\Psi^{(j)}_{x_{j+1}})+\frac{\sigma^{2}}{4}(\Psi^{(j)}_{x_{j}x_{j}x_{j+1}}+\Psi^{(j)}_{x_{j+1}x_{j+1}x_{j+1}})\\\
-\frac{1}{2}(\Psi^{(j)}_{x_{j}x_{j+1}}\Psi^{(j)}_{x_{j}}+\Psi^{(j)}_{x_{j}x_{j}}\Psi^{(j)}_{x_{j+1}})-\frac{\sigma^{2}}{4}(\Psi^{(j)}_{x_{j}x_{j+1}x_{j+1}}+\Psi^{(j)}_{x_{j}x_{j}x_{j}})\\\
0\\\ \vdots\\\ 0\\\ \end{pmatrix}.$

And the first $j-1$ components and the last $d-j-1$ components of
$\tilde{\bm{v}}_{j}^{S}$ are all $0$. Then a similar analysis can be done for
the structure-preserving scheme.

$\square$

#### Proof of Theorem 4.2

###### Proof

Consider for the structure-preserving scheme when diffusion coefficient
$\sigma^{2}$ has the same scale as a fixed $\Delta t$, we can obtain that

$\displaystyle(\bm{v}+$ $\displaystyle\Delta
t\tilde{\bm{v}}_{S})\cdot\hat{\bm{\mu}}_{S}=\bm{v}\cdot\int_{0}^{\infty}\bm{v}(\bm{y}(s))ds+\Delta
t(\tilde{\bm{v}}_{S}\cdot\int_{0}^{\infty}\bm{v}(\bm{y}(s))ds+$
$\displaystyle\bm{v}\cdot\int_{0}^{\infty}\tilde{\bm{v}}_{S}(\bm{y}(s))ds)+\Delta
t^{2}\tilde{\bm{v}}_{S}\cdot\int_{0}^{\infty}\tilde{\bm{v}}_{S}(\bm{y}(s))ds,\
\text{where}\ \bm{y}(s)=\bm{x}_{S}(t+s).$ (92)

Here we use $\bm{x}^{S}_{t}$ to denote $\bm{x}_{S}(t)$. Then we can obtain
that,

$\displaystyle\mathbb{E}$
$\displaystyle[\Phi(t)]=\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\bm{v}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds+\Delta
t\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds\ $ $\displaystyle+\Delta
t\int_{0}^{t}\int_{0}^{\infty}\langle\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s})\cdot\bm{v}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds+\Delta
t^{2}\int_{0}^{t}\int_{0}^{\infty}\langle\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds$
$\displaystyle-\langle[\bm{x}^{S}\cdot\hat{\bm{\mu}}_{S}]_{0}^{t}\rangle+\sigma^{2}t+O(\sigma^{2}\Delta
t^{2}t).$ (93)

Then for the energy spectrum $E(k)$ with $\sigma\to 0$, the corresponding
terms for the structure-preserving scheme can be represented by

$\displaystyle|\int_{0}^{t}$
$\displaystyle\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds|$ $\displaystyle<$ $\displaystyle\
\frac{1}{2}|\int_{0}^{t}\int_{0}^{\infty}\langle\Psi_{x_{1}x_{2}}(\bm{v}(\bm{x}^{S}_{s})\cdot\bm{v}(\bm{x}^{S}_{s+\tau}))\rangle
d\tau ds|$
$\displaystyle+\frac{r^{2}}{2}|\int_{0}^{t}\int_{0}^{\infty}\langle\Psi\
(v_{1}(\bm{x}^{S}_{s})v_{2}(\bm{x}^{S}_{s+\tau}))\rangle d\tau ds|$
$\displaystyle+\frac{\sigma^{2}r^{2}}{4}|\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds|$ $\displaystyle\leq$ $\displaystyle\
(\frac{C_{1}}{2}+\frac{\sigma^{2}r^{2}}{4})|\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\bm{v}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds|$ $\displaystyle=$ $\displaystyle\
(\frac{C_{1}}{2}+\frac{\sigma^{2}r^{2}}{4})\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}^{S}_{s})\cdot\hat{\bm{\mu}}(\bm{x}^{S}_{s})ds].$
(94)

Here
$\hat{\bm{\mu}}(\bm{x}^{S}_{s})=\int_{0}^{\infty}\bm{v}(\bm{x}^{S}_{\tau+s})d\tau$.
Then we consider the following error estimate

$||\int_{0}^{t}\bm{v}(\bm{x}^{S}_{s})\cdot\hat{\bm{\mu}}(\bm{x}^{S}_{s})ds-\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds||_{L^{2}}\leq
C_{2}t(\Delta t)^{\frac{2c_{1}}{2c_{1}+c_{3}}}.$ (95)

Here $C_{1}$ and $C_{2}$ are two positive constants. Now we substitute (95)
back (94) and obtain that

$\displaystyle|\Delta
t\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds|\leq$ $\displaystyle(\frac{C_{1}}{2}+\frac{r^{2}\sigma^{2}}{4})\Delta
t\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]$
$\displaystyle+(\frac{C_{1}C_{2}}{2}+\frac{C_{2}r^{2}\sigma^{2}}{4})t\Delta
t^{\varrho}.$ (96)

here $\varrho=1+\frac{2c_{1}}{2c_{1}+c_{3}}>1$. A similar analysis can be done
for the third term on the right-hand side of (93). The fourth term on the
right-hand side of (93) is a high-order term of $\Delta t$. For d-dimensional
cases using the structure-preserving scheme, the term
$\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds$ can be transformed into

$\displaystyle\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}_{S}(\bm{x}^{S}_{s+\tau})\rangle
d\tau
ds=\int_{0}^{t}\int_{0}^{\infty}\langle(\sum_{i=1}^{d-1}\bm{v}_{i}(\bm{x}_{s}^{S}))\cdot(\sum_{j=1}^{d-1}\tilde{\bm{v}}^{S}_{j}(\bm{x}^{S}_{s+\tau}))\rangle
d\tau ds$
$\displaystyle=\int_{0}^{t}\int_{0}^{\infty}\Big{(}\sum_{i=1}^{d-1}\langle\bm{v}_{i}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}^{S}_{i}(\bm{x}^{S}_{s+\tau})\rangle+\sum_{|i-j|=1}\langle\bm{v}_{i}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}^{S}_{j}(\bm{x}^{S}_{s+\tau})\rangle\Big{)}d\tau
ds.$ (97)

The analysis of the terms
$\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}_{i}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}^{S}_{i}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds$ can just follow the idea of (94) and (95). Then for the other terms
in (97), we have

$\displaystyle\int_{0}^{t}\int_{0}^{\infty}\langle\bm{v}_{i}(\bm{x}^{S}_{s})\cdot\tilde{\bm{v}}^{S}_{j}(\bm{x}^{S}_{s+\tau})\rangle
d\tau ds$
$\displaystyle=\int_{0}^{t}\int_{0}^{\infty}\Big{\langle}\frac{1}{2}\Psi^{(i)}_{x_{i}}(\bm{x}^{S}_{s})\Big{(}\Psi^{(j)}_{x_{j+1}x_{j+1}}(\bm{x}^{S}_{s+\tau})\Psi^{(j)}_{x_{j}}(\bm{x}^{S}_{s+\tau})+\Psi^{(j)}_{x_{j}x_{j+1}}(\bm{x}^{S}_{s+\tau})\cdot$
$\displaystyle\
\Psi^{(j)}_{x_{j+1}}(\bm{x}^{S}_{s+\tau})\Big{)}+\frac{\sigma^{2}}{4}\Psi^{(i)}_{x_{i}}(\bm{x}^{S}_{s})\Big{(}\Psi^{(j)}_{x_{j}x_{j}x_{j+1}}(\bm{x}^{S}_{s+\tau})+\Psi^{(j)}_{x_{j+1}x_{j+1}x_{j+1}}(\bm{x}^{S}_{s+\tau})\Big{)}\Big{\rangle}d\tau
ds$ $\displaystyle=0.$ (98)

This is due to the random variables in $\Psi^{(i)}$ and $\Psi^{(j)}$ are
independent and mean zero. So equation (96) also holds in $d$-dimensional
cases. Then we can obtain that for the structure-preserving scheme

$\displaystyle\frac{\mathbb{E}[\Phi(t)]}{2t}=$
$\displaystyle\frac{\sigma^{2}}{2}+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}+\rho_{1}(\Delta
t)+O(\Delta
t\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t})$
$\displaystyle+O(\Delta t^{\varrho})+O(\sigma^{2}\Delta
t^{2})+O(\frac{1}{\sqrt{t\Delta t}}).$ (99)

For the structure-preserving scheme, the original error estimate term
$\rho(\Delta t)$ can be separated into $\rho_{1}^{S}(\Delta t)+O(\Delta
t\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t})+O(\Delta
t^{\varrho})$. Here $\rho_{1}^{S}(\Delta t)$ is just the error estimate for
(95), $O(\Delta t^{\frac{2c_{1}}{2c_{1}+c_{3}}})$, which is caused by the
original random velocity field under different paths. The other terms caused
by the modified velocity field can be represented by $O(\Delta
t\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t})+O(\Delta
t^{\varrho})$. This suggests that the two additional errors introduced in the
structure-preserving scheme are either higher-order terms of $\Delta t$ or
correlated with the original convection-enhanced behavior. Then when
$\sigma^{2}$ has the same scale as a fixed $\Delta t$, the other two error
terms of the structure-preserving method for the original system with sub-
diffusion behavior when $\sigma^{2}$ has the same scale as a fixed $\Delta t$
will have the following analysis

$\displaystyle\lim_{t\to\infty}\frac{\rho_{2}^{S}(\Delta
t)}{\frac{\sigma^{2}}{2}+\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}}$
$\displaystyle=\lim_{t\to\infty}\frac{h_{1}\Delta
t\frac{\mathbb{E}[\int_{0}^{t}\bm{v}(\bm{x}_{s})\cdot\hat{\bm{\mu}}(s)ds]}{2t}+h_{2}\Delta
t^{\varrho}}{h_{3}\Delta t}$ $\displaystyle=0+O(\Delta t^{\varrho-1})<\infty.$
(100)

Here $h_{1},h_{2}$, and $h_{3}$ are all finite constants. While for the Euler-
Maruyama scheme, it does not hold this behavior. For $\rho_{2}^{E}(\Delta t)$
in a 2-dimensional case, it will become that

$\displaystyle\int_{0}^{t}$
$\displaystyle\int_{0}^{\infty}\langle\bm{v}(\bm{x}^{E}_{s})\cdot\tilde{\bm{v}}_{E}(\bm{x}^{E}_{s+\tau})\rangle
d\tau
ds=\frac{1}{2}\int_{0}^{t}\int_{0}^{\infty}\langle\Psi_{x_{1}x_{2}}(\bm{v}(\bm{x}^{E}_{s})\cdot\mathbf{D}_{2}\bm{v}(\bm{x}^{E}_{s+\tau}))\rangle
d\tau ds$
$\displaystyle+\frac{1}{2}\int_{0}^{t}\int_{0}^{\infty}\langle(\Psi_{x_{1}x_{1}}-\Psi_{x_{2}x_{2}})v_{1}(\bm{x}^{E}_{s})v_{2}(\bm{x}^{E}_{s+\tau}))\rangle
d\tau ds.$ (101)

Here

$\mathbf{D}_{2}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}.$ (102)

It can be seen here that the errors introduced by the modified system
corresponding to the Euler-Maruyama scheme cannot be transformed into a form
that is consistent with the original random velocity field’s convection-
enhanced behavior. $\square$

#### Proof of Theorem 4.3

###### Proof

We first compute the stream function for the structure-preserving scheme’s
modified equation system

$\Psi^{s}=\Psi-\frac{\Delta
t}{2}\Psi_{x_{1}}\Psi_{x_{2}}-\frac{\sigma^{2}}{4}\Delta\Psi.$ (103)

Then we consider the correlation function of the modified system,

$R_{ij}^{s}(\bm{r})=\langle\Psi_{x_{1}}(\bm{r})\Psi_{x_{1}}(0)+\Psi_{x_{2}}(\bm{r})\Psi_{x_{2}}(0)\rangle+\frac{\sigma^{4}\Delta
t^{2}}{16}\tilde{R_{1}}(\bm{r})+\frac{\Delta
t^{2}}{4}\tilde{R_{2}}(\bm{r})-\frac{\sigma^{2}\Delta
t}{4}\tilde{R_{3}}(\bm{r}).$ (104)

Here

$\displaystyle\tilde{R_{1}}(\bm{r})=\Big{\langle}$ $\displaystyle\sum_{i\neq
j}\Big{(}\big{(}\Psi_{x_{i}x_{i}x_{j}}(\bm{r})+\Psi_{x_{j}x_{j}x_{j}}(\bm{r})\big{)}\big{(}\Psi_{x_{i}x_{i}x_{j}}(0)+\Psi_{x_{j}x_{j}x_{j}}(0)\big{)}\Big{)}\Big{\rangle}.$
$\displaystyle\tilde{R_{2}}(\bm{r})=\Big{\langle}$ $\displaystyle\sum_{i\neq
j}\Big{(}\big{(}\Psi_{x_{i}}(\bm{r})\Psi_{x_{j}x_{j}}(\bm{r})+\Psi_{x_{j}}(\bm{r})\Psi_{x_{i}x_{j}}(\bm{r})\big{)}\cdot$
$\displaystyle\quad\quad\quad\big{(}\Psi_{x_{i}}(0)\Psi_{x_{j}x_{j}}(0)+\Psi_{x_{j}}(0)\Psi_{x_{i}x_{j}}(0)\big{)}\Big{)}\Big{\rangle}.$
$\displaystyle\tilde{R_{3}}(\bm{r})=\Big{\langle}$ $\displaystyle\sum_{i\neq
j}\Big{(}\Psi_{x_{i}}(0)\Psi_{x_{i}x_{i}x_{i}}(\bm{r})+\Psi_{x_{i}}(\bm{r})\Psi_{x_{i}x_{i}x_{i}}(0)$
$\displaystyle\quad\quad+\Psi_{x_{i}}(0)\Psi_{x_{i}x_{j}x_{j}}(\bm{r})+\Psi_{x_{i}}(\bm{r})\Psi_{x_{i}x_{j}x_{j}}(0)\Big{)}\Big{\rangle}.$
(107)

Here $i=1,2$ and $j=1,2$. Then we can compute that

$\displaystyle\frac{\sigma^{4}\Delta t^{2}}{16}\tilde{R_{1}}(\bm{r})$
$\displaystyle=\frac{\sigma^{4}\Delta
t^{2}}{16}\frac{1}{N}\sum_{i=1}^{N}\frac{1}{|\bm{k}_{i}|^{2}}(k_{i,1}^{6}+3k_{i,1}^{4}k_{i,2}^{2}+3k_{i,1}^{2}k_{i,2}^{4}+k_{i,2}^{6})\cos(\bm{k}_{i}\cdot\bm{r})$
$\displaystyle=\frac{\sigma^{4}\Delta
t^{2}}{16N}\sum_{i=1}^{N}|\bm{k}_{i}|^{4}\cos(\bm{k}_{i}\cdot\bm{r}).$ (108)

Due to the analysis in Majda1999 , as $N\to\infty$, this summation (108)
converges, more precisely

$\displaystyle\lim_{N\to\infty}\frac{\sigma^{4}\Delta
t^{2}}{16}\tilde{R_{1}}(\bm{r})=$ $\displaystyle\frac{\sigma^{4}\Delta
t^{2}}{16}\int_{0}^{L}\int_{0}^{2\pi}k^{2(3-\alpha)}cos(k|\bm{r}|\cos(\theta))d\theta
dk$ (109) $\displaystyle=$ $\displaystyle\frac{\sigma^{4}\Delta
t^{2}}{16}\int_{0}^{L}k^{2(3-\alpha)}J_{0}(k|\bm{r}|)dk.$

Here $J_{0}(\cdot)$ is the Bessel function of the first kind. Similar analysis
can be done for $\tilde{R_{2}}(\bm{r})$ and $\tilde{R_{3}}(\bm{r})$ and we can
obtain that

$\displaystyle\lim_{N\to\infty}\frac{\Delta t^{2}}{4}\tilde{R_{2}}(\bm{r})$
$\displaystyle\sim\lim_{N\to\infty}\frac{\Delta
t^{2}}{4N}\sum_{i=1}^{N}|\bm{k}_{i}|^{4}\cos(\bm{k}_{i}\cdot\bm{r})$
$\displaystyle=\frac{\Delta
t^{2}}{4}\int_{0}^{L}\int_{0}^{2\pi}k^{2(3-\alpha)}cos(k|\bm{r}|\cos(\theta))d\theta
dk$ (110) $\displaystyle=\frac{\Delta
t^{2}}{4}\int_{0}^{L}k^{2(3-\alpha)}J_{0}(k|\bm{r}|)dk.$
$\displaystyle\lim_{N\to\infty}\frac{\sigma^{2}\Delta
t}{4}\tilde{R_{3}}(\bm{r})$
$\displaystyle=\lim_{N\to\infty}\frac{\sigma^{2}\Delta
t}{2N}\sum_{i=1}^{N}|\bm{k}_{i}|^{2}\cos(\bm{k}_{i}\cdot\bm{r})$
$\displaystyle=\frac{\sigma^{2}\Delta
t}{2}\int_{0}^{L}\int_{0}^{2\pi}k^{2(2-\alpha)}cos(k|\bm{r}|\cos(\theta))d\theta
dk$ $\displaystyle=\frac{\sigma^{2}\Delta
t}{2}\int_{0}^{L}k^{2(2-\alpha)}J_{0}(k|\bm{r}|)dk.$

Then we can observe that for further analysis we need to estimate the moments
of the Bessel function of the first kind. Following the analysis in
Wang2016NumericalAF  and MILNE-THOMSON1945 , we can obtain that

$\displaystyle\int_{0}^{L}k^{\beta}J_{0}(k|\bm{r}|)dk=$
$\displaystyle\frac{L}{|\bm{r}|^{\beta}}\big{[}(\beta-1)J_{0}(L|\bm{r}|)s^{(2)}_{\beta-1,-1}(L|\bm{r}|)+J_{1}(L|\bm{r}|)s^{(2)}_{\beta,0}(L|\bm{r}|)\big{]}$
$\displaystyle+\frac{2^{\beta}\Gamma(\frac{\beta+1}{2})}{|\bm{r}|^{\beta+1}\Gamma(\frac{1-\beta}{2})}.$
(111)

For $x$ large enough, the Lommel function $s^{(2)}_{n,m}(x)$ has an asymptotic
expansion MILNE-THOMSON1945 ,

$s^{(2)}_{n,m}(x)=x^{n-1}\big{(}1-\frac{(n-1)^{2}+m^{2}}{x^{2}}+\frac{\big{(}(n-1)^{2}-m^{2}\big{)}\big{(}(n-3)^{2}-m^{2}\big{)}}{x^{4}}-\cdots\big{)}.$
(112)

In this case, the Bessel function also has an asymptotic function for $x$
large enough

$J_{\alpha}(x)=\sqrt{\frac{2}{\pi
x}}(\cos(x-\frac{\alpha\pi}{2}-\frac{\pi}{4})+O(\frac{1}{x})).$ (113)

Substitute (112) and (113) back into (111), and we can obtain that as
$|\bm{r}|\to\infty$, the leading term in (111) will become

$\int_{0}^{L}k^{\beta}J_{0}(k|\bm{r}|)dk\sim\left\\{\begin{aligned}
&\frac{1}{|\bm{r}|^{\beta+1}},0<\beta<\frac{1}{2}\\\
&\frac{1}{|\bm{r}|^{\frac{3}{2}}},\quad\beta\geq\frac{1}{2}\end{aligned}\right.$
(114)

Then we go back to the analysis of
$\tilde{R_{1}}(\bm{r}),\tilde{R_{2}}(\bm{r})$ and $\tilde{R_{3}}(\bm{r})$, due
to $\frac{1}{4}<\alpha<1$, so we can obtain that for $|\bm{r}|$ large enough

$\frac{\sigma^{4}\Delta t^{2}}{16}\tilde{R_{1}}(\bm{r})+\frac{\Delta
t^{2}}{4}\tilde{R_{2}}(\bm{r})-\frac{\sigma^{2}\Delta
t}{4}\tilde{R_{3}}(\bm{r})\sim\frac{O(\sigma^{2}\Delta t+\Delta
t^{2})}{|\bm{r}|^{\frac{3}{2}}}.$ (115)

Together with the original correlation function $R_{i,j}(\bm{r})$, we can
obtain that when $\frac{1}{4}\leq\alpha<1$, the asymptotic performance of
$R^{s}_{i,j}(\bm{r})$ is consistent with the original system, i.e.

$R^{s}_{i,j}(\bm{r})\sim\frac{1}{|\bm{r}|^{2(1-\alpha)}},\
\frac{1}{4}\leq\alpha<1.$ (116)

Then from Proposition 5, we can obtain that for $t$ large enough,

$\mathbf{E}[(\bm{x}_{t}^{S})^{2}]\sim t^{\frac{2}{2-\alpha}},\
\frac{1}{4}\leq\alpha<1.$ (117)

$\square$

##  5 Numerical experiments

In this section, we present several numerical results computed by the
structure-preserving scheme. To start with, we first discuss how to generate
random velocity fields $\bm{v}(\bm{x},t)$. Notice that we have suppressed the
dependence of the velocity on $\omega$ in $\bm{v}(\bm{x},t)$ for notation
simplicity here. In our numerical experiments, we will consider the particles
moving in an isotropic velocity field satisfies Assumption 2 and 2.1. And we
have mentioned that the covariance functions of $\bm{v}(\bm{x},t)$ have the
forms like (10) and (9). We take the time-correlation function
$D(t)=\exp(-\frac{1}{2}\theta_{0}^{2}t^{2})$, where $\theta_{0}\geq 0$
represents the strength of the decorrelation in time $t$. As for the energy
spectrum function $E(k)$, we will consider the spectrums with the following
forms:

$\displaystyle E_{1}(k)=\delta(k-k_{0}),\ \text{2D flows},$ (118)
$\displaystyle
E_{2}(k)=\frac{9}{2}k^{3}k_{0}^{-4}\exp(-\frac{3}{2}k^{2}/k_{0}^{2}),\
\text{2D flows},$ (119) $\displaystyle E_{3}(k)=\frac{3}{2}\delta(k-k_{0}),\
\text{3D flows},$ (120) $\displaystyle
E_{4}(k)=16\sqrt{\frac{2}{\pi}}k^{4}k_{0}^{-5}\exp(-2k^{2}/k_{0}^{2}),\
\text{3D flows},$ (121)

with a fixed $k_{0}$. Here $E_{1}$ and $E_{3}$ describe excitation confined to
a thin shell, in two- and three- dimensional spaces, respectively, while
$E_{2}$ and $E_{4}$ describe less concentrated excitations, in two- and three-
dimensional spaces, respectively; see Kraichnan1970 .

For the energy spectrum functions $E_{1}$ and $E_{3}$, we choose $\bm{k}_{n}$
isotropically distributed on the surface of a sphere and a circle of radius
$k_{0}$, respectively. For the energy spectrum function $E_{2}$, each
component of $\bm{k}_{n}$ is chosen from a Gaussian distribution of standard
deviation $k_{0}/\sqrt{3}$, while for the energy spectrum function $E_{4}$,
each component of $\bm{k}_{n}$ is chosen from a Gaussian distribution of
standard deviation $k_{0}/2$. Finally, $\theta_{n}$ in (11) is chosen from a
Gaussian distribution with standard deviation $\theta_{0}$. In the remaining
part of this paper, we will use $D_{0}$ as a parameter in our numerical
experiments, which is equivalent to $\sigma^{2}/2$.

###  5.1 Comparison between structure-preserving and Euler-Maruyama

In this subsection, we compare the effective diffusivity computed by the
structure-preserving scheme and Euler-Maruyama scheme, where the velocity
field is generated by (11) with the setting (2.2). Here $E(k)=E_{1}(k)$
defined in (118). Theorem 1 of Fannjiang1996  has explained why $\theta_{0}=0$
and $D_{0}=0$ leads to anomalous diffusion phenomena. The occurrence of the
trapping phenomena results in the random velocity field with this setting
shall exhibit the sub-diffusion performance.

Figure 1: Comparison of $D^{11}_{E}$ by two numerical schemes with different
time step $\Delta t$ in 2D random field with $E_{1}(k)$ and
$\theta_{0}=0,D_{0}=0$ Figure 2: Comparison of $D^{11}_{E}$ by two numerical
schemes with different $D_{0}$ in 2D random field with $E_{1}(k)$ and
$\theta_{0}=0,\Delta t=0.1$ Figure 3: Comparison for dispersion
$\mathbb{E}[\bm{x}_{t}^{2}]$, by Log-log plot in 2D random field with
$E_{1}(k)$ and $\theta_{0}=0,D_{0}=0$

From Figure 1 we can see that, the structure-preserving scheme is not
sensitive to the time step $\Delta t$ and their results are very close. they
both converge to zero and show the trapping phenomena. While the Euler-
Maruyama scheme is sensitive to the time step $\Delta t$ and as we decrease
$\Delta t$, it gets closer to the reference result but still converges to a
non-zero limit and does not show the sub-diffusion performance. To get a
relatively close result, the Euler-Maruyama scheme needs a time step which is
much smaller than the structure-preserving scheme. Then we also test for the
performance of the Euler-Maruyama scheme under a long time scale and at this
time, the Euler-Maruyama scheme will no longer perform sub-diffusion
phenomena. This is due to the numerical diffusion phenomena caused by the
Euler-Maruyama scheme itself. This is also well illustrated in Figure 2. We
can see that when $\Delta t$ is fixed, the effective diffusivity calculated by
the Euler-Maruyama method does not decrease with decreasing $D_{0}$. This
phenomenon does not occur in the structure-preserving scheme. This explains
that the Euler-Maruyama scheme is very sensitive to the scale relationship
between $\Delta t$ and $D_{0}$. When $D_{0}\ll\Delta t$, this scheme
drastically changes the behavior of the effective diffusivity. This also
supports our Theorem 4.2. To more conveniently tell the diffusion behavior of
the random flows, we can observe the dispersion $\mathbb{E}[\bm{x}_{t}^{2}]$
at time $t$. Then Figure 3 illustrates that even for a relatively small
$\Delta t$, the Euler-Maruyama scheme will approach a diffusion phenomenon in
the long time behavior, while the structure-preserving scheme can keep the
sub-diffusion phenomenon which is consistent with the theoretical result of
the original system.

###  5.2 Computing effective diffusivity of 2D and 3D random flows

In this subsection, we investigate the convection-enhanced diffusion phenomena
and compute the effective diffusivity of two-dimensional incompressible random
flows. Here we choose a time step $\Delta t=0.05$, number of Monte Carlo
samples $N_{mc}=100,000$ and number of Fourier modes $N=200$ in approximating
the random velocity field. For the velocity field generated by (11) according
to $E_{1}(k)$ defined in (118) and $E_{2}(k)$ defined in (119), we compute
their effective diffusivity with different $\theta_{0}$ and $D_{0}$ settings.

(a) $\theta_{0}=0$

(b) $\theta_{0}=1$

Figure 4: Effective diffusivity $D_{E}^{11}$ v.s. time $t$ in 2D random
velocity field with $E_{1}(k)$

From Figure 4 , we can see that the effective diffusivity all converges to a
non-zero limit except in the case when $\theta_{0}$ and $D_{0}$ are both equal
to $0$. This implies that for $E_{1}(k)$ if $D_{0}>0$, its corresponding
velocity field will perform normal diffusion phenomena. The energy spectrum
density $E_{2}(k)$ has a similar performance as $E_{1}(k)$.

(a) $\theta_{0}=0$

(b) $\theta_{0}=1$

Figure 5: Effective diffusivity $D_{E}^{11}$ v.s. time $t$ in 3D random
velocity field with $E_{3}(k)$

From Figure 5, we can observe that the effective diffusivity all converges to
a non-zero limit in this 3D case. This implies that for $E_{3}(k)$, its
corresponding velocity field will perform normal diffusion phenomena. The
energy spectrum $E_{4}(k)$ performs a similar result to $E_{3}(k)$.

###  5.3 Anomalous diffusion phenomena in 2D random flows

In this subsection, we investigate and compute some 2D random velocity fields
which will perform anomalous diffusion phenomena. With the condition (14), we
can find a $E(k)$ series that does not satisfy this condition.

$\displaystyle
E_{5}(k)=\sqrt{\frac{6}{\pi}}k_{0}^{-1}\cdot\exp(-\frac{3k^{2}}{2k_{0}^{2}}),$
(122) $\displaystyle
E_{6}(k)=\frac{\sqrt[4]{54}}{\Gamma(\frac{3}{4})}k_{0}^{-1.5}k^{0.5}\cdot\exp(-\frac{3k^{2}}{2k_{0}^{2}}).$
(123)

Then we plot their dispersion computed by the structure-preserving scheme.
Here we set $\theta_{0}=0$ and other parameters remain the same.

Figure 6: Log-log plot of Dispersion $\mathbb{E}[\bm{x}_{t}^{2}]$ v.s. time
$t$ in 2D random velocity field with $E_{5}(k)$ Figure 7: Log-log plot of
Dispersion $\mathbb{E}[\bm{x}_{t}^{2}]$ v.s. time $t$ in 2D random velocity
field with $E_{6}(k)$

From Figure 6, and 7, we can see that they both exhibit anomalous diffusion
phenomena, specifically super-diffusion phenomena. This is due to the
diffusivity is always enhanced in incompressible flows and the computation
results are consistent with the theoretical sharp condition (14). Details can
be found in Fannjiang1998 .

Then, as mentioned in Fannjiang1998 , we calculate the power law between space
and time, $\mu$, of the velocity fields with the correlation functions (73) by
two different numerical schemes and compare them with the theoretical results.
From Figure 8, we can see that the results of both numerical schemes match the
theoretical results well when alpha is relatively large and close to $1$.
Meanwhile, when $0.25<\alpha<0.5$, the power law computed by the structure-
preserving scheme is more in line with the theoretical results. This also
coincides with our proof of Theorem 4.3.

Figure 8: $\mu$ v.s. $\alpha$, comparison for the power law computed by two
numerical schemes, $\mathbb{E}[\bm{x}_{t}^{2}]\sim t^{\mu}$, in a 2D random
field with spectrum $\hat{R}_{i,j}(\bm{k})$ in (73).

###  5.4 Anomalous diffusion phenomena in 3D random flows

In this subsection, we investigate and compute some 3D random velocity fields
which will perform anomalous diffusion phenomena. With the sharp condition
(16), we can find the $E_{7}(k)$ that does not satisfy this condition as a 3D
case.

$E_{7}(k)=6k_{0}^{-2}k\cdot\exp(-\frac{2k^{2}}{k_{0}^{2}}).$ (124)

We can compute its corresponding velocity field’s effective diffusivity and
dispersion with the same parameters.

Figure 9: Log-log plot of Dispersion $\mathbb{E}[\bm{x}_{t}^{2}]$ v.s. time
$t$ in 3D random velocity field with $E_{7}(k)$

From Figure 9, we can see that the 3D random velocity field with energy
spectrum $E_{7}(k)$ performs super-diffusion phenomena which is consistent
with the theoretical condition. We also conduct experiments to calculate the
power law, $\mu$, of the random field mentioned in Fannjiang1998  in both
numerical schemes in the three-dimensional cases. Now the random velocity
field will have the following spectrum

$\hat{R}_{i,j}(\bm{k})\sim\frac{1}{|\bm{k}|^{2\alpha+1}}(\delta_{i,j}-\frac{k_{i}k_{j}}{|\bm{k}|^{2}}).$
(125)

Figure 10: $\mu$ v.s. $\alpha$, comparison for the power law computed by two
numerical schemes, $\mathbb{E}[\bm{x}_{t}^{2}]\sim t^{\mu}$, in a 3D random
field with spectrum $\hat{R}_{i,j}(\bm{k})$ in (125).

From Figure 10, we can see that, when $\alpha$ is not so small, i.e.
$\alpha>0.25$, the power law $\mu$ computed by the structure-preserving scheme
matches the theoretical result better than the Euler-Maruyama scheme.

##  6 Conclusions

In this paper, we developed an efficient structure-preserving scheme to
compute the effective diffusivity $D^{11}_{E}$ in 2D and 3D stochastic flows,
with convergence analysis based on the Markov semigroup and modified flow
system. We also compare the structure-preserving scheme with the Euler-
Maruyama scheme and prove that the structure-preserving scheme can keep the
original diffusion behavior of the system unchanged in the long time scale
while the Euler-Maruyama method cannot. We also conducted numerical
experiments to address this point, and the experimental results are consistent
with it. For a 2D random field which should exhibit sub-diffusion phenomena
when diffusion vanishes, the structure-preserving scheme maintains this
property, while the Euler-Maruyama scheme exhibits diffusion phenomena in the
long time scale. We also computed the effective diffusivity for some random
fields with diffusion phenomena in 2D and 3D. Then based on the criterion in
Fannjiang1996 , we derived and discovered a series of 2D and 3D random fields
that would perform anomalous diffusion phenomena. The corresponding numerical
experiments also yielded consistent results. Meanwhile, we also prove that for
random flows that follow the spectrums (73), (125) mentioned in Fannjiang1998
, the structure-preserving scheme can keep the original super-diffusive
behavior. Compared to the Euler-Maruyama scheme and other Eulerian methods,
our structure-preserving particle method has the following advantages: (1) It
is scalable in spatial dimension and can handle stochastic 3D problems at
ease; (2) It can keep the diffusion behavior of the original system in the
long time scale and is not very sensitive to the time step.

###### Acknowledgements.

JX was partially supported by NSF grant DMS-2309520. ZZ was supported by Hong
Kong RGC grant (Projects 17300318 and 17307921), National Natural Science
Foundation of China (Project 12171406), Seed Funding Programme for Basic
Research (HKU) and Seed Funding for Strategic Interdisciplinary Research
Scheme 2021/22 (HKU).

##  Appendix A Some formulations and spectral gap conditions in random flows

We first define a function space that holds the stationary and ergodic
properties in $\mathcal{X}$. Then let $(\mathcal{X},\mathcal{H},P_{0})$ be a
probability space. Here $\mathcal{H}$ is a $\sigma$-algebra over
$\mathcal{X}$, and $P_{0}$ is a probability measure. Then let $\tau_{\bm{x}}$,
$\bm{x}\in\mathbb{R}^{d}$ be an almost surely continuous, jointly measurable
group of measure-preserving transformation on $\mathcal{X}$ with the following
properties:

  1. P1.

$\tau_{\bm{0}}=Id_{\mathcal{X}}$, and
$\tau_{\bm{x}+\bm{y}}=\tau_{\bm{x}}\tau_{\bm{y}},\forall\
\bm{x},\bm{y}\in\mathbb{R}^{d}$.

  2. P2.

The mapping $(\tilde{\bm{x}},\bm{x})\to\tau_{\bm{x}}\tilde{\bm{x}}$ is jointly
measurable.

  3. P3.

$P_{0}(\tau_{\bm{x}}A)=P_{0}(A),\forall\bm{x}\in\mathbb{R}^{d},A\in\mathcal{H}$.

  4. P4.

$\lim_{\bm{x}\to\bm{0}}P_{0}(\tilde{\bm{x}}:|f\circ\tau_{\bm{x}}(\tilde{\bm{x}})-f(\tilde{\bm{x}})|\geq\epsilon)=0,\
\forall f\in L^{2}(\mathcal{X}),\epsilon>0$.

  5. P5.

If $P_{0}(A\Delta\tau_{\bm{x}}A)=0,\forall\bm{x}\in\mathbb{R}^{d}$, then $A$
is a trivial event.

Then $\tau_{\bm{x}}$ induces a strongly continuous group of unitary mapping
$U^{\bm{x}}$ on $L^{2}(\mathcal{X})$, which satisfies

$U^{\bm{x}}f(\tilde{\bm{x}})=f(\tau_{\bm{x}}(\tilde{\bm{x}})),\ f\in
L^{2}(\mathcal{X}),\ \bm{x}\in\mathbb{R}^{d}.$ (126)

Next, we include the time variable $t$ and investigate the Markovian property.
Let $\Omega$ be the space of $\mathcal{X}$-valued continuous function
$C([0,\infty);\mathcal{X})$ and let $\mathcal{B}$ be its Borel
$\sigma$-algebra. Let $P^{t},t\geq 0$, be a strongly continuous Markov
semigroup on $L^{2}(\mathcal{X})$, which satisfies the following properties.

  1. Q1.

$P^{t}\textbf{1}=\textbf{1}$ and $P^{t}f\geq 0$ if $f\geq 0$.

  2. Q2.

$\int P^{t}fdP_{0}=\int fdP_{0}$, for all $f\in L^{2}(\mathcal{X}),t\geq 0$.

  3. Q3.

$\textbf{E}_{\tilde{\bm{x}}}[f(\theta_{t+h}(\omega))|\mathcal{B}_{\leq
t}]=P^{h}F(\omega(t))$, where
$F(\tilde{\bm{x}}):=\textbf{E}_{\tilde{\bm{x}}}f$ for any $f\in
L^{1}(\Omega)$, $t,h\geq 0,\tilde{\bm{x}}\in\mathcal{X}$.

Here $\mathcal{B}_{\leq t}$ are the $\sigma$-algebras generated by measurable
events preceding time $t$, and
$\theta_{t}(\omega(\cdot))=\omega(\cdot+t),t\geq 0$ is an standard shift
operator on the path space $(\Omega,\mathcal{B})$. Then a stationary measure
$P$ can be defined on the path space $(\Omega,\mathcal{B})$ by

$P(B)=\int P_{\tilde{\bm{x}}}(B)P_{0}(d\tilde{\bm{x}}),B\in\mathcal{B},$

and we use $\textbf{E},\ \textbf{E}_{\tilde{\bm{x}}}$ to denote the
corresponding expectation operators taken over the measure $P,\
P_{\tilde{\bm{x}}}$. Here $P_{\tilde{\bm{x}}}$ can be considered as the
conditional probability for all events in $\mathcal{B}$ condition with the
initial point lies on $\tilde{\bm{x}}$. We denote $\theta_{t}$ as a standard
shift operator on the path space $(\Omega,\mathcal{B})$.

###### Proposition 6

$P$ is invariant under the action of $\theta_{t}$ and $\tau_{\bm{x}}$ for any
$(t,\bm{x})\in\mathbb{R}^{+}\times\mathbb{R}^{d}$.

Let $L:\mathcal{D}(L)\to L^{2}(\mathcal{X})$ be the generator of the semigroup
$P^{t}$. We assume that the generator $L$ satisfies the following spectral gap
condition

$-(Lf,f)_{L^{2}(\mathcal{X})}\geq c_{1}||f||^{2}_{L^{2}(\mathcal{X})},\
c_{1}>0.$ (127)

The spectral gap condition (127) guarantees the exponential decay property

$||P^{t}f||_{L^{2}(\mathcal{X})}\leq\exp(-2c_{1}t)||f||_{L^{2}(\mathcal{X})},\
c_{1}>0,\ f\in L_{0}^{2}(\mathcal{X}).$ (128)

See Doukhan1994 , Rosen2012 . The exponential decay property (128) for $P^{t}$
is equivalent to $\varrho$-mixing of the process $X(t),t\geq 0$ and
contributes significantly to proving the existence of the effective
diffusivity.

##  Appendix B The behavior of the stream function $\Psi$ in two numerical
schemes

###### Lemma 3

Let $\Psi$ defined in (38) be the stream function of a 2D divergence-free
random velocity field. We can treat the mean value of the Hamiltonian $\Psi$
as a function of time $t$. Then it can be obtained that

$\frac{d\Psi}{dt}=\frac{\sigma^{2}}{2}\Delta\Psi+\sigma\nabla\Psi\cdot\frac{d\bm{w}_{t}}{dt}.$
(129)

###### Proof

We apply Ito’s formula to the stream function $\Psi$ for $\bm{x}$ satisfies
the original SDEs (1), to find

$d\Psi=[(\nabla\Psi)\cdot\bm{v}+\frac{1}{2}\text{Tr}((\sigma\mathbf{I}_{2})(\mathbf{H}_{\Psi})(\sigma\mathbf{I}_{2}))]dt+(\nabla\Psi)\cdot(\sigma\mathbf{I}_{2})d\bm{w}_{t},$
(130)

here $\mathbf{H}_{\Psi}$ is the Hessian matrix of $\Psi$ with respect to
$\bm{x}$. Then due to $\bm{v}=\nabla^{\perp}\Psi$ is divergence-free, we can
compute that

$(\nabla\Psi)\cdot\bm{v}=(\nabla\Psi)\cdot(\nabla^{\perp}\Psi)=0,\quad\frac{1}{2}\text{Tr}((\sigma\mathbf{I}_{2})(\mathbf{H}_{\Psi})(\sigma\mathbf{I}_{2}))=\frac{\sigma^{2}}{2}\Delta\Psi.$
(131)

Then we can obtain that

$\frac{d\Psi}{dt}=\frac{\sigma^{2}}{2}\Delta\Psi+\sigma\nabla\Psi\cdot\frac{d\bm{w}_{t}}{dt}.$
(132)

The integral of $\sigma\nabla\Psi\cdot\frac{d\bm{w}_{t}}{dt}$ is a mean zero
martingale, i.e. $\mathbf{E}[\int_{0}^{t}\sigma\nabla\Psi\cdot
d\bm{w}_{t}]=0$.

$\square$

###### Lemma 4

Let $\Psi^{E}$ and $\Psi^{S}$ be the stream function of a 2D divergence-free
random velocity field satisfying the modified equation systems of Euler-
Maruyama scheme and structure-preserving scheme correspondingly. The Euler-
Maruyama method will introduce an extra numerical diffusion term in the
behavior of the mean value of $\Psi^{E}$, i.e.

$\displaystyle\frac{d\Psi^{E}}{dt}$ $\displaystyle=\frac{\Delta
t}{2N\sqrt{N}}\sum_{n=1}^{N}M_{n}\Delta\Psi_{n}+\frac{\sigma^{2}}{2}M_{\sigma}\Delta\Psi^{E}+\nabla\Psi^{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt}.$
(133)

Here

$\displaystyle M_{n}$
$\displaystyle=\frac{\big{(}\sum_{i=1}^{N}(k_{n,1}k_{i,2}-k_{n,2}k_{i,1})(\xi_{i}\cos(\bm{k}_{i}\cdot\bm{x})+\zeta_{i}\sin(\bm{k}_{i}\cdot\bm{x}))\big{)}^{2}}{k_{n,1}^{2}+k_{n,2}^{2}}\geq
0,$ $\displaystyle M_{\sigma}$ $\displaystyle=1+\frac{3\Delta
t^{2}}{4}\Psi_{x_{1}x_{2}}^{2}+\frac{\Delta
t^{2}}{4}(\Psi_{x_{1}x_{1}}^{2}-\Psi_{x_{1}x_{1}}\Psi_{x_{2}x_{2}}+\Psi_{x_{2}x_{2}}^{2})\geq
0,$

and the integral of $\nabla\Psi^{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt}$
is a mean zero martingale, i.e.
$\mathbb{E}[\int_{0}^{t}\nabla\Psi^{E}\cdot\mathbf{G}_{E}d\bm{w}_{t}]$ $=0$.
Compared with (129) in Lemma 3, the Euler-Maruyama scheme will introduce an
additional numerical diffusion term. While the structure-preserving scheme
will not introduce this error term, i.e.

$\displaystyle\frac{d\Psi^{S}}{dt}$
$\displaystyle=\frac{\sigma^{2}}{2}M_{\sigma}\Delta\Psi^{S}+\nabla\Psi^{S}\cdot\mathbf{G}_{S}\frac{d\bm{w}_{t}}{dt}.$
(134)

###### Proof

Inspired by the idea in PAVLIOTIS20091030 , we apply Ito’s formula to the
stream function $\Psi^{E}(\bm{x})$ for the random flow for $\bm{x}$ satisfies
the modified equation (54) of the Euler-Maruyama method. We can find that

$d\Psi^{E}=[(\nabla\Psi^{E})\cdot(\bm{v}+\Delta
t\tilde{\bm{v}}_{E})+\frac{\sigma^{2}}{2}\text{Tr}(\mathbf{G}_{E}^{T}(\mathbf{H}_{\Psi})\mathbf{G}_{E})]dt+(\nabla\Psi^{E})\cdot\mathbf{G}_{E}d\bm{w}_{t}.$
(135)

Then we can compute that,

$\displaystyle(\nabla\Psi^{E})\cdot($ $\displaystyle\bm{v}+\Delta
t\tilde{\bm{v}}_{E})=\frac{\Delta
t}{2}(\Psi_{x_{1}}^{2}\Psi_{x_{2}x_{2}}+\Psi_{x_{2}}^{2}\Psi_{x_{1}x_{1}}-2\Psi_{x_{1}}\Psi_{x_{2}}\Psi_{x_{1}x_{2}})$
$\displaystyle\quad+\frac{\sigma^{2}\Delta
t}{4}(\Psi_{x_{1}}\Psi_{x_{1}x_{1}x_{2}}+\Psi_{x_{1}}\Psi_{x_{2}x_{2}x_{2}}-\Psi_{x_{2}}\Psi_{x_{1}x_{2}x_{2}}-\Psi_{x_{2}}\Psi_{x_{1}x_{1}x_{1}}).$
(136) $\displaystyle\frac{\sigma^{2}}{2}\text{Tr}(\mathbf{G}_{E}^{T}$
$\displaystyle(\mathbf{H}_{\Psi})\mathbf{G}_{E})=\frac{\sigma^{2}}{2}M_{2}\Delta\Psi^{E}=-\frac{\sigma^{2}k_{0}^{2}}{2}M_{2}\Psi^{E}.$
(137)

Here

$M_{2}=1+\frac{3\Delta t^{2}}{4}\Psi_{x_{1}x_{2}}^{2}+\frac{\Delta
t^{2}}{4}(\Psi_{x_{1}x_{1}}^{2}-\Psi_{x_{1}x_{1}}\Psi_{x_{2}x_{2}}+\Psi_{x_{2}x_{2}}^{2})\geq
1.$

And for the two terms in (136) for the case that the flow with the stream
function (38), we can obtain that

$\displaystyle\frac{\Delta
t}{2}(\Psi_{x_{1}}^{2}\Psi_{x_{2}x_{2}}+\Psi_{x_{2}}^{2}\Psi_{x_{1}x_{1}}-2\Psi_{x_{1}}\Psi_{x_{2}}\Psi_{x_{1}x_{2}})=$
$\displaystyle\frac{-\Delta
t}{2N\sqrt{N}}\sum_{n=1}^{N}\big{(}\sum_{i=1}^{N}(k_{n,1}k_{i,2}-k_{n,2}k_{i,1})(\xi_{i}\cos(\bm{k}_{i}\cdot\bm{x})+\zeta_{i}\sin(\bm{k}_{i}\cdot\bm{x}))\big{)}^{2}\Psi_{n}.$
(138)

And we can easily verify that

$\frac{\sigma^{2}\Delta
t}{4}(\Psi_{x_{1}}\Psi_{x_{1}x_{1}x_{2}}+\Psi_{x_{1}}\Psi_{x_{2}x_{2}x_{2}}-\Psi_{x_{2}}\Psi_{x_{1}x_{2}x_{2}}-\Psi_{x_{2}}\Psi_{x_{1}x_{1}x_{1}})=0.$
(139)

Then together with (137), (138) and (139), we can obtain that

$\displaystyle\frac{d\Psi^{E}}{dt}$ $\displaystyle=\frac{\Delta
t}{2N\sqrt{N}}\sum_{n=1}^{N}M_{n}\Delta\Psi_{n}+\frac{\sigma^{2}}{2}M_{2}\Delta\Psi^{E}+\nabla\Psi^{E}\cdot\mathbf{G}_{E}\frac{d\bm{w}_{t}}{dt}$
(140)

where $\mathbf{E}[\int_{0}^{t}\nabla\Psi^{E}\cdot\mathbf{G}_{E}d\bm{w}_{t}]=0$
and

$M_{n}=\big{(}\sum_{i=1}^{N}(k_{n,1}k_{i,2}-k_{n,2}k_{i,1})(\xi_{i}\cos(\bm{k}_{i}\cdot\bm{x})+\zeta_{i}\sin(\bm{k}_{i}\cdot\bm{x}))\big{)}^{2}\Big{/}(k_{n,1}^{2}+k_{n,2}^{2})\geq
0$ (142)

A similar analysis can be done for the structure-preserving scheme, and then
we can finally obtain that

$\displaystyle\frac{d\Psi^{S}}{dt}$
$\displaystyle=\frac{\sigma^{2}}{2}M_{2}\Delta\Psi^{S}+\nabla\Psi^{S}\cdot\mathbf{G}_{S}\frac{d\bm{w}_{t}}{dt}$
(143)

$\square$

###### Remark 2

The first term on the right-hand side of (133) does not depend on the
diffusion coefficient $\sigma$ but on the time step $\Delta t$. When
$\sigma^{2}\ll\Delta t$, it drastically changes the behavior of
$\mathbf{E}[\Psi(t)]$ compared with (129).

###### Remark 3

For energy spectral densities that are distributed on the same unit circle as
$E_{1}(k)=\delta(k-k_{0})$, we have $k_{n,1}^{2}+k_{n,2}^{2}=k_{0}^{2}$. Then
we can further obtain that

$\frac{d\Psi}{dt}=-\frac{\sigma^{2}k_{0}^{2}}{2}\Psi+\sigma\nabla\Psi\cdot\frac{d\bm{w}_{t}}{dt},$
(144)

which means the mean value of the stream function, $\mathbf{E}[\Psi(t)]$ will
decay like $e^{-2t/(\sigma^{2}k_{0}^{2})}$, i.e. $\mathbf{E}[\Psi(t)]\approx
e^{-2t/(\sigma^{2}k_{0}^{2})}\mathbf{E}[\Psi(0)]$.

## References

  * (1) Bensoussan, A., Lions, J.L., Papanicolaou, G.: Asymptotic analysis for periodic structures, vol. 374.  American Mathematical Soc. (2011)
  * (2) Carmona, R.A., Xu, L.: Homogenization for time-dependent two-dimensional incompressible Gaussian flows.  The Annals of Applied Probability 7(1), 265 – 279 (1997)
  * (3) Doukhan, P.: Mixing: Properties and Examples.  Springer Berlin, Heidelberg (1994)
  * (4) Fannjiang, A.: Anomalous Diffusion in Random Flows, pp. 81–99.  Springer New York, New York, NY (1998)
  * (5) Fannjiang, A., Komorowski, T.: Turbulent diffusion in markovian flows.  The Annals of Applied Probability 9(3), 591–610 (1999)
  * (6) Fannjiang, A., Papanicolaou, G.: Diffusion in turbulence.  Probability Theory and Related Fields 105(3), 279 – 334 (1996)
  * (7) Feng, K., Shang, Z.: Volume-preserving algorithms for source-free dynamical systems.  Numerische Mathematik 71(4), 451–463 (1995)
  * (8) Feng, K., Wu, H., Qin, M., Wang, D.: Construction of canonical difference schemes for hamiltonian formalism via generating functions.  Journal of Computational Mathematics 7(1), 71–96 (1989)
  * (9) Garnier, J.: Homogenization in a periodic and time-dependent potential.  SIAM Journal on Applied Mathematics 57(1), 95–111 (1997)
  * (10) Hairer, E., Wanner, G., Lubich, C.: Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations.  Springer Berlin, Heidelberg (2006)
  * (11) Jikov, V.V., Kozlov, S.M., Oleinik, O.A.: Homogenization of Differential Operators and Integral Functionals.  Springer Berlin, Heidelberg (1994)
  * (12) Komorowski, T., Olla, S.: On homogenization of time-dependent random flows.  Probability Theory and Related Fields 121(1), 98–116 (2001)
  * (13) Kraichnan, R.H.: Diffusion by a Random Velocity Field.  The Physics of Fluids 13(1), 22–31 (1970)
  * (14) Landim, C., Olla, S., Yau, H.T.: Convection–diffusion equation with space–time ergodic random flow.  Probability Theory and Related Fields 112(2), 203–220 (1998)
  * (15) Lyu, J., Wang, Z., Xin, J., Zhang, Z.: Convergence analysis of stochastic structure-preserving schemes for computing effective diffusivity in random flows.  SIAM Journal on Numerical Analysis 58(5), 3040–3067 (2020)
  * (16) Majda, A., Kramer, P.: Simplified models for turbulent diffusion: Theory, numerical modelling, and physical phenomena.  Physics Reports 314(4), 237–574 (1999)
  * (17) McLachlan, R.I., Quispel, G.R.W.: Splitting methods.  Acta Numerica 11, 341–434 (2002)
  * (18) Milne-Thomson, L.: A treatise on the theory of bessel functions.  Nature 156(3955), 190–191 (1945)
  * (19) Pavliotis, G.A., Stuart, A.M.: Multiscale Methods: Averaging and Homogenization.  Springer New York, NY (2008)
  * (20) Pavliotis, G.A., Stuart, A.M., Zygalakis, K.C.: Calculating effective diffusivities in the limit of vanishing molecular diffusion.  Journal of Computational Physics 228(4), 1030–1055 (2009)
  * (21) Peter E. Kloeden, E.P.: Numerical Solution of Stochastic Differential Equations.  Springer Berlin, Heidelberg (1992)
  * (22) Rosenblatt, M.: Markov Processes, Structure and Asymptotic Behavior.  Springer Berlin, Heidelberg (2012)
  * (23) Taylor, G.I.: Diffusion by continuous movements.  Proceedings of the London Mathematical Society s2-20(1), 196–212 (1922)
  * (24) Wang, Y., Li, Y., Luo, J.: Numerical analysis for the moments of bessel functions and bessel-trigonometric functions.  arXiv: Numerical Analysis (2016)