# A convergent interacting particle method and computation of KPP front speeds
in chaotic flows

Junlong Lyu<EMAIL_ADDRESS>Zhongjian Wang
<EMAIL_ADDRESS>Jack Xin<EMAIL_ADDRESS>Zhiwen Zhang
<EMAIL_ADDRESS>Department of Mathematics, The University of Hong Kong,
Pokfulam Road, Hong Kong SAR, China.  Department of Statistics, The University
of Chicago, Chicago, IL 60637, USA.  Department of Mathematics, University of
California at Irvine, Irvine, CA 92697, USA.

###### Abstract

In this paper, we study the propagation speeds of reaction-diffusion-advection
(RDA) fronts in time-periodic cellular and chaotic flows with Kolmogorov-
Petrovsky-Piskunov (KPP) nonlinearity. We first apply the variational
principle to reduce the computation of KPP front speeds to a principal
eigenvalue problem of a linear advection-diffusion operator with space-time
periodic coefficients on a periodic domain. To this end, we develop efficient
Lagrangian particle methods to compute the principal eigenvalue through the
Feynman-Kac formula. By estimating the convergence rate of Feynman-Kac
semigroups and the operator splitting methods for approximating the linear
advection-diffusion solution operators, we obtain convergence analysis for the
proposed numerical methods. Finally, we present numerical results to
demonstrate the accuracy and efficiency of the proposed method in computing
KPP front speeds in time-periodic cellular and chaotic flows, especially the
time-dependent Arnold-Beltrami-Childress (ABC) flow and time-dependent
Kolmogorov flow in three-dimensional space.
AMS subject classification: 35K57, 47D08, 65C35, 65L20, 65N25.

###### keywords:

KPP front speeds; cellular and chaotic flows; Feynman-Kac semigroups;
interacting particle method; operator splitting methods; eigenvalue problems;
convergence analysis.

##  1 Introduction

Front propagation in complex fluid flows arises in many scientific areas such
as chemical kinetics, combustion, biology, transport in porous media, and
industrial deposition processes (see [46] for a review). A fundamental problem
is to analyze and compute large-scale front speeds in complex flows. An
extensively studied model problem is the reaction-diffusion-advection (RDA)
equation with Kolmogorov-Petrovsky-Piskunov (KPP) nonlinearity [22]. To be
specific, the KPP equation is

$\displaystyle
u_{t}=\kappa\Delta_{\textbf{x}}u+(\textbf{v}\cdot\nabla_{\textbf{x}})u+\tau^{-1}f(u),\quad
t\in\mathbb{R}^{+},\quad\textbf{x}=(x_{1},...,x_{d})^{T}\in\mathbb{R}^{d},$
(1)

where $\kappa$ is the diffusion constant, $\tau$ is the time scale of reaction
rate, v is an incompressible velocity field (its precise definition will be
discussed later), $u$ is the concentration of reactant or population, and the
KPP reaction $f(u)=u(1-u)$ satisfying $f(u)\leq uf^{\prime}(0)$. In our
analysis and numerical examples, we will keep $\tau$ and $\kappa$ fixed, while
change the magnitude of the velocity field v, which equivalently means
changing the Péclet number.

Since the pioneering work of Kolmogorov, Petrovsky, and Piskunov [22] and
Fisher [13] on traveling fronts of the reaction-diffusion equations, this
field has gone through enormous growth and development. Reaction-diffusion
front propagation in fluid flows has been an active research topic for
decades; see e.g.[16, 44, 24, 45, 3, 32, 33, 46, 29] and references therein.
Significant amounts of mathematical analysis and numerical works in this
direction have been accomplished when the streamlines of fluid flow are either
well-structured (regular motion) or fully random (ergodic motion). Yet, the
often encountered less studied case is when the streamlines consist of both
regular and irregular motions, while neither one takes up the entire phase
space, such as the chaotic Arnold-Beltrami-Childress (ABC) flow [10, 4] and
Kolmogorov flows [15, 7].

In recent years, much progress has been made in finite element computation of
the KPP front propagation in time-periodic cellular and chaotic flows based on
a linearized corrector equation. If the velocity field
$\textbf{v}=\textbf{v}(\textbf{x})$ in the KPP equation (1) is time-
independent, the minimal front speed in direction e is given by the
variational formula [16]:
$c^{*}(\textbf{e})=\inf_{\lambda>0}\mu(\lambda)/\lambda$, where $\mu(\lambda)$
is the principal eigenvalue of the elliptic operator,
$\mathcal{A}^{\lambda}_{1}$, namely,

$\displaystyle\mathcal{A}^{\lambda}_{1}\Phi\equiv\kappa\Delta_{\textbf{x}}\Phi+(2\lambda\textbf{e}+\textbf{v})\cdot\nabla_{\textbf{x}}\Phi+\big{(}\kappa\lambda^{2}+\lambda\textbf{v}\cdot\textbf{e}+\tau^{-1}f^{\prime}(0)\big{)}\Phi=\mu(\lambda)\Phi.$
(2)

In Eq.(2), $\Phi\in L^{2}(\mathbb{T}^{d})$, $\mathbb{T}=\mathbb{R}/\mathbb{Z}$
is the one-dimensional torus, and v is period $1$ in all direction
$x_{i},1\leq i\leq d$. Accurate estimation of $c^{*}(\textbf{e})$ boils down
to computing the principal eigenvalue of the operator
$\mathcal{A}^{\lambda}_{1}$ in (2). Adaptive finite element methods (FEM) were
successfully applied to solve (2) in [40, 39]. If the velocity field
$\textbf{v}=\textbf{v}(t,\textbf{x})$ in the KPP equation (1) is periodic in
time $t$, then the variational formula
$c^{*}(\textbf{e})=\inf_{\lambda>0}\mu(\lambda)/\lambda$ still holds [30],
where $\mu(\lambda)$ is the principal eigenvalue [17] of the time-periodic
parabolic operator, $\mathcal{A}^{\lambda}_{2}$, namely,

$\displaystyle\mathcal{A}^{\lambda}_{2}\Phi\equiv\kappa\Delta_{\textbf{x}}\Phi+(2\lambda\textbf{e}+\textbf{v})\cdot\nabla_{\textbf{x}}\Phi+\big{(}\kappa\lambda^{2}+\lambda\textbf{v}\cdot\textbf{e}+\tau^{-1}f^{\prime}(0)\big{)}\Phi-\Phi_{t}=\mu(\lambda)\Phi,$
(3)

on the space-time domain $\mathbb{T}^{d}\times[0,T]$ ($T$ is the period of v
in $t$), subject to the same boundary condition in x as (1) and periodic in
$t$. An edge-averaged FEM with algebraic multigrid acceleration was developed
in [47] to study KPP front speeds in two-dimensional time-periodic cellular
flows with chaotic streamlines. Adaptive FEM methods provide an efficient way
to investigate the KPP front speeds in time-periodic cellular and chaotic
flows. However, when the magnitude of velocity field is large and/or the
dimension of spatial variables is big (e.g. $d=3$), it is extremely expensive
to compute KPP front speeds by using the FEM.

Recently, we have made significant progress in developing Lagrangian particle
methods for computing effective diffusivities in chaotic and random flows [43,
42, 23]. This motivates us to develop interacting particle methods to compute
KPP front propagation in time-periodic cellular and chaotic flows in this
paper, especially in three-dimensional flows.

In this paper, we first apply operator splitting methods to approximate the
solution operator of the linear advection-diffusion operator (see Eq.(4)),
which is a non-autonomous evolution equation and corresponding to the
linearization of the KPP equation. Then, we develop Lagrangian numerical
schemes to compute the KPP front propagation through the Feynman-Kac formula,
which establishes a link between parabolic PDEs and SDEs. Direct approximation
of the Feynman-Kac formula is unstable, since the main contribution to the
expectation comes from sample paths that visit maximal points of the
potential; see Eq.(7). Alternatively, we study a normalized version, i.e., the
Feynman-Kac semigroup. Specifically, the principal eigenvalue of
$\mathcal{A}_{1}^{\lambda}$ and $\mathcal{A}_{2}^{\lambda}$ can be obtained by
studying the convergence of Feynman-Kac semigroups for SDEs associated with
operators $\mathcal{A}_{1}^{\lambda}$ and $\mathcal{A}_{2}^{\lambda}$ [8, 12].
We approximate the evolution of probability measures by an interacting
particle system and use the resampling technique to reduce the variance.
Moreover, we estimate the approximation of semigroups associated with the
solution operators of non-autonomous evolution equations and obtain
convergence analysis for our method in computing the KPP front speeds.

We point out that using Feynman-Kac semigroups to estimate the principal
eigenvalue of differential operators has a long history. It was developed in
large deviation theory, where Feynman-Kac semigroups were used to calculate
cumulant generating functions [9]. They were also used in important practical
applications, such as the diffusion Monte Carlo (DMC) method [14]. In the case
when the flow is autonomous, [12] gave a rigorous proof of convergence and
error analysis using a backward error analysis approach. When the flow becomes
non-autonomous, their method cannot be directly applied. There are several
novelties in our paper. Firstly, we analyze the solution operator by an
operator splitting method and estimate the error in the $L_{2}$ operator norm.
Secondly, we prove the convergence of estimating principal eigenvalues by the
Feynman-Kac semigroups for non-autonomous periodic systems. Furthermore, we
develop the $N$-interacting particle system ($N$-IPS) method to calculate the
principal eigenvalue, where several important 3D chaotic flows are
investigated. Notice that when the magnitude of the velocity field is large
and/or the dimension of spatial variables is three, it is extremely expensive
to calculate the principal eigenvalue using the FEM and spectral method.

Finally, we carry out numerical experiments to demonstrate the accuracy and
efficiency of the proposed method in computing KPP front speeds for time-
periodic cellular and chaotic flows. Most importantly, we investigate the
dependence of KPP front speeds on the chaos (disorder) and flow intensities.
Let $A$ denote the magnitude of the velocity field. For space-time-periodic
shear flow, the speed $c^{*}(A)$ obeys a quadratic enhancement law:
$c^{*}(A)=c_{0}(1+\alpha A^{2})+O(A^{3})$, $A\ll 1$, where $c_{0}$ is the KPP
front speed in homogeneous media ($A=0$) and $\alpha>0$ depends only on flow v
[31]. The study for complicated flows, e.g. 3D flows remains largely open. At
large $A$, the solution of the principal eigenvalue problem (2) develops
internal layers, which brings difficulties for the FEM and spectral method. We
will study this issue in Section 4.3. Numerical results show that our
interacting particle method is still very efficient when the magnitude of
velocity field $A$ is large and computational cost linearly depends on the
dimension $d$ of spatial variables in the KPP equation (1). Thus, we are able
to compute the KPP front speeds for time-dependent cellular and chaotic flows
of physical interests, including the ABC flows and the Kolmogorov flows in
three-dimensional space. We numerically verify that the relationship between
the KPP front speed $c^{*}(A)$ and the effective diffusivity $D^{E}(A)$, i.e.
$c^{*}(A)=O(\sqrt{D^{E}(A)})$, is true in 2D steady cellular flows. To the
best of our knowledge, our work appears to be the first one in the literature
to develop numerical methods to compute KPP front speeds in 3D time-dependent
flows.

The rest of the paper is organized as follows. In Section 2, we propose
Lagrangian interactive particle methods in computing KPP front speeds in time-
periodic cellular and chaotic flows. In Section 3, we estimate the
approximation of semigroups associated with the solution operators of non-
autonomous evolution equations and obtain convergence analysis for our method.
In Section 4, we present numerical results to demonstrate the accuracy and
efficiency of our method. In addition, we investigate the dependence of KPP
front speeds on the chaos (disorder) and flow intensities, especially in 3D
time-dependent chaotic flows. Concluding remarks are made in Section 5.
Finally, we collect several fundamental results for abstract linear evolution
equations by semigroup theory in the Appendix.

##  2 Efficient Lagrangian methods in computing KPP front speeds

###  2.1 Computing principal eigenvalue via the Feynman-Kac formula

In this section, we develop Lagrangian interacting particle methods to compute
KPP front propagation via the Feynman-Kac formula. We consider the linearized
corrector equation of the KPP equation (1), where the velocity field
$\textbf{v}(t,\textbf{x})$ is space-time periodic, mean zero, and divergence-
free. To compute the KPP front speed $c^{*}(\textbf{e})$ along direction e,
let $w$ solve a linearized equation parameterized by $\lambda>0$:

$w_{t}=\mathcal{A}w:=\kappa\Delta
w+(2\lambda\textbf{e}+\textbf{v})\cdot\nabla_{\textbf{x}}w+\big{(}\kappa\lambda^{2}+\lambda\textbf{v}\cdot\textbf{e}+\tau^{-1}f^{\prime}(0)\big{)}w,$
(4)

with initial condition $w(\textbf{x},0)=1$. Then, the principal eigenvalue
$\mu(\lambda)$ is given by

$\mu(\lambda)=\lim_{t\to\infty}\frac{1}{t}\ln\int_{\mathbb{T}^{d}}w(t,\textbf{x})d\textbf{x}.$
(5)

The number $\mu(\lambda)$ is also the principal Lyapunov exponent of the
parabolic equation (4), which is convex and superlinear in large $\lambda$
[30, 47]. Finally, we compute the KPP front speed using the variational
formula $c^{*}(\textbf{e})=\inf_{\lambda>0}\mu(\lambda)/\lambda$.

To design Lagrangian particle methods, we decompose the operator $\mathcal{A}$
in (4) into $\mathcal{A}=\mathcal{L}+\mathcal{C}$, where
$\mathcal{L}:=\kappa\Delta+(2\lambda\textbf{e}+\textbf{v})\cdot\nabla_{\textbf{x}}$
and
$\mathcal{C}:=c(t,\textbf{x})=\big{(}\kappa\lambda^{2}+\lambda\textbf{v}\cdot\textbf{e}+\tau^{-1}f^{\prime}(0)\big{)}$.
To approximate the operator $\mathcal{L}$, we define a SDE system as follows

${d\textbf{X}^{t,\textbf{x}}=\textbf{b}(t,\textbf{X}^{t,\textbf{x}})dt+\sqrt{2\kappa}d\textbf{w}(t),\quad\textbf{X}^{0,\textbf{x}}=\textbf{x}},$
(6)

where the drift term $\textbf{b}=2\lambda\textbf{e}+\textbf{v}$ is determined
by the advection field in the operator $\mathcal{L}$ and $\textbf{w}(t)$ is a
$d$-dimensional Brownian motion. The principal eigenvalue $\mu(\lambda)$ of
(4) can be represented via the Feynman-Kac formula as

$\mu(\lambda)=\lim_{t\to\infty}\frac{1}{t}\ln\mathbb{E}\Big{(}\exp\big{(}\int_{0}^{t}c(s,\textbf{X}^{s,\textbf{x}})ds\big{)}\Big{)},$
(7)

where the expectation $\mathbb{E}[\cdot]$ is over randomness induced by the
Brownian motion $\textbf{w}(t)$.

If we apply the formula (5) to compute the principal eigenvalue
$\mu(\lambda)$, we need to solve a parabolic-type PDE (4) using numerical
methods, such as FEM and spectral method. When the magnitude of velocity field
is large and/or the dimension of spatial variables $d$ is big (say $d=3$), the
FEM and spectral method become extremely expensive. The Feynman-Kac formula
(7) provides an alternative strategy to design Lagrangian methods to compute
the principal eigenvalue $\mu(\lambda)$, and thus allows us to compute the KPP
front speeds. As we will demonstrate in Section 4, the proposed Lagrangian
methods are efficient for computing KPP front speeds in 3D time-dependent
chaotic flows.

######  Remark 2.1.

When the velocity field in the KPP equation (1) is time-independent,
construction of the Lagrangian method for computing KPP front speeds is
straightforward. We simply replace the drift term b in (6) and the potential
$c$ in (7) by their time-independent counterparts.

###  2.2 Feynman-Kac semigroups

Directly using the Feynman-Kac formula (7) and Monte Carlo method to compute
the principal eigenvalue $\mu(\lambda)$ is unstable as the main contribution
to
$\mathbb{E}\Big{(}\exp\big{(}\int_{0}^{t}c(s,\textbf{X}^{s,\textbf{x}})ds\big{)}\Big{)}$
comes from sample paths that visit maximal or minimal points of the potential
function $c$, which leads to inaccurate or even divergent results.

Accurate principal eigenvalue $\mu(\lambda)$ can be obtained by studying the
convergence of the Feynman-Kac semigroup associated with the SDE system (6)
and the potential $c$. Specifically, let $\mathcal{P}(\mathbb{T}^{d})$ denote
the set of probability measures over $\mathbb{T}^{d}$ and
$S=\mathcal{C}^{\infty}(\mathbb{T}^{d})$. We define the evolution operator
associated with the process $(\textbf{X}^{t,\textbf{x}})_{t\geq 0}$ in (6) as

$(\nu)(P_{t}\phi)=\mathbb{E}_{\nu}[\phi(\textbf{X}^{t,\textbf{x}})],\quad\forall\nu\in\mathcal{P}(\mathbb{T}^{d}),~{}\phi\in
S.$ (8)

Similarly, we define its weighted counterpart as

$(\nu)(P^{c}_{t}\phi)=\mathbb{E}_{\nu}\big{[}\phi(\textbf{X}^{t,\textbf{x}})\exp\big{(}\int_{0}^{t}c(s,\textbf{X}^{s,\textbf{x}})\big{)}ds\big{]},\quad\forall\nu\in\mathcal{P}(\mathbb{T}^{d}),~{}\phi\in
S.$ (9)

In other words, infinitesimal generators of $P_{t}$ and $P^{c}_{t}$ are
$\mathcal{L}$ and $\mathcal{A}=\mathcal{L}+\mathcal{C}$, respectively.

Equipped with the definitions of the evolution operators $P_{t}$ and
$P^{c}_{t}$, we can define the Feynman-Kac semigroup $\Phi^{c}_{t}$ as follows

$\Phi^{c}_{t}(\nu)(\phi):=\frac{(\nu)(P^{c}_{t}\phi)}{(\nu)(P^{c}_{t}1)}=\frac{\mathbb{E}_{\nu}\Big{(}\phi(\textbf{X}^{t,\textbf{x}})\exp\big{(}\int_{0}^{t}c(s,\textbf{X}^{s,\textbf{x}})\big{)}ds\Big{)}}{\mathbb{E}_{\nu}\Big{(}\exp\big{(}\int_{0}^{t}c(s,\textbf{X}^{s,\textbf{x}})\big{)}ds\Big{)}}.$
(10)

One can easily verify that for all $\nu\in\mathcal{P}(\mathbb{T}^{d})$ and
$t_{1}$, $t_{2}\in\mathbb{R}_{+}$,
$\Phi^{c}_{t_{1}}(\Phi^{c}_{t_{2}}(\nu))=\Phi^{c}_{t_{1}+t_{2}}(\nu)$. Thus,
the family of maps $\\{\Phi^{c}_{t}\\}_{t\geq 0}$ is a measure-valued
semigroup.

Notice that we use $T$ to denote the period of the velocity in time.
Therefore, we consider the Feynman-Kac semigroup for $t=nT,n\in\mathbb{N}$.
Namely, we consider $\Phi^{c}_{nT}=(\Phi^{c}_{T})^{n}$. One can easily verify
the Feynman-Kac semigroup $\Phi^{c}_{nT}$ satisfies the following property,
where the proof is a direct conclusion of Theorem 3.9 and Theorem 3.10.

######  Proposition 2.2.

For any $\nu\in\mathcal{P}(\mathbb{T}^{d})$ and $\phi\in S$, there exists
$C>0$ such that

$\Big{|}\Phi^{c}_{nT}(\nu)(\phi)-\int_{\Omega}\phi d\nu_{c}\Big{|}\leq
C||\phi||\exp(-\delta_{c}nT),$ (11)

where
$\delta_{c}=\inf\\{\mu(\lambda)-\Re(z):z\in\sigma(\mathcal{A})\setminus\\{\mu(\lambda)\\}\\}>0$
is the spectral gap of the operator $\mathcal{A}$.

The exponential-decay property stated above ensures us to obtain an invariant
measure $\nu_{c}$ for $\Phi^{c}_{T}$ from any initial measure $\nu$. From the
definition of $\nu_{c}$, we know that $\Phi_{T}^{c}(\nu_{c})=\nu_{c}$, which
means that for any $\phi\in S$

$\int_{\mathbb{T}^{d}}\phi
d\nu_{c}=(\int_{\mathbb{T}^{d}}P_{T}^{c}1d\nu_{c})^{-1}\int_{\mathbb{T}^{d}}P_{T}^{c}\phi
d\nu_{c}.$ (12)

Therefore, we can find that the principal eigenvalue of $P_{T}^{c}$ is just
$\int_{\mathbb{T}^{d}}P_{T}^{c}1d\nu_{c}$, which provides a feasible way to
compute the principal eigenvalue.

###  2.3 Numerical discretization and resampling techniques

Let $M$ be the number of time discretization interval for each period and
$\Delta t=T/M$. We use the Euler-Maruyama scheme to discretize the SDE (6) and
obtain

${\bm{X}}_{i+1}=\bm{X}_{i}+\textbf{b}(t_{i},\bm{X}_{i})\Delta
t+\sqrt{2\kappa\Delta t}\bm{\omega}_{i},$ (13)

where $t_{i}=i\Delta t$ and $\bm{\omega}_{i}$’s are i.i.d. $d$-dimensional
standard Gaussian random variables. The numerical scheme (13) defines an
evolution operator $P^{\Delta t}_{i}$ (also known as transition operator) as
follows

$P^{\Delta
t}_{i}\phi(\bm{x})=\mathbb{E}(\phi({\bm{X}}_{i+1})|\bm{X}_{i}=\bm{x}).$ (14)

The evolution operator $P^{\Delta t}_{i}$ describes how the values of a given
function evolve in $L_{2}$ sense over one time step $\Delta t$. One can easily
verify that

$\big{|}\big{|}P^{\Delta t}_{i}-e^{\Delta
t\mathcal{L}(t_{i})}\big{|}\big{|}_{L^{2}}\leq C(\Delta t)^{2},$ (15)

where $C$ is a positive constant [26]. Specially, when $\textbf{b}=0$,
$P^{\Delta t}_{t_{i}}=e^{\Delta t\mathcal{L}(t_{i})}$ for all $i$. Therefore,
solving the SDE system (6) by the numerical scheme (13) provides a good
approximation to the evolution operator $e^{\Delta t\mathcal{L}(t_{i})}$,
which plays an important role in the error estiamte of our Lagrangian methods
in Section 3.

In addition, we can define the approximation operator for $P^{c}_{t}$ in (9).
For instance, if we choose the left-point rectangular integration, we obtain

$(\nu)(P_{i}^{\Delta t}e^{\Delta
t\mathcal{C}(t_{i})}\phi)=\mathbb{E}\Big{[}\phi(\bm{X}_{i+1})\exp\big{(}c(t_{i},\bm{X}_{i+1})\Delta
t\big{)}\big{|}\bm{X}_{i}\sim\nu\Big{]},\quad i=1,2,....,M.$ (16)

The time discretization for Feynman-Kac semigroup (10) reads:

$\Phi_{i}^{\mathcal{C},\Delta t}(\nu)(\phi)=\frac{(\nu)(P_{i}^{\Delta
t}e^{\Delta t\mathcal{C}(t_{i})}\phi)}{(\nu)(P_{i}^{\Delta t}e^{\Delta
t\mathcal{C}(t_{i})}1)},\quad i=0,1,2,....,M-1.$ (17)

It is difficult to obtain a closed-form solution to the evolution of
probability measure in (17). Therefore, we approximate the evolution of
probability measure in (17) by an $N$-interacting particle system ($N$-IPS)
[28]. Let us introduce the notation $\mathcal{K}^{\Delta
t}=\mathcal{K}^{\Delta t,M-1}\mathcal{K}^{\Delta
t,M-2}\cdots\mathcal{K}^{\Delta t,0}$, where $\mathcal{K}^{\Delta
t,i}=P_{i}^{\Delta t}e^{\Delta t\mathcal{C}(t_{i})}$, $\Delta t=T/M$, and $T$
is time period. We denote

$\displaystyle\Phi^{\mathcal{K}^{\Delta
t,i}}(\nu)(\phi)=\frac{(\nu)(\mathcal{K}^{\Delta
t,i}\phi)}{(\nu)(\mathcal{K}^{\Delta t,i}1)}$ (18)

the Feynman-Kac semigroup associated with the operator $\mathcal{K}^{\Delta
t,i}$. Then, according to Lemma 3.8, it satisfies

$\Phi^{\mathcal{K}^{\Delta t}}=\prod_{i=0}^{M-1}\Phi^{\mathcal{K}^{\Delta
t,i}}=\Phi^{\mathcal{K}^{\Delta t,0}}\Phi^{\mathcal{K}^{\Delta
t,1}}\cdots\Phi^{\mathcal{K}^{\Delta t,M-1}}.$ (19)

Suppose the Markov process $({\Theta},(\mathcal{F}_{n})_{n\geq
0},(\xi^{n})_{n\geq 0},\mathbb{P})$ is defined in the product space
$(\mathbb{T}^{d})^{N}$. For any initial probability measure $\pi_{0}=\nu$, we
approximate it by a $N$-particle system as

$P(\bm{\xi}^{0}\in d\bm{z})=\prod_{p=1}^{N}\pi_{0}(dz^{p}).$ (20)

Then, we evolve the $N$-particle system according to

$P(\bm{\xi}^{n}\in
d\bm{z}|\bm{\xi}^{n-1}=\bm{x})=\prod_{p=1}^{N}\Phi^{\mathcal{K}^{\Delta
t}}(\frac{1}{N}\sum_{i=1}^{N}\delta_{x^{i}})(dz^{p})=\prod_{p=1}^{N}(\prod_{i=0}^{M-1}\Phi^{\mathcal{K}^{\Delta
t,i}})(\frac{1}{N}\sum_{i=1}^{N}\delta_{x^{i}})(dz^{p}),$ (21)

where $\bm{x}=(x^{1},...,x^{N})^{T}$ and $n$ denotes the iteration number in
the evolution of probability measure by the Feynman-Kac semigroup (17).

From Eq.(20), we can compute the evolution from $\bm{\xi}^{n-1}$ to
$\bm{\xi}^{n}$. It will be divided into $M$ small steps. Let us denote
$\bm{\xi}^{n}_{0}=\bm{\xi}^{n}$ for all $n$. Within each iteration stage
$n-1$, we evolve the particles from $t=0$ to $t=T$ by the evolution operator
$P^{\Delta t}_{i}$ and resample these particles according to weights
determined by the potential function. Specifically, at $t_{i}=i\Delta t$,
$i=0,...,M-1$, we evolve the particles $\bm{\xi}_{i}^{n-1}$ by the numerical
scheme (13) and get $\widetilde{\bm{\xi}}_{i+1}^{n-1}$. Each particle is
updated by

$\widetilde{\xi}_{i+1}^{p,n-1}=\xi_{i}^{p,n-1}+\textbf{b}(t_{i},\xi_{i}^{p,n-1})\Delta
t+\sqrt{2\kappa\Delta t}\bm{\omega}_{i}^{p,n-1},\quad p=1,2,...,N,$ (22)

where $\bm{\omega}_{i}^{p,n-1}$ are i.i.d. $d$-dimensional standard Gaussian
random variables.

Then, we resample the particles in $\widetilde{\bm{\xi}}_{i+1}^{n-1}$
according to the multinomial distribution with the weights

${w}_{i}^{p,n-1}=\frac{\exp{\big{(}c(t_{i},\widetilde{\xi}_{i}^{p,n-1})\Delta
t\big{)}}}{\sum_{p=1}^{N}\exp{\big{(}c(t_{i},\widetilde{\xi}_{i}^{p,n-1})\Delta
t\big{)}}},\quad p=1,...,N,$ (23)

and obtain $\bm{\xi}_{i+1}^{n-1}$.

The evolution of $N$-IPS from $(n-1)T$ to $nT$ can be represented as follows

$\displaystyle\bm{\xi}_{0}^{n-1}=(\xi_{0}^{1,n-1},\cdots,\xi_{0}^{N,n-1})\longrightarrow\bm{\xi}_{1}^{n-1}=(\xi_{1}^{1,n-1},\cdots,\xi_{1}^{N,n-1})\longrightarrow$
$\displaystyle\cdots\longrightarrow\bm{\xi}_{M}^{n-1}=(\xi_{M}^{1,n-1},\cdots,\xi_{M}^{N,n-1})=\bm{\xi}_{0}^{n}=(\xi_{0}^{1,n},\cdots,\xi_{0}^{N,n}).$
(24)

We know that the empirical distribution of the particles $\bm{\xi}_{0}^{n}$
will weakly converge to the distribution $\Phi_{n}^{\mathcal{K}^{\Delta
t}}(\pi_{0})$ as $N\to\infty$.

After obtaining the empirical distribution of the particles
$\bm{\xi}_{0}^{n}$, we can compute the principal eigenvalues. At the iteration
stage $n$, we first define the change of the mass as follows

$e_{i,n}^{N}=N^{-1}\sum_{p=1}^{N}\exp(c(t_{i},\widetilde{\xi}_{i}^{p,n-1})\Delta
t).$ (25)

Then, we compute the approximation of the principal eigenvalue by

$\mu^{n}_{\Delta t}(\lambda)=(M\Delta
t)^{-1}\sum_{i=1}^{M}\log\Big{(}N^{-1}\sum_{p=1}^{N}\exp(c(t_{i},\widetilde{\xi}_{i}^{p,n-1})\Delta
t)\Big{)}.$ (26)

Finally, we give the complete algorithm in Algorithm 1. The performance of our
method will be demonstrated in Section 4.

Algorithm 1 Algorithm for computing the principal eigenvalues of parabolic
equations

1:velocity field $\textbf{v}(\textbf{x},t)$, potential $c(\textbf{x},t)$,
number of $N$-IPS system (i.e., $N$), initial probability measure $\nu_{0}$,
iteration number $n$, time period $T$, and time step $\Delta t=T/M$.

2:Generate $N$ i.i.d. $\nu_{0}$-distributed random variables on $[0,1]^{d}$:
$\bm{\xi}_{0}^{1}=(\xi_{0}^{1,1},\cdots,\xi_{0}^{N,1})$.

3:for $k=1:n$ do

4: for $i=0:M-1$ do

5: Generate standard normal variables $\omega_{i}^{k}$ and compute
$\widetilde{\bm{\xi}}_{i+1}^{k-1}$ according to ${\bm{\xi}}_{i}^{k-1}$ by
(22).

6: Compute the pointwise value $\bm{S}=(e^{C^{1}},\cdots,e^{C^{N}})$, where
$C^{i}=c(i\Delta t,\widetilde{\bm{\xi}}_{i}^{k-1})\Delta t$.

7: Compute weights $\bm{w}=(w^{1},\cdots,w^{N})=\bm{S}/\text{SUM}(\bm{S})$ and
${\bf E}_{k,i}=\frac{1}{\Delta t}\text{log}(\text{mean}(\bm{S}))$.

8: Resample $\widetilde{\bm{\xi}}_{i+1}^{k-1}$ according to multinomial
distribution with weight $\bm{w}$, and get $\bm{\xi}_{i+1}^{k-1}$.

9: end for

10: Compute ${\bf E}_{k}=M^{-1}\sum_{i=0}^{M-1}({\bf E}_{k,i})$ and define
$\bm{\xi}_{M}^{k-1}=\bm{\xi}_{0}^{k}$.

11:end for

12:The approximate invariant distribution $\Phi_{n}^{\mathcal{K}^{\Delta
t}}(\nu)$-distributed variable $\bm{\xi}_{n}^{0}$ and approximate eigenvalue
$\widetilde{\lambda}_{n,\Delta t}=n^{-1}\sum_{k=1}^{n}{\bf E}_{k}$;

######  Remark 2.1.

When the flow is time-independent, we can view it as a periodic flow with any
given period $T$. The principal eigenvalue calculated by (5) will be the same
with any $T>0$. Hence the numerical schemes and the convergence analysis
proposed in time-dependent flow can be applied by assigning $T=\Delta t$ and
$M=1$.

##  3 Convergence analysis of the Lagrangian particle method

In this section, we will prove the convergence of the Lagrangian particle
method in computing KPP front speed. We divide the analysis into two parts.
The first part studies the approximation of the evolution of parabolic
operators by using an operator splitting method. The second part studies the
error estimate of the Lagrangian particle method in computing the principal
eigenvalue of parabolic operators.

###  3.1 Approximation of the evolution of parabolic operators

We first rewrite the linearized corrector equation of the KPP equation (4)
into the following non-autonomous parabolic equation

$w_{t}=\kappa\Delta_{\textbf{x}}w+\textbf{b}(t,\textbf{x})\cdot\nabla_{\textbf{x}}w+c(t,\textbf{x})w,\quad\textbf{x}=(x_{1},...,x_{d})^{T}\in\mathbb{T}^{d}=[0,1]^{d},\quad
t\in[0,T],$ (27)

where the initial condition $w(0,\textbf{x})=w_{0}$,
$\textbf{b}(t,\textbf{x})=2\lambda\textbf{e}+\textbf{v}$,
$c(t,\textbf{x})=\kappa\lambda^{2}+\lambda\textbf{v}\cdot\textbf{e}+\tau^{-1}f^{\prime}(0)$,
and $T$ is final computational time. Since the velocity
$\textbf{v}=\textbf{v}(t,\textbf{x})$ is space-time periodic, so do
$\textbf{b}(t,\textbf{x})$ and $c(t,\textbf{x})$. We assume the period of
$\textbf{b}(t,\textbf{x})$ and $c(t,\textbf{x})$ is one in each dimension and
they are smooth functions. For notational simplicity, we define

$\mathcal{A}(t)=\mathcal{L}(t)+\mathcal{C}(t),$ (28)

where
$\mathcal{L}(t):=\kappa\Delta_{\textbf{x}}+\textbf{b}(t,\textbf{x})\cdot\nabla_{\textbf{x}}$
and $\mathcal{C}(t)=c(t,\textbf{x})$. The operator $\mathcal{A}(t)$ has a real
isolated principal eigenvalue $\mu(\lambda)$ [17]. We aim to obtain error
estimates of our Lagrangian method in approximating the principal eigenvalue
$\mu(\lambda)$. To this end, we study the approximation the solution operator
for the parabolic equation (27) by using an operator splitting method.

We define the solution operator $\mathcal{U}(t,s)$ corresponding to the
parabolic equation (27), which satisfies the following properties:

  1. 1.

$\mathcal{U}(s,s)=Id$, for any $s\geq 0$;

  2. 2.

$\mathcal{U}(t,r)\circ\mathcal{U}(r,s)=\mathcal{U}(t,s)$, for any $t\geq r\geq
s\geq 0$;

  3. 3.

$\frac{d}{dt}\mathcal{U}(t,s)w_{0}=\mathcal{A}(t)\mathcal{U}(t,s)w_{0}$, for
any $t\geq s\geq 0,w_{0}\in L^{2}([0,1]^{d})$.

The solution operator $\mathcal{U}(t,s)$ enables us to study the evolution of
parabolic operator in (27), e.g., the principal eigenvalue of
$\mathcal{U}(T,0)$ gives the principal eigenvalue of the parabolic operator
$\mathcal{A}(t)$. It has been proven that the principal eigenvalue of
$\mathcal{U}(T,0)$ exists and is real [17]. It is difficult to obtain a
closed-form for the solution operator $\mathcal{U}(T,0)$. Therefore, we
approximate the solution operator $\mathcal{U}(T,0)$ by using an operator
splitting method.

We set $t_{i}=i\Delta t$ with $\Delta t=\frac{T}{M}$ and consider the
following parabolic equation with freezing time coefficients.

$w_{t}=\kappa\Delta_{\textbf{x}}w+\textbf{b}(t_{i},\textbf{x})\cdot\nabla_{\textbf{x}}w+c(t_{i},\textbf{x})w,\quad
t_{i}<t\leq t_{i+1},\quad i\geq 0.$ (29)

The corresponding solution operator can be formally represented as

$w(t)=e^{(t-t_{i})(\mathcal{L}+\mathcal{C})(t_{i})}\prod_{k=0}^{i-1}e^{\Delta
t(\mathcal{L}+\mathcal{C})(t_{k})}w_{0},\quad t_{i}\leq t<t_{i+1}.$ (30)

Furthermore, we can apply the first-order Lie-Trotter operator splitting
method to approximate the solution operator defined in (30) and obtain

$w(t)=e^{(t-t_{i})\mathcal{L}(t_{i})}e^{(t-t_{i})\mathcal{C}(t_{i})}\prod_{k=0}^{i-1}e^{\Delta
t\mathcal{L}(t_{j})}e^{\Delta t\mathcal{C}(t_{j})}w_{0},\quad t_{i}\leq
t<t_{i+1}.$ (31)

We will prove the solution operator $\prod_{j=0}^{M-1}e^{\Delta
t\mathcal{L}(t_{j})}e^{\Delta t\mathcal{C}(t_{j})}$ obtained by the Lie-
Trotter operator splitting method converges to the solution operator
$\mathcal{U}(T,0)$ in certain operator norm as $\Delta t$ approaches zero. As
a consequence of this convergence result, we can further prove the convergence
of the principal eigenvalue associated with these two solution operators.

To make our paper self-contained, we collect several fundamental results for
abstract linear evolution equations by semigroup theory in A. We begin with
following lemma, which is as a special case of Theorem 1 in [41].

######  Lemma 3.3.

Let $t$ be fixed. If $\textbf{b}(t,{\bf x})$ and $c(t,{\bf x})$ are smooth and
bounded, then the operator $\mathcal{A}(t)$ defined in (28) is a strongly
elliptic operator on $\mathbb{T}^{d}$. Moreover, $\mathcal{A}(t)$ generates an
analytic semigroup $e^{\cdot\mathcal{A}(t)}$ in $L^{p}(\mathcal{D})$, for all
$1\leq p\leq\infty$.

We will prove that, in our non-autonomous parabolic equation setting, the
assumptions made in A are all satisfied, so we can obtain the error of the
operator splitting method in approximation the non-autonomous parabolic
operator.

We first prove that the operator $\mathcal{A}$ defined in (28) satisfies a
Hölder continuous condition.

######  Lemma 3.4.

Suppose $\textbf{b}(t,\textbf{x})$ and $c(t,\textbf{x})$ in the operator
$\mathcal{A}(t)$ are bounded, smooth and periodic in each component of x, and
uniformly Hölder continuous in $t$, i.e., for any $t,s\in\mathbb{R}^{+}$,

$\big{|}\big{|}\textbf{b}(t,\textbf{x})-\textbf{b}(s,\textbf{x})\big{|}\big{|}\leq
C_{1}|t-s|^{\beta},\quad\big{|}c(t,\textbf{x})-c(s,\textbf{x})\big{|}\leq
C_{1}|t-s|^{\beta},$ (32)

for some positive $C_{1}$ and $\beta$. Let
$v\in\mathcal{D}(\mathcal{A}(\cdot))=H^{2}(\mathbb{T}^{d})$ is periodic. Then,
for any $0<s\leq\tau$, there exists ${\gamma_{1}}>0$, such that

$\big{|}\big{|}\mathcal{A}(\tau)v-\mathcal{A}(s)v\big{|}\big{|}_{L^{2}}\leq
C_{2}(\tau-s)^{\beta}\big{|}\big{|}(\mathcal{A}(t)-{\gamma_{1}})v\big{|}\big{|}_{L^{2}}^{1/2}\big{|}\big{|}v\big{|}\big{|}_{L^{2}}^{1/2},$
(33)

for any $t\in\mathbb{R}^{+}$. Specifically, if $\textbf{b}(t,\textbf{x})=0$,
then

$\big{|}\big{|}\mathcal{A}(\tau)v-\mathcal{A}(s)v\big{|}\big{|}_{L^{2}}\leq
C_{3}(\tau-s)^{\beta}\big{|}\big{|}v\big{|}\big{|}_{L^{2}}.$ (34)

###### Proof.

We first consider the case when $\textbf{b}(t,\textbf{x})\neq 0$. By using the
uniformly Hölder continuous conditions for $b(t,\textbf{x})$ and
$c(t,\textbf{x})$, we have

$\displaystyle\big{|}\big{|}\mathcal{A}(\tau)v-\mathcal{A}(s)v\big{|}\big{|}_{L^{2}}=$
$\displaystyle\big{|}\big{|}(\textbf{b}(\tau,\textbf{x})-\textbf{b}(s,\textbf{x}))\cdot\nabla_{\textbf{x}}v+(c(\tau,\textbf{x})-c(s,\textbf{x}))v\big{|}\big{|}_{L^{2}}$
$\displaystyle\leq$ $\displaystyle
C_{1}(t-s)^{\beta}(||\nabla_{\textbf{x}}v||_{L^{2}}+||v||_{L^{2}}).$ (35)

For the operator $\mathcal{A}(t)$, we claim that there exists ${\gamma_{1}}>0$
such that,

$\big{|}\big{|}(\mathcal{A}(t)-{\gamma_{1}})v\big{|}\big{|}_{L^{2}}>C(\kappa,\textbf{b},c)\big{(}||\Delta_{\textbf{x}}v||_{L^{2}}+||v||_{L^{2}}\big{)},\quad\forall
v\in\mathcal{D}(\mathcal{A}(\cdot)),$ (36)

where the constant $C(\kappa,\textbf{b},c)$ depends on $\kappa$,
$\textbf{b}(t,\textbf{x})$ and $c(t,\textbf{x})$.

We prove the statement in (36) before move to the main results. Let
$c_{\gamma_{1}}=c-{\gamma_{1}}$ and assume
$\big{|}\big{|}\textbf{b}(t,\textbf{x})\big{|}\big{|}\leq M_{1}$,
$|c(t,\textbf{x})|\leq M_{2}$,
$\big{|}\big{|}\nabla_{\textbf{x}}c(t,\textbf{x})\big{|}\big{|}\leq M_{3}$. We
know that

$\displaystyle\big{|}\big{|}(\mathcal{A}(t)-{\gamma_{1}})v\big{|}\big{|}_{L^{2}}$
$\displaystyle=\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+\textbf{b}(t,\textbf{x})\cdot\nabla_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}$
$\displaystyle\geq\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}-\big{|}\big{|}b(t,\textbf{x})\cdot\nabla_{\textbf{x}}v\big{|}\big{|}_{L^{2}}.$
(37)

For the term
$\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}$,
the periodic condition of $v$ implies that

$\displaystyle\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}^{2}$
$\displaystyle=$
$\displaystyle||\kappa\Delta_{\textbf{x}}v||_{L^{2}}^{2}+\big{|}\big{|}c_{\gamma_{1}}(t,\textbf{x})v\big{|}\big{|}_{L^{2}}^{2}-2\langle\kappa\nabla_{\textbf{x}}v,c_{\gamma_{1}}(t,\textbf{x})\nabla_{\textbf{x}}v\rangle_{L^{2}}-2\langle\kappa\nabla_{\textbf{x}}v,v\nabla_{\textbf{x}}c(t,\textbf{x})\rangle_{L^{2}}$
(38)

Notice that if we choose ${\gamma_{1}}=\frac{2M_{1}^{2}}{\kappa}+M_{2}$, then
we obtain

$\displaystyle-2\langle\kappa\nabla_{\textbf{x}}v,c_{\gamma_{1}}(t,\textbf{x})\nabla_{\textbf{x}}v\rangle_{L^{2}}\geq
4\kappa({\gamma_{1}}-M_{2})||\nabla_{\textbf{x}}v||_{L^{2}}\geq
4\big{|}\big{|}\textbf{b}(t,\textbf{x})\cdot\nabla_{\textbf{x}}v\big{|}\big{|}_{L^{2}}.$
(39)

In addition, we have

$\displaystyle
2\langle\kappa\nabla_{\textbf{x}}v,v\nabla_{\textbf{x}}c(t,\textbf{x})\rangle_{L^{2}}\leq
2\kappa M_{3}||\nabla_{\textbf{x}}v||_{L^{2}}||v||_{L^{2}}\leq 2\kappa
M_{3}C||\Delta_{\textbf{x}}v||_{L^{2}}^{\frac{1}{2}}||v||_{L^{2}}^{\frac{3}{2}}.$
(40)

Here, we use the fact that $||\nabla_{\textbf{x}}v||_{L^{2}}\leq
C||\Delta_{\textbf{x}}v||_{L^{2}}^{\frac{1}{2}}||v||_{L^{2}}^{\frac{1}{2}}$,
which is the moment inequality in interpolation theory; see Theorem 5.34 of
[11]. If we take ${\gamma_{1}}$ large enough such that
$4(\frac{{\gamma_{1}}-M_{2}}{3})^{\frac{3}{4}}\kappa^{\frac{1}{4}}\geq 2\kappa
M_{3}C$, we get that

$||\kappa\Delta_{\textbf{x}}v||_{L^{2}}^{2}+\big{|}\big{|}c_{\gamma_{1}}(t,\textbf{x})v\big{|}\big{|}_{L^{2}}^{2}\geq
2\kappa
M_{3}C||\Delta_{\textbf{x}}v||_{L^{2}}^{\frac{1}{2}}||v||_{L^{2}}^{\frac{3}{2}}\geq
2\langle\kappa\nabla_{\textbf{x}}v,v\nabla_{\textbf{x}}c(t,\textbf{x})\rangle_{L^{2}}.$
(41)

Substititing the estiamtes (39)-(41) into (38), we obtain

$\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}\geq
2\big{|}\big{|}\textbf{b}(t,\textbf{x})\cdot\nabla_{\textbf{x}}v\big{|}\big{|}_{L^{2}}.$
(42)

Thus, from (37) we get that

$\big{|}\big{|}(\mathcal{A}(t)-{\gamma_{1}})v\big{|}\big{|}_{L^{2}}\geq\frac{1}{2}\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}.$
(43)

Using the same argument, we can prove that for ${\gamma_{1}}$ large enough,

$\big{|}\big{|}(\kappa\Delta_{\textbf{x}}+c_{\gamma_{1}}(t,\textbf{x}))v\big{|}\big{|}_{L^{2}}\geq\hat{C}(||\Delta_{\textbf{x}}v||_{L^{2}}+||v||_{L^{2}}).$
(44)

Finally, using the moment inequality we prove the statement in (33).

The case when $\textbf{b}(t,\textbf{x})=0$ is simple since we have

$\big{|}\big{|}\mathcal{A}(\tau)v-\mathcal{A}(s)v\big{|}\big{|}_{L^{2}}=\big{|}\big{|}(c(\tau,\textbf{x})-c(s,\textbf{x}))v\big{|}\big{|}_{L^{2}}\leq
C_{3}(t-s)^{\beta}||v||_{L^{2}}.$ (45)

∎

We then verify the operators $\mathcal{L}(t)$ and $\mathcal{C}(t)$ defined in
(28) satisfy the assumption A.24. Given $\tau\geq 0$, we assume the bounded
conditions as follows

$||e^{\tau\mathcal{L}(t)}||_{L^{2}}\leq
1,\quad||e^{\tau\mathcal{C}(t)}||_{L^{2}}\leq
1,\quad||e^{\tau(\mathcal{L}(t)+\mathcal{C}(t))}||_{L^{2}}\leq 1.$ (46)

######  Lemma 3.5.

Suppose $\textbf{b}(t,\textbf{x})$ and $c(t,\textbf{x})$ in the operator
$\mathcal{A}$ satisfy the same assumption as that in Lemma 3.4. Then, there
exists ${\gamma_{2}}>0$ such that, for any periodic $v\in
L^{2}(\mathbb{T}^{d})$, commutator of $\mathcal{L}$ and $\mathcal{C}$ acting
on $v$ follows,

$\big{|}\big{|}[\mathcal{L}(t),\mathcal{C}(t)]v\big{|}\big{|}_{L^{2}}\leq
C_{1}\big{|}\big{|}(\mathcal{L}(t)-{\gamma_{2}})v\big{|}\big{|}^{\frac{1}{2}}_{L^{2}}||v||_{L^{2}}^{\frac{1}{2}},\quad\forall
t\geq 0.$ (47)

###### Proof.

We first observe that, for any $v$ periodic in $L^{2}(\mathbb{T}^{d})$,

$\displaystyle\big{|}\big{|}[\mathcal{L}(t),\mathcal{C}(t)]v\big{|}\big{|}_{L^{2}}=\big{|}\big{|}\mathcal{L}(t)(\mathcal{C}(t)v)-\mathcal{C}(t)(\mathcal{L}(t)v)\big{|}\big{|}_{L^{2}},$
$\displaystyle=$
$\displaystyle\big{|}\big{|}\big{(}\kappa\Delta_{\textbf{x}}c(t,\textbf{x})+\textbf{b}(t,\textbf{x})\cdot\nabla_{\textbf{x}}c(t,\textbf{x})\big{)}v+2\kappa\nabla_{\textbf{x}}c(t,\textbf{x})\cdot\nabla
v\big{|}\big{|}_{L^{2}},$ $\displaystyle\leq$ $\displaystyle(\kappa
M_{4}+M_{1}M_{3})||v||_{L^{2}}+2\kappa M_{3}||\nabla_{\textbf{x}}v||_{L^{2}},$
(48)

where $||b(t,\textbf{x})||\leq M_{1}$,
$||\nabla_{\textbf{x}}c(t,\textbf{x})||\leq M_{3}$, and
$|\Delta_{\textbf{x}}c(t,\textbf{x})|\leq M_{4}$.

Following the same procedure as in the proof of Lemma 3.4, we have

$\displaystyle\big{|}\big{|}(\mathcal{L}(t)-{\gamma_{1}})v\big{|}\big{|}_{L^{2}}\geq
C(\kappa,\textbf{b})(||\Delta_{\textbf{x}}v||_{L^{2}}+||v||_{L^{2}}).$ (49)

Using the fact that $||\nabla_{\textbf{x}}v||_{L^{2}}\leq
C||\Delta_{\textbf{x}}v||_{L^{2}}^{\frac{1}{2}}||v||_{L^{2}}^{\frac{1}{2}}$,
we finally prove the assertion in (47). ∎

######  Remark 3.1.

If the bounded conditions (46) for $\mathcal{L}(t)$ and $\mathcal{C}(t)$ do
not hold, we can shift the operators by a constant so that the shifted
operators satisfy the bounded condition. Shift the operator by a constant will
not affect the commutator in (47).

Now we are in the position to present the main result in approximating the
solution operator $\mathcal{U}(t,s)$ for the parabolic equation (27).

######  Theorem 3.6.

The solution operator (30) has the following error in approximating the
solution operator $\mathcal{U}(T,0)$ in $L^{2}$ operator norm

$\displaystyle||\mathcal{U}(T,0)-\prod_{k=0}^{M-1}e^{\Delta
t\mathcal{A}(k\Delta t)}||_{L^{2}(\mathbb{T}^{d})}\leq C_{1}(T)(\Delta
t)^{\beta-\frac{1}{2}},$ (50)

where $T>0$, $M$ is an integer, and $\Delta t=\frac{T}{M}$. In addition, the
Lie-Trotter operator splitting method has the following error in approximaing
the solution operator (30)

$\displaystyle\Big{|}\Big{|}\prod_{k=0}^{M-1}e^{\Delta t\mathcal{A}(k\Delta
t)}-\prod_{k=0}^{M-1}e^{\Delta t\mathcal{L}(k\Delta t)}e^{\Delta
t\mathcal{C}(k\Delta t)}\Big{|}\Big{|}_{L^{2}(\mathbb{T}^{d})}\leq
C_{2}(T)(\Delta t)^{\frac{1}{2}}.$ (51)

###### Proof.

We take $\gamma=\max(\gamma_{1},\gamma_{2})$, where $\gamma_{1}$ and
$\gamma_{2}$ are defined in Lemma 3.4 and Lemma 3.5 respectively. Let
$\mathcal{U}_{\gamma}(t,s)=e^{-{\gamma}(t-s)}\mathcal{U}(t,s)$ be the solution
operator that corresponds to the parabolic equation (27) with
$\mathcal{A}_{\gamma}(t)=\mathcal{A}(t)-{\gamma}$ and
$\mathcal{L}_{\gamma}(t)=\mathcal{L}(t)-{\gamma}$. Then, we have

$\displaystyle\mathcal{U}(T,0)-\prod_{k=0}^{M-1}e^{\Delta t\mathcal{A}(k\Delta
t)}=e^{{\gamma}T}(U_{\gamma}(T,0)-\prod_{k=0}^{M-1}e^{\Delta
tA_{\gamma}(k\Delta t)})$ (52)

The statement in (50) is proved accoeding to Theorem A.22.

For the Lie-Trotter operator splitting method, we know that

$\displaystyle\prod_{k=0}^{M-1}e^{\Delta t\mathcal{A}(k\Delta
t)}-\prod_{k=0}^{M-1}e^{\Delta t\mathcal{L}(k\Delta t)}e^{h\mathcal{C}(k\Delta
t)}=e^{{\gamma}T}(\prod_{k=0}^{M-1}e^{\Delta t\mathcal{A}_{\gamma}(k\Delta
t)}-\prod_{k=0}^{M-1}e^{\Delta t\mathcal{L}_{\gamma}(k\Delta t)}e^{\Delta
t\mathcal{C}(k\Delta t)})$ (53)

Now according to the Lemma 3.2 and 3.3, $\mathcal{A}_{\gamma}(k\Delta
t)=\mathcal{L}_{\gamma}(k\Delta t)+\mathcal{C}(k\Delta t)$ and
$\mathcal{L}_{\gamma}$ and $\mathcal{C}$ satisfy the assumptions A.23 and
A.24. Thus, applying Theorem A.25 and Theorem A.26, we can prove the estimate
(51). ∎

The convergence of $\mathcal{K}^{\Delta t}$ in the operator norm
$\mathcal{L}(L^{2},H^{1})$ has been proved in [2]. In Theorem 3.6 we obtain
the convergence of $\mathcal{K}^{\Delta t}$ in the operator norm
$\mathcal{L}(L^{2})$. Finally, we can obtain the error estimate for the
principal eigenvalue.

######  Theorem 3.7.

Let $e^{\mu(\lambda)T}$ and $e^{\mu_{\Delta t}(\lambda)T}$ denote the
principal eigenvalue of the solution operator $\mathcal{U}(T,0)$ and the
approximated solution operator $\prod_{k=0}^{M-1}e^{\Delta
t\mathcal{L}(k\Delta t)}e^{\Delta t\mathcal{C}(k\Delta t)}$, respectively.
Then, we have the error estimate as follows:

$\displaystyle\big{|}e^{\mu(\lambda)T}-e^{\mu_{\Delta t}(\lambda)T}\big{|}\leq
C_{1}(T)(\Delta t)^{\beta-\frac{1}{2}}+C_{2}(T)(\Delta t)^{\frac{1}{2}}.$ (54)

Hence, $|\mu(\lambda)-\mu_{\Delta t}(\lambda)|={O}((\Delta
t)^{\min(\beta-\frac{1}{2},\frac{1}{2})})$.

###### Proof.

According to the standard spectral theorem [21], the principal eigenvalue
$e^{\mu(\lambda)}$ for the solution operator $\mathcal{U}(T,0)$ and the
principal eigenvalue $e^{\mu_{\Delta t}(\lambda)}$ for the approximated
solution operator $\prod_{k=0}^{M-1}e^{\Delta t\mathcal{L}(k\Delta
t)}e^{\Delta t\mathcal{C}(k\Delta t)}$ satisfy

$\displaystyle\big{|}e^{\mu(\lambda)T}-e^{\mu_{\Delta t}(\lambda)T}\big{|}\leq
C_{3}\big{|}\big{|}U(T,0)-\prod_{k=0}^{N-1}e^{\Delta t\mathcal{L}(k\Delta
t)}e^{\Delta t\mathcal{C}(k\Delta t)}\big{|}\big{|}_{L^{2}(\mathbb{T}^{d})}.$
(55)

By using triangle inequality for the right hand side of (55) and the estimated
results from Theorem 3.6, we can get the error estimate (54). ∎

If $b(t,\textbf{x})$ and $c(t,\textbf{x})$ in the operator $\mathcal{A}$ are
uniformly Lipschitz, then the error of the principal eigenvalue obtained by
the Lie-Trotter operator splitting method is $O((\Delta t)^{\frac{1}{2}})$.

###  3.2 Analysis of the Lagrangian particle method

We consider the Feynman-Kac semigroup $\Phi^{\mathcal{A}}$ associated with an
arbitary operaor $\mathcal{A}$. The action of the Feynman-Kac semigroup
$\Phi^{\mathcal{A}}$ on a probability measure $\nu$ is defined by

$\Phi^{\mathcal{A}}(\nu)(\phi)=\frac{(\nu)(\mathcal{A}\phi)}{(\nu)(\mathcal{A}1)},\quad\forall\phi\in
L^{2}(\mathbb{T}^{d}).$ (56)

Moreover, we denote $\Phi^{\mathcal{A}}_{n}=(\Phi^{\mathcal{A}})^{n}$. The
Feynman-Kac semigroup operation satisfies the following property.

######  Lemma 3.8.

For any operaors $\mathcal{A}$, $\mathcal{B}$ in
$\mathcal{L}(L^{2}(\mathbb{T}^{d}))$,
$\Phi^{\mathcal{A}\mathcal{B}}=\Phi^{\mathcal{B}}\Phi^{\mathcal{A}}.$

###### Proof.

Let $\nu$ be a probability measure and $\phi$ be a function in
$L^{2}(\mathbb{T}^{d})$. Then, we can easily verify that

$\displaystyle\Phi^{\mathcal{A}\mathcal{B}}(\nu)(\phi)$
$\displaystyle=\frac{(\nu)(\mathcal{A}\mathcal{B}\phi)}{(\nu)(\mathcal{A}\mathcal{B}1)}=\frac{(\nu)(\mathcal{A}\mathcal{B}\phi)}{(\nu)(\mathcal{A}1)}\frac{(\nu)(\mathcal{A}1)}{(\nu)(\mathcal{A}\mathcal{B}1)},$
$\displaystyle=\frac{\Phi^{\mathcal{A}}(\nu)(\mathcal{B}\phi)}{\Phi^{\mathcal{A}}(\nu)(\mathcal{B}1)}=\Phi^{\mathcal{B}}\Phi^{\mathcal{A}}(\nu)(\phi).$
(57)

∎

Recall the operator $\Phi_{n}^{\mathcal{K}^{\Delta t}}$ defined in (19) is a
compostion of the Feynman-Kac semigroup $\Phi^{\mathcal{K}^{\Delta t,i}}$
associated with the operator $\mathcal{K}^{\Delta t,i}$; see (18). In the
sequel, we will prove the operator $\Phi_{n}^{\mathcal{K}^{\Delta t}}$
satisfies the uniform minorization and boundedness condition with respect to
$\Delta t$, which guarantees the existence of an invariant measure.

######  Theorem 3.9.

There exists a probability measure $\eta$ so that the operator
$\mathcal{K}^{\Delta t}$ satisfies a uniform minorization and boundedness
condition as follows

$\epsilon\eta(\phi)\leq\mathcal{K}^{\Delta
t}(\phi)(\textbf{x})\leq\gamma\eta(\phi),\quad\forall\textbf{x}\in\mathbb{T}^{d},\forall\phi\in
L^{2}(\mathbb{T}^{d}),$ (58)

where $0<\epsilon<\gamma$ is independent with $\Delta t$. Moreover, when
$\Delta t\to 0$ the limit operator is just the exact solution operator
$\mathcal{U}(T,0)$, which also satisfies this uniform minorization and
boundedness condition.

###### Proof.

We first define an operator $P^{\Delta t}=\prod_{i=0}^{M-1}P_{t_{i}}^{\Delta
t}$, which corresponds to the case when $c(t,\textbf{x})=0$ in Eq.(27). Since
$c(t,\textbf{x})$ is bounded (i.e. $c_{1}\leq c(t,\textbf{x})\leq c_{2}$), one
can easily obtain the following estimate based on the Feynman-Kac formula

$P^{\Delta t}(\phi)e^{c_{1}T}\leq\mathcal{K}^{\Delta t}(\phi)\leq P^{\Delta
t}(\phi)e^{c_{2}T}.$ (59)

Thus, to estimate the bounds for $\mathcal{K}^{\Delta t}$, we only need to
study the operator $P^{\Delta t}$. Moreover, it is sufficient to prove that
there exist a probability measure $\eta$ and a constant $\epsilon>0$ so that
for any indicator function of a Borel set $S\subset\mathbb{T}^{d}$ the
following result holds

$\mathbb{P}(\bm{X}_{M}\in S|\bm{X}_{0}=\textbf{x})\geq\epsilon\eta(S),$ (60)

where $\bm{X}_{i}$ are defined in Eq.(13) as the numerical solution to the SDE
(6). The idea of the proof is to explicitly rewrite $\bm{X}_{M}$ as a
perturbation of the reference evolution corresponding to $\bm{b}=0$. According
to the SDE scheme (13), we have

$\bm{X}_{M}=\bm{X}_{0}+\bm{G}_{M}+\bm{F}_{M},$ (61)

where

$\bm{G}_{M}=\sqrt{2\kappa\Delta
t}\sum_{i=0}^{M-1}\bm{\omega}_{i},\quad\text{and},\quad\bm{F}_{M}=\Delta
t\sum_{i=0}^{M-1}\bm{b}(i\Delta t,\bm{X}_{i}).$ (62)

We know that $|\bm{F}_{m}|\leq T||\bm{b}||_{L^{\infty}}$ and $\bm{G}_{M}$ is a
Gaussian random variable with covariance matrix $2\kappa T\text{Id}_{d}$.
Therefore

$\displaystyle\mathbb{P}(\bm{X}_{M}\in S|\bm{X}_{0}=\bm{x})$
$\displaystyle\geq\mathbb{P}(\bm{G}_{M}\in S-\bm{x}-\bm{F}_{M})$
$\displaystyle=(\frac{1}{2\pi\kappa
T})^{d/2}\int_{S-\bm{x}-\bm{F}_{m}}\exp(-\frac{|\bm{y}|^{2}}{2\kappa
T})d\bm{y}.$ (63)

Since the state space $\mathbb{T}^{d}$ is compact, we can find $R>0$ such that
$|\bm{x}+F_{M}|\leq R$ for all $\bm{x}\in\mathbb{T}^{d}$. Thus, we define the
probability measure $\eta$ as

$\displaystyle\eta(S)=Z_{R}^{-1}\inf_{|Q|\leq
R}\int_{S+Q}\exp(-\frac{|\bm{y}|^{2}}{2\kappa T})d\bm{y},\quad\forall
S\subset\mathbb{T}^{d},$ (64)

where $Z_{R}$ is a normalization constant. Setting $\epsilon=Z_{R}(4\pi\kappa
T)^{-d/2}$, we can easily verify that $\eta(S)\geq
Z_{R}^{-1}\exp(-\frac{|R+1|^{2}}{2\kappa T})|S|$, which satisfies a uniform
minorization condition.

The uniform boundedness condition is automatically satisfied since $\eta$ has
a positive density with respect to Lebesgue measure.

The situation when the exact solution operator is considered can be samely
proved by changing Eq.(61) into an Ito integration form

$\displaystyle\textbf{X}^{T,\textbf{x}}=\textbf{X}^{0,\textbf{x}}+\int_{0}^{T}\textbf{b}(t,\textbf{X}^{t,\textbf{x}})dt+\int_{0}^{T}\sqrt{2\kappa}d\textbf{w}(t)$
(65)

and then go through the same procedure. ∎

We now represent an important result that ensures the existence of the
limiting messure for the discretized Feymann-Kac dynamics. The detailed proof
of Theorem 3.10 can be found in [25] or Corollary 2.5 in [27].

######  Theorem 3.10.

Suppose the minorization and boundedness conditions (58) hold true. Then, the
$\Phi_{n}^{\mathcal{K}^{\Delta t}}$ admits an invariant measure $\nu_{\Delta
t}$, whose density function is the eigenfunction of the operator
$(\mathcal{K}^{\Delta t})^{\star}$, adjoint operator of the solution operator
$\mathcal{K}^{\Delta t}$. Moreover, for any initial distribution
$\nu_{0}\in\mathcal{P}(\mathbb{T}^{d})$

$\big{|}\big{|}\Phi_{n}^{\mathcal{K}^{\Delta t}}(\nu_{0})-\nu_{\Delta
t}\big{|}\big{|}_{TV}\leq 2(1-\frac{\epsilon}{\gamma})^{n},$ (66)

where $||\cdot||_{TV}$ is the total variation norm and $0<\epsilon<\gamma$ are
from the minorization and boundedness conditions introduced in (58). This is
also true when changing $\mathcal{K}^{\Delta t}$ to the exact solution
operator $\mathcal{U}(T,0)$.

######  Corollary 3.11.

The principal eigenvalue $\mu_{\Delta t}$ of $\mathcal{K}^{\Delta t}$
satisfies the following relation

$e^{\mu_{\Delta t}(\lambda)T}=\nu_{\Delta t}\mathcal{K}^{\Delta
t}1=\Phi_{n}^{\mathcal{K}^{\Delta t}}(\nu_{0})\mathcal{K}^{\Delta
t}1+\rho_{n},$ (67)

where $\nu_{0}$ is any bounded non-negative initial probability measure, $T$
is the period of the time parameter, and
$\rho_{n}=O(1-\frac{\epsilon}{\gamma})^{n}$.

###### Proof.

Lemma 3.10 implies that for any bounded non-negative measure $\nu_{0}$, the
measure $\Phi_{n}^{\mathcal{K}^{\Delta t}}(\nu_{0})$ converges to an invariant
measure $\nu_{\Delta t}$ in the weak sense, that is

$\nu_{\Delta t}\phi:=\int_{\mathbb{T}^{d}}\phi d\nu_{\Delta
t}=\Phi_{n}^{\mathcal{K}^{\Delta
t}}(\nu_{0})(\phi)+O(1-\frac{\epsilon}{\gamma})^{n},$ (68)

for any bounded non-negative measurable function $\phi$.

Then, we take $\phi=\mathcal{K}^{\Delta t}1$. From the fact that the density
function of $\nu_{\Delta t}$ is the eigenfunction of the operator
$(\mathcal{K}^{\Delta t})^{\star}$, we get that

$\nu_{\Delta t}(\mathcal{K}^{\Delta t}1)=((\mathcal{K}^{\Delta
t})^{\star}\nu_{\Delta t})1=e^{\mu_{\Delta t}(\lambda)T}(\nu_{\Delta
t}1)=e^{\mu_{\Delta t}(\lambda)T}.$ (69)

Thus we finish the proof. ∎

Now we compute the principal eigenvalue $\mu_{\Delta t}(\lambda)$.

######  Lemma 3.12.

Denote $\nu_{\Delta t}^{k}=\prod_{i=0}^{k-1}\Phi^{\mathcal{K}^{\Delta
t,i}}\nu_{\Delta t}$. Let $e_{k}=(\nu_{\Delta t}^{k})(\mathcal{K}^{\Delta
t,k}1)$ denote the changing of mass. We have

$e^{\mu_{\Delta
t}(\lambda)T}=\prod_{k=0}^{M-1}e_{k},\quad\text{and}\quad\mu_{\Delta
t}(\lambda)=\frac{1}{M\Delta t}\sum_{k=0}^{M-1}\log(e_{k}).$ (70)

###### Proof.

It is easy to verify that

$\displaystyle\nu_{\Delta t}^{M}=\nu_{\Delta t}^{0}=\nu_{\Delta
t},\quad(\mathcal{K}^{\Delta t,k})^{*}\nu_{\Delta t}^{k}=e_{k}\nu_{\Delta
t}^{k+1},$ (71)

for some positive numbers $e_{k}$’s. These $e_{k}$’s are referred to the
changing of the mass for each small step $\mathcal{K}^{\Delta t,k}$. Thus, we
have $(\mathcal{K}^{\Delta t})^{*}\nu_{\Delta
t}^{0}=(\prod_{k=0}^{M-1}e_{k})\nu_{\Delta t}^{0}$, which means
$e^{\mu_{\Delta t}(\lambda)T}=\prod_{k=0}^{M-1}e_{k}$. ∎

Finally, we obtain the error estimate of the Lagrangian particle method in
computing the principal eigenvalue of parabolic operators as follows.

######  Theorem 3.13.

Let $\xi_{k}^{p,n-1}$, $k=1,\cdots,M$, $p=1,\cdots,N$, $n\in\mathbb{Z}^{+}$,
and $\Delta t$ are defined in the $N$-IPS system Algrithm 1. Then, we have the
following convergence result

$\lim_{N\to\infty}(M\Delta
t)^{-1}\sum_{i=1}^{M}\log\Big{(}N^{-1}\sum_{p=1}^{N}\exp(c(t_{i},\widetilde{\xi}_{i}^{p,n-1})\Delta
t)\Big{)}=\mu(\lambda)+O((1-\frac{\epsilon}{\gamma})^{n})+O((\Delta
t)^{\frac{1}{2}}),$ (72)

where $\epsilon$ and $\gamma$ are defined in Theorem 3.10.

###### Proof.

By the theory of $N$-IPS system, $\\{\xi_{k}^{p,n-1}\\}_{p=1,\cdots,N}$ will
approximately distributed as $\prod_{i=0}^{k-1}\Phi^{\mathcal{K}^{\Delta
t,i}}\Phi_{n-1}^{\mathcal{K}^{\Delta t}}\nu_{0}$ when $N\to\infty$. Let
$e_{k,n}^{N}=N^{-1}\sum_{p=1}^{N}\text{exp}(c(t_{k},\widetilde{\xi}_{k}^{p,n-1})\Delta
t)$ denote the increasing of the mass for each small step $\mathcal{K}^{\Delta
t,k}$. Then, we can get that

$\displaystyle\lim_{N\to\infty}e_{k,n}^{N}=(\prod_{i=0}^{k-1}\Phi^{\mathcal{K}^{\Delta
t,i}}\Phi_{n-1}^{\mathcal{K}^{\Delta t}}\nu_{0})(\mathcal{K}^{\Delta t,k}1).$
(73)

According to Theorem 3.10, we have that $\Phi_{n-1}^{\mathcal{K}^{\Delta
t}}\nu_{0}=\nu_{\Delta t}+\delta_{n}$, where $||\delta_{n}||_{TV}\leq
2(1-\frac{\epsilon}{\gamma})^{n}$. This implies

$\displaystyle\lim_{N\to\infty}e_{k,n}^{N}=(\prod_{i=0}^{k-1}\Phi^{\mathcal{K}^{\Delta
t,i}}\nu_{\Delta t})(\mathcal{K}^{\Delta
t,k}1)+O\big{(}(1-\frac{\epsilon}{\gamma})^{n}\big{)}.$ (74)

Combining Lemma 3.12, we conclude that

$\displaystyle\lim_{N\to\infty}(M\Delta
t)^{-1}\sum_{i=1}^{M}\log(e_{k,n}^{N})$ $\displaystyle=(M\Delta
t)^{-1}\sum_{k=0}^{M-1}\log(e_{k})+O\big{(}(1-\frac{\epsilon}{\gamma})^{n}\big{)}$
$\displaystyle=\mu_{\Delta
t}(\lambda)+O\big{(}(1-\frac{\epsilon}{\gamma})^{n}\big{)}.$ (75)

From Theorem 3.7, we know that $\big{|}\mu(\lambda)-\mu_{\Delta
t}(\lambda)\big{|}={O}\big{(}(\Delta t)^{\frac{1}{2}}\big{)}$, which concludes
the proof. ∎

##  4 Numerical results

In this section, we first present numerical examples to verify the convergence
analysis of the proposed method in computing eigenvalues. Then, we compute the
KPP front speeds in 2D and 3D chaotic flows. In addition, we investigate the
dependence of the KPP front speed on the magnitude of velocity field and the
evolution of the empirical distribution of the $N$-IPS. To be consistent with
the setting of numerical experiments in the literature, e.g., [40, 39], we
choose the torus space $\mathbb{T}^{d}=[0,2\pi]^{d}$, $d=2,3$.

###  4.1 Convergence tests in computing principal eigenvalue

We first verify the convergence of the operator splitting method in
approximating solution operator. Let $\textbf{x}=(x_{1},x_{2})^{T}$. We
consider a two-dimensional non-autonomous equation on $[0,2\pi]^{2}$ as
follows:

$u_{t}=\mathcal{L}(t)u+\mathcal{C}(t)u,$ (76)

where $\mathcal{L}(t)=\Delta_{\textbf{x}}+(\sin(x_{2})\cos(2\pi
t),\sin(x_{1})\cos(2\pi t))\cdot\nabla_{\textbf{x}}$ and
$\mathcal{C}(t)=\big{(}\sin(x_{1}+x_{2})+\cos(x_{1}+x_{2})\big{)}\sin(2\pi
t)$.

We use spectral method to discretize Eq.(76), in order to obtain an accurate
approximation in the physical space of the solution operator of Eq.(76).
Speficially, let $V_{H}=\text{span}\\{e^{i(k_{1}x_{1}+k_{2}x_{2})}:-H\leq
k_{1},k_{2}\leq H\\}$ denote a finite dimensional space spanned by fourier
basis functions, where $H$ is a positive integer. First, we compute the
approximations of the operators $\mathcal{L}(t)$ and $\mathcal{C}(t)$ in the
space $V_{H}$. Let matrices $L^{H}(t)$ and
$C^{H}(t)\in\mathbb{C}^{(2H+1)^{2}\times(2H+1)^{2}}$ denote the approximations
of $\mathcal{L}(t)$ and $\mathcal{C}(t)$, respectively [38]. Then, we use the
matrix exponential functions $e^{\Delta tL^{H}(t)}$ and $e^{\Delta tC^{H}(t)}$
to approximate $e^{\Delta t\mathcal{L}(t)}$ and $e^{\Delta t\mathcal{C}(t)}$,
respectively. Thus, we get an approximation formula for $\mathcal{K}^{\Delta
t}$ as

$K^{H,\Delta t}=\prod_{j=0}^{T/\Delta t-1}e^{\Delta tL^{H}(t_{j})}e^{\Delta
tM^{H}(t_{j})}.$ (77)

For the reference solution, we choose a much finer time step $\Delta t_{ref}$
and compute the approximation formula

$\widetilde{K}^{H,\Delta t_{ref}}=\prod_{j=0}^{T/\Delta t_{ref}-1}e^{\Delta
t_{ref}(L^{H}(t_{j})+M^{H}(t_{j}))}.$ (78)

In this experiment, we choose $H=24$, $\Delta t=2^{-1},2^{-2},\cdots,2^{-9}$,
and $\Delta t_{ref}=2^{-12}$. Then, we compute $||K^{H,\Delta
t}-\widetilde{K}^{H,\Delta t_{ref}}||_{L^{2}}$ to verify our result. Figure 1
shows the convergence results for the splitting method. The convergence rate
is $(\Delta t)^{1.05}$. This numerical result suggests that the convergence
analysis in Theorem 3.7 is not sharp. More studies on the convergence analysis
of our method will be reported in our future work.

Figure 1: Numerical errors for $||K^{H,\Delta t}-\widetilde{K}^{H,\Delta
t_{ref}}||_{L^{2}}$.

Then, we test the convergence of the Lagrangian method, i.e., Algorithm 1, in
computing principal eigenvalues of parabolic-type equations. We still consider
the problem (76) with the same $\mathcal{L}(t)$ and $\mathcal{C}(t)$. In this
experiment, we choose $\Delta t=2^{-1},2^{-2},2^{-3},2^{-4},2^{-5}$,
$N=200,000$ in the $N$-IPS system, and iteration number $n=200$ and $n=400$ in
the Feynman-Kac semigroup iteration method. Figure 2 shows the convergence of
principal eigenvalues with respect to $\Delta t$ by spectral method and our
Lagrangian method, where the reference solution is computed from spectral
method with a finer grid $\Delta t_{ref}=2^{-10}$. So given sufficient large
$N$ and $n$, the error in calculating principal eigenvalues of linearized KPP
operator $\mathcal{A}$ via our proposed Lagrangian approach only comes from
the error of operator splitting. Also as the Lagrangian method will eventually
converges to some invariant measure approximating the ground truth invariant
measure, there is no error accumulation for long time integration.

Figure 2: In the Lagrangian method, iteration number $n=200$ and $n=400$. The
reference solution is obtained by the spectral method.

###  4.2 Computing KPP front speeds in different flows

We first compute the KPP front speeds in two different time-independent flows,
i.e., a 2D steady cellular flow and a 3D ABC flow. Let
$\textbf{x}=(x_{1},...,x_{d})^{T}\in[0,2\pi]^{d}$ with $d=2,3$. We use the
Lagrangian method to compute the following principal eigenvalue problem with
periodic boundary condition

$\displaystyle\kappa\Delta_{\textbf{x}}\Phi+(2\kappa{\lambda}{\bf
e+v})\cdot\nabla_{\textbf{x}}\Phi+\big{(}\kappa{\lambda}^{2}+{\lambda}{\bf
e\cdot v}+\tau^{-1}f^{\prime}(0)\big{)}\Phi=\mu(\lambda)\Phi,$ (79)

$f(u)=u(1-u)$, and $(\mu(\lambda),\Phi)$ are principal eigenvalue of (79) and
its associated eigenfunction, respectively. The velocity field ${\bf v}=(-\sin
x_{1}\cos x_{2},\cos x_{1}\sin x_{2})$ in the 2D steady cellular flow and
${\bf v}=(\sin x_{3}+\cos x_{2},\sin x_{1}+\cos x_{3},\sin x_{2}+\cos x_{1})$
in the 3D ABC flow, respectively.

We choose the parameters $\kappa=1$ and $\tau=1$ in (79). We use the spectral
method to obtain an accurate reference solution for the principal eigenvalue
of (79). Figure 3 shows the convergence results of the Lagrangian method in
computing the principal eigenvalue, where $\lambda=0.35$ for the 2D cellular
flow and $\lambda=0.55$ for the 3D ABC flow. We find the convergence rate of
the Lagrangian method is $(\Delta t)^{1.51}$ for the 2D steady cellular flow,
and $(\Delta t)^{1.70}$ for the 3D ABC flow. Thus, we can use the Lagrangian
method to compute the KPP front speeds in both 2D and 3D flows.

(a) 2D convergence test, fitted slope $\approx$ 1.51

(b) 3D convergence test, fitted slope $\approx$ 1.60.

Figure 3: Errors of the principal eigenvalue computed by using different time
steps.

After getting the principal eigenvalue, we compute the KPP front speed $c^{*}$
through the formula $c^{*}=\inf_{{\lambda}>0}\frac{\mu(\lambda)}{\lambda}$. We
only show the numerical results for the 3D ABC flow here since the results for
the 2D steady cellular flow is quantitatively similar. We choose the velocity
field ${\bf v}=A(\sin x_{3}+\cos x_{2},\sin x_{1}+\cos x_{3},\sin x_{2}+\cos
x_{1})$, where $A$ is the strength of the convection. In Figure 4, we show the
results of $\frac{\mu(\lambda)}{\lambda}$ for ABC flows with $A=1$ and $A=10$.
The amplitude of the principal eigenvalue increases fast and the convergence
speed becomes slower. Notice that in this case, the flow becomes very unstable
since the convection becomes dominant comparing to the diffusion. This issue
will be studied in subsection 4.3.

(a) A = 1

(b) A = 10

Figure 4: Numerical results of $\frac{\mu({\lambda})}{{\lambda}}$ for
different $\lambda$’s in the ABC flow.

Next, we compute the KPP front speed in a 2D unsteady (time-dependent)
cellular flow. Let $\textbf{x}=(x_{1},x_{2})^{T}$. We use the Lagrangian
method to compute the following principal eigenvalue problem with periodic
boundary condition

$\displaystyle\kappa\Delta_{\textbf{x}}\Phi+(2\kappa\lambda{\bf e}+{\bf
v})\cdot\nabla_{\textbf{x}}\Phi+\big{(}\kappa{\lambda}^{2}+{\lambda}{\bf
e\cdot v}+\tau^{-1}f^{\prime}(0)\big{)}\Phi-\Phi_{t}=\mu(\lambda)\Phi,$ (80)

where $(t,\textbf{x})\in[0,T]\times[0,2\pi]^{2}$, $T$ is the period of v in
$t$, $f(u)=u(1-u)$, and $(\mu(\lambda),\Phi)$ are principal eigenvalue of (80)
and its associated eigenfunction, respectively. The velocity field of the 2D
unsteady cellular flow is ${\bf v}=\big{(}-\sin x_{1}\cos x_{2}(1+\delta\cos
2\pi t),\cos x_{1}\sin x_{2}(1+\delta\cos 2\pi t)\big{)}$, where $\delta>0$ is
a parameter.

We choose the parameters $\kappa=1$ and $\tau=1$ in (80) and $\delta=0.5$ in
the velocity field ${\bf v}$. We use the spectral method to obtain an accurate
reference solution for the principal eigenvalue of (80). For figure 5(a), we
choose $\lambda=0.57$. Figure 5(a) shows the convergence results of the
Lagrangian method in computing the principal eigenvalue, where the convergence
rate is $(\Delta t)^{1.31}$. Figure 5(b) shows the numerical results of
$\frac{\mu(\lambda)}{\lambda}$ for different $\lambda$’s, from which we can
compute the KPP front speed in the 2D unsteady cellular flow. We can see that
$\frac{\mu(\lambda)}{\lambda}$ is convex within the computational domain of
$\lambda$. Thus, we can compute the KPP front speed by finding the minimizer
of $\frac{\mu(\lambda)}{\lambda}$.

(a) Convergence test for different $\Delta t$’s. The fitted slope is $\approx
1.31$.

(b) Numerical results of $\frac{\mu({\lambda})}{{\lambda}}$ for different
$\lambda$’s.

Figure 5: Numerical results for a 2D unsteady cellular flow.

###  4.3 Investigate the dependence of front speed on the strength of the
flows

To further test the performance of the Lagrangian method, we study the
dependence of the KPP front speeds on the strength of different flows.
Moreover, we study the relationship between the KPP front speeds in the
chaotic flows and the effective diffusivity of the passive tracer model in the
same chaotic flows. We refer the interested reader to [43, 42, 23] for the
recent development in computing effective diffusivities in chaotic and random
flows. We set the diffusion constant $\kappa=1$ and the time scale of reaction
rate $\tau=1$.

Let us first consider this issue in KPP front speeds of time-independent
flows. If we scale ${\bf v}\rightarrow A{\bf v}$, Eq.(79) can be rewritten as
the following form

$\displaystyle\Delta_{\textbf{x}}\Phi+(2\lambda{\bf e}+A{\bf
v})\cdot\nabla_{\textbf{x}}\Phi+\big{(}{\lambda}^{2}+{\lambda}{\bf e}\cdot
A{\bf v}+f^{\prime}(0)\big{)}\Phi=\mu(\lambda)\Phi.$ (81)

The KPP front speed is
$c^{*}=\inf_{{\lambda}>0}\frac{\mu(\lambda)}{{\lambda}}$. Notice that the KPP
front speed $c^{*}$ depends on $A$, i.e., $c^{*}=c^{*}(A)$. Therefore, we
consider the equivalent equation

$\displaystyle A^{-1}\Delta_{\textbf{x}}\Phi+(2A^{-1}\lambda{\bf e}+{\bf
v})\cdot\nabla_{\textbf{x}}\Phi+\big{(}A^{-1}{\lambda}^{2}+{\lambda}{\bf
e}\cdot{\bf v}+A^{-1}f^{\prime}(0)\big{)}\Phi=\widetilde{\mu}(\lambda)\Phi,$
(82)

where $\widetilde{\mu}(\lambda)=A^{-1}\mu(\lambda)$. Let $\widetilde{c}^{*}$
denote the KPP front speed of the rescaled equation (82). We have that

$\widetilde{c}^{*}=\inf_{\lambda>0}\frac{\widetilde{\mu}(\lambda)}{\lambda}=\frac{c^{*}}{A}.$
(83)

We denote $\sigma=A^{-1}$. For the 2D steady cellular flow ${\bf v}=(-\sin
x_{1}\cos x_{2},\cos x_{1}\sin x_{2})$, it has been proved that
$c^{*}(A)=O(A^{1/4})$ [1, 34]. Let $D^{E}(A)$ denote the effective diffusivity
corresponding of the passive tracer model in the same 2D steady cellular flow
${\bf v}$. It has been proved by a boundary layer analysis that
$D^{E}(A)=O(A^{1/2})$ in [1, 6]. By scaling analysis, we obtain that for the
2D steady cellular flow the following result holds

$c^{*}(A)=O(\sqrt{D^{E}(A)}).$ (84)

To the best of our knowledge, the above relationship between the KPP front
speeds and the effective diffusivity was only proved in 2D steady cellular
flows; see [34, 36]. The result (84) implies that
$\widetilde{c}^{*}(\sigma)=\sigma O(\sigma^{-1/4})=O(\sigma^{3/4})$, which
provides a theoretical guidence for our numerical experiments. Figure 6(a)
shows the numerical results of $\widetilde{c}^{*}(\sigma)$ in the 2D steady
cellular flow obtained by our method. From the numerical results, we compute
regression and obtain $\widetilde{c}^{*}(\sigma)=O(\sigma^{0.74})$, which
agrees with the theoretical result (84).

For other flows, such as unsteady flows and 3D chaotic flows, the
understanding of $c^{*}(A)$ for large $A$’s (or $\widetilde{c}^{*}(\sigma)$
for small $\sigma$’s) remains open. We will study these flows here. In our
previous work [42], we computed the effective diffusivity of the passive
tracer model in the 3D Kolmogorov flow, where ${\bf v}=(\sin x_{1},\sin
x_{2},\sin x_{3})$, and obtained that $D^{E}(A)=O(A^{1.13})$. Notice that in
[42] the effective diffusivity is represented in terms of the diffusion and we
have converted the result in terms of the strength of the flows here, which
are equivalent. The result (84) implies that $\widetilde{c}^{*}(\sigma)=\sigma
O(\sigma^{-0.56})=O(\sigma^{0.44})$. Using our method, we compute
$\widetilde{c}^{*}(\sigma)$ for $\sigma$ in 3D Kolmogorov flow and show the
numerical results in Figure 6(b). We obtain that
$\widetilde{c}^{*}(\sigma)=O(\sigma^{0.43})$, which means that the result (84)
also holds in the 3D Kolmogorov flow. We conjecture that the result (84) also
holds true in other 3D chaotic flows. We will study this issue in future
works.

(a) Numerical results of $\widetilde{c}^{*}(\sigma)$ in 2D cellular flow. The
fitted slope is $\approx 0.74$.

(b) Numerical results of $\widetilde{c}^{*}(\sigma)$ in 3D Kolmogorov flow.
The fitted slope is $\approx 0.43$.

Figure 6: Numerical results of $\widetilde{c}^{*}(\sigma)$ in different flows.

Next, we study the dependence of the KPP front speeds on the strength of time-
dependent flows. Specifically, we will consider two 3D flows. The first one is
a time-dependent Kolmogorov flow with ${\bf
v}=\big{(}\sin(x_{3}+\theta\sin(2\pi t)),\sin(x_{1}+\theta\sin(2\pi
t)),\sin(x_{2}+\theta\sin(2\pi t))\big{)}$, and the second one is a time-
dependent ABC flow with ${\bf v}=\big{(}\sin(x_{3}+\sin(2\pi\Omega
t))+\cos(x_{2}+\sin(2\pi\Omega t)),\sin(x_{1}+\sin(2\pi\Omega
t))+\cos(x_{3}+\sin(2\pi\Omega t)),\sin(x_{2}+\sin(2\pi\Omega
t))+\cos(x_{1}+\sin(2\pi\Omega t))\big{)}$.

For the 3D time-dependent Kolmogorov flow, we choose iteration time $n=256$,
time step $\Delta t=2^{-9}$ and particle number $N=400,000$. Figure 7 shows
the result of $\widetilde{c}^{*}(\sigma)$ for small $\sigma$’s and different
$\theta$’s. Again, we find the KPP front speed $\widetilde{c}^{*}(\sigma)$ is
not very sensitive to the paramater $\theta$. When $\theta=1$, we obtain that
$\widetilde{c}^{*}(\sigma)=O(\sigma^{0.39})$.

Figure 7: Numerical results of $\widetilde{c}^{*}(\sigma)$ in a 3D time-
dependent Kolmogorov flow.

In Figure 8, we plot out procedure searching for the $\lambda$ when the
minimal in Eq.(83) was reached. We use $a\lambda+b\lambda^{-1}+c$ to fit a
curve, then find the minimum of the curve. When $\sigma$ is large, the
relative flutuation is small and the minimun is easily to be find. When
$\sigma$ is small, the relative flutuation becomes strong enough, so we decide
to fit the curve, then find the minimun point.

Figure 8: Numerical results of $\frac{\mu({\lambda})}{{\lambda}}$ for
different $\lambda$’s and $\sigma$’s in the 3D time-dependent Kolmogorov flow.
The red dash curve is fitted by $a\lambda+b\lambda^{-1}+c$.

For the 3D time-dependent ABC flow, we choose the iteration time $n=2048$
(since the ABC flow is more chaotic), time step $\Delta t=2^{-9}$, and
particle number $N=400,000$. Figure 9(a) shows the KPP front speeds
$\widetilde{c}^{*}(\sigma)$ for different $\Omega$’s, where $\Omega$ ranges
from $2^{-7}$ to $2^{0}$. Figure 9(b) shows the slope of each approximation
line for each $\Omega$ in Figure 9(a). If we assume
$\widetilde{c}^{*}(\sigma)=O(\sigma^{\alpha})$ is true, the slope values in
Figure 9(b) give the power value $\alpha$’s for different $\Omega$’s. We find
that when $\Omega$ is near $0.1$, the power value $\alpha$ is large. When
$\Omega$ is away from $0.1$, say $\Omega<2^{-4}$ or $\Omega>2^{-2}$, the power
value $\alpha$ is small. A similar sensitive dependence on the frequency of
time-dependent ABC flows was reported in [4], where the Lyapunov exponent of
the deterministic time-dependent ABC flow problem (i.e., $\Omega=0$) was
studied as the indicator of the extent of chaos; see Figure 2 and Figure 3 of
[4].

(a) Numerical results of $\widetilde{c}^{*}(\sigma)$ for different $\Omega$’s
and different $\sigma$’s.

(b) Power values of $\widetilde{c}^{*}(\sigma)=O(\sigma^{\alpha})$ for
different $\Omega$’s.

Figure 9: Numerical results for the time-dependent ABC flows.

We compare the computational time of the interacting particle method and the
spectral method in the 2D cellular flow example. The numerical experiments are
carried out on the same core of the HPC2015 system at HKU with 10-core Intel
Xeon E5-2600 v3(Haswell) processors and 96 GB physical memory. We compute the
front speed using the spectral method mentioned in Section 4.1. We set the
Fourier modes $H=2^{k}$ and $k$ is a positive integer. When $\sigma=2^{-2}$,
for the spectral method, $H=2^{3}$ is enough and it spends $1.13$ seconds to
calculate the front speed; while for our interacting particle method, the
computational time is about $45.01$ seconds. When $\sigma=2^{-5}$, for the
spectral method, $H=2^{4}$ is enough and it spends $42.35$ seconds to
calculate the front speed, and the interacting particle method costs $172.76$
seconds. When $\sigma=2^{-8}$, for the spectral method, $H=2^{4}$ is needed
and it costs $1203.12$ seconds to calculate the front speed; on the other
hand, our interacting particle method costs $676.23$ seconds. When $\sigma$
becomes extremely small, the spectral method becomes very expensive, however,
our interacting particle method is still very efficient. For instance, when
$\sigma=2^{-12}$, the spectral method may need several days to calculate the
front speed, but our interacting particle method only costs $5378.24$ seconds.
We remark that the spectral method becomes very expensive in computing front
speeds for 3D chaotic flows. However, the computational time of the
interacting particle method only weakly depends on the dimension of the
physical space. Thus, we can compute KPP front speeds in 3D chaotic flows.

###  4.4 Evolution of the empirical distribution of the particles

As stated in Theorem 3.10, the empirical distribution converges to the
invariant measure of Feynman-Kac semigroup as $n$ approaches infinity. Our
Lagrangian method can not only calculate the principal eigenvalue but also
compute the evolution of the distribution. In this subsection, we will study
the empirical distribution of the $N$-IPS system moduled to torus space. In
this subsection, the particle number for all tests is $N=200,000$.

Figure 10 shows the invariant distribution generated by the $N$-IPS system in
the 2D steady cellular flow, where ${\bf v}=(-\sin x_{1}\cos x_{2},\cos
x_{1}\sin x_{2})$. The parameter $\sigma$ varies from $2^{0}$ to $2^{-5}$. The
strength of the convection is then proportion to $1/\sigma$. Then we can see
from the graph when increasing the strength, the main part of invariant
measure gets into a smaller domain and the gradient becomes sharper, which is
a common phenomenon in fluid dynamics. In addition, by comparing to the
pattern at the boundary of the plot, one can find that the invariant measure
is periodic in physical space.

Figure 10: Empirical distributions for the 2D steady cellular flow with
$\sigma$ varies from $2^{0}$ to $2^{-5}$. First row from left to right:
$\sigma=2^{0}$, $\sigma=2^{-1}$, and $\sigma=2^{-2}$. Second row from left to
right: $\sigma=2^{-3}$, $\sigma=2^{-4}$, and $\sigma=2^{-5}$.

Next, we study the evolution of invariant distribution generated by our
$N$-IPS system in a 2D time-periodic mixing flow, where ${\bf v}=(-\cos
x_{2}-\theta\sin x_{1}\cos(2\pi t),\cos x_{1}+\theta\sin x_{2}\cos(2\pi t))$.
Figure 11 shows the empirical distribution of the $N$-IPS system at different
times within one period when the iteration time $n=400$. From these numerical
results, we can see the invariant distribution varies at different times of
one period. But as the starting time of different periods, the first and last
subfigure are identical. These are consistent with our expectation obtained in
Lemma 3.12, where we proved that the invariant measure changes periodically
with the same period as the flow.

Figure 11: Empirical distributions for the 2D time-periodic mixing flow with
$\theta=1$, $\sigma=1$, in different phase of one period: $t$ varies from $0$
to $1$ with time interval equal to $1/9$.

Finally, we let the parameter $\sigma$ vary from $2^{0}$ to $2^{-5}$ and study
the evolution of invariant distribution generated by our $N$-IPS system in the
2D time-periodic mixing flow. Figure 12 shows that with the increasing of the
strength of the convection, the invariant measure becomes compactly supported
with a sharp gradient.

From these numerical results, we get two conclusions. First, the invariant
measure of the Feynman-Kac semigroup associated with the KPP operator is no
longer uniform distribution. This is due to the effect from the potential
function $c(t,\textbf{x})$. Second, the invariant measure converges to a
limiting measure as $\sigma\to 0$. Notice that when $\sigma$ is small, the
invariant measure develops sharp gradients, which requires more particles to
compute. Moreover, it may take more iteration time steps to converges.
Developing effective sampling method to compute invariant measure for the KPP
operator with small diffusion constant will be studied in our future works.

Figure 12: Empirical distributions for the 2D time-periodic mixing flow with
$\theta=1$, $\sigma$ varies from $2^{0}$ to $2^{-5}$. First row from left to
right: $\sigma=2^{0}$, $\sigma=2^{-1}$, and $\sigma=2^{-2}$. Second row from
left to right: $\sigma=2^{-3}$, $\sigma=2^{-4}$, and $\sigma=2^{-5}$.

##  5 Conclusion

In this paper, we developed efficient Lagrangian particle methods to compute
the KPP front speeds in time-periodic cellular and chaotic flows and provided
rigorous convergence analysis for the numerical schemes. In the convergence
analysis, we first obtained the error of the operator splitting methods in
approximating the solution operator corresponding to the linearized KPP
equation. Then, we proved the convergence of the Lagrangian particle method in
computing the principal eigenvalue based on the Feynman-Kac semigroup theory.
Finally, we presented numerical results to verify the convergence rate of the
proposed method for computing the principal eigenvalues. In addition, we
computed the KPP front speeds in several typical chaotic flow problems of
physical interests, including the Arnold-Beltrami-Childress (ABC) flow and the
Kolmogorov flow. It has been proved that the KPP front speed and the effective
diffusivity satisfies the relation $c^{*}(A)=O(\sqrt{D^{E}(A)})$ in 2D
cellular flows [34, 36]. We numerically verified this relation and found that
this relation still holds in 3D Kolmogorov flows and ABC flows.

There are three directions we plan to explore in our future work. First, we
will extend the Lagrangian particle method to compute KPP front speeds in
time-stochastic and space-periodic flows. Second, we will develop Lagrangian
particle methods to compute KPP fronts speeds in more complex fluid flows,
where the computational domain is not compact. This type of problem is more
challenging both analytically and numerically. As stated in the introduction
part, there is limited literature on studying the existence of KPP front
speeds in complex flows. In the aspect of numerical computation, our current
method cannot be adapted to non-compact domains. We shall adopt some
relaxation techniques to address this problem. In addition, we shall develop
adaptive sampling methods for our Lagrangian particle methods in order to
resolve the sharp gradients in the invariant measure when the magnitude of the
velocity field is very large.

## Acknowledgement

The research of J. Lyu and Z. Wang are 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 the HKU ITS research computing
facilities that are supported in part by the Hong Kong UGC Special Equipment
Grant (SEG HKU09).

##  Appendix A Error bounds for exponential operator splitting in non-
autonomous evolution equations

###  A.1 Euler methods for non-autonomous evolution equations

In this section, we review the fundamental rsults for abstract linear
evolution equations by semigroup theory; see e.g. [11, 5] for more details. We
consider the non-autonomous Cauchy problem (NCP) as follows

$\left\\{\begin{aligned} \frac{d}{dt}u(t)&=\mathcal{A}(t)u(t),\quad t\geq
s\in\mathbb{R}\\\ u(s)&=x\in X,\end{aligned}\right.$ (85)

where $X$ is a Banach space and
$(\mathcal{A}(t),\mathcal{D}(\mathcal{A}(t)))_{t\in\mathbb{R}}$ is a family of
linear operators on $X$.

######  Definition A.14.

A continuous function $u:[s,\infty)\longrightarrow X$ is called a classical
solution of (85) if $u\in
C^{1}([s,\infty);X),u(t)\in\mathcal{D}(\mathcal{A}(t))$ for all $t\geq s$,
$u(s)=x$, and $\frac{d}{dt}u(t)=\mathcal{A}(t)u(t)$ for all $t\geq s$.

######  Definition A.15.

For a family $(\mathcal{A}(t),\mathcal{D}(\mathcal{A}(t)))_{t\in\mathbb{R}}$
of linear operators on a Banach space $X$, the NCP (85) is well-posed with
regularity subspace $(Y_{s})_{s\in\mathbb{R}}$ and exponentially bounded
solutions, if

  1. 1.

(Existence) For all $s\in R$ the subspace

$Y_{s}=\\{y\in X:\text{there exists a classical solution for the NCP
\eqref{NCPE-appendix}}\\}\subset\mathcal{D}(\mathcal{A}(s))$ (86)

is dense in X.

  2. 2.

(Uniqueness) For every $y\in Y_{s}$, the solution $u_{s}(\cdot,y)$ is unique.

  3. 3.

(Continuous dependence) The solution continuously depends on $s$ and $y$,
i.e., if $s_{n}\to s\in\mathbb{R},||y_{n}-y||_{X}\to 0$ with $y_{n}\in
Y_{y_{n}}$, then we have $||\hat{u}_{s_{n}}(t,y_{n})-\hat{u}_{s}(t,y)||_{X}\to
0$ uniformly for $t$ in compact subsets of $\mathbb{R}$, where

$\hat{u}_{s}(t,y)=\left\\{\begin{aligned} &u_{r}(t,y)&if~{}&r\leq t,\\\
&y&if~{}&r>t.\end{aligned}\right.$

  4. 4.

(Exponential boundedness) There exists a constant $\omega\in\mathbb{R}$ such
that

$||u_{s}(t,y)||_{X}\leq e^{\omega(t-s)}||y||_{X}$

for all $y\in Y_{s}$ and $t\geq s$.

######  Definition A.16.

A family $\\{\mathcal{U}(t,s),t\geq s\\}$ of linear, bounded solution
operators on Banach space $X$ is called an exponentially bounded evolution
family if

  1. 1.

$\mathcal{U}(t,r)\mathcal{U}(r,s)=\mathcal{U}(t,s)$ and $\mathcal{U}(t,t)=Id$
hold for all $t\geq r\geq s\in\mathbb{R}$,

  2. 2.

the mapping $(t,s)\to\mathcal{U}(t,s)$ is strongly continuous,

  3. 3.

$||\mathcal{U}(t,s)||_{X}\leq e^{\omega(t-s)}$ for some $\omega\in\mathbb{R}$
and all $t\geq s\in\mathbb{R}$.

In contrast to the behavior of $\text{C}_{0}$-semigroups, the algebraic
proposition of an evolution family do not imply any differentiability on a
dense subspace. Therefore, we need extra assumptions in order to solve an NCP.

######  Definition A.17.

An evolution family $\\{\mathcal{U}(t,s),t\geq s\\}$ is called evolution
family solving NCP (85) if for every $s\in\mathbb{R}$ the regularity space

$Y_{s}=\\{y\in X:[s,\infty)\ni t\mapsto\mathcal{U}(t,s)y\text{ solves NCP
\eqref{NCPE-appendix}}\\}$

is dense in X.

In this case, the unique classical solution of the NCP (85) is given by
$u(t)=\mathcal{U}(t,s)x$. The well-posedness of the NCP (85) can now be
characterized by the existence of solving an evolution family
$\\{\mathcal{U}(t,s),t\geq s\\}$.

######  Proposition A.18.

Let $X$ be a Banach space and
$(\mathcal{A}(t),\mathcal{D}(\mathcal{A}(t)))_{t\in\mathbb{R}}$ be a family of
linear operators on $X$. The following assertions are equivalent [11].

  1. 1.

The NCP (85) is well-posed.

  2. 2.

There exits a unique evolution family $\\{\mathcal{U}(t,s),t\geq s\\}$ solving
the NCP (85).

In addition, if $||e^{\tau\mathcal{A}(t)}||_{X}\leq e^{\omega\tau}$ for any
$\tau\geq 0,t\in\mathbb{R}$, then we have $||\mathcal{U}(t,s)||_{X}\leq
e^{\omega(t-s)}$.

The well-posedness of non-autonomous evolution equations is complicated and
there is no general theory describing it. Conditions implying well-posedness
are generally divided into parabolic-type assumptions and hyperbolic-type
ones. Due to the property of the KPP equation, we only study the parabolic-
type conditions in this paper, where the domain $(\mathcal{D}(\mathcal{A}(t))$
is independent of $t\in\mathbb{R}$. We refer the interested reader to [37] for
more general cases.

######  Assumption A.19.

(Parabolic-type conditions)

  1. 1.

The domain $\mathcal{D}=\mathcal{D}(\mathcal{A}(t))$ is independent of
$t\in\mathbb{R}$.

  2. 2.

For each $t\in R$ the operator $\mathcal{A}(t)$ is sectorial and generates an
analytic semigroup $e^{\cdot\mathcal{A}(t)}$. For all $t\in\mathbb{R}$, the
resolvent $\mathcal{R}({\gamma_{1}},\mathcal{A}(t))$ exists for all
${\gamma_{1}}\in\mathbb{C}$ with $\text{Real}{\gamma_{1}}\geq 0$ and there is
a constant $M\geq 1$ such that

$\big{|}\big{|}R({\gamma_{1}},\mathcal{A}(t))\big{|}\big{|}_{X}\leq\frac{M}{|{\gamma_{1}}|+1}$
(87)

for $\text{Real}{\gamma_{1}}\geq 0$ and $t\in\mathbb{R}$. The semigroups
$e^{\cdot\mathcal{A}(t)}$ satisfy $||e^{\tau\mathcal{A}(t)}||_{X}\leq
e^{\omega\tau}$ for some constant $\omega\in\mathbb{R}$.

  3. 3.

There exist constants $L\geq 0$ and $0<\theta\leq 1$ such that

$\big{|}\big{|}(\mathcal{A}(t)-\mathcal{A}(s))\mathcal{A}(0)^{-1}\big{|}\big{|}_{X}\leq
L|t-s|^{\theta},\text{ for all }t,s\in\mathbb{R}.$ (88)

To obtain a convergence estimate for the operator in certain norm, we need an
additional assumption on $\mathcal{A}(t)$ as follows.

######  Assumption A.20.

The operator $\mathcal{A}(t)$ satisfies a Hölder continuous condition. Namely,
there exists $0\leq\alpha<\beta$ such that for any
$x\in\mathcal{D}(\mathcal{A})$,

$\big{|}\big{|}(\mathcal{A}(t)-\mathcal{A}(s))x\big{|}\big{|}_{X}\leq
C|t-s|^{\beta}||\mathcal{A}(\tau)x||^{\alpha}_{X}||x||_{X}^{1-\alpha},$ (89)

for any $s\leq\tau\leq t$.

For forward Euler type discretization, Assumption A.20 can be relaxed to
$\tau=s$ only. The backword Euler type discretization needs $\tau=t$, and
other discretization methods need different $\tau$’s instead. For analytic
semigroups, the following estimate holds true [11, 35].

######  Lemma A.21.

Let $e^{t\mathcal{A}}$ be an anlytical semigroup on $X$. Let $\mathcal{A}$ be
the infinitesimal generator. There is a constant $C\geq 0$ such that

$||\mathcal{A}e^{t\mathcal{A}}||_{X}\leq\frac{C}{t},\quad t>0,\quad
0\leq\alpha\leq 1.$ (90)

Now we state the first result, which gives the approximation error of the
freezing time coefficients methods for solving the NCP (85).

######  Theorem A.22.

Suppose assumptions A.19 and A.20 hold true. Let $\mathcal{U}(T,0)$ be the
solution operator associated with the NCP (85). Then the solution operator
obtained by the freezing time coefficients methods has the following
approximation error to $\mathcal{U}(T,0)$

$\big{|}\big{|}\mathcal{U}(T,0)-\prod_{k=0}^{M-1}e^{\Delta
t\mathcal{A}(k\Delta t)}\big{|}\big{|}_{X}\leq C(T)(\Delta t)^{\beta-\alpha},$
(91)

where $T>0$, $M$ is an integer, and $\Delta t=\frac{T}{M}$.

###### Proof.

First we refer to [37] for the abstract version of the method of freezing
coefficients,

$\mathcal{U}(t,s)=e^{(t-s)\mathcal{A}(s)}+\int_{s}^{t}\mathcal{U}(t,\tau)(\mathcal{A}(\tau)-\mathcal{A}(s))e^{(\tau-s)\mathcal{A}(s)}d\tau$
(92)

which immediately gives us that, for every $x\in X$,

$\displaystyle\big{|}\big{|}(\mathcal{U}(t,s)-e^{(t-s)\mathcal{A}(s)})x\big{|}\big{|}_{X}$
$\displaystyle=$
$\displaystyle\big{|}\big{|}\int_{s}^{t}\mathcal{U}(t,\tau)(\mathcal{A}(\tau)-\mathcal{A}(s))e^{(\tau-s)\mathcal{A}(s)}xd\tau\big{|}\big{|}_{X}$
$\displaystyle\leq$
$\displaystyle\int_{s}^{t}\big{|}\big{|}\mathcal{U}(t,\tau)\big{|}\big{|}_{X}(\tau-s)^{\beta}\big{|}\big{|}\mathcal{A}(s)e^{(\tau-s)\mathcal{A}(s)}x\big{|}\big{|}_{X}^{\alpha}\big{|}\big{|}e^{(\tau-s)\mathcal{A}(s)}x\big{|}\big{|}^{1-\alpha}_{X}d\tau.$
(93)

In (93), we have used the fact that
$e^{(\tau-s)\mathcal{A}(s)}x\in\mathcal{D}(\mathcal{A})$ for any $x\in X$.
Notice that $\mathcal{A}(s)$ generates an analytic semigroup
$e^{\cdot\mathcal{A}(s)}$, according to (A.21) we have the following estimate

$\big{|}\big{|}\mathcal{A}(s)e^{(\tau-s)\mathcal{A}(s)}\big{|}\big{|}^{\alpha}_{X}\leq
C(\tau-s)^{-\alpha}e^{\omega\alpha(\tau-s)}.$ (94)

Substituting (94) into (93), we obtain that,

$\displaystyle\big{|}\big{|}(\mathcal{U}(t,s)-e^{(t-s)\mathcal{A}(s)})x\big{|}\big{|}_{X}$
$\displaystyle\leq$
$\displaystyle\int_{s}^{t}Ce^{\omega(t-\tau)}(\tau-s)^{\beta-\alpha}e^{\omega(\tau-s)}d\tau||x||_{X}=\frac{C}{1+\beta-\alpha}e^{\omega(t-s)}(t-s)^{1+\beta-\alpha}||x||_{X}.$
(95)

Thus, we get the estimate for the operator in the norm $||\cdot||_{X}$

$\big{|}\big{|}\mathcal{U}(t,s)-e^{(t-s)\mathcal{A}(s)}\big{|}\big{|}_{X}\leq\frac{C}{1+\beta-\alpha}e^{\omega(t-s)}(t-s)^{1+\beta-\alpha}.$
(96)

We denote $\mathcal{U}(T,0)=\prod_{k=0}^{M-1}\mathcal{U}((k+1)\Delta t,k\Delta
t)$. Using the telescoping sum argument, we obtain

$\displaystyle\big{|}\big{|}\mathcal{U}(T,0)-\prod_{k=0}^{M-1}e^{\Delta
t\mathcal{A}(k\Delta t)}\big{|}\big{|}_{X}$ (97) $\displaystyle=$
$\displaystyle\Big{|}\Big{|}\sum_{j=0}^{M-1}\prod_{k=j+1}^{M-1}U((k+1)\Delta
t,k\Delta t)\big{(}\mathcal{U}((j+1)\Delta t,j\Delta t)-e^{\Delta
t\mathcal{A}(j\Delta t)}\big{)}\prod_{l=0}^{j-1}e^{\Delta t\mathcal{A}(l\Delta
t)}\Big{|}\Big{|}_{X}$ (98) $\displaystyle\leq$
$\displaystyle\sum_{j=0}^{M-1}e^{\omega(N-j-1)\Delta
t}\frac{C}{1+\beta-\alpha}e^{\omega\Delta t}(\Delta
t)^{1+\beta-\alpha}e^{\omega j\Delta t}=\frac{Ce^{\omega
T}}{1+\beta-\alpha}(\Delta t)^{\beta-\alpha}.$ (99)

The statement in (91) is proved. ∎

For higher order operator splitting methods, in some specific situation the
higher order convergence has been proved in [18, 19]. In their works, the
assumption A.20 was largely strengthen, both for the operator $\mathcal{A}(t)$
and initial condition, and the convergence was largely depends on the graph
norm $||v||_{\alpha}:=||\mathcal{A}(t)^{\alpha}v||_{X}$. The convergence in
norm $||\cdot||_{X}$ is still open and will be our future reasearch plan.

###  A.2 Operator splitting methods for solving non-autonomous evolution
equations

We study the approximation error of operator splitting methods in solving non-
autonomous evolution equations. To be specific, we consider an abstract NCP as
follows

$\left\\{\begin{aligned}
\frac{d}{dt}u(t)&=(\mathcal{A}(t)+\mathcal{B}(t))u(t),\quad t\geq
s\in\mathbb{R}\\\ u(s)&=x\in X\end{aligned}\right.$ (100)

on a Banach space $X$, where $\mathcal{A}(t)$ and $\mathcal{B}(t)$ are linear
operators, $\mathcal{D}(\mathcal{A}(t))$ is independent of $t$ and dense in
$X$, and for each $t\in\mathbb{R}$, $\mathcal{A}(t)$, $\mathcal{B}(t)$ and
$\mathcal{A}(t)+\mathcal{B}(t)$ generate strongly continuous semigroups
$e^{\cdot\mathcal{A}(t)}$, $e^{\cdot\mathcal{B}(t)}$ and
$e^{\cdot(\mathcal{A}(t)+\mathcal{B}(t))}$, respectively. Furthermore, due to
the property of evolution equation, solving $u(t)$ and solving
$e^{{\gamma_{1}}t}u(t)$ is equivalent, we assume
$\big{|}\big{|}e^{\tau\mathcal{A}(t)}\big{|}\big{|}_{X}\leq
1$,$\big{|}\big{|}e^{\tau\mathcal{B}(t)}\big{|}\big{|}_{X}\leq
1$,$\big{|}\big{|}e^{\tau(\mathcal{A}(t)+\mathcal{B}(t))}\big{|}\big{|}_{X}\leq
1$.

We will study the NCP (100) based on the perturbation theory. We assume
$\mathcal{A}(t)$ is a sectorial operator, which generates an analytical
semigroups $e^{\cdot\mathcal{A}(t)}$, and assume $\mathcal{B}(t)$ is bounded,
thus $\mathcal{A}(t)+\mathcal{B}(t)$ is also sectorial and generates an
analytical semigroups $e^{\cdot(\mathcal{A}(t)+\mathcal{B}(t))}$, where
$\mathcal{D}(\mathcal{A}(t)+\mathcal{B}(t))=\mathcal{D}(\mathcal{A}(t))$. In
addition, we assume that the operator $\mathcal{A}(t)+\mathcal{B}(t)$
satisfies assumptions A.19 and A.20. Therefore, the corresponding evolution
family $\mathcal{U}(t,s)$ solves the NCP problem 100 and admits an Euler-type
approximation, i.e.,

$\big{|}\big{|}\mathcal{U}(T,0)-\prod_{k=0}^{M-1}e^{\Delta
t(\mathcal{A}+\mathcal{B})(k\Delta t)}\big{|}\big{|}_{X}\leq C(T)(\Delta
t)^{\beta-\alpha},$ (101)

where $T=M\Delta t$, $\alpha,\beta$ are constants defined in assumptions A.19
and A.20.

In the sequel, we analyze the error between $\prod_{k=0}^{M-1}e^{\Delta
t(\mathcal{A}+\mathcal{B})(k\Delta t)}$ and $\prod_{k=0}^{M-1}e^{\Delta
t\mathcal{A}(k\Delta t)}e^{\Delta t\mathcal{B}(k\Delta t)}$.

First, we list all the assumptions as follows:

######  Assumption A.23.

  1. 1.

${\mathcal{A}(t)}_{t\geq 0}$ and ${\mathcal{B}(t)}_{t\geq 0}$ are all linear
operators (may be unbounded) on X,

  2. 2.

$\mathcal{D}(\mathcal{A}(t))$ are the same for all $t$ and dense in $X$,

  3. 3.

$||\mathcal{B}(t)||_{X}<C$ for all $t\geq 0$,

  4. 4.

$\mathcal{A}(t)$ satisfies A.19 and $\mathcal{A}(t)+\mathcal{B}(t)$ satisfies
A.19 and A.20,

  5. 5.

$||e^{\tau\mathcal{A}(t)}||_{X}\leq 1$,$||e^{\tau\mathcal{B}(t)}||_{X}\leq
1$,$||e^{\tau(\mathcal{A}(t)+\mathcal{B}(t))}||_{X}\leq 1$ for all $\tau\geq
0$.

To obtain a convergence theorem, we need an extra assumption in $\mathcal{A}$
and $\mathcal{B}$:

######  Assumption A.24.

For the commutator
$[\mathcal{A}(t),\mathcal{B}(t)]=\mathcal{A}(t)\mathcal{B}(t)-\mathcal{B}(t)\mathcal{A}(t)$,
we assume that there is a non-negative $\gamma$ with

$\big{|}\big{|}[\mathcal{A}(t),\mathcal{B}(t)]x\big{|}\big{|}_{X}\leq
c_{1}\big{|}\big{|}\mathcal{A}(t)x\big{|}\big{|}^{\gamma}_{X}||x||_{X}^{1-\gamma},~{}\forall~{}x\in\mathcal{D}(\mathcal{A}).$
(102)

Next is a standard result from [20], and we prove it here.

######  Theorem A.25.

Suppose assumptions A.23 and A.24 are satisfied. We have the following error
estimate for the operator splitting method,

$\big{|}\big{|}(e^{\tau\mathcal{A}(t)}e^{\tau\mathcal{B}(t)}-e^{\tau(\mathcal{A}(t)+\mathcal{B}(t))})x\big{|}\big{|}_{X}\leq
C_{1}\tau^{2-\gamma}||x||_{X},~{}\forall x\in X,$ (103)

where $C_{1}$ depends only on $c_{1}$, $\gamma$ and $||\mathcal{B}||_{X}$.

###### Proof.

We use the freezing coefficient formula and obtain

$e^{\tau(\mathcal{A}(t)+\mathcal{B}(t))}x=e^{\tau(\mathcal{A}(t))}x+\int_{0}^{\tau}e^{s\mathcal{A}(t)}\mathcal{B}(t)e^{(\tau-s)(\mathcal{A}(t)+\mathcal{B}(t))}xds.$
(104)

Expressing the term $e^{(\tau-s)(\mathcal{A}(t)+\mathcal{B}(t))}$ using the
integral form (104), we have

$e^{\tau(\mathcal{A}(t)+\mathcal{B}(t))}x=e^{\tau(\mathcal{A}(t))}x+\int_{0}^{\tau}e^{s\mathcal{A}(t)}\mathcal{B}(t)e^{(\tau-s)\mathcal{A}(t)}xds+R_{1}x,$
(105)

where

$R_{1}=\int_{0}^{\tau}e^{s\mathcal{A}(t)}\mathcal{B}(t)\int_{0}^{\tau-s}e^{\sigma\mathcal{A}(t)}\mathcal{B}(t)e^{(\tau-s-\sigma)(\mathcal{A}(t)+\mathcal{B}(t))}d\sigma
ds.$ (106)

We can easily verify that the term $R_{1}$ is bounded, i.e.,
$||R_{1}||_{X}\leq\frac{1}{2}\tau^{2}||\mathcal{B}(t)||_{X}^{2}$.

On the other hand side, we express the term $e^{\tau\mathcal{B}(t)}$ into
exponential series and obtain

$e^{\tau\mathcal{A}(t)}e^{\tau\mathcal{B}(t)}x=e^{\tau\mathcal{A}(t)}x+\tau
e^{\tau\mathcal{A}(t)}\mathcal{B}(t)x+R_{2}x,$ (107)

where $||R_{2}||_{X}\leq\frac{1}{2}\tau^{2}||\mathcal{B}(t)||_{X}^{2}$.

Denoted by
$f(s)=e^{s\mathcal{A}(t)}\mathcal{B}(t)e^{(\tau-s)\mathcal{A}(t)}x$, we have

$e^{\tau\mathcal{A}(t)}e^{\tau\mathcal{B}(t)}x-e^{\tau(\mathcal{A}(t)+\mathcal{B}(t))}x=\tau
f(\tau)-\int_{0}^{\tau}f(s)ds+r=d+r,$ (108)

where $d=\tau f(\tau)-\int_{0}^{\tau}f(s)ds=\tau^{2}\int_{0}^{1}\theta
f^{\prime}(\theta\tau)d\theta$ and $r=R_{2}x-R_{1}x$.

Since
$f^{\prime}(s)=e^{s\mathcal{A}(t)}[\mathcal{A}(t),\mathcal{B}(t)]e^{(\tau-s)\mathcal{A}(t)}x$,
assumption A.24 implies

$\displaystyle\big{|}\big{|}e^{s\mathcal{A}(t)}[\mathcal{A}(t),\mathcal{B}(t)]e^{(\tau-s)\mathcal{A}(t)}x\big{|}\big{|}_{X}\leq
c_{1}||e^{s\mathcal{A}(t)}||_{X}||\mathcal{A}(t)e^{(\tau-s)\mathcal{A}(t)}x\big{|}\big{|}^{\gamma}_{X}\big{|}\big{|}e^{(\tau-s)\mathcal{A}(t)}x\big{|}\big{|}_{X}^{1-\gamma}.$
(109)

By using the property of analytic semigroup A.21, we know that

$\big{|}\big{|}\mathcal{A}(t)e^{(\tau-s)\mathcal{A}(t)}x\big{|}\big{|}_{X}\leq
C(\tau-s)^{-1}||x||_{X}.$ (110)

Thus, we have

$\displaystyle||d||_{X}=\big{|}\big{|}\tau^{2}\int_{0}^{1}\theta
f^{\prime}(\theta\tau)d\theta\big{|}\big{|}_{X}\leq$
$\displaystyle|\tau^{2}\int_{0}^{1}C\theta(\tau-\theta\tau)^{-\gamma}d\theta|||v||_{X}$
$\displaystyle=$
$\displaystyle\frac{C}{(1-\gamma)(2-\gamma)}\tau^{2-\gamma}||v||_{X}.$ (111)

Notice that $||r||_{X}\leq\tau^{2}||\mathcal{B}||_{X}^{2}$. We finish the
proof. ∎

Using the one step estimate obtained in Theorem A.25, we finally obtain the
error estimate for the operator splitting method.

######  Theorem A.26.

Suppose assumptions A.23 and A.24 hold true. We have the following error
estimate for the operator splitting method in solving the NCP (100).

$\big{|}\big{|}\prod_{k=0}^{N-1}e^{\Delta t(\mathcal{A}+\mathcal{B})(k\Delta
t)}-\prod_{k=0}^{N-1}e^{h\mathcal{A}(k\Delta t)}e^{\Delta t\mathcal{B}(k\Delta
t)}\big{|}\big{|}_{X}\leq C_{1}(\Delta t)^{1-\gamma},$ (112)

where $C_{1}$ is a constant independent of $\gamma$.

###### Proof.

We take $t=j\Delta t$ and $s=(j-1)\Delta t$ for $j=1,\cdots,M-1$ in Theorem
A.25, and by using the telescoping sum argument, we obtain that for any $x\in
X$,

$\displaystyle\big{|}\big{|}\prod_{k=0}^{M-1}e^{\Delta
t(\mathcal{A}+\mathcal{B})(k\Delta t)}x-\prod_{k=0}^{M-1}e^{\Delta
t\mathcal{A}(k\Delta t)}e^{\Delta t\mathcal{B}(k\Delta t)}x\big{|}\big{|}_{X}$
$\displaystyle=$
$\displaystyle\Big{|}\Big{|}\sum_{j=0}^{M-1}\prod_{k=j+1}^{M-1}e^{\Delta
t(\mathcal{A}+\mathcal{B})(k\Delta t)}\big{(}e^{\Delta
t(\mathcal{A}+\mathcal{B})(j\Delta t)}-e^{\Delta t\mathcal{A}(j\Delta
t)}e^{\Delta t\mathcal{B}(j\Delta t)}\big{)}\prod_{l=0}^{j-1}e^{\Delta
t\mathcal{A}(l\Delta t)}e^{\Delta t\mathcal{B}(l\Delta t)}x\Big{|}\Big{|}_{X}$
$\displaystyle\leq$ $\displaystyle\sum_{j=0}^{M-1}C_{1}(\Delta
t)^{2-\gamma}\big{|}\big{|}\prod_{l=0}^{j-1}e^{\Delta t\mathcal{A}(l\Delta
t)}e^{\Delta t\mathcal{B}(l\Delta
t)}x\big{|}\big{|}_{X}\leq\sum_{j=0}^{M-1}C_{1}(\Delta
t)^{2-\gamma}||x||_{X}=C_{1}(\Delta t)^{1-\gamma}||x||_{X}.$ (113)

Thus, we finish the proof. ∎

## References

  * [1] B. Audoly, H. Berestycki, and Y. Pomeau, Réaction diffusion en écoulement stationnaire rapide, Comptes Rendus de l’Académie des Sciences-Series IIB-Mechanics-Physics-Astronomy, 328 (2000), pp. 255–262.
  * [2] A. Bátkai, P. Csomós, B. Farkas, and G. Nickel, Operator splitting for non-autonomous evolution equations, Journal of Functional Analysis, 260 (2011), pp. 2163–2190.
  * [3] H. Berestycki, F. Hamel, and N. Nadirashvili, The speed of propagation for KPP type problems. I: Periodic framework, Journal of The European Mathematical Society, 7 (2005), pp. 173–213.
  * [4] N. Brummell, F. Cattaneo, and S. Tobias, Linear and nonlinear dynamo properties of time-dependent ABC flows, Fluid Dynamics Research, 28 (2001), p. 237.
  * [5] C. Chicone and Y. Latushkin, Evolution semigroups in dynamical systems and differential equations, no. 70, American Mathematical Soc., 1999\.
  * [6] S. Childress, Alpha-effect in flux ropes and sheets, Physics of the Earth and Planetary Interiors, 20 (1979), pp. 172–180.
  * [7] S. Childress and A. D. Gilbert, Stretch, twist, fold: the fast dynamo, vol. 37, Springer Science & Business Media, 1995.
  * [8] P. Del Moral, Feynman-Kac formulae, in Feynman-Kac Formulae, Springer, 2004, pp. 47–93.
  * [9] F. Den Hollander, Large deviations, vol. 14, American Mathematical Soc., 2008.
  * [10] T. Dombre, U. Frisch, J. M. Greene, M. Henon, A. Mehr, and M. Soward, Chaotic streamlines in the ABC flows, J. Fluid Mech., 167 (1986), pp. 353–391.
  * [11] K. Engel and R. Nagel, One-parameter semigroups for linear evolution equations, vol. 194, Springer Science & Business Media, 1999.
  * [12] G. Ferré and G. Stoltz, Error estimates on ergodic properties of discretized Feynman–Kac semigroups, Numerische Mathematik, 143 (2019), pp. 261–313.
  * [13] R. Fisher, The wave of advance of advantageous genes, Annals of eugenics, 7 (1937), pp. 355–369.
  * [14] W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Quantum Monte Carlo simulations of solids, Reviews of Modern Physics, 73 (2001), p. 33.
  * [15] D. Galloway and M. Proctor, Numerical calculations of fast dynamos in smooth velocity fields with realistic diffusion, Nature, 356 (1992), p. 691.
  * [16] J. Gärtner and M. Freidlin, On the propagation of concentration waves in periodic and random media, in Doklady Akademii Nauk, vol. 249, Russian Academy of Sciences, 1979, pp. 521–525.
  * [17] P. Hess, Periodic-parabolic boundary value problems and positivity, Longman, 1991.
  * [18] M. Hochbruck and C. Lubich, On Magnus integrators for time-dependent Schrödinger equations, SIAM Journal on Numerical Analysis, 41 (2003), pp. 945–963.
  * [19] M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numerica, 19 (2010), pp. 209–286.
  * [20] T. Jahnke and C. Lubich, Error bounds for exponential operator splittings, BIT Numerical Mathematics, 40 (2000), pp. 735–744.
  * [21] T. Kato, Perturbation theory for linear operators, vol. 132, Springer Science & Business Media, 2013.
  * [22] A. Kolmogorov, I. Petrovsky, and N. Piskunov, Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problem, Bull. Moscow State Univ. Ser. A: Math. Mech, 1 (1937), pp. 1–25.
  * [23] J. Lyu, Z. Wang, J. Xin, and Z. Zhang, Convergence analysis of stochastic structure-preserving schemes for computing effective diffusivity in random flows, SIAM Journal on Numerical Analysis, 58 (2020), pp. 3040–3067.
  * [24] A. Majda and P. Souganidis, Large scale front dynamics for turbulent reaction-diffusion equations with separated velocity scales, Nonlinearity, 7 (1994), p. 1.
  * [25] S. Meyn and R. L. Tweedie, Stochastic Stability of Markov Chains, Springer, New York, 1992.
  * [26] G. Milstein, G. John, and S. Vladimir, Transition density estimation for stochastic differential equations via forward-reverse representations, Bernoulli, 10, pp. 281–312.
  * [27] P. D. Moral and A. Guionnet, On the stability of interacting processes with applications to filtering and genetic algorithms, Annales De L Institut Henri Poincare-probabilites Et Statistiques, 37 (2001), pp. 155–194.
  * [28] P. D. Moral and L. Miclo, Branching and interacting particle systems approximations of Feynman-Kac formulae with applications to non-linear filtering, in Seminaire de probabilites XXXIV, Springer, 2000, pp. 1–145.
  * [29] J. Nolen, J. Roquejoffre, L. Ryzhik, and A. Zlatoš, Existence and non-existence of Fisher-KPP transition fronts, Archive for Rational Mechanics and Analysis, 203 (2012), pp. 217–246.
  * [30] J. Nolen, M. Rudd, and J. Xin, Existence of KPP fronts in spatially-temporally periodic advection and variational principle for propagation speeds, Dynamics of PDEs, 2 (2005), pp. 1–24.
  * [31] J. Nolen and J. Xin, Reaction-diffusion front speeds in spatially-temporally periodic shear flows, Multiscale Modeling & Simulation, 1 (2003), pp. 554–570.
  * [32] , Computing reactive front speeds in random flows by variational principle, Physica D: Nonlinear Phenomena, 237 (2008), pp. 3172–3177.
  * [33] , Asymptotic spreading of KPP reactive fronts in incompressible space-time random flows, Ann Inst. H. Poincare, Analyse Non Lineaire, 26 (2009), pp. 815–839.
  * [34] A. Novikov and L. Ryzhik, Boundary layers and KPP fronts in a cellular flow, Archive for rational mechanics and analysis, 184 (2007), pp. 23–48.
  * [35] A. Pazy, Semigroups of linear operators and applications to Partial Differential Equations, vol. 44, Springer Science & Business Media, 2012.
  * [36] L. Ryzhik and A. Zlatoš, KPP pulsating front speed-up by flows, Communications in Mathematical Sciences, 5 (2007), pp. 575–593.
  * [37] R. Schnaubelt, Well-posedness and asymptotic behaviour of non-autonomous linear evolution equations, in Evolution equations, semigroups and functional analysis, Springer, 2002, pp. 311–338.
  * [38] J. Shen, T. Tang, and L. Wang, Spectral methods: algorithms, analysis and applications, vol. 41, Springer Science & Business Media, 2011\.
  * [39] L. Shen, J. Xin, and A. Zhou, Finite element computation of KPP front speeds in 3D cellular and ABC flows, Mathematical Modelling of Natural Phenomena, 8 (2013), pp. 182–197.
  * [40] , Finite element computation of KPP front speeds in cellular and cat’s eye flows, Journal of Scientific Computing, 55 (2013), pp. 455–470.
  * [41] H. Stewart, Generation of analytic semigroups by strongly elliptic operators, Transactions of the American Mathematical Society, 199 (1974), pp. 141–162.
  * [42] 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.
  * [43] , Computing effective diffusivity of chaotic and stochastic flows using structure-preserving schemes, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2322–2344.
  * [44] J. Xin, Existence of planar flame fronts in convective-diffusive periodic media, Archive for rational mechanics and analysis, 121 (1992), pp. 205–233.
  * [45] , Front propagation in heterogeneous media, SIAM review, 42 (2000), pp. 161–230.
  * [46] , An introduction to fronts in random media, vol. 5, Springer Science & Business Media, 2009.
  * [47] P. Zu, L. Chen, and J. Xin, A computational study of residual KPP front speeds in time-periodic cellular flows in the small diffusion limit, Physica D: Nonlinear Phenomena, 311 (2015), pp. 37–44.