Convergence analysis of a Lagrangian numerical scheme in computing effective diffusivity of 3D time-dependent flows

Zhongjian Wang zhongjian@statistics.uchicago.edu Jack Xin jxin@math.uci.edu Zhiwen Zhang zhangzw@hku.hk Department of Statistics, The University of Chicago, Chicago, IL 60637, USA. Department of Mathematics, University of California at Irvine, Irvine, CA 92697, USA. Department of Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong SAR, China.

Abstract

In this paper, we study the convergence analysis for a robust stochastic structure-preserving Lagrangian numerical scheme in computing effective diffusivity of time-dependent chaotic flows, which are modeled by stochastic differential equations (SDEs). Our numerical scheme is based on a splitting method to solve the corresponding SDEs in which the deterministic subproblem is discretized using structure-preserving schemes while the random subproblem is discretized using the Euler-Maruyama scheme. We obtain a sharp and uniform-in-time convergence analysis for the proposed numerical scheme that allows us to accurately compute long-time solutions of the SDEs. As such, we can compute the effective diffusivity for time-dependent flows. Finally, we present numerical results to demonstrate the accuracy and efficiency of the proposed method in computing effective diffusivity for the time-dependent Arnold-Beltrami-Childress (ABC) flow and Kolmogorov flow in three-dimensional space.
AMS subject classification: 35B27, 37A30, 60H35, 65M12, 65M75.

keywords:

Convection-enhanced diffusion; time-dependent chaotic flows; effective diffusivity; structure-preserving scheme; ergodic theory; convergence analysis.

1 Introduction

In this paper, we study the convection-enhanced diffusion phenomenon for particles moving in time-dependent chaotic flows, which is defined by the following passive tracer model, i.e., a stochastic differential equation (SDE),

\displaystyle\mathrm{d}\textbf{X}(t)=\textbf{v}(t,\textbf{X})\mathrm{d}t+\sigma\mathrm{d}\textbf{w}(t),\quad\textbf{X}\in\mathbb{R}^{d},

(1)

where $\textbf{X}(t)=(x_{1}(t),...,x_{d}(t))^{T}\in\mathbb{R}^{d}$ is the position of the particle, $\sigma>0$ is the molecular diffusivity, and $\{{\bf w}(t)\}_{t\geq 0}$ is the standard $d$ -dimensional Browinian motion. The velocity field $\textbf{v}(t,\textbf{x})$ is time-dependent and divergence free, i.e. $\nabla_{\textbf{x}}\cdot\textbf{v}(t,\textbf{x})=0$ , for all $t\geq 0$ . In order to guarantee the existence of the solution to the SDE (1), we also assume that $\textbf{v}(t,\textbf{x})$ is Lipschitz in x. The passive tracer model (1) has many applications in physical and engineering sciences, including atmosphere science, ocean science, chemical engineering, and combustion.

We will study the long-time large-scale behavior of the particle $\textbf{X}(t)$ in the passive tracer model (1), i.e., whether the motion of the particle ${\bf X}(t)$ has a long-time diffusive limit. Let $\textbf{X}^{\epsilon}(t)\equiv\epsilon\textbf{X}(t/\epsilon^{2})$ denote the rescaled process of (1). We aim to investigate whether $\textbf{X}^{\epsilon}(t)$ converge in law to a Brownian motion with a covariance matrix $D^{E}\in\mathbb{R}^{d\times d}$ as $\epsilon\rightarrow 0$ , where $D^{E}$ is called the effective diffusivity matrix. The dependence of $D^{E}$ on the velocity field of the passive tracer model is complicated and highly nontrivial. There are many theoretical works, where the homogenization theory was applied to study the effective diffusivity matrix $D^{E}$ of the passive tracer model with spatial periodic velocity fields or random velocity fields with short-range correlations; see e.g. [3, 12, 16, 27] and references therein.

For many complicated velocity fields of physical interests, one cannot apply the homogenization theory to compute the corresponding effective diffusivity matrix $D^{E}$ , or even determine its existence. Therefore, many numerical methods were developed to compute $D^{E}$ . These results include, among others, for time-independent Taylor-Green flows, the authors of [28] proposed a stochastic splitting method and calculated the effective diffusivity in the limit of vanishing molecular diffusion. For time-dependent chaotic flows, an efficient model reduction method based on the spectral method was developed to compute $D^{E}$ using the Eulerian framework [21]. The reader is referred to [22] for an extensive review of many existing mathematical theories and numerical simulations for the passive tracer model with different velocity fields.

Recently, we developed a robust structure-preserving Lagrangian scheme to compute the effective diffusivity for chaotic and stochastic flows in [32]. We also obtained a rigorous error estimate for the numerical scheme in [32]. Specifically, let $D^{E,num}$ denote the numerical effective diffusivity obtained by our method. We got the error estimate, $||D^{E,num}-D^{E}||\leq C\Delta t+C(T)(\Delta t)^{2}$ , where the computational time $T$ should be greater than the diffusion time (also known as mixing time). This error estimate is not sharp in the sense that the pre-factor $C(T)$ may grow fast with respect to $T$ , since the error estimation is based on a Gronwall inequality technique. Later, we obtained a sharp convergence rate for our numerical scheme and got rid of the term $C(T)$ in the error estimate in [31]. However, this result only holds for the passive tracer models with time-independent flows.

In this paper, we aim to obtain a sharp convergence analysis for our numerical scheme in computing the effective diffusivity of passive tracer models in spatial-temporal periodic velocity fields. These types of flow fields are well-known for exhibiting chaotic streamlines and have many applications in turbulent diffusion [22]. Since in this case the velocity field depends on the temporal variable, the generator associated with the stochastic process, i.e. the solution $\textbf{X}(t)$ in Eq.(1) becomes non-autonomous. The generator is now a parabolic-type operator (see Eq.(3)), instead of an elliptic-type operator that was studied in [31] when the flows are time-independent. The uniqueness only happens regarding each time period. Hence the extension of the analysis developed in [31] to time-dependent flows is not straightforward. We will develop new techniques to overcome this difficulty; see Theorem 4.3 and Lemma 4.5 in Section 4. We also emphasize that when the flows are time-independent, we can construct the ballistic orbits of the ABC flows and Kolmogorov flows and study their dynamic behaviors. When the flows are time-dependent, however, their streamlines are complicated, which makes it is difficult to construct the orbits of these flows.

Though there are several prior works on structure-preserving schemes for ODEs and SDEs, see e.g. [14, 15, 30, 1, 20] and references therein, our work has several novel contributions. The first novelty is the convergence analysis, where we develop new techniques to deal with time-dependent flows. To handle the parabolic-type generator, we pile up snapshots of each time step within a single time period together. By viewing the numerical solutions as a Markov process and exploring the ergodicity of the solution process, we succeed in obtaining a sharp convergence analysis for our method in computing the effective diffusivity, where the error estimate does not depend on the computational time. Therefore, we can compute the long-time solutions of passive tracer models without losing accuracy; see Fig.2 and Fig.4. If we choose the Gronwall inequality in the error estimate, we cannot get rid of the exponential growth pre-factor in the error term, which makes the convergence analysis not sharp. Most importantly, our convergence result reveals the equivalence of the definition of effective diffusivity using the Eulerian framework and the Lagrangian framework; see Theorem 4.8, which is fundamental and important. For 3D time-dependent flow problems, the Eulerian framework has good theoretical value but the Lagrangian framework is computationally accessible.

Another novelty is that the stochastic structure-preserving Lagrangian scheme is robust and quite cheap in computing the long-time solutions of the passive tracer model (1), especially for problems in three-dimensional space. If one adopts the Eulerian framework to compute the effective diffusivity of the passive tracer model (1), one needs to solve a convection-diffusion-type cell problem; see Eq.(5). When the molecular diffusivity $\sigma$ is small and/or the dimension of spatial variables is big, say $d=3$ , it is exorbitantly expensive to solve the cell problems. As indicated in Eq.(6), the effective diffusivities depend on the integration of the gradient of the solution to the cell problem. In many cases, e.g. time-dependent ABC flow, the effective diffusivities grows up rapidly as $\sigma$ decreases; see Fig.5. In our Lagrangian approach, we can overcome the difficulties of long-time integration of the SDEs (raised as $\sigma$ decreases) by using robust structure-preserving schemes. However, for the Eulerian approach, one needs to solve a four-dimensional PDE (three variables in space dimension and one variable in the time dimension) and solutions have sharp gradients as the diffusivity decreases, which makes the Eulerian approach for computing effective diffusivities expensive.

Numerical results show that our Lagrangian scheme is insensitive to the molecular diffusivity $\sigma$ and computational cost linearly depends on the dimension $d$ of spatial variables in the passive tracer models (1). Thus, we are able to investigate the convection-enhanced diffusion phenomenon for several typical time-dependent chaotic flows of physical interests, including the time-dependent ABC flow and the time-dependent Kolmogorov flow in three-dimensional space. We discover that the maximal enhancement is achieved in the former case, while a submaximal enhancement is observed in the latter case; see Fig.5 and Fig.3, respectively. In addition, we find that the level of chaos and the strength of diffusion enhancement seem to compete with each other in the time-dependent ABC flow; see Fig.6. To the best of our knowledge, our work appears to be the first one in the literature to develop numerical methods to study the convection-enhanced diffusion phenomenon in 3D time-dependent flows.

The rest of the paper is organized as follows. In Section 2, we give the definition of the effective diffusivity matrix using the Eulerian framework and the Lagrangian framework. In Section 3, we propose the stochastic structure-preserving Lagrangian scheme in computing effective diffusivity for the passive tracer model (1). In Section 4, we provide a sharp convergence analysis for the proposed method based on a probabilistic approach. In addition, we shall show that our method can be used to solve high-dimensional flow problems and the error estimate can be obtained naturally. In Section 5, we present numerical results to demonstrate the accuracy and efficiency of our method. We also investigate the convection-enhanced diffusion phenomenon for time-dependent chaotic flows. Concluding remarks are made in Section 6.

2 Effective diffusivity of the passive tracer models

There are two frameworks to define the effective diffusivity of the passive tracer models. We first discuss the Eulerian framework. One natural way to study the expectation of the paths for the SDE given by the Eq.(1) is to consider its associated backward Kolmogorov equation [26]. Due to the time-dependence nature of the velocity field, we need to deal with a space-time ergodic random flow. Specifically, given a sufficiently smooth function $\phi(\tau,\textbf{x})$ in $\mathbb{R}\times\mathbb{R}^{d}$ , let $u(t,\tau,\textbf{x})=\mathbb{E}\big{[}\phi(t+\tau,\textbf{X}_{t+\tau})|\textbf{X}_{\tau}=\textbf{x}\big{]}$ and $\textbf{X}(t)$ be the solution to Eq.(1). Then, $u(t,\tau,\textbf{x})$ satisfies the following backward Kolmogorov equation

\displaystyle u_{t}=\mathcal{L}u,\quad u(0,\tau,\textbf{x})=\phi(\tau,\textbf{x}).

(2)

In Eq.(2), the generator $\mathcal{L}$ is defined as

\displaystyle\mathcal{L}u=\partial_{\tau}{\color[rgb]{0,0,0}u}+\textbf{v}\cdot\nabla_{\textbf{x}}u+D_{0}\Delta_{\textbf{x}}u,

(3)

where $D_{0}=\sigma^{2}/2$ is the diffusion coefficient, v is the velocity field, and $\nabla_{\textbf{x}}$ and $\Delta_{\textbf{x}}$ denote the gradient operator and Laplace operator, respectively.

Remark 2.1.

Let $\rho(t,\tau,\textbf{x})$ denote the density function of the particle $(t+\tau,\textbf{X}(t+\tau))$ of Eq.(1). One can define the adjoint operator $\mathcal{L}^{*}$ as $\mathcal{L}^{*}\rho=-\partial_{\tau}{\color[rgb]{0,0,0}\rho}-\nabla\cdot(\textbf{v}\rho)+D_{0}\Delta\rho$ . Then, $\rho$ satisfies the Fokker-Planck equation $\rho_{t}=\mathcal{L}^{*}\rho$ with the initial density $\rho(t,\tau,\textbf{x})=\rho_{0}(\tau,\textbf{x})$ .

When v is incompressible (i.e. $\nabla_{\textbf{x}}\cdot\textbf{v}(t,\textbf{x})=0,\ \forall t$ ), deterministic and space-time periodic in $O(1)$ scale, where we assume the period of v is $1$ in each dimension of the physical and temporal space, the formula for the effective diffusivity matrix is [3, 27]

\displaystyle D^{E}=D_{0}I+\big{\langle}\textbf{v}(t,\textbf{x})\otimes\boldsymbol{\chi}(t,\textbf{x})\big{\rangle}_{p},

(4)

where we have assumed that the fluid velocity $\textbf{v}(t,\textbf{x})$ is smooth and the (vector) corrector field $\boldsymbol{\chi}$ satisfies the cell problem,

\displaystyle\mathcal{L}\boldsymbol{\chi}=-\textbf{v}(t,\textbf{y}),\quad(t,\textbf{y})\in\mathbb{T}\times\mathbb{T}^{d},

(5)

and $\langle\cdot\rangle_{p}$ denotes temporal and spatial average over $\mathbb{T}\times\mathbb{T}^{d}$ . Since v is incompressible, the solution $\boldsymbol{\chi}$ to the cell problem (5) is unique up to an additive constant by the Fredholm alternative. By multiplying $\boldsymbol{\chi}$ to Eq.(5), integrating the corresponding result in $\mathbb{T}\times\mathbb{T}^{d}$ and using the periodic conditions of $\boldsymbol{\chi}$ and v, we get an equivalent formula for the effective diffusivity as follows

\displaystyle D^{E}=D_{0}I+\big{\langle}\nabla\boldsymbol{\chi}(t,\textbf{x})\otimes\nabla\boldsymbol{\chi}(t,\textbf{x})\big{\rangle}_{p}.

(6)

The correction to $D_{0}I$ in Eq.(6) is nonnegative definite. We can see that $\textbf{e}^{T}D^{E}\textbf{e}\geq D_{0}$ for all unit column vectors $\textbf{e}\in\mathbb{R}^{d}$ , which is called convection-enhanced diffusion. By using a variational principle for time-periodic velocity flows, one can find a upper bound for the effective diffusivity, i.e., there exists a nonzero unit column vector $\textbf{e}\in\mathbb{R}^{d}$ , such that

\textbf{e}^{T}D^{E}\textbf{e}\sim\frac{1}{D_{0}},\quad\text{as }D_{0}\to 0,

(7)

which is known as the maximal enhancement. More details of the derivation can be found in [4, 24, 9]. We point out that many theoretical results were built upon the passive tracer models with time-independent flows. We are interested in studying the convection-enhanced diffusion phenomenon for time-dependent chaotic flows in this paper. Especially, whether the time-dependent chaotic flows still have the maximal enhancement.

In practice, the cell problem (5) can be solved using numerical methods, such as finite element methods and spectral methods. However, when $D_{0}$ becomes extremely small, the solutions of the cell problem (5) develop sharp gradients and demand a large number of finite element basis or Fourier basis to resolve, which makes the Eulerian framework expensive. In addition, when the dimension of spatial variables is big, say $d=3$ , the Eulerian framework becomes expensive too.

Alternatively, one can use the Lagrangian framework to compute the effective diffusivity matrix, which is defined as follows,

\displaystyle D_{ij}^{E}=\lim_{t\rightarrow\infty}\frac{\Big{\langle}\big{(}x_{i}(t)-x_{i}(0))(x_{j}(t)-x_{j}(0)\big{)}\Big{\rangle}}{2t},\quad 1\leq i,j\leq d,

(8)

where $\textbf{X}(t)=(x_{1}(t),...,x_{d}(t))^{T}$ is the position of a particle tracer at time $t$ and the average $\langle\cdot\rangle$ is taken over an ensemble of particles. If the above limit exists, that means the transport of the particle is a standard diffusion process, at least on a long-time scale. For example, when the velocity field is the Taylor-Green velocity field [9, 28], the long-time and large-scale behavior of the passive tracer model is a diffusion process. However, there are cases showing that the spreading of particles does not grow linearly with time but has a power-law $t^{\gamma}$ , where $\gamma>1$ and $\gamma<1$ correspond to super-diffusive and sub-diffusive behaviors, respectively; see e.g. [4, 22, 2].

We shall use the Lagrangian framework in this paper. The Lagrangian framework has the advantages that: (1) it is easy to implement; (2) it does not directly suffer from a small molecular diffusion coefficient $\sigma$ during the computation; and (3) its computational cost only linearly depends on the dimension of spatial variables in the passive tracer models. However, the major difficulty in solving Eq.(1) is that the computational time should be long enough to approach the diffusion time scale. To address this challenge, we shall develop robust numerical schemes, which are structure-preserving and accurate for long-time integration. Moreover, we aim to develop the convergence analysis of the proposed numerical schemes. Finally, we shall investigate the relationship between parameters of the time-dependent chaotic flows and the corresponding effective diffusivity.

3 Stochastic structure-preserving schemes

3.1 Derivation of numerical schemes

To demonstrate the main idea, we first construct stochastic structure-preserving schemes for a two-dimensional passive tracer model. The derivation of the numerical schemes for high-dimensional passive tracer models will be discussed in Section 4.5. Specifically, let $\textbf{X}=(x_{1},x_{2})^{T}$ denote the position of the particle, then the model can be written as

\begin{cases}\mathrm{d}x_{1}=v_{1}\mathrm{d}t+\sigma\mathrm{d}w_{1,t},\quad x_{1}(0)=x_{1,0},\\ \mathrm{d}x_{2}=v_{2}\mathrm{d}t+\sigma\mathrm{d}w_{2,t},\quad x_{2}(0)=x_{2,0},\end{cases}

(9)

where $w_{i,t}$ , $i=1,2$ are independent Brownian motions. We assume that $\textbf{v}=(v_{1},v_{2})^{T}$ is divergence-free and mean-zero at any time $t$ , i.e.,

\nabla\cdot\textbf{v}:=\partial_{x_{1}}v_{1}+\partial_{x_{2}}v_{2}=0\quad\forall t,

(10)

and

\begin{cases}\int_{\mathbb{T}}v_{1}(t,x_{1},x_{2})\mathrm{d}x_{2}=0\quad\forall x_{1},~t,\\ \int_{\mathbb{T}}v_{2}(t,x_{1},x_{2})\mathrm{d}x_{1}=0\quad\forall x_{2},~t,\end{cases}

(11)

where $\mathbb{T}=[0,1]$ . We also assume that v is smooth and its first-order derivatives $v_{i}(t,x_{1},x_{2})$ , $i=1,2$ are bounded. These conditions are necessary to guarantee the existence and uniqueness of solutions of the SDE (9); see [26]. Moreover, we assume that the diagonal of the Jacobian of the velocity field $\textbf{v}=(v_{1},v_{2})^{T}$ are all zeros. A typical example is a Hamiltonian system with a separable Hamiltonian, i.e., there exists $H(t,x_{1},x_{2})=H_{1}(t,x_{1})+H_{2}(t,x_{2})$ such that,

v_{1}=-\partial_{x_{2}}H,\quad v_{2}=\partial_{x_{1}}H.

(12)

In this paper, we denote with slightly abuse of notation that $v_{1}(t,x_{2})=v_{1}(t,x_{1},x_{2})$ and $v_{2}(t,x_{1})=v_{2}(t,x_{1},x_{2})$ . These notations simplify our derivation. Whenever a statement corresponds to $v_{1}(t,x_{2})$ (or $v_{2}(t,x_{1})$ ) is made, it is equivalent to that for $v_{1}(t,x_{1},x_{2})$ (or $v_{2}(t,x_{1},x_{2})$ ).

In [32], we proposed a stochastic structure-preserving scheme based on a Lie-Trotter splitting scheme to solve the SDE (9). Specifically, we split the problem (9) into a deterministic subproblem,

\begin{cases}\mathrm{d}x_{1}=v_{1}(t,x_{2})\mathrm{d}t,\\ \mathrm{d}x_{2}=v_{2}(t,x_{1})\mathrm{d}t,\end{cases}

(13)

which is solved by using a symplectic-preserving scheme (e.g., the symplectic Euler scheme for deterministic equations with frozen time), and a stochastic subproblem,

\begin{cases}\mathrm{d}x_{1}=\sigma\mathrm{d}w_{1,t},\\ \mathrm{d}x_{2}=\sigma\mathrm{d}w_{2,t},\end{cases}

(14)

which is solved by using the Euler-Maruyama scheme [26]. Notice that when $\sigma$ is a constant in (14), the Euler-Maruyama scheme exactly solves Eq.(14)

Now we discuss how to discretize Eq.(9). From time $t=t_{n}$ to time $t=t_{n+1}$ , where $t_{n+1}=t_{n}+\Delta t$ , $t_{0}=0$ , and $\Delta t$ is the time step, we assume the numerical solution $\textbf{X}^{n}=(x_{1}^{n},x_{2}^{n})^{T}$ is given, which approximates the exact solution $\textbf{X}(t_{n})$ to the SDE (9) at time $t_{n}=n\Delta t$ . Then, we apply the Lie-Trotter splitting method to solve the SDE (9) and obtain,

\begin{cases}x_{1}^{n+1}=x_{1}^{n}+v_{1}(t_{n+\frac{1}{2}},x_{2}^{n})\Delta t+\sigma N^{n}_{1},\\ x_{2}^{n+1}=x_{2}^{n}+v_{2}\big{(}t_{n+\frac{1}{2}},x_{1}^{n}+v_{1}(t_{n+\frac{1}{2}},x_{2}^{n})\Delta t\big{)}\Delta t+\sigma N^{n}_{2},\end{cases}

(15)

where $t_{n+\frac{1}{2}}=t_{n}+\frac{\Delta t}{2}$ , $N^{n}_{1}=\sqrt{\Delta t}\xi_{1}$ , $N^{n}_{2}=\sqrt{\Delta t}\xi_{2}$ , and $\xi_{1}$ , $\xi_{2}\sim\mathcal{N}(0,1)$ are i.i.d. normal random variables. In this paper, we view the solution sequence $\textbf{X}^{n}=(x_{1}^{n},x_{2}^{n})^{T}$ , $n=1,2,3,...$ , generated by the scheme (15) as a discrete Markov stochastic process, which enables us to use techniques from stochastic process to obtain a sharp convergence analysis for the numerical solutions; see Section 4.

In a 2D Hamiltonian system, when the system contains an additive temporal noise, the noise itself is considered to be symplectic pathwise [25]. Therefore, we state that the scheme (15) is stochastic symplectic-preserving since it preserves symplecticity. Specifically, the scheme (15) can be viewed as a composition of two symplectic transforms. In addition, we know that the numerical solution converges to the exact one as the time step $\Delta t$ approaches zero. In high-dimensional systems, a structure-preserving scheme refers to a volume-preserving scheme; see Section 4.5.

3.2 The backward Kolmogorov equation and related results

We first define the backward Kolmogorov equation associated with the Eq.(9) as

\displaystyle u_{t}=\mathcal{L}u,\quad u(0,\tau,\textbf{x})=\phi(\tau,\textbf{x}),

(16)

where the generator $\mathcal{L}$ associated with the Markov process in Eq.(9) is given by

\displaystyle\mathcal{L}=\partial_{\tau}+v_{1}(\tau,x_{2})\partial_{x_{1}}+v_{2}(\tau,x_{1})\partial_{x_{2}}+\frac{\sigma^{2}}{2}(\partial_{x_{1}x_{1}}+\partial_{x_{2}x_{2}}).

(17)

Recall that the solution $u(t,\tau,\textbf{x})$ to the Eq.(16) satisfies, $u(t,\tau,\textbf{x})=\mathbb{E}\big{[}\phi(t+\tau,\textbf{X}_{t+\tau})|\textbf{X}_{\tau}=\textbf{x}\big{]}$ where $\textbf{X}_{t}$ is the solution to Eq.(9) and $\phi$ is a smooth function in $\mathbb{R}^{1}\times\mathbb{R}^{2}$ . In other words, $u(t,\tau,\textbf{x})$ is the flow generated by the original SDE (9).

Similarly, we can study the flow generated by the stochastic structure-preserving scheme (15). According to the splitting method used in the derivation of the scheme in Section 3.1, we respectively define $\mathcal{L}_{1}=\partial_{\tau}$ , $\mathcal{L}_{2}=v_{1}\partial_{x_{1}}$ , $\mathcal{L}_{3}=v_{2}\partial_{x_{2}}$ , and $\mathcal{L}_{4}=\frac{\sigma^{2}}{2}(\partial_{x_{1}x_{1}}+\partial_{x_{2}x_{2}})$ . Starting from $u(0,\cdot,\cdot)$ , during one time step $\Delta t$ , we compute

\begin{cases}\partial_{t}u^{1}&=\mathcal{L}_{1}u^{1},\quad u^{1}(0,\cdot,\cdot)=u(0,\cdot,\cdot),\\ \partial_{t}u^{2}&=\mathcal{L}_{2}u^{2},\quad u^{2}(0,\cdot,\cdot)=u^{1}(\frac{\Delta t}{2},\cdot,\cdot),\\ \partial_{t}u^{3}&=\mathcal{L}_{3}u^{3},\quad u^{3}(0,\cdot,\cdot)=u^{2}(\Delta t,\cdot,\cdot),\\ \partial_{t}u^{4}&=\mathcal{L}_{1}u^{4},\quad u^{4}(0,\cdot,\cdot)=u^{3}(\Delta t,\cdot,\cdot),\\ \partial_{t}u^{5}&=\mathcal{L}_{4}u^{5},\quad u^{5}(0,\cdot,\cdot)=u^{4}(\frac{\Delta t}{2},\cdot,\cdot).\\ \end{cases}

(18)

Then, $u^{5}(\Delta t,\cdot,\cdot)$ will be the flow at time $t=\Delta t$ generated by our stochastic structure-preserving scheme (15) and it approximates the solution $u(\Delta t,\cdot,\cdot)$ to the Eq.(16) well when $\Delta t$ is small. It is also worth mentioning that, $u^{3}(\Delta t,\cdot,\cdot)$ is the exact flow generated by the deterministic symplectic Euler scheme in solving Eq.(13). We repeat this process to compute the flow equations of our scheme at other time steps, which approximate the solution $u(n\Delta t,\cdot,\cdot),n=2,3,...$ to the Eq.(16) at different time steps.

Remark 3.1.

Given the operators $\mathcal{L}_{i}$ , $i=1,2,3,4$ , there are many possible choices in setting the coefficients for each operator $\mathcal{L}_{i}$ and designing the splitting method; see Section 2.5 of [14]. Eq.(18) is a simple choice that was used in this paper.

To analyze the error between the flow operator in Eq.(16) and the composition of operators in Eq.(18), we shall resort to the Baker-Campbell-Hausdorff (BCH) formula, which is widely used in non-commutative algebra [13]. For example, in the matrix theory,

\exp(tA)\exp(tB)=\exp\bigg{(}t(A+B)+t^{2}\frac{[A,B]}{2}+\frac{t^{3}}{12}\Big{(}\big{[}A,[A,B]\big{]}+\big{[}B,[B,A]\big{]}\Big{)}+\cdots\bigg{)},

(19)

where $t$ is a scalar, $A$ and $B$ are two square matrices of the same size, $[,]$ is the Lie-Bracket, and the remaining terms on the right hand side are all nested Lie-brackets.

In our analysis, we replace the matrices in Eq.(19) by differential operators and the BCH formula yields critical insights into the particular structure of the splitting error. Let $I_{\Delta t}$ denote the composite flow operator associated with Eq.(18), i.e.,

I_{\Delta t}u(0,\cdot,\cdot):=\exp(\Delta t\mathcal{L}_{4})\exp(\frac{\Delta t}{2}\mathcal{L}_{1})\exp(\Delta t\mathcal{L}_{3})\exp(\Delta t\mathcal{L}_{2})\exp(\frac{\Delta t}{2}\mathcal{L}_{1})u(0,\cdot,\cdot).

(20)

Notice that after propagating time $t=\Delta t$ , the exact solution to the Eq.(16) started at any $\tau$ can be represented as

u(\Delta t,\cdot,\cdot)=\exp(\Delta t\mathcal{L})u(0,\cdot,\cdot)=\exp\big{(}\Delta t(\mathcal{L}_{1}+\mathcal{L}_{2}+\mathcal{L}_{3}+\mathcal{L}_{4})\big{)}u(0,\cdot,\cdot).

(21)

Therefore, we can apply the BCH formula to analyze the error between the original flow and the approximated flow. Moreover, we find that computing the $k$ -th order modified equation associated with Eq.(9) in the backward error analysis (BEA) [29, 7] is equivalent to computing the terms of BCH formula up to order $(\Delta t)^{k}$ in the Eq.(20). To show that the solution generated by Eq.(15) follows a perturbed Hamiltonian system (with divergence-free velocity and additive noise) at any order $k$ , we only need to consider the $(k+1)$ -nested Lie bracket consisting of $\big{\{}\partial_{\tau},v_{1}\partial_{x_{1}},\ v_{2}\partial_{x_{2}},\ \frac{\sigma^{2}}{2}(\partial_{x_{1}x_{1}}+\partial_{x_{2}x_{2}})\big{\}}$ and we can easily see that they generate divergence-free fields.

We remark that given any explicit splitting scheme for deterministic systems, by adding additive noise we shall obtain a similar form of flow propagation. And we shall see in later proof that, such operator formulation is very effective in analyzing the order of convergence and volume-preserving property.

4 Convergence analysis

In this section, we prove the convergence rate of our stochastic structure-preserving schemes in computing effective diffusivity based on a probabilistic approach, which allows us to get rid of the exponential growth factor in the error estimate. We first limit our analysis to 2D separable Hamiltonian velocity fields. Then, in Section 4.5 we will show that all the derivations can be generalized to high-dimensional cases.

4.1 Convergence to an invariant measure

To compute the effective diffusivity of a passive tracer model using a Lagrangian numerical scheme is closely related to study the limit of a solution sequence (a stochastic process) generated by the numerical scheme. Therefore, we can apply the results from ergodic theory to study the convergence behaviors of the solution.

Let $(S,\Sigma)$ be a probability space, on which a family $P(\textbf{x},E)$ , $\textbf{x}\in S$ , $E\in\Sigma$ , of probability measure is defined. We assume $\textbf{x}\to P(\textbf{x},E)$ is measurable, $\forall E\in\Sigma$ . This corresponds to a linear bounded operator on $\mathcal{B}(S)$ , which is the space of bounded measurable functions on $S$ . This operator, denoted by $P$ , is defined by,

P\phi(\textbf{x})=\int_{S}P(\textbf{x},\mathrm{d}\textbf{z})\phi(\textbf{z}),\quad\forall\phi\in\mathcal{B}(S).

(22)

Clearly $||P||\leq 1$ . One of the main objectives of ergodic theory is to study the limit of the operator sequence $P^{n}$ as $n\rightarrow+\infty$ . The result can be summarized into the following proposition, which plays a fundamental role in our convergence analysis.

Proposition 4.1 (Theorem 3.3.1 of [3]).

We assume that,

1.

$S$ is a compact metric space and $\Sigma$ is the Borel $\sigma$ -algebra;
2.

there exists a probability measure $\mu$ on $(S,\Sigma)$ such that $P(\textbf{x},E)=\int_{E}p(\textbf{x},\textbf{y})\mu(\mathrm{d}\textbf{y})$ ;
3.

$p(\textbf{x},\textbf{y}):S\times S\to\mathbb{R}^{+}$ is continuous;
4.

there exists a ball $U_{0}$ such that $\mu(U_{0})>0$ and a positive number $\delta>0$ (depending on $U_{0}$ ) such that $p(\textbf{x},\textbf{y})\geq\delta$ , $\textbf{x}\in S$ , $\forall\textbf{y}\in U_{0}$ .

Then, there exists one and only one invariant probability measure $\pi$ on $(S,\Sigma)$ and one has,

\sup_{\textbf{x}\in S}\Big{|}P^{n}\phi(\textbf{x})-\int\phi\pi(\mathrm{d}\textbf{x})\Big{|}\leq C||\phi||e^{-\rho n},\ \forall\phi\in\mathcal{B}(S),

(23)

where $\rho=\log\frac{1}{1-\delta\mu(U_{0})}>0$ and $C=\frac{2}{1-\delta\mu(U_{0})}>0$ are independent of $\phi$ .

Now we study the convergence behaviors of the solution generated by our stochastic structure-preserving scheme (15). We first prove a lemma as follows.

Lemma 4.2.

Let $\tilde{Y}=\mathbb{R}^{2}/\mathbb{Z}^{2}$ denote the physical torus space and $\mathbb{T}$ be the time periodic space. Let $I_{\tau,1+\tau}^{*}$ denote the transform of the density on $\tilde{Y}$ during $[\tau,1+\tau]$ (time period is $1$ ) using the numerical scheme (15). In addition, let $I_{\tau,1+\tau}$ denote the adjoint operator (i.e., the flow operator) of $I_{\tau,1+\tau}^{*}$ in the space of $\mathcal{B}(\tilde{Y})$ , which is the set of bounded measurable functions on $\tilde{Y}$ . Then, there exists one and only one invariant probability measure on $(\tilde{Y},\Sigma)$ , denoted by $\pi_{\tau}$ , satisfying,

\sup_{\textbf{x}\in\tilde{Y}}\Big{|}\big{(}(I_{\tau,1+\tau})^{n}\phi\big{)}(\textbf{x})-\int\phi(\textbf{x}^{\prime})\pi_{\tau}(\mathrm{d}\textbf{x}^{\prime})\Big{|}\leq C||\phi||_{L_{\infty}}e^{-\rho n},\quad\forall\phi\in\mathcal{B}(\tilde{Y}),

(24)

where $\rho>0$ , $C>0$ are independent of $\phi(\cdot)$ . Moreover, the kernel space of $(I_{d}-I_{\tau,1+\tau})$ is the constant functions in $\tilde{Y}$ , where $I_{d}$ is the identity operator.

Proof.

We shall verify that the transition kernel associated with the numerical scheme (15) satisfies the assumptions required by Prop. 4.1. First notice that in the space $\mathbb{R}^{2}$ , the integration process associated with the numerical scheme (15) can be expressed as a Markov process with the transition kernel,

		$\displaystyle K_{t}\big{(}\textbf{X}^{n},\textbf{X}^{n+1}\big{)}=\frac{1}{2\pi\sigma^{2}\Delta t}\cdot$
		$\displaystyle\leavevmode\resizebox{422.77661pt}{}{$\exp\Bigg{(}-\frac{\Big{(}x_{1}^{n+1}-x_{1}^{n}-v_{1}(t+\frac{\Delta t}{2},x_{2}^{n})\Delta t\Big{)}^{2}+\Big{(}x_{2}^{n+1}-x_{2}^{n}-v_{2}\big{(}t+\frac{\Delta t}{2},x^{n+1}_{1}-x^{n}_{1}-v_{1}(t+\frac{\Delta t}{2},x^{n}_{2})\Delta t\big{)}\Delta t\Big{)}^{2}}{2\sigma^{2}\Delta t}\Bigg{)}$},$		(25)

where $\textbf{X}^{n}=(x^{n}_{1},x^{n}_{2})^{T}$ and $\textbf{X}^{n+1}=(x^{n+1}_{1},x^{n+1}_{2})^{T}$ are the numerical solutions at time $t=t_{n}$ and $t=t_{n+1}$ , respectively.

Then, using the periodicity of v, we directly extend Eq.(4.2) to the torus space $\tilde{Y}$ as

		$\displaystyle\tilde{K}_{\tau}\big{(}\textbf{X}^{n},\textbf{X}^{n+1}\big{)}=\sum_{i,j\in\mathbb{Z}}\frac{1}{2\pi\sigma^{2}\Delta t}\cdot$
		$\displaystyle\leavevmode\resizebox{422.77661pt}{}{$\exp\Bigg{(}-\frac{\Big{(}x^{n+1}_{1}+i-x^{n}_{1}-v_{1}(\tau+\frac{\Delta t}{2},x^{n}_{2})\Delta t\Big{)}^{2}+\Big{(}x^{n+1}_{2}+j-x^{n}_{2}-v_{2}\big{(}\tau+\frac{\Delta t}{2},x^{n+1}_{1}-x^{n}_{1}-v_{1}(\tau+\frac{\Delta t}{2},x^{n}_{2})\Delta t\big{)}\Delta t\Big{)}^{2}}{2\sigma^{2}\Delta t}\Bigg{)}$}.$		(26)

Let $\tilde{\textbf{K}}_{\tau,\tau+k\Delta t}$ denote the kernel from $\tau$ to $\tau+k\Delta t$ , which is the density of the transition kernel associated with applying our scheme starting from time $\tau$ for $k$ steps. Then, we have

\tilde{\textbf{K}}_{\tau,\tau+k\Delta t}(\textbf{X}^{0},\textbf{X}^{k})=\int_{(\tilde{Y})^{k-1}}\prod_{m=0}^{k-1}\tilde{K}_{\tau+m\Delta t}(\textbf{X}^{m},\textbf{X}^{m+1})\mathrm{d}\textbf{X}^{1}\mathrm{d}\textbf{X}^{2}\cdots\mathrm{d}\textbf{X}^{k-1}.

(27)

We choose $k=\frac{1}{\Delta t}$ and obtain $\tilde{\textbf{K}}_{\tau,\tau+1}$ . One can see that the kernel $\tilde{\textbf{K}}_{\tau,\tau+1}$ is essentially bounded above zero since $\tilde{K}_{\tau+m\Delta t}$ in (27) are all positive. Moreover, if $0<\Delta t\ll 1$ , $\tilde{\textbf{K}}_{\tau,\tau+1}$ is a continuous function on the domain $\tilde{Y}\times\tilde{Y}$ . Then by noticing that the domain $\tilde{Y}\times\tilde{Y}$ is compact, the kernel $\tilde{\textbf{K}}_{\tau,\tau+1}$ is strictly positive. Namely, there exists $\delta_{\tau}>0$ such that $\tilde{\textbf{K}}_{\tau,\tau+1}(\textbf{X}^{0},\textbf{X}^{k})>\delta_{\tau},\ \forall(\textbf{X}^{0},\textbf{X}^{k})\in\tilde{Y}\times\tilde{Y}$ . If we apply Prop.4.1 to $I_{\tau,1+\tau}$ (whose kernel is $\tilde{\textbf{K}}_{\tau,\tau+1}$ ), we prove the statement in (24).

Finally, we know that the operator $I_{\tau,1+\tau}$ is compact since it is an integral operator with a continuous kernel. By using the Fredholm alternative, we know that $\dim\ker(I_{d}-I_{\tau,1+\tau})=\dim\ker(I_{d}-I_{\tau,1+\tau}^{*})=1$ . Therefore, it is easy to verify that the constant functions are in the kernel of $(I_{d}-I_{\tau,1+\tau})$ . ∎

Equipped with the Lemma 4.2, we study the convergence rate of the space-time transition kernel associated with our numerical scheme (15).

Theorem 4.3.

Let $\Delta t=\frac{1}{N}$ , $N$ is a positive integer. We have the following properties hold:

(a)

Given $\Delta t$ , there exists $C>0$ and $\rho>0$ , such that,

\sup_{\tau,\textbf{x}}\Big{|}\big{(}I_{\Delta t}^{N}\big{)}^{n}\phi(\tau,\textbf{x})-\int\phi(\tau,\textbf{x}^{\prime})\pi_{\tau}(\mathrm{d}\textbf{x}^{\prime})\Big{|}\leq C||\phi||_{L_{\infty}}e^{-\rho n},\quad\forall\phi\in\mathcal{B}(\mathbb{T}\times\tilde{Y}),

(28)

where $C$ and $\rho$ do not depend on $\phi$ and $\tau$ .

(b)

If $\int_{\tilde{Y}}\phi\pi_{\tau}=0$ , then we get

$\lim_{n\to\infty}\sum_{i=1}^{n}\mathbb{E}\phi(\tau,\textbf{X}^{N\tau+i})<\infty,\quad\forall\tau\in\mathbb{T}.$ (29)
(c)

The kernel space of $(I_{d}-I_{\Delta t}^{N})$ is $\big{\{}c(\tau)~|~c(\tau)\text{ is a periodic function in}~\mathbb{T}~\text{with period 1}\big{\}}$ .

Proof.

By definition of $I_{\Delta t}$ and $I_{\tau,1+\tau}$ in Eq.(20) and Lemma 4.2, we have $(I_{\Delta t})^{N}\phi(\tau,\cdot)\equiv I_{\tau,1+\tau}\phi(\tau,\cdot)$ . To prove the property (a), we need to show that the lower bound of the kernel $\tilde{\textbf{K}}_{\tau,\tau+1}$ , which is defined in the proof of Lemma 4.2, does not depend on $\tau$ . For all $\tau\in\mathbb{T}$ , $\textbf{X}^{n}=(x^{n}_{1},x^{n}_{2})^{T}\in\mathbb{T}^{2}$ and $\textbf{X}^{n+1}=(x^{n+1}_{1},x^{n+1}_{2})^{T}\in\mathbb{T}^{2}$ , we pick $i_{0}=\lfloor-x^{n+1}_{1}+x^{n}_{1}+v_{1}(\tau+\frac{\Delta t}{2},x^{n}_{2})\Delta t\rfloor$ and $j_{0}=\lfloor-x^{n+1}_{2}+x^{n}_{2}+v_{2}\big{(}\tau+\frac{\Delta t}{2},x^{n+1}_{1}-x^{n}_{1}-v_{1}(\tau+\frac{\Delta t}{2},x^{n}_{2})\Delta t\big{)}\Delta t\rfloor$ , where $\lfloor a\rfloor$ denotes the largest integer not greater than $a$ . Applying to Eq.(4.2), we can see that

	$\displaystyle\tilde{K}_{\tau}\big{(}\textbf{X}^{n},\textbf{X}^{n+1}\big{)}\geq\frac{1}{2\pi\sigma^{2}\Delta t}\cdot$

	$\displaystyle\geq\frac{1}{2\pi\sigma^{2}\Delta t}\exp\big{(}-\frac{1}{\sigma^{2}\Delta t}\big{)}>0.$		(30)

According to the definition of the kernel $\tilde{\textbf{K}}_{\tau,\tau+1}$ ; see Eq.(27), we know the minimal value of $\tilde{\textbf{K}}_{\tau,\tau+1}$ is above zero and is independent of $\tau$ . Now, we apply this observation to Lemma 4.2 and conclude the proof of the property (a). The property (b) is a simple conclusion of the exponential decay property proved in (a). For the property (c), we consider the equation $I_{\Delta t}^{N}w=w$ . Then, for a given time $\tau$ , we have $I_{\tau,1+\tau}w(\tau,\cdot)=w(\tau,\cdot)$ . Notice the fact that in Lemma 4.2 the invariant space of $I_{\tau,1+\tau}$ is constant in the spacial variable. Thus, we obtain $w=w(\tau)$ . ∎

Before we close this subsection, we provide a convergence result for the inverse of operator sequences, which will be useful in our convergence analysis.

Proposition 4.4.

Let $X,Y$ denote two Banach spaces. Assume $T_{n}$ , $T$ are bounded linear operators from $X$ to $Y$ , satisfying $\lim_{n\to\infty}||T_{n}-T||_{\mathcal{B}(X,Y)}=0$ , and $T^{-1}\in\mathcal{B}(Y,X)$ . Given $f\in Y$ , if $T_{n}^{-1}f$ , $n=1,2,...$ uniquely exist, then we have a convergence estimate as follows,

\displaystyle\lim_{n\to\infty}\big{|}\big{|}(T_{n}^{-1}-T^{-1})f\big{|}\big{|}=0.

(31)

The proof is quite standard. It can also be viewed as a modification of Theorem 1.16 in Section IV of [17].

4.2 A discrete-type cell problem

In the Eulerian framework, the periodic solution of the cell problem (5) and the corresponding formula for the effective diffusivity (4) play a key role in studying the behaviors of chaotic and stochastic flows. In the Lagrangian framework, we shall define a discrete analogue of cell problem that enables us to compute the effective diffusivity. Let $\textbf{X}^{0}=(x_{1}^{0},x_{2}^{0})^{T}$ be the initial data and $\textbf{X}^{n}=(x_{1}^{n},x_{2}^{n})^{T}$ denote the numerical solution generated by the scheme (15) at $t_{n}=n\Delta t$ , i.e.

\begin{cases}x_{1}^{n}=x_{1}^{n-1}+v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})\Delta t+\sigma N^{n-1}_{1},\\ x_{2}^{n}=x_{2}^{n-1}+v_{2}\big{(}t_{n-\frac{1}{2}},x_{1}^{n-1}+v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})\Delta t\big{)}\Delta t+\sigma N^{n-1}_{2},\end{cases}

(32)

where $N^{n-1}_{1}=\sqrt{\Delta t}\xi_{1}$ , $N^{n-1}_{2}=\sqrt{\Delta t}\xi_{2}$ , and $\xi_{1}$ , $\xi_{2}\sim\mathcal{N}(0,1)$ are i.i.d. normal random variables. For convenience we have replaced $n+1$ by $n$ .

First of all, we show that the solutions $x_{1}^{n}$ and $x^{n}_{2}$ obtained by the scheme (32) have bounded expectations if the initial values are bounded. Taking expectation of the first equation of Eq.(32) on both sides, we obtain

\displaystyle\mathbb{E}x_{1}^{n}=\mathbb{E}x_{1}^{n-1}+\Delta t\mathbb{E}v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})=\mathbb{E}x_{1}^{0}+\Delta t\sum_{k=0}^{n-1}\mathbb{E}v_{1}(t_{k+\frac{1}{2}},x_{2}^{k}).

(33)

As a symplectic scheme in 2D, the numerical scheme (32) admits the uniform measure as its invariant measure. Applying the results (a) and (b) of Theorem 4.3 and using the fact that v is a periodic function with zero mean, we know that,

\sup_{\textbf{X}^{0}\in\tilde{Y}}|\mathbb{E}v_{1}(t_{k+\frac{1}{2}},\textbf{X}^{k})|\leq e^{-\rho k}C_{N}\sup_{m=1,2,...,N,\ \textbf{x}\in\mathbb{T}^{2}}||v_{1}(t_{m+\frac{1}{2}},\textbf{x})||_{\infty}.

(34)

Notice that $v_{1}(t_{k+\frac{1}{2}},\textbf{X}^{k})$ is equivalent to $v_{1}(t_{k+\frac{1}{2}},x_{2}^{k})$ , since $v_{1}$ is indpentdent of $x_{1}^{k}$ . By applying triangle inequalities in Eq.(33) and using the result in Eq.(34), we arrive at,

|\mathbb{E}x_{1}^{n}|\leq|\mathbb{E}x_{1}^{0}|+C_{1}||v_{1}||_{\infty},

(35)

where $C_{1}$ does not depend on $n$ . Using the same approach, we know that expectation of the second component $\mathbb{E}x^{n}_{2}$ is also bounded.

Now, we are in the position to define the discrete-type cell problem. Starting at time $\tau$ with time step $\Delta t=\frac{1}{N}$ , we denote the starting time index to be $N\tau$ . Then, we define

\hat{v}_{1,N}(\tau,\textbf{x})=\Delta t\sum_{i=0}^{\infty}\mathbb{E}\big{[}v_{1}(t_{i+\frac{1}{2}}+\tau,\textbf{X}^{N\tau+i})|\textbf{X}^{N\tau}=\textbf{x}\big{]},

(36)

where the summation is well defined due to the fact stated in Eq.(34). We will show that $\hat{v}_{1,N}(\tau,\textbf{x})$ satisfies the following properties. Namely, $\hat{v}_{1,N}(\tau,\textbf{x})$ is the solution of the discrete-type cell problem defined in Eq.(37).

Lemma 4.5.

According to our assumption on v, we know that $v_{1}$ is a periodic function with zero mean on $\tilde{Y}$ , $\forall\tau$ , i.e., $\int_{\tilde{Y}}v_{1}=0$ . Therefore, $\hat{v}_{1,N}(\tau,\textbf{x})$ is the unique solution in $\mathcal{B}_{0}(\mathbb{T}\times\tilde{Y})$ such that

\hat{v}_{1,N}(\tau,\textbf{x})=(I_{\Delta t}\hat{v}_{1,N})(\tau,\textbf{x})+\Delta tv_{1}(\tau+\frac{\Delta t}{2},\textbf{x}),

(37)

where $\Delta t=\frac{1}{N}$ and the operator $I_{\Delta t}$ is defined in (20). Moreover, $\hat{v}_{1,N}(\tau,\textbf{x})$ is smooth.

Proof.

Throughout the proof, we shall use the fact that if $X$ , $Y$ are random processes and $Y$ is measurable under a filtration $\mathcal{F}$ , then with appropriate integrability assumption, we have

\mathbb{E}[XY]=\mathbb{E}\Big{[}\mathbb{E}[XY|\mathcal{F}]\Big{]}=\mathbb{E}\Big{[}\mathbb{E}[X|\mathcal{F}]Y\Big{]}.

(38)

Some simple calculations will give that

$\displaystyle\hat{v}_{1,N}(\tau,\textbf{x})-\Delta tv_{1}(\tau+\frac{\Delta t}{2},\textbf{x})=$	$\displaystyle\Delta t\sum_{i=1}^{\infty}\mathbb{E}\big{[}v_{1}(t_{i+\frac{1}{2}}+\tau,\textbf{X}^{N\tau+i})\|\textbf{X}^{N\tau}=\textbf{x}\big{]}$
$\displaystyle=$	$\displaystyle\mathbb{E}\Big{[}\Delta t\sum_{i=1}^{\infty}\mathbb{E}\big{[}v_{1}(t_{i+\frac{1}{2}}+\tau,\textbf{X}^{N\tau+i})\|\textbf{X}^{N\tau+1}\big{]}\|\textbf{X}^{N\tau}=\textbf{x}\Big{]}$
$\displaystyle=$	$\displaystyle\mathbb{E}\big{[}\hat{v}_{1,N}(\tau+\Delta t,\textbf{X}^{N\tau+1})\|\textbf{X}^{N\tau}=\textbf{x}\big{]}.$	(39)

Recall the definition of the operator $I_{\Delta t}$ in (20), Eq.(39) implies that

\hat{v}_{1,N}(\tau,\textbf{x})-\Delta tv_{1}(\tau+\frac{\Delta t}{2},\textbf{x})=(I_{\Delta t}\hat{v}_{1,N})(\tau,\textbf{x}).

(40)

Suppose we have that $I_{\Delta t}w=w$ . Then, we get $(I_{\Delta t})^{N}w=w$ . According to Theorem 4.3, we know that $w=0$ if $\int_{\tilde{Y}}w\mathrm{d}\textbf{x}=0$ , $\forall t$ . So $\ker(I_{\Delta t}-I_{d})=\{0\}$ and $\hat{v}_{1,N}$ is unique. Finally, by the definition of $\hat{v}_{1,N}$ , we obtain that

	$\displaystyle\hat{v}_{1,N}(\tau,\textbf{x})=$	$\displaystyle\Delta t\sum_{i=0}^{\infty}\mathbb{E}\big{[}v_{1}(t_{i+\frac{1}{2}}+\tau,\textbf{X}^{N\tau+i})\|\textbf{X}^{N\tau}=\textbf{x}\big{]}$
	$\displaystyle=$	$\displaystyle\Delta t\sum_{i=0}^{\infty}\int_{\tilde{Y}}v_{1}(t_{i+\frac{1}{2}}+\tau,\textbf{y})\tilde{K}_{\tau,\tau+i\Delta t}(\textbf{x},\textbf{y})\mathrm{d}\textbf{y},$		(41)

which indicates that $\hat{v}_{1,N}$ has the same regularity as $v_{1}$ does. Notice the kernel $\tilde{K}_{\tau,\tau+i\Delta t}(\textbf{x},\textbf{y})$ has a fast decay property, which guarantees the order of the differentiation and summation is interchangeable. ∎

Remark 4.1.

Notice that $v_{1}$ and $\hat{v}_{1,N}$ only depend on the second component of the numerical solution $\textbf{X}^{n}=(x_{1}^{n},x_{2}^{n})^{T}$ . However, we will write $v_{1}$ and $\hat{v}_{1,N}$ as functions of $\textbf{X}^{n}$ when we view $\textbf{X}^{n}$ as a Markov process in the convergence analysis.

Remark 4.2.

When the flow is time-independent, we obtain

\mathbb{E}[\hat{v}_{1,N}(\textbf{X}^{n+1})|\textbf{X}^{n}]-\hat{v}_{1,N}(\textbf{X}^{n})=-\Delta tv_{1}(\textbf{X}^{n}),\quad a.s.\quad\forall n\in\mathbb{N}.

(42)

Therefore, the discrete-type cell problem defined in (37) is a generalization of the discrete-type cell problem for time-independent flow problems studied in our previous work [31], although technically it is more involved.

In the remaining part of this paper, we only need the result that $\hat{v}_{1,N}(\tau,\textbf{x})$ is unique in an Hölder space $\mathbb{C}^{p_{1},p_{2},\alpha}_{0}(\mathbb{T}\times\tilde{Y})\subsetneq\mathcal{B}(\mathbb{T}\times\tilde{Y})$ . To be precise, given a smooth drift function $v_{1}$ , $\hat{v}_{1,N}(\tau,\textbf{x})$ will be in $\mathbb{C}^{p_{1},p_{2},\alpha}_{0}(\tilde{Y})$ , where $p_{1}\geq 2,p_{2}\geq 6,0<\alpha<1$ and the subscript index $0$ indicates that it is a subspace with zero-mean functions.

4.3 Convergence estimate of the discrete-type cell problem

In this section, we shall prove that the solution $\hat{v}_{1,N}(\tau,\textbf{x})$ of the discrete-type cell problem (i.e., Eq.(37)) converges to the solution of a continuous cell problem in certain subspace. Here, we choose the space $\mathbb{C}^{2,6,\alpha}_{0}(\mathbb{T}^{1}\times\tilde{Y})$ to carry out our analysis. However, there is no requirement that we have to choose this one. In fact, any space that has certain regularity (belongs to the domain of the operator $\mathcal{L}$ ) will work. Notice that the continuous cell problem (5) is defined for a vector function, where the first component satisfies

\mathcal{L}\chi_{1}=-v_{1}.

(43)

For the two-dimensional problem, the operator $\mathcal{L}$ is defined in Eq.(17). Given the fact that $v_{1}$ is a smooth function defined on $\mathbb{T}^{1}\times\tilde{Y}$ , which satisfies $\int_{\tilde{Y}}v_{1}(\tau,\textbf{x})\mathrm{d}\textbf{x}=0,\ \forall\tau\in\mathbb{T}^{1}$ . Then, Eq.(43) admits a unique solution $\chi_{1}$ in $\mathbb{C}_{0}^{2,6,\alpha}(\mathbb{T}^{1}\times\tilde{Y})$ . This is a standard result of parabolic PDEs in Hölder space (see, e.g., the Theorem 8.7.3 in [18]). The following theorem states that under certain conditions the solution of the discrete-type cell problem converges to the solution of the continuous one.

Theorem 4.6.

Assume $v_{1}$ is a smooth function defined on $\mathbb{T}^{1}\times\tilde{Y}$ , satisfying $\int_{\tilde{Y}}v_{1}(\tau,\textbf{x})\mathrm{d}\textbf{x}=0,\ \forall\tau\in\mathbb{T}^{1}$ . Let $\hat{v}_{1}$ and $\chi_{1}$ be the solutions of the discrete-type cell problem (37) and continuous cell problem (43), respectively. Then, we have the following convergence estimate holds

||\chi_{1}-\hat{v}_{1}||=\mathcal{O}(\Delta t),

(44)

where $||\cdot||$ is a function norm associated with the space $\mathbb{C}_{0}^{2,6,\alpha}(\mathbb{T}^{1}\times\tilde{Y})$ .

Proof.

Using Prop. 4.4, one can easily verify that $\mathcal{L}$ is a bijection between two Banach spaces $\mathbb{C}_{0}^{2,6,\alpha}(\mathbb{T}^{1}\times\tilde{Y})$ and $\mathbb{C}_{0}^{1,4,\alpha}(\mathbb{T}^{1}\times\tilde{Y})$ and its inverse is bounded. Integrating Eq.(43) along time gives,

\displaystyle\exp(\Delta t\mathcal{L})\chi_{1}-\chi_{1}=-v_{1}\Delta t+\mathcal{O}\big{(}(\Delta t)^{2}\big{)}\equiv-\Delta t\bar{v}_{1},

(45)

where $\bar{v}_{1}=v_{1}+O(\Delta t)$ . Combining Eqns.(40) and (45), we obtain

\exp(\Delta t\mathcal{L})\chi_{1}-I_{\Delta t}\hat{v}_{1}-(\chi_{1}-\hat{v}_{1})=\Delta t(v_{1}-\bar{v}_{1}).

(46)

Notice that Eq.(46) shows the connection between $\chi_{1}$ and $\hat{v}_{1}$ . After some simple calculations, we get that

\mathcal{L}(\chi_{1}-\hat{v}_{1})=(\mathcal{L}-\tilde{L}_{1})(\chi_{1}-\hat{v}_{1})+\tilde{L}_{2}\hat{v}_{1}+(v_{1}-\bar{v}_{1}),

(47)

where

\tilde{L}_{1}=\frac{\exp(\Delta t\mathcal{L})-I_{d}}{\Delta t},\quad\text{and}\quad\tilde{L}_{2}=\frac{I_{\Delta t}-\exp(\Delta t\mathcal{L})}{\Delta t}.

(48)

Moreover, we can verify that in the space of bounded linear operators from $\mathbb{C}_{0}^{2,6,\alpha}(\tilde{Y})$ to $\mathbb{C}_{0}^{1,4,\alpha}(\tilde{Y})$ , there is a strong convergence in the operator norm $||\cdot||$ ,

||\mathcal{L}-\tilde{L}_{1}||=\mathcal{O}(\Delta t)\quad\text{as }\Delta t\to 0.

(49)

For the operator $\tilde{L}_{2}$ , noticing that $\mathcal{L}=\mathcal{L}_{1}+\mathcal{L}_{2}+\mathcal{L}_{3}+\mathcal{L}_{4}$ and operator $I_{\Delta t}$ is defined in (20), we can use the BCH formula and obtain

	$\displaystyle\tilde{L}_{2}=$	$\displaystyle\frac{\exp\Big{(}\frac{(\Delta t)^{2}}{2}\big{(}[\mathcal{L}_{4},\mathcal{L}_{3}]+[\mathcal{L}_{4},\mathcal{L}_{2}]+[\mathcal{L}_{4},\mathcal{L}_{1}]+[\mathcal{L}_{3},\mathcal{L}_{2}]+[\mathcal{L}_{2},\mathcal{L}_{1}]+[\mathcal{L}_{3},\mathcal{L}_{1}]\big{)}+\mathcal{O}(\Delta t)^{3}\Big{)}-I_{d}}{\Delta t}\cdot\exp(\Delta t\mathcal{L})$
	$\displaystyle\to$	$\displaystyle\frac{\Delta t}{2}\big{(}([\mathcal{L}_{4},\mathcal{L}_{3}]+[\mathcal{L}_{4},\mathcal{L}_{2}]+[\mathcal{L}_{4},\mathcal{L}_{1}]+[\mathcal{L}_{3},\mathcal{L}_{2}]+[\mathcal{L}_{2},\mathcal{L}_{1}]+[\mathcal{L}_{3},\mathcal{L}_{1}]\big{)}+\mathcal{O}\big{(}(\Delta t)^{2}\big{)}.$		(50)

Denoting $\tilde{L}_{3}\equiv\tilde{L}_{1}+\tilde{L}_{2}=\frac{I_{\Delta t}-I_{d}}{\Delta t}$ , we have $\tilde{L}_{3}\to\mathcal{L}$ in $\mathcal{B}\big{(}\mathbb{C}_{0}^{2,6,\alpha}(\mathbb{T}^{1}\times\tilde{Y}),\mathbb{C}_{0}^{1,4,\alpha}(\mathbb{T}^{1}\times\tilde{Y})\big{)}$ as $\Delta t$ approaches zero. Then, applying the Prop. 4.4, we get,

\lim_{\Delta t\to 0}\hat{v}_{1}=\lim_{\Delta t\to 0}\tilde{L}_{3}^{-1}(-v_{1})=\mathcal{L}^{-1}(-v_{1})=\chi_{1}.

(51)

In addition, combining the results of the Eqns.(45), (49), (4.3) and (51) for the right hand side of Eq.(47), we know that when $\Delta t$ is small enough, the assertion in (44) is proved. The constant in the $\mathcal{O}(\Delta t)$ of (44) does not depend on the total computational time $T$ , but may depend on the regularities of $v_{1}$ , $v_{2}$ and the constant $\sigma$ . ∎

4.4 Convergence analysis for the effective diffusivity

This section contains the main results of our convergence analysis. We first prove that the second-order moment of the solution obtained by using our numerical scheme has an (at most) linear growth rate. Secondly, we provide the convergence rate of our numerical method in computing the effective diffusivity.

Theorem 4.7.

Let $\textbf{X}^{n}=(x^{n}_{1},x^{n}_{2})^{T}$ denote the solution of the two-dimensional passive tracer model (9) obtained by using our numerical scheme (32) with time step $\Delta t$ . If the Hamiltonian function $H(t,x_{1},x_{2})$ is separable, periodic and smooth (in order to guarantee the existence and uniqueness of the solution to the SDE (9)), then we can prove that the second-order moment of the solution $\textbf{X}^{n}$ (which can be viewed as a discrete Markov process) is at most linear growth, i.e.,

\max_{n}\big{\{}\mathbb{E}\frac{||\textbf{X}^{n}||^{2}}{n}\big{\}}\ \text{is bounded.}

(52)

Proof.

We first estimate the second-order moment of the first component of $\textbf{X}^{n}=(x^{n}_{1},x^{n}_{2})^{T}$ , since the other one can be estimated in the same manner. Simple calculations show that

	$\displaystyle\mathbb{E}[(x^{n}_{1})^{2}\|(x_{1}^{n-1},x_{2}^{n-1})]$	$\displaystyle=\mathbb{E}\big{(}x_{1}^{n-1}+v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})\Delta t+\sigma N_{1}^{n-1}\big{)}^{2}$
		$\displaystyle=\mathbb{E}(x_{1}^{n-1})^{2}+\Delta t\big{(}\sigma^{2}+2\mathbb{E}[x_{1}^{n-1}v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})]\big{)}+(\Delta t)^{2}\mathbb{E}v_{1}^{2}(t_{n-\frac{1}{2}},x_{2}^{n-1}).$		(53)

The term $\mathbb{E}[x_{1}^{n-1}v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})]$ corresponds to the strength of the convection-enhanced diffusion. Our goal here is to prove that it is bounded over $n$ , though it may depend on $v_{1}$ , $v_{2}$ and $\sigma$ . Notice that we are calculating the expectation of $(x_{1}^{n})^{2}$ , which is not defined in the torus space. But in the following derivation, we will show that it can be decomposed into sums of periodic functions acting on $\textbf{X}^{n}=(x_{1}^{n},x_{2}^{n})^{T}$ . Hence after the decomposition (see Eq.(56)) we can still apply the previous analysis on torus space.

We now directly compute the contribution of the term $\mathbb{E}[x_{1}^{n-1}v_{1}(t_{n-\frac{1}{2}},x_{2}^{n-1})]$ to the effective diffusivity with the help of Eq.(39),

\displaystyle\Delta t\sum_{i=0}^{n-1}\mathbb{E}[x_{1}^{i}v_{1}(t_{i+\frac{1}{2}},x_{2}^{i})]=\sum_{i=0}^{n-1}\mathbb{E}\big{[}x_{1}^{i}\big{(}\hat{v}_{1}(t_{i},\textbf{X}^{i})-\mathbb{E}[\hat{v}_{1}(t_{i+1},\textbf{X}^{i+1})|\textbf{X}^{i}]\big{)}\big{]}.

(54)

Let $\mathcal{F}_{i}$ denote the filtration generated by the solution process until $\textbf{X}^{i}$ . Notice that $x_{1}^{i}\in\mathcal{F}_{i}$ . For the Eq.(54), we have

RHS	$\displaystyle=\sum_{i=0}^{n-1}\mathbb{E}\big{[}x_{1}^{i}\big{(}\hat{v}_{1}(t_{i},\textbf{X}^{i})-\hat{v}_{1}(t_{i+1},\textbf{X}^{i+1})\big{)}\big{]}$
	$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})(x_{1}^{i}-x_{1}^{i-1})\big{]}+\hat{v}_{1}(t_{0},\textbf{X}^{0})x_{1}^{0}-\mathbb{E}[\hat{v}_{1}(t_{n},\textbf{X}^{n})x^{n}_{1}]$
	$\displaystyle=\sum_{i=1}^{n}\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}\big{]}+\hat{v}_{1}(t_{0},\textbf{X}^{0})x_{1}^{0}-\mathbb{E}[\hat{v}_{1}(t_{n},\textbf{X}^{n})x^{n}_{1}].$	(55)

Hence, we obtain the following result

$\displaystyle\frac{1}{n}\mathbb{E}\big{[}(x^{n}_{1})^{2}\|(x_{1}^{0},x_{2}^{0})\big{]}=$	$\displaystyle\frac{1}{n}(x_{1}^{0})^{2}+\Delta t\sigma^{2}+2\Delta t\frac{1}{n}\sum_{i=0}^{n-1}\mathbb{E}[x_{1}^{i}v_{1}(t_{i+\frac{1}{2}},x_{2}^{i})]+(\Delta t)^{2}\frac{1}{n}\sum_{i=0}^{n-1}\mathbb{E}v_{1}^{2}(t_{i+\frac{1}{2}},x_{2}^{i})$
$\displaystyle=$	$\displaystyle\frac{1}{n}(x_{1}^{0})^{2}+\Delta t\sigma^{2}+(\Delta t)^{2}\frac{1}{n}\sum_{i=0}^{n-1}\mathbb{E}v_{1}^{2}(t_{i+\frac{1}{2}},x_{2}^{i})$
	$\displaystyle+\frac{2}{n}\sum_{i=1}^{n}\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}\big{]}$
	$\displaystyle+\frac{2}{n}\big{(}\hat{v}_{1}(t_{0},\textbf{X}^{0})x_{1}^{0}-\mathbb{E}[\hat{v}_{1}(t_{n},\textbf{X}^{n})x^{n}_{1}]\big{)}.$	(56)

By using the Cauchy-Schwarz inequality, we know the term

	$\displaystyle\frac{2}{n}\sum_{i=1}^{n}\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}\big{]}$
$\displaystyle\leq$	$\displaystyle\frac{2}{n}\sum_{i=1}^{n}\mathbb{E}\big{[}2(\hat{v}_{1}(t_{i},\textbf{X}^{i}))^{2}+\big{(}(v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t)^{2}+(\sigma N^{i-1}_{1})^{2}\big{)}\big{]}$
$\displaystyle=$	$\displaystyle\frac{2}{n}\sum_{i=1}^{n}\mathbb{E}\big{[}2(\hat{v}_{1}(t_{i},\textbf{X}^{i}))^{2}+(v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1}))^{2}(\Delta t)^{2}+\sigma^{2}\Delta t\big{]}.$	(57)

Notice that if $v_{1}$ and $\hat{v}_{1}$ are bounded in sup norm, right-hand-side of Eq.(57) is bounded for any $n$ . Other terms on the right-hand side of Eq.(56) are also bounded, which can be checked easily. Therefore, we can prove that $\frac{1}{n}\mathbb{E}\big{[}(x^{n}_{1})^{2}|(x_{1}^{0},x_{2}^{0})\big{]}$ is bounded. Repeat the same trick, we know that $\frac{1}{n}\mathbb{E}\big{[}(x^{n}_{2})^{2}|(x_{1}^{0},x_{2}^{0})\big{]}$ is also bounded. Thus, the assertion in Eq.(52) is proved. ∎

In practice, we shall first choose a time step $\Delta t$ and run our numerical scheme (15) to compute the effective diffusivity until the result converges to a constant, which may depend on $\Delta t$ . As such, we shall prove that the limit of the constant converges to the exact effective diffusivity of the original passive tracer model as $\Delta t$ approaches zero. Namely, we shall prove that our numerical scheme is robust in computing effective diffusivity. More details on the numerical results will be given in Section 5.

Theorem 4.8.

Let $x^{n}_{1}$ , $n=0,1,....$ be the first component of the numerical solution obtained by using the scheme (15) and $\Delta t$ denote the time step. We have the convergence estimate of the effective diffusivity as

\displaystyle\lim_{n\to\infty}\frac{\mathbb{E}(x^{n}_{1})^{2}}{n\Delta t}=\sigma^{2}+2\int_{\mathbb{T}^{2}}\chi_{1}v_{1}+\mathcal{O}(\Delta t),

(58)

where the constant in $\mathcal{O}(\Delta t)$ may depend on the regularity of $v_{1}$ , $v_{2}$ and the constant $\sigma$ , but does not depend on the computational time $T$ .

Proof.

We will prove the statement by direct computation. We divide both sides of the Eq.(56) by $\Delta t$ and obtain

$\displaystyle\frac{1}{n\Delta t}\mathbb{E}\big{[}(x^{n}_{1})^{2}\|(x_{1}^{0},x_{2}^{0})\big{]}=$	$\displaystyle\frac{1}{n\Delta t}(x_{1}^{0})^{2}+\sigma^{2}+\frac{\Delta t}{n}\sum_{i=0}^{n-1}\mathbb{E}v_{1}^{2}(t_{i+\frac{1}{2}},x_{2}^{i})$
	$\displaystyle+\frac{2}{n\Delta t}\sum_{i=1}^{n}\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}\big{]}$
	$\displaystyle+\frac{2}{n\Delta t}\big{(}\hat{v}_{1}(t_{0},\textbf{X}^{0})x_{1}^{0}-\mathbb{E}[\hat{v}_{1}(t_{n},\textbf{X}^{n})x^{n}_{1}]\big{)}.$	(59)

First, we notice that for a fixed $\Delta t$ , the terms $\frac{1}{n\Delta t}(x_{1}^{0})^{2}$ and $\frac{2}{n\Delta t}\hat{v}_{1}(t_{0},\textbf{X}^{0})x_{1}^{0}$ converge to zero as $n\to\infty$ , where we have used the fact that $\hat{v}_{1}(t_{0},\textbf{X}^{0})$ is bounded. Also notice that the term $\frac{\Delta t}{n}\sum_{i=0}^{n-1}\mathbb{E}v_{1}^{2}(t_{i+\frac{1}{2}},x_{2}^{i})$ is $\mathcal{O}(\Delta t)$ , due to the term $(v_{1})^{2}$ is bounded. Then, for a fixed $\Delta t$ , we have

\lim_{n\to\infty}\frac{2}{n\Delta t}\big{|}\mathbb{E}[\hat{v}_{1}(\textbf{X}^{n})x^{n}_{1}]\big{|}\leq\lim_{n\to\infty}\frac{2}{\sqrt{n}\Delta t}||\hat{v}_{1}||_{\infty}\mathbb{E}|\frac{x^{n}_{1}}{\sqrt{n}}|\leq\lim_{n\to\infty}\frac{1}{\sqrt{n}\Delta t}||\hat{v}_{1}||_{\infty}\mathbb{E}\big{[}\frac{(x^{n}_{1})^{2}}{n}+1\big{]}=0,

(60)

where the term $\mathbb{E}[\frac{(x^{n}_{1})^{2}}{n}]$ is bounded due to the Theorem 4.7 and $||\hat{v}_{1}||_{\infty}\to||\chi_{1}||_{\infty}<\infty$ due to the Theorem 4.6.

Therefore, we only need to focus on the estimate of terms in the second line of Eq.(4.8), which corresponds to the convection-enhanced diffusion effect. Notice that $\hat{v}_{1}\in\mathbb{C}^{2,6,\alpha}$ , we compute the Ito-Taylor series approximation of $\hat{v}_{1}(t_{i},\textbf{X}^{i})$ ,

$\displaystyle\hat{v}_{1}(t_{i},\textbf{X}^{i})$	$\displaystyle=\hat{v}_{1}(t_{i-1},\textbf{X}^{i-1})+\hat{v}_{1,x_{1}}(t_{i-1},\textbf{X}^{i-1})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}$
	$\displaystyle+\hat{v}_{1,x_{2}}(t_{i-1},\textbf{X}^{i-1})\big{(}v_{2}\big{(}t_{i-\frac{1}{2}},x_{1}^{i-1}\big{)}\Delta t+\sigma N^{i-1}_{2}\big{)}$
	$\displaystyle+\frac{1}{2}\big{(}\hat{v}_{1,x_{1}x_{1}}(t_{i-1},\textbf{X}^{i-1})+\hat{v}_{1,x_{2}x_{2}}(t_{i-1},\textbf{X}^{i-1})\big{)}\sigma^{2}\Delta t+\mathcal{O}\big{(}(\Delta t)^{2}\big{)},$	(61)

where we have used the fact that $v_{2}\big{(}t_{i-\frac{1}{2}},x_{1}^{i-1}+v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t\big{)}=v_{2}\big{(}t_{i-\frac{1}{2}},x_{1}^{i-1}\big{)}+\mathcal{O}(\Delta t)$ , when $\Delta t$ is small and $v_{2}$ is smooth. Since $\hat{v}_{1}\to\chi_{1}$ in $\mathbb{C}_{0}^{2,6,\alpha}$ , the truncated term $\mathcal{O}\big{(}(\Delta t)^{2}\big{)}$ in Eq.(4.8) is uniformly bounded when $\Delta t$ is small enough. Substituting the Taylor expansion of $\hat{v}_{1}(t_{i},\textbf{X}^{i})$ in Eq.(4.8) into the target term of our estimate (i.e., terms in the second line of Eq.(4.8)), we get

$\displaystyle\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})$	$\displaystyle(v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1})\big{]}=\mathbb{E}\Big{[}\Big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\Big{)}\cdot$
	$\displaystyle\Big{(}\hat{v}_{1}(t_{i-1},\textbf{X}^{i-1})+\hat{v}_{1,x_{1}}(t_{i-1},\textbf{X}^{i-1})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}$
	$\displaystyle+\hat{v}_{1,x_{2}}(t_{i-1},\textbf{X}^{i-1})\big{(}v_{2}(t_{i-\frac{1}{2}},x_{1}^{i-1})\Delta t+\sigma N^{i-1}_{2}\big{)}$
	$\displaystyle+\frac{1}{2}\big{(}\hat{v}_{1,x_{1}x_{1}}(t_{i-1},\textbf{X}^{i-1})+\hat{v}_{1,x_{2}x_{2}}(t_{i-1},\textbf{X}^{i-1})\big{)}\sigma^{2}\Delta t+\mathcal{O}((\Delta t)^{2})\Big{)}\Big{]}.$	(62)

Combining the terms with the same order of $\Delta t$ , we obtain

		$\displaystyle\mathbb{E}\big{[}\hat{v}_{1}(t_{i},\textbf{X}^{i})\big{(}v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})\Delta t+\sigma N^{i-1}_{1}\big{)}\big{]}$
	$\displaystyle=$	$\displaystyle\Delta t\mathbb{E}\big{[}\hat{v}_{1}(t_{i-1},\textbf{X}^{i-1})v_{1}(t_{i-\frac{1}{2}},x_{2}^{i-1})+\sigma^{2}\hat{v}_{1,x_{1}}(t_{i-1},\textbf{X}^{i-1})\big{]}+\mathcal{O}((\Delta t)^{2}),$		(63)

where we have used the facts that: (1) $\textbf{X}^{i-1}$ is independent of $N^{i-1}_{1}$ and $N^{i-1}_{2}$ so the expectations of the corresponding terms vanish; (2) $N^{i-1}_{1}$ and $N^{i-1}_{2}$ are independent so $\mathbb{E}(N^{i-1}_{1}N^{i-1}_{2})=0$ ; and (3) $\mathbb{E}(N^{i-1}_{1})^{2}=\Delta t$ .

Finally, by using the Theorem 4.3 and noticing the invariant measure is the uniform measure, we obtain from Eq.(4.8) that

\displaystyle\lim_{n\to\infty}\frac{1}{n\Delta t}\mathbb{E}\big{[}(x^{n}_{1})^{2}|(x_{1}^{0},x_{2}^{0})\big{]}=\sigma^{2}+2\int(\hat{v_{1}}v_{1}+\sigma^{2}\hat{v}_{1,x_{1}})+\mathcal{O}(\Delta t).

(64)

Thus, our statement in the Eq.(58) is proved using the facts that $\hat{v}_{1}$ converges to $\chi_{1}$ (see Theorem 4.6) and $\int\hat{v}_{1,x_{1}}=0$ . ∎

Remark 4.3.

If we divide two on both sides of the Eq.(58), we can find that our result recovers the definition of the effective diffusivity $D^{E}_{11}$ defined in the Eq.(4). Recall that $D_{0}=\sigma^{2}/2$ . Theorem 4.8 reveals the connection of the definition of the effective diffusivity using the Eulerian framework and Lagrangian framework; see Eq.(4) and Eq.(8), which is fundamental in this context. For 3D time-dependent flow problems, the Eulerian framework has good theoretical values but the Lagrangian framework is computationally accessible.

Remark 4.4.

For the second component of the numerical solution, i.e., $x^{n}_{2}$ , $n=0,1,...$ , we can obtain the similar convergence result in computing the effective diffusivity. First we consider $\tilde{v}_{2}(t,x_{1},x_{2}):=v_{2}(t,x_{1}+v_{1}(t,x_{2})\Delta t)$ and notice that $\tilde{v}_{2}-v_{2}=\mathcal{O}(\Delta t)$ and $\int_{\mathbb{T}}\tilde{v}_{2}\mathrm{d}x_{1}=0$ . The remaining part of the proof is essentially the same as the results obtained in Sections 4.2, 4.3 and 4.4, so we skip the details here.

4.5 Generalizations to high-dimensional cases

To show the essential idea of our probabilistic approach in proving the convergence rate of the numerical schemes, we have carried out our convergence analysis based on a two-dimensional model problem (9). In fact, the extension of our approach to higher-dimensional problems is straightforward. Now we consider a high-dimensional problem as follow,

\mathrm{d}\textbf{X}=\textbf{v}(t,\textbf{X})\mathrm{d}t+\Sigma\mathrm{d}\textbf{w}(t),

(65)

where $\textbf{X}=(x_{1},x_{2},\cdots,x_{d})^{T}\in\mathbb{R}^{d}$ is the position of a particle, $\textbf{v}=(v_{1},v_{2},\cdots,v_{d})^{T}\in\mathbb{R}^{d}$ is the Eulerian velocity field at position X, $\Sigma$ is a $d\times d$ constant non-singular matrix, and $\textbf{w}(t)$ is a $d$ -dimension Brownian motion vector. In particular, we assume the component $v_{i}$ does not depend on $x_{i}$ , $i=1,...,d$ . Thus, the incompressible condition for $\textbf{v}(t,\textbf{X})$ (i.e. $\nabla_{\textbf{X}}\cdot\textbf{v}(t,\textbf{X})=0$ ) is easily guaranteed.

For a deterministic and divergence-free dynamical system, Feng et. al. proposed a volume-preserving method [10], which splits a $d$ -dimensional problem into $d-1$ subproblems with each of them being a two-dimensional problem and thus being volume-preserving. We shall modify Feng’s method (first-order case) by including the randomness as the last subproblem to take into account the additive noise, i.e.,

{\color[rgb]{0,0,0}\begin{cases}x_{1}^{*}=x_{1}^{n-1}+\Delta tv_{1}(t_{n}+\frac{\Delta t}{2},x_{2}^{n-1},x_{3}^{n-1},x_{4}^{n-1},\cdots,x_{d-1}^{n-1},x_{d}^{n-1}),\\ x_{2}^{*}=x_{2}^{n-1}+\Delta tv_{2}(t_{n}+\frac{\Delta t}{2},x_{1}^{*},x_{3}^{n-1},x_{4}^{n-1},\cdots,x_{d-1}^{n-1},x_{d}^{n-1}),\\ x_{3}^{*}=x_{3}^{n-1}+\Delta tv_{3}(t_{n}+\frac{\Delta t}{2},x_{1}^{*},x_{2}^{*},x_{4}^{n-1},\cdots,x_{d-1}^{n-1},x_{d}^{n-1}),\\ \cdots,\\ x_{d}^{*}=x_{d}^{n-1}+\Delta tv_{d}(t_{n}+\frac{\Delta t}{2},x_{1}^{*},x_{2}^{*},x_{3}^{*},x_{4}^{*},\cdots,x_{d-1}^{*}),\\ \textbf{X}^{n}=\textbf{X}^{*}+\Sigma(\textbf{W}^{n}-\textbf{W}^{n-1}),\end{cases}}

(66)

where $\textbf{X}^{*}=(x_{1}^{*},x_{2}^{*},\cdots,x_{d}^{*})^{T}$ , $\textbf{W}^{n}-\textbf{W}^{n-1}$ is a $d$ -dimensional independent random vector with each component of the form $\sqrt{\Delta t}\xi_{i}$ , $\xi_{i}\sim\mathcal{N}(0,1)$ , and $\textbf{X}^{n}=(x_{1}^{n},x_{2}^{n},\cdot\cdot\cdot,x_{d}^{n})^{T}$ is the numerical approximation to the exact solution $\textbf{X}(t_{n})$ to the SDE (65) at time $t_{n}=n\Delta t$ .

The techniques of the convergence analysis for the two-dimensional problem can be applied to high-dimensional problems without much difficulty. For the high-dimensional problem (65), the smoothness and strict positivity of the transition kernel in the discrete process can be guaranteed if one assumes that the covariance matrix $\Sigma$ is non-singular and the scheme (66) is explicit. According to our assumption for the velocity field, the scheme (66) is volume-preserving for each step. Thus, the solution to the first-order modified equation is divergence-free and the invariant measure on the torus (defined by $\mathbb{R}^{d}/\mathbb{Z}^{d}$ , when the period is $1$ ) remains uniform for all $t$ . Finally, the convergence of the cell problem can be studied by using the BCH formula (19) with $d+2$ differential operators. Recall that in the Eq.(20) we have four differential operators when we study the two-dimensional problem.

Therefore, our numerical methods are robust in computing effective diffusivity for high-dimensional problems, which will be demonstrated through time-dependent chaotic flow problems in three-dimensional space in Section 5.

5 Numerical results

In this section, we will present numerical examples to verify the convergence analysis of the proposed method in computing effective diffusivity for time-dependent chaotic flows. In addition, we will investigate the convection-enhanced diffusion phenomenon in 3D time-dependent flow, i.e., the time-dependent ABC flow and the time-dependent Kolmogorov flow. Without loss of generality, we compute the quantity $\frac{\mathbb{E}(x_{1}(T))^{2}}{2T}$ , which is used to approximate $D^{E}_{11}$ in the effective diffusivity matrix (4).

5.1 Verification of the convergence rate

We first consider a two-dimensional passive tracer model. Let $(x_{1},x_{2})^{T}\in\mathbb{R}^{2}$ denote the position of a particle. Its motion is described by the following SDE,

\begin{cases}\mathrm{d}x_{1}=\sin\big{(}4x_{2}+1+\sin(2\pi t)\big{)}\exp\big{(}\cos\big{(}4x_{2}+1+\sin(2\pi t)\big{)}\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{1,t},\\ \mathrm{d}x_{2}=\cos\big{(}2x_{1}+\sin(2\pi t)\big{)}\exp\big{(}\sin\big{(}2x_{1}+\sin(2\pi t)\big{)}\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{2,t},\end{cases}

(67)

where $\sigma=\sqrt{2\times 0.1}$ , $w_{i,t}$ , $i=1,2$ are independent Brownian motions, and the initial data $(x_{1}^{0},x_{2}^{0})^{T}$ follows uniform distributions in $[-0.5,0.5]^{2}$ . One can easily verify the velocity field in (67) is time-dependent and divergence-free.

In our numerical experiments, we use Monte Carlo samples to discretize the Brownian motions $w_{1,t}$ and $w_{2,t}$ . The sample number is denoted by $N_{mc}$ . We choose $\Delta t_{ref}=\frac{1}{2^{12}}$ and $N_{mc}=3,200,000$ to solve the SDE (67) to compute the reference solution, i.e., the “exact” effective diffusivity, where the final computational time is $T=3000$ so that the calculated effective diffusivity converges to a constant. In fact, we find that the passive tracer model will enter a mixing stage if the computational time is bigger than $T=1000$ . It takes about $17$ hours to compute the reference solution on a 80-core server (HPC2015 System at HKU). The reference solution for the effective diffusivity is $D^{E}_{11}=0.219$ .

In Fig.1, we plot the convergence results of the effective diffusivity using our method (i.e., $\frac{\mathbb{E}(x_{1}(T))^{2}}{2T}$ ) with respective to different time-step $\Delta t$ at $T=3000$ . In addition, we show a fitted straight line with the slope $1.04$ , i.e., the convergence rate is about $(\Delta t)^{1.04}$ . This numerical result verifies the convergence analysis in Theorem 4.8.

Refer to caption — Figure 1: Error of $D^{E}_{11}$ for two time-dependent flows with different time-steps. (a) 2D time-dependent chaotic flow, fitted slope $\approx 1.04$ ; (b) 3D time-dependent Kolmogorov flow, fitted slope $\approx 1.22$ .

To further study the accuracy and robustness of our method for long-time integration, we consider a 3D time-dependent Kolmogorov flow problem. Let $(x_{1},x_{2},x_{3})^{T}\in\mathbb{R}^{3}$ denote the position of a particle. The motion of a particle moving in the 3D time-dependent Kolmogorov flow is described by the following SDE,

\begin{cases}\mathrm{d}x_{1}=\sin\big{(}x_{3}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{1,t},\\ \mathrm{d}x_{2}=\sin\big{(}x_{1}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{2,t},\\ \mathrm{d}x_{3}=\sin\big{(}x_{2}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{3,t}.\end{cases}

(68)

where $w_{1,t}$ , $w_{2,t}$ and $w_{3,t}$ are independent Brownian motions. When $\epsilon=0$ , the velocity field in (68) corresponds to the Kolmogorov flow [11]. The Kolmogorov flow possesses very chaotic behaviors [6], which brings challenges for our method.

In our numerical experiment, we choose $\epsilon=10^{-1}$ and $\sigma=\sqrt{2\times 10^{-3}}$ in the Eq.(68). We choose $\Delta t_{ref}=\frac{1}{2048}$ and $N_{mc}=480,000$ to compute the reference solution for the SDE (68), i.e., the “exact” effective diffusivity. In our numerical tests, we find that the passive tracer model will enter a mixing stage if the computational time is bigger than $T=2000$ . To show the accuracy and robustness of our numerical scheme, we set $T=10^{5}$ here. It takes about 59 hours to compute the reference solution on the server and the reference solution for the effective diffusivity is $D^{E}_{11}=0.693$ .

In Fig.1, we plot the convergence results of the effective diffusivity using our method with respect to different time-step $\Delta t$ . In addition, we show a fitted straight line with the slope $1.22$ , i.e., the convergence rate is about $(\Delta t)^{1.22}$ . This numerical result again agrees with our error analysis.

5.2 Investigation of the convection-enhanced diffusion phenomenon

As we have already demonstrated in Section 5.1, our method is very accurate and robust for long-time integration. Here, we will study the dependence of the effective diffusivity $D^{E}_{11}$ on different parameters in the time-dependent flows. First of all, we solve Eq.(68) and carry out the test for the 3D time-dependent Kolmogorov flow.

In Fig.2, we show the time evolution of $\frac{\mathbb{E}(x_{1}(t))^{2}}{2t}$ for different $D_{0}$ ’s (here $D_{0}=\sigma^{2}/2$ ) and for four different $\epsilon$ ’s, where the result in Fig.2 corresponding to the time-independent Kolmogorov flow (see Figure 6 of [31]). Notice that in Eq.(68) the parameter $\epsilon$ controls the strength of the time dependence. For each $D_{0}$ and $\epsilon$ , we use $N_{mc}=240,000$ particles to solve the SDE (68). We find that for each given $D_{0}$ , the time evolution of $\frac{E(x_{1}(t))^{2}}{2t}$ converges as $\epsilon$ approaches zero. This can be rigorously justified through analysis; see A. In addition, we observe a certain amount of enhanced diffusion when $D_{0}$ decreases. However, the dependency of $D^{E}_{11}$ on $D_{0}$ is quite different from the pattern of the time-dependent ABC flow, which is known as the maximal enhancement and will be discussed later; see Fig.5.

To study the dependence of $D^{E}_{11}$ on $D_{0}$ and $\epsilon$ , we choose different $\epsilon$ ’s and $D_{0}$ ’s and compute the corresponding effective diffusivity $D^{E}_{11}$ . In this experiment, we use $\Delta t=2^{-7}$ and $N_{mc}=240,000$ particles to compute. The final computational time is $T=10^{5}$ so that the particles are fully mixed. We show the numerical results in Fig.3.

We find that for each given $D_{0}$ as $\epsilon$ decreases the corresponding effective diffusivity $D^{E}_{11}$ converges to the effective diffusivity $D^{E}_{11}$ associated with $\epsilon=0$ . This means the time dependency of $\epsilon$ improves the chaotic property of Kolmorogov flow though, it does not change the pattern of convection-enhanced diffusion in the Kolmogorov flow. When $\epsilon\leq 1$ the fitted slope within $D_{0}\in[10^{-5},10^{-3}]$ is $-0.2$ , which indicates that $D^{E}_{11}\sim\mathcal{O}(1/D_{0}^{0.2})$ . We call this behavior a sub-maximal enhancement, which may be explained by the fact that the Kolmogorov flow is more chaotic than the ABC flow [11]. The chaotic trajectories in Kolmogorov flow enhance diffusion much less than channel-like structures such as the ballistic orbits of ABC flows [23, 33].

Next, we use our stochastic structure-preserving scheme to solve time-dependent ABC flow problems. Let $(x_{1},x_{2},x_{3})^{T}\in\mathbb{R}^{3}$ denote the position of a particle in the 3D Cartesian coordinate system. The motion of a particle moving in the 3D time-dependent ABC flow is described by the following SDE,

\begin{cases}\mathrm{d}x_{1}=A\sin\big{(}x_{3}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+C\cos\big{(}x_{2}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{1,t},\\ \mathrm{d}x_{2}=B\sin\big{(}x_{1}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+A\cos\big{(}x_{3}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{2,t},\\ \mathrm{d}x_{3}=C\sin\big{(}x_{2}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+B\cos\big{(}x_{1}+\epsilon\sin(2\pi t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{3,t},\end{cases}

(69)

where $w_{1,t}$ , $w_{2,t}$ and $w_{3,t}$ are independent Brownian motions. For $\epsilon=0$ and $\sigma=0$ , the velocity field in (69) corresponds to the standard ABC flow [8]. The ABC flow is a three-dimensional incompressible velocity field which is an exact solution to the Euler’s equation. It is notable as a simple example of a fluid flow that can have chaotic trajectories. In our numerical experiments, we set $A=B=C=1$ .

In Fig.4, we show the time evolution of the $\frac{\mathbb{E}(x_{1}(t))^{2}}{2t}$ for different $D_{0}$ ’s (here $D_{0}=\sigma^{2}/2$ ) and for four different $\epsilon$ ’s, where the result in Fig.4 corresponding to the time-independent ABC flow (see Figure 3 of [31]). Again the parameter $\epsilon$ controls the strength of the time dependence. For each $D_{0}$ and $\epsilon$ , we use $N_{mc}=240,000$ particles to solve the SDE (69). We find that for each given $D_{0}$ , the time evolution of the $\frac{\mathbb{E}(x_{1}(t))^{2}}{2t}$ converges when $\epsilon$ converges to zero. However, we observe two different patterns compared with the results shown in Fig.2. First, when we decrease $D_{0}$ , it takes a longer time for the system to enter a mixing stage. Second, we observe a large amount of enhanced diffusion when $D_{0}$ decreases.

To further investigate the dependence of $D^{E}_{11}$ on $D_{0}$ and $\epsilon$ , we choose different $\epsilon$ ’s and $D_{0}$ ’s and compute the corresponding effective diffusivity $D^{E}_{11}$ . In this experiment, we use $\Delta t=2^{-7}$ and $N_{mc}=240,000$ particles to compute. The final computational time is $T=10^{5}$ so that the particles are fully mixed.

In Fig.5, we show the numerical results. We find that for each given $D_{0}$ , as $\epsilon$ decreases the corresponding effective diffusivity $D^{E}_{11}$ converges to the effective diffusivity $D^{E}_{11}$ associated with $\epsilon=0$ . Thus, the time-dependent ABC flow has a similar convection-enhanced diffusion behavior as the time-independent ABC flow. The fitted slope within $D_{0}\in[10^{-5},10^{-1}]$ is about $-1.0$ , which indicates that $D^{E}_{11}\sim\mathcal{O}(1/D_{0}^{1})$ . This result indicates that the $D^{E}_{11}$ of the time-dependent ABC flow achieves the upper-bound of Eq.(7), i.e. the maximal enhancement. This maximal enhancement phenomenon may be attributed to the ballistic orbits of the ABC flow, where the time-independent case was discussed in [23, 33].

Moreover, our result for $D_{0}\in[10^{-3},10^{-1}]$ and $\epsilon=0$ recovers the same phenomenon as the Fig.2 in [4], which was obtained by using the Eulerian framework, i.e., solving a cell problem. In Fig.5, our method can be easily used to compute the effective diffusivity when $D_{0}\in[10^{-5},10^{-4}]$ . It will be, however, extremely expensive for the Eulerian framework since one needs to solve a convection-dominated PDE (5) in 3D space, whose Péclet number is proportion to $\frac{1}{D_{0}}$ .

Finally, we investigate the dependence of $D^{E}_{11}$ on the frequency of the time-dependent ABC flow. Specifically, we solve the following SDE,

\begin{cases}\mathrm{d}x_{1}=A\sin\big{(}x_{3}+\sin(\Omega t)\big{)}\mathrm{d}t+C\cos\big{(}x_{2}+\sin(\Omega t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{1,t},\\ \mathrm{d}x_{2}=B\sin\big{(}x_{1}+\sin(\Omega t)\big{)}\mathrm{d}t+A\cos\big{(}x_{3}+\sin(\Omega t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{2,t},\\ \mathrm{d}x_{3}=C\sin\big{(}x_{2}+\sin(\Omega t)\big{)}\mathrm{d}t+B\cos\big{(}x_{1}+\sin(\Omega t)\big{)}\mathrm{d}t+\sigma\mathrm{d}w_{3,t},\end{cases},

(70)

where $A=B=C=1$ and $\Omega$ is the frequency. Here we first choose $\Delta t=2^{-7}$ , $N_{mc}=240,000$ and $T=10^{5}$ . Then, we choose different $\Omega$ and compute the corresponding effective diffusivity $D^{E}_{11}$ .

In Fig.6, we show the numerical results. We find that when $\Omega$ is near $0.1$ the diffusion enhancement is weak. When $\Omega$ is away from $0.1$ , say $\Omega<0.05$ or $\Omega>0.2$ , we observe the maximal enhancement phenomenon. A similar sensitive dependence on the frequency of time-dependent ABC flows was reported in [5], where the Lyapunov exponent of the deterministic time-dependent ABC flow problem (i.e., $\sigma=0$ in Eq. (69)) was studied as the indicator of the extent of chaos; see Fig.2 and Fig.3 of [5].

When $\Omega=0$ , the flow of (70) is the same as that for $\epsilon=0$ case in (69), which will give the maximal enhancement phenomenon. When $\Omega$ is positive, the flow becomes time-dependent and the regions of chaos expand until the extent of chaos (i.e. the Lyapunov exponent) appears to reach a maximum, which is corresponding to $\Omega=0.1$ . It seems that the diffusion enhancement is significantly weakened in this range of $\Omega$ . When $\Omega$ continues to grow, the islands of the integrability regrow and the chaotic regions have shrunk significantly. We again observe the maximal enhancement phenomenon in this range of $\Omega$ . Our numerical results suggest that the level of chaos and the strength of diffusion enhancement seem to compete with each other. More intensive theoretic and numerical studies will be reported in our future work.

6 Conclusion

In this paper, we developed a robust stochastic structure-preserving Lagrangian scheme in computing effective diffusivity of passive tracer models and provided a sharp convergence analysis on the proposed numerical scheme. Our convergence analysis is based on a probabilistic approach, which interprets the solution process generated by our numerical scheme as a Markov process. By exploring the ergodicity of the solution process, we gave a sharp and uniform-in-time error estimate for our numerical scheme, which allows us to compute the effective diffusivity over infinite time. Numerical results verify that the proposed method is robust and accurate in computing effective diffusivity of time-dependent chaotic flows. We observed the maximal enhancement phenomenon in time-dependent ABC flows and the sub-maximal enhancement phenomenon in time-dependent Kolmogorov flows, respectively. Moreover, we found that the time dependency in the velocity field improves the chaotic property of ABC flow and Kolmorogov flow though, it does not change the pattern of convection-enhanced diffusion in both flows.

There are two directions we plan to explore in our future work. First, we intend to study the convection-enhanced diffusion phenomenon and provide a sharp convergence analysis for general time-dependent chaotic flows, where the flows have a quasi-periodic property in the time domain. In addition, we shall investigate the convection-enhanced diffusion phenomenon for general spatial-temporal stochastic flows [19, 22] and develop convergence analysis for the corresponding numerical methods.

Acknowledgement

The research of Z. Wang is partially supported by the Hong Kong PhD Fellowship Scheme. The research of J. Xin is partially supported by NSF grants DMS-1924548 and DMS-1952644. The research of Z. Zhang is supported by Hong Kong RGC grants (Projects 17300817 and 17300318), Seed Funding Programme for Basic Research (HKU), and Basic Research Programme (JCYJ20180307151603959) of The Science, Technology and Innovation Commission of Shenzhen Municipality. The computations were performed using research computing facilities offered by Information Technology Services, the University of Hong Kong.

Appendix A Limit of the parameter $\epsilon$ in a time-dependent chaotic flow

We shall prove that when $\epsilon$ approaches zero, the effective diffusivity corresponding to the time-dependent chaotic flow, e.g. the flow in (68) will converge to the one corresponding to the time-independent one, e.g. $\epsilon=0$ in the flow of (68). For notational simplicity, let $\textbf{v}=\textbf{v}^{\epsilon}$ denote the velocity field in (68) and $\textbf{v}=\textbf{v}^{0}$ denote the velocity field when $\epsilon=0$ in $\textbf{v}=\textbf{v}^{\epsilon}$ . Moreover, we denote $\mathcal{L}^{\epsilon}(\cdot)=\textbf{v}^{\epsilon}\cdot\nabla_{x}(\cdot)+D_{0}\Delta_{x}(\cdot)$ . Now, the vector corrector field $\boldsymbol{\chi}^{\epsilon}$ associated with the velocity field $\textbf{v}^{\epsilon}$ satisfies the following cell problem,

(\partial_{\tau}+\mathcal{L}^{\epsilon})\boldsymbol{\chi}^{\epsilon}=-\textbf{v}^{\epsilon}.

(71)

Let $\boldsymbol{\chi}^{\epsilon}_{0}$ denote the solution of the following equation

(\partial_{\tau}+\mathcal{L}^{\epsilon})\boldsymbol{\chi}_{0}^{\epsilon}=-\textbf{v}^{0}.

(72)

We aim to prove $\boldsymbol{\chi}^{\epsilon}$ converges to $\boldsymbol{\chi}^{\epsilon}_{0}$ as $\epsilon$ approaches zero. At the same time we know the vector corrector field $\boldsymbol{\chi}_{0}$ associated with the velocity field $\textbf{v}^{0}$ satisfies the following cell problem,

\mathcal{L}^{0}\boldsymbol{\chi}_{0}=-\textbf{v}^{0},

(73)

where $\mathcal{L}^{0}(\cdot)=\textbf{v}^{0}\cdot\nabla_{x}(\cdot)+D_{0}\Delta_{x}(\cdot)$ . Now we consider, $\boldsymbol{\chi}_{0}^{0}(t,x)=\boldsymbol{\chi}_{0}(x)$ , which solves,

(\partial_{\tau}+\mathcal{L}^{0})\boldsymbol{\chi}_{0}^{0}=-\textbf{v}^{0},

(74)

since $\partial_{\tau}\boldsymbol{\chi}_{0}^{0}=0$ . Comparing Eqns.(72) and (74) and using Prop.4.4, we know that $\boldsymbol{\chi}^{\epsilon}_{0}$ converges to $\boldsymbol{\chi}_{0}^{0}$ when $\epsilon$ approaches zero. Finally, comparing Eqns.(71) and (72), we know $\boldsymbol{\chi}^{\epsilon}$ converges to $\boldsymbol{\chi}_{0}^{\epsilon}$ when $\epsilon$ approaches zero. Therefore, we prove $\boldsymbol{\chi}^{\epsilon}$ converges to $\boldsymbol{\chi}_{0}$ when $\epsilon$ approaches zero.

References

[1] B. Afkham and J. Hesthaven, Structure preserving model reduction of parametric hamiltonian systems, SIAM Journal on Scientific Computing, 39 (2017), pp. A2616–A2644.
[2] G. Ben Arous and H. Owhadi, Multiscale homogenization with bounded ratios and anomalous slow diffusion, Communications on Pure and Applied Mathematics, 56 (2003), pp. 80–113.
[3] A. Bensoussan, J. L. Lions, and G. Papanicolaou, Asymptotic analysis for periodic structures, vol. 374, American Mathematical Soc., 2011.
[4] L. Biferale, A. Crisanti, M. Vergassola, and A. Vulpiani, Eddy diffusivities in scalar transport, Phys. Fluids, 7 (1995), pp. 2725–2734.
[5] N. Brummell, F. Cattaneo, and S. Tobias, Linear and nonlinear dynamo properties of time-dependent ABC flows, Fluid Dynamics Research, 28 (2001), p. 237.
[6] S. Childress and A. D. Gilbert, Stretch, twist, fold: the fast dynamo, vol. 37, Springer Science & Business Media, 1995.
[7] A. Debussche and E. Faou, Weak backward error analysis for SDEs, SIAM Journal on Numerical Analysis, 50 (2012), pp. 1735–1752.
[8] T. Dombre, U. Frisch, J. Greene, M. Hénon, A. Mehr, and A. Soward, Chaotic streamlines in the ABC flows, Journal of Fluid Mechanics, 167 (1986), pp. 353–391.
[9] A. Fannjiang and G. Papanicolaou, Convection-enhanced diffusion for periodic flows, SIAM J Appl. Math., 54 (1994), pp. 333–408.
[10] K. Feng and Z. Shang, Volume-preserving algorithms for source-free dynamical systems, Numerische Mathematik, 71 (1995), pp. 451–463.
[11] D. Galloway and M. Proctor, Numerical calculations of fast dynamos in smooth velocity fields with realistic diffusion, Nature, 356 (1992), p. 691.
[12] J. Garnier, Homogenization in a periodic and time-dependent potential, SIAM Journal on Applied Mathematics, 57(1) (1997), pp. 95–111.
[13] R. Gilmore, Baker-Campbell-Hausdorff formulas, Journal of Mathematical Physics, 15 (1974), pp. 2090–2092.
[14] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, Springer Science and Business Media, 2006.
[15] J. Hong, H. Liu, and G. Sun, The multi-symplecticity of partitioned Runge-Kutta methods for Hamiltonian PDEs, Mathematics of computation, 75 (2006), pp. 167–181.
[16] V. V. Jikov, S. Kozlov, and O. A. Oleinik, Homogenization of Differential Operators and Integral Functionals, Springer, Berlin, 1994.
[17] T. Kato, Perturbation theory for linear operators, vol. 132, Springer Science & Business Media, 2013.
[18] N. V. Krylov, Lectures on elliptic and parabolic equations in Hölder spaces, Graduate studies in mathematics.
[19] C. Landim, S. Olla, and H. T. Yau, Convection–diffusion equation with space–time ergodic random flow, Probability theory and related fields, 112 (1998), pp. 203–220.
[20] T. Lelièvre and G. Stoltz, Partial differential equations and stochastic methods in molecular dynamics, Acta Numerica, 25 (2016), pp. 681–880.
[21] J. Lyu, J. Xin, and Y. Yu, Computing residual diffusivity by adaptive basis learning via spectral method, Numerical Mathematics: Theory, Methods and Applications, 10(2) (2017), pp. 351–372.
[22] A. J. Majda and P. R. Kramer, Simplified models for turbulent diffusion: theory, numerical modelling, and physical phenomena, Phys. Rep., 314 (1999), pp. 237–574.
[23] T. McMillen, J. Xin, Y. F. Yu, and A. Zlatos, Ballistic orbits and front speed enhancement for ABC flows, SIAM Journal on Applied Dynamical Systems, 15 (2016), pp. 1753–1782.
[24] I. Mezić, J. F. Brady, and S. Wiggins, Maximal effective diffusivity for time-periodic incompressible fluid flows, SIAM Journal on Applied Mathematics, 56 (1996), pp. 40–56.
[25] G. Milstein, Y. Repin, and M. Tretyakov, Symplectic integration of Hamiltonian systems with additive noise, SIAM J. Numer. Anal, 39 (2002), pp. 2066–2088.
[26] B. Oksendal, Stochastic Differential Equations: an introduction with applications., Springer Science and Business Media, 2013.
[27] G. Pavliotis and A. Stuart, Multiscale methods: averaging and homogenization, Springer Science and Business Media, 2008.
[28] G. Pavliotis, A. Stuart, and K. Zygalakis, Calculating effective diffusivities in the limit of vanishing molecular diffusion, J. Comput. Phys., 228 (2009), pp. 1030–1055.
[29] S. Reich, Backward error analysis for numerical integrators, SIAM J. Numer. Anal., 36 (1999), pp. 1549–1570.
[30] M. Tao, H. Owhadi, and J. Marsden, Nonintrusive and structure preserving multiscale integration of stiff ODEs, SDEs, and Hamiltonian systems with hidden slow dynamics via flow averaging, Multiscale Modeling & Simulation, 8 (2010), pp. 1269–1324.
[31] Z. Wang, J. Xin, and Z. Zhang, Sharp uniform in time error estimate on a stochastic structure-preserving Lagrangian method and computation of effective diffusivity in 3D chaotic flows, to appear in SIAM Multiscale Model. Simul., arXiv:1808.06309, (2018).
[32] Z. J. Wang, J. Xin, and Z. W. Zhang, Computing effective diffusivity of chaotic and stochastic flows using structure-preserving schemes, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2322–2344.
[33] J. Xin, Y. Yu, and A. Zlatos, Periodic orbits of the ABC flow with ${A}={B}={C}=1$ , SIAM Journal on Mathematical Analysis, 48 (2016), pp. 4087–4093.

Convergence analysis of a Lagrangian numerical scheme in computing effective diffusivity of 3D time-dependent flows

Abstract

keywords:

1 Introduction

2 Effective diffusivity of the passive tracer models

Remark 2.1.

3 Stochastic structure-preserving schemes

3.1 Derivation of numerical schemes

3.2 The backward Kolmogorov equation and related results

Remark 3.1.

4 Convergence analysis

4.1 Convergence to an invariant measure

Proposition 4.1 (Theorem 3.3.1 of [3]).

Lemma 4.2.

Proof.

Theorem 4.3.

Proof.

Proposition 4.4.

4.2 A discrete-type cell problem

Lemma 4.5.

Proof.

Remark 4.1.

Remark 4.2.

4.3 Convergence estimate of the discrete-type cell problem

Theorem 4.6.

Proof.

4.4 Convergence analysis for the effective diffusivity

Theorem 4.7.

Proof.

Theorem 4.8.

Proof.

Remark 4.3.

Remark 4.4.

4.5 Generalizations to high-dimensional cases

5 Numerical results

5.1 Verification of the convergence rate

5.2 Investigation of the convection-enhanced diffusion phenomenon

6 Conclusion

Acknowledgement

Appendix A Limit of the parameter ϵ\epsilon in a time-dependent chaotic flow

References

Appendix A Limit of the parameter $\epsilon$ in a time-dependent chaotic flow