# Spectral measures and iterative bounds for effective diffusivity
of steady and space-time periodic flows

N. B. Murphy, D. Hallman, E. Cherkaev, J. Xin, and K. M. Golden

###### Abstract

Over three decades ago the advection-diffusion equation for a steady fluid
velocity field was homogenized, leading to a Stieltjes integral representation
for the effective diffusivity, which is given in terms of a spectral measure
of a compact self-adjoint operator and the Péclet number of the fluid flow.
This result was recently extended to space-time periodic flows, instead
involving an unbounded self-adjoint operator. Padé approximants provide
rigorous upper and lower bounds for Stieltjes functions in terms of the
moments of the spectral measure. However, with the lack of a method for
calculating the moments of the spectral measure for general fluid velocity
fields, the utility of this powerful mathematical framework for calculating
bounds for the effective diffusivity has not been fully realized. Here we
significantly expand the applicability of this framework by providing an
iterative method that enables an arbitrary number of moments, hence bounds, to
be calculated analytically in closed form for both spatially and space-time
periodic flows. The method is demonstrated for periodic flows in two spatial
dimensions. The known asymptotic behavior of the effective diffusivity for a
steady flow is accurately captured by high order upper and lower bounds,
demonstrating the ability of the method to provide accurate estimates for the
effective diffusivity for a broad range of parameter values.

##  1 Introduction

The long time large scale transport of passive tracers by an incompressible
fluid velocity field is equivalent to a diffusion process [24] involving an
effective diffusivity matrix $\mathfrak{D}^{*}$ [12]. A Stieltjes integral
representation for $\mathfrak{D}^{*}$ was developed for steady flows,
involving a spectral measure $\mu$ of a compact self-adjoint operator [1, 2,
16]. Recently this result was extended to space-time periodic flows, involving
a spectral measure for an unbounded self-adjoint operator [15]. This integral
representation _separates_ the Péclet number from the geometry and dynamics of
the fluid velocity field, which is encoded in the spectral measure through its
moments $\mu^{n}$. The theory of Padé approximants for Stieltjes functions
provide a nested sequence of rigorous bounds for the diagonal components of
the matrix $\mathfrak{D}^{*}$, given in terms of the measure moments [3]. The
bounds get tighter as more moments are incorporated and can converge to the
true value of $\mathfrak{D}^{*}$ for certain Péclet number values [3].
However, for over 3 decades, the lack of a way to calculate the measure
moments for general fluid velocity fields has hindered the progress of
providing rigorous bounds for $\mathfrak{D}^{*}$.

Here, we develop an iterative method that, in principle, enables an arbitrary
number of measure moments to be calculated analytically in closed form for any
spatially periodic or space-time periodic fluid velocity field represented as
a finite trigonometric Fourier series. This, in turn, enables an arbitrary
number of nested bounds to be calculated for such flows. This iterative method
is implemented into a numerical algorithm using Maple’s symbolic math toolbox
which can be used to calculate measure moments in closed form for such flows
up to a given order only limited by computational resources. Moreover, we
extend this numerical algorithm to Matlab, which enables hundreds of moments
to be computed using floating point arithmetic. We incorporate the moment
values into an existing numerical algorithm padeapprox [10] which computes
Padé approximants in a robust, stable way and compute several nested bounds
for the diagonal components of the matrix $\mathfrak{D}^{*}$ for some model
steady and space-time periodic flows. High order bounds accurately capture the
known [6, 7] asymptotic behavior of the effective diffusivity for a steady
cell-flow as a function of Péclet number in the advection dominated regime
[15, 16]. Adding a space-time periodic term to the fluid velocity field of
this steady flow results in an appreciable additional enhancement of
$\mathfrak{D}^{*}$, shown both in the Padé approximant bounds here and in
numerical results involving direct computation of the spectral measure $\mu$
[15].

The organization of the paper is as follows. In Section 2 the homogenization
problem for the advection-diffusion equation is reviewed [12, 6, 17, 11]. In
Section 2.1 an abstract Hilbert space framework is reviewed [15] which is used
in Section 2.2 to provide Stieltjes integral representations for the
components $\mathfrak{D}^{*}_{jk}$, $j,k=1,\ldots,d$, of the effective
diffusivity matrix $\mathfrak{D}^{*}$ involving spectral measures $\mu_{jk}$
of a self-adjoint operator [15], where $d$ is the spatial dimension of the
system. This abstract framework is utilized in Section 3 to develop an
iterative method for calculating the moments $\mu_{jk}^{n}$, $n=0,1,2,\ldots$,
of $\mu_{jk}$ for spatially and space-time periodic fluid velocity fields
$\boldsymbol{u}$. The moments for a spatially and a space-time periodic
$\boldsymbol{u}$ are calculated in Section 4, which are incorporated into Padé
approximant bounds for the diagonal components of $\mathfrak{D}^{*}$ [3] in
Section 5. Concluding remarks are given in Section 6.

##  2 Effective transport by advection-diffusion

The density $\phi$ of a cloud of passive tracer particles diffusing along with
molecular diffusivity $\varepsilon$ and being advected by an incompressible
velocity field $\boldsymbol{u}$ satisfies the advection-diffusion equation

$\displaystyle\partial_{t}\phi(t,\boldsymbol{x})=\boldsymbol{u}(t,\boldsymbol{x})\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\phi(t,\boldsymbol{x})+\varepsilon\Delta\phi(t,\boldsymbol{x}),\quad\phi(0,\boldsymbol{x})=\phi_{0}(\boldsymbol{x}),$
(1)

for $t>0$ and $\boldsymbol{x}\in\mathbb{R}^{d}$. Here, the initial density
$\phi_{0}(\boldsymbol{x})$ and the fluid velocity field $\boldsymbol{u}$ are
assumed given, and $\boldsymbol{u}$ satisfies
$\mbox{\boldmath${\nabla}$}\mbox{\boldmath${\cdot}$}\boldsymbol{u}=0$. In
equation (1), $\varepsilon>0$ is the molecular diffusion constant,
$\partial_{t}$ denotes partial differentiation with respect to time $t$, and
$\Delta=\mbox{\boldmath${\nabla}$}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}=\nabla^{2}$
is the Laplacian. Moreover,
$\boldsymbol{\psi}\mbox{\boldmath${\cdot}$}\boldsymbol{\varphi}=\boldsymbol{\psi}^{\,T}\overline{\boldsymbol{\varphi}}$,
where $\boldsymbol{\psi}^{\,T}$ denotes transposition of the vector
$\boldsymbol{\psi}$ and $\overline{\boldsymbol{\varphi}}$ denotes component-
wise complex conjugation, with
$\boldsymbol{\psi}\mbox{\boldmath${\cdot}$}\boldsymbol{\psi}=|\boldsymbol{\psi}|^{2}$.
Later, we will extensively use this form of the dot product over complex
fields, with built in complex conjugation. However, we emphasize that all
quantities considered in this section are _real-valued_.

In our analysis of the effective diffusivity matrix $\mathfrak{D}^{*}$, it is
beneficial to use non-dimensional parameters. We therefore assume that
equation (1) has been non-dimensionalized as follows. Let $\ell$ and $\tau$ be
typical length and time scales associated with the problem of interest.
Mapping to the non-dimensional variables $t\mapsto t/\tau$ and
$\boldsymbol{x}\mapsto\boldsymbol{x}/\ell$, one finds that $\phi$ satisfies
the advection diffusion equation in (1) with a non-dimensional molecular
diffusivity and fluid velocity field,

$\displaystyle\varepsilon\mapsto\tau\varepsilon/\ell^{\,2},\quad\boldsymbol{u}\mapsto\tau\,\boldsymbol{u}/\ell.$
(2)

This non-dimensionalization demonstrates that the fluid velocity field
$\boldsymbol{u}$ is divided by a quantity with dimensions of velocity and the
molecular diffusivity is divided by a quantity with dimensions of velocity
multiplied by spatial length. A detailed discussion of various non-
dimensionalizations involving the Strouhal number, the Péclet number, and the
periodic Péclet number is given in [13, 11]. It is convenient to choose the
rescaled $\boldsymbol{u}$ and $\varepsilon$ in a way that captures information
about the fluid velocity field. However, it is also convenient to choose these
rescaled variables in a way that _separates_ the rescaled $\varepsilon$ from
the _geometry and dynamics_ of $\boldsymbol{u}$; this leads to mathematically
and physically meaningful properties of rigorous bounds for $\mathfrak{D}^{*}$
which follow from the analytic structure of Stieltjes integral representations
for $\mathfrak{D}^{*}$ [2, 3] — discussed in Section 2.2 below.

We accomplish both of these goals as follows. Define the dimensional fluid
velocity field by $\boldsymbol{u}=u_{0}\boldsymbol{v}$, where the parameter
$u_{0}$ has dimensions of velocity and represents the “ _flow strength_ ” of
$\boldsymbol{u}$ which is independent of the geometry and dynamics of
$\boldsymbol{u}$ which, in turn, is encapsulated in the non-dimensional vector
field $\boldsymbol{v}$. With these definitions, we choose reference scales
$\tau$ and $\ell$ in equation (2) to satisfy $u_{0}=\ell/\tau$ so that
$\boldsymbol{u}\mapsto\boldsymbol{v}$ and
$\varepsilon\mapsto\varepsilon/u_{0}\,\ell\,$. For example, for $BC$-flow [5],
we define the dimensional fluid velocity field by
$\boldsymbol{u}=u_{0}\,(C\cos{y},B\cos{x})$, where the flow strength
$u_{0}\in(0,\infty)$ is chosen to be independent of the non-dimensional
parameters $B,C\in[0,1]$ which determine the streamline geometry of
$\boldsymbol{u}$ in $\boldsymbol{v}=(C\cos{y},B\cos{x})$.

An example of a non-dimensional parameter that compares the rate of scalar
advection to the rate of diffusion is the Péclet number. We define it by the
ratio $\mbox{{Pe}}=\ell u_{0}/\varepsilon$, although other definitions have
been used [13, 11, 11]. Therefore, our choice of the rescaled $\varepsilon$
satisfies $\mbox{{Pe}}=1/\varepsilon$. The advection and diffusion dominated
regimes are characterized by $\mbox{{Pe}}\gg 1$ and $\mbox{{Pe}}\ll 1$,
respectively.

The _parameter separation_ between Pe and the geometry of the flow is
important for rigorous upper and lower Padé approximant bounds for
$\mathfrak{D}^{*}$ [2] discussed in Section 5. Padé approximants of
$\mathfrak{D}^{*}$ are given in terms of ratios of polynomials [3]
$P(z)/Q(z)$, where $z=\mbox{{Pe}}^{2}$, $0<z<\infty$, and the coefficients of
these polynomials depend on the moments of a _spectral measure_ that, in turn,
depend on the fluid velocity field $\boldsymbol{u}$ [2]. For example, when
$\boldsymbol{u}$ is given by $BC$-flow the moments of the measure depend on
the parameters $B$ and $C$. Our numerical investigations have shown if the
non-dimensionalization of equation (1) is chosen in a way that the variable
$z$ also depends on the flow geometry through the ratio $B/C$, then this gives
rise to _positive real_ roots for the polynomials $P(z)$ and $Q(z)$. This, in
turn, gives rise to positive real roots and poles in the (rigorous) Padé
approximant bounds for $\mathfrak{D}^{*}$, which is not physically or
mathematically consistent with the known behavior of $\mathfrak{D}^{*}$ [6,
18, 5, 11]. This demonstrates the importance of _parameter separation_ between
$z$ and the flow geometry for Padé approximant bounds for $\mathfrak{D}^{*}$.

This way of non-dimensionalizing equation (1) is also convenient in the case
of a time-dependent fluid velocity field [15], where the parameter $u_{0}$
again represents the flow strength and the vector field $\boldsymbol{v}$
encapsulates the _geometric and dynamical_ properties of the flow. For
example, the space-time periodic flow with velocity field
$\boldsymbol{u}=u_{0}(\,(C\cos{y},B\cos{x})+\theta\cos{t}\,(\sin{y},\sin{x})\,)$
has dynamical behavior exhibiting Lagrangian chaos [5, 15]. Here, the flow
strength $u_{0}\in(0,\infty)$ is independent of the parameters
$B,C,\theta\in[0,1]$ which determine the geometric and dynamical properties of
$\boldsymbol{u}$. This choice of non-dimensionalization gives a clearer
interpretation of the advection and diffusion dominated regimes in terms of
$\mbox{{Pe}}=1/\varepsilon$ than the non-dimensionalization given in [15].

We now discuss the effective transport properties of advection enhanced
diffusion, as described by the advection diffusion equation in (1). We will
assume in this manuscript that the fluid velocity field $\boldsymbol{u}$ is
mean-zero in space for steady $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x})$
and mean-zero in space-time when
$\boldsymbol{u}=\boldsymbol{u}(t,\boldsymbol{x})$ is space-time dependent
(also see [19]). The long time, large scale dispersion of diffusing tracers,
such as heat or pollutants, being advected by an incompressible fluid velocity
field is equivalent to an enhanced diffusion process [25] with an effective
diffusivity matrix $\mathfrak{D}^{*}$. In recent decades, methods of
homogenization theory [12, 6, 17, 11] have been used to provide an explicit
representation for $\mathfrak{D}^{*}$. In particular, these methods have
demonstrated that the averaged or _homogenized_ behavior of the advection-
diffusion equation in (1), with space-time periodic velocity field
$\boldsymbol{u}$, is determined by a diffusion equation involving an averaged
scalar density $\bar{\phi}$ and an effective diffusivity tensor
$\mathfrak{D}^{*}$ [11]

$\displaystyle\partial_{t}\bar{\phi}(t,\boldsymbol{x})=\mbox{\boldmath${\nabla}$}\mbox{\boldmath${\cdot}$}[\mathfrak{D}^{*}\mbox{\boldmath${\nabla}$}\bar{\phi}(t,\boldsymbol{x})],\quad\bar{\phi}(0,\boldsymbol{x})=\phi_{0}(\boldsymbol{x}).$
(3)

Equation (3) follows from the assumption that the initial tracer density
$\phi_{0}$ varies slowly relative to the variations of the fluid velocity
field $\boldsymbol{u}$ [12, 7, 11]. This information is incorporated into
equation (1) by introducing a small dimensionless parameter $\delta\ll 1$ and
writing [12, 7, 11]

$\displaystyle\phi(0,\boldsymbol{x})=\phi_{0}(\delta\boldsymbol{x})\,.$ (4)

Anticipating that $\phi$ will have diffusive dynamics as $t\to\infty$, space
and time are rescaled according to the standard diffusive relation

$\displaystyle\boldsymbol{\xi}=\boldsymbol{x}/\delta,\quad\tau=t/\delta^{2}.$
(5)

The rescaled form of equation (1) is given by [11]

$\displaystyle\partial_{t}\phi^{\delta}(t,\boldsymbol{x})=\delta^{-1}\boldsymbol{u}(t/\delta^{2},\boldsymbol{x}/\delta)\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\phi^{\delta}(t,\boldsymbol{x})+\varepsilon\Delta\phi^{\delta}(t,\boldsymbol{x}),\quad\phi^{\delta}(0,\boldsymbol{x})=\phi_{0}(\boldsymbol{x}),$
(6)

where we have denoted
$\phi^{\delta}(t,\boldsymbol{x})=\phi(t/\delta^{2},\boldsymbol{x}/\delta)$.
The convergence of $\phi^{\delta}$ to $\bar{\phi}$ can be rigorously
established in the following sense [11]

$\displaystyle\lim_{\delta\to 0}\;\sup_{0\leq t\leq
t_{0}}\,\sup_{\boldsymbol{x}\in\mathbb{R}^{d}}|\phi^{\delta}(t,\boldsymbol{x})-\bar{\phi}(t,\boldsymbol{x})|=0,$
(7)

for every finite $t_{0}>0$, provided that $\phi_{0}$ and $\boldsymbol{u}$ obey
some mild smoothness and boundedness conditions.

An explicit representation of the effective diffusivity tensor
$\mathfrak{D}^{*}$ is given in terms of the (unique) mean zero, space-time
periodic solution $\chi_{j}$ of the following _cell problem_ [5, 11],

$\displaystyle\partial_{\tau}\chi_{j}(\tau,\boldsymbol{\xi})-\varepsilon\Delta_{\xi}\chi_{j}(\tau,\boldsymbol{\xi})-\boldsymbol{u}(\tau,\boldsymbol{\xi})\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}_{\xi}\chi_{j}(\tau,\boldsymbol{\xi})=u_{j}(\tau,\boldsymbol{\xi}),$
(8)

where the subscript $\xi$ in $\Delta_{\xi}$ and
$\mbox{\boldmath${\nabla}$}_{\xi}$ indicates that differentiation is with
respect to the fast variable $\boldsymbol{\xi}$ defined in equation (5). The
components $\mathfrak{D}^{*}_{jk}$, $j,k=1,\ldots,d$, of the matrix
$\mathfrak{D}^{*}$ are given by [12, 6, 17, 11]

$\displaystyle\mathfrak{D}^{*}_{jk}=\varepsilon\delta_{jk}+\langle
u_{j}\chi_{k}\rangle,$ (9)

where $\delta_{jk}$ is the Kronecker delta and $u_{j}$ is the $j$th component
of the vector $\boldsymbol{u}$. The averaging $\langle\cdot\rangle$ in (9) is
with respect to the fast variables defined in equation (5). The averaging is
over the bounded sets $\mathcal{T}\subset\mathbb{R}$ and
$\mathcal{V}\subset\mathbb{R}^{d}$, with $\tau\in\mathcal{T}$ and
$\boldsymbol{\xi}\in\mathcal{V}$, which define the space-time period cell
($(d+1)$–torus) $\mathcal{T}\times\mathcal{V}$.

In general, the effective diffusivity tensor $\mathfrak{D}^{*}$ has a
symmetric $\mathfrak{S}^{*}$ and antisymmetric $\mathfrak{A}^{*}$ part defined
by

$\displaystyle\mathfrak{D}^{*}=\mathfrak{S}^{*}+\mathfrak{A}^{*},\qquad\mathfrak{S}^{*}=\frac{1}{2}\left(\mathfrak{D}^{*}+[\mathfrak{D}^{*}]^{\,T}\right),\quad\mathfrak{A}^{*}=\frac{1}{2}\left(\mathfrak{D}^{*}-[\mathfrak{D}^{*}]^{\,T}\right),$
(10)

where $[\mathfrak{D}^{*}]^{\,T}$ denotes transposition of the matrix
$\mathfrak{D}^{*}$. Denote by $\mathfrak{S}^{*}_{jk}$ and
$\mathfrak{A}^{*}_{jk}$, $j,k=1,\ldots,d$, the components of
$\mathfrak{S}^{*}$ and $\mathfrak{A}^{*}$ in (10). When the fluid velocity
field is mean-zero and divergence-free, as discussed above, then equation (7)
holds and the effective diffusivity tensor $\mathfrak{D}^{*}$ defined in (9)
is constant [11]. Consequently, only the symmetric part of $\mathfrak{D}^{*}$
plays a role in the effective transport equation shown in (3), as the
antisymmetric part of $\mathfrak{D}^{*}$ cancels out in the sum
$\sum_{ij}\mathfrak{D}^{*}_{ij}\partial_{i}\partial_{j}\bar{\phi}$, where
$\partial_{i}$ denotes differentiation in the $i$th spatial direction [18].

In Section 4.3 we consider the fluid velocity field $\boldsymbol{u}$

$\displaystyle\boldsymbol{u}(t,\boldsymbol{x})=(C\cos{y},B\cos{x})+\theta\,\cos{t}\;(\sin{y},\sin{x}),\quad\theta\in[0,1].$
(11)

with temporal periodicity $\mathcal{T}=[0,2\pi]$ and spatial periodicity
$\mathcal{V}=[0,2\pi]^{d}$, with $d=2$. In the case of a time-dependent fluid
velocity field, $\langle\cdot\rangle$ denotes space-time averaging over
$\mathcal{T}\times\mathcal{V}$. In the special case of a time-independent
fluid velocity field, the function $\chi_{j}$ is time-independent and
satisfies equation (8) with $\partial_{\tau}\chi_{j}\equiv 0$, and
$\langle\cdot\rangle$ in (9) denotes spatial averaging over $\mathcal{V}$ [6,
17, 11].

###  2.1 Hilbert space

In this section we provide an abstract Hilbert space formulation of the
effective parameter problem for advection-diffusion that was proposed in [18],
based on [4], and generalized to the setting of a space-time periodic fluid
velocity field in [15]. To fix ideas, consider the following sets
$\mathcal{T}=[0,T]$ and $\mathcal{V}=\otimes_{j=1}^{d}[0,L]$ which define the
space-time period cell $\mathcal{T}\times\mathcal{V}$. Now consider the
Hilbert spaces $L^{2}(\mathcal{T})$ and $L^{2}(\mathcal{V})$ of Lebesgue
measurable scalar functions over the complex field $\mathbb{C}$ that are also
square integrable [9]. Define the associated Hilbert spaces
$\mathscr{H}_{\mathcal{T}}$, $\mathscr{H}_{\mathcal{V}}$, and
$\mathscr{H}_{\mathcal{T}\mathcal{V}}=\mathscr{H}_{\mathcal{T}}\otimes\mathscr{H}_{\mathcal{V}}$
of periodic functions, where

$\displaystyle\mathscr{H}_{\mathcal{T}}$ $\displaystyle=\big{\\{}\psi\in
L^{2}(\mathcal{T})\,|\,\psi(t)=\psi(t+T)\big{\\}},$ (12)
$\displaystyle\mathscr{H}_{\mathcal{V}}$ $\displaystyle=\big{\\{}\psi\in
L^{2}(\mathcal{V})\,|\,\psi(\boldsymbol{x})=\psi(\boldsymbol{x}+L\boldsymbol{e}_{j}),\
j=1,\ldots,d\big{\\}},$

and the $\boldsymbol{e}_{j}$ are standard basis vectors.

More specifically, denote time average over $\mathcal{T}$ by
$\langle\cdot\rangle_{\mathcal{T}}$, space average over $\mathcal{V}$ by
$\langle\cdot\rangle_{\mathcal{V}}$, and space-time average over
$\mathcal{T}\times\mathcal{V}$ by $\langle\cdot\rangle$. The space-time
average $\langle\cdot\rangle$, induces a sesquilinear inner-product
$\langle\cdot,\cdot\rangle$ given by
$\langle\psi,\varphi\rangle=\langle\psi\;\overline{\varphi}\rangle$, with
$\langle\varphi,\psi\rangle=\overline{\langle\psi,\varphi\rangle}$. This
$\mathscr{H}_{\mathcal{T}\mathcal{V}}$–inner-product, in turn, induces a norm
$\|\cdot\|$ given by $\|\psi\|=\langle\psi,\psi\rangle^{1/2}$ [9]. The set of
space-time periodic Lebesgue measurable functions
$\mathscr{H}_{\mathcal{T}\mathcal{V}}$ satisfying $\|f\|<\infty$ is a
(complete) Hilbert space [9]. Similarly, the space and time averages,
$\langle\cdot\rangle_{\mathcal{V}}$ and $\langle\cdot\rangle_{\mathcal{T}}$,
induce sesquilinear inner-products, $\langle\cdot,\cdot\rangle_{\mathcal{V}}$
and $\langle\cdot,\cdot\rangle_{\mathcal{T}}$, that induce norms,
$\|\cdot\|_{\mathcal{V}}$ and $\|\cdot\|_{\mathcal{T}}$, associated with the
Hilbert spaces $\mathscr{H}_{\mathcal{V}}$ and $\mathscr{H}_{\mathcal{T}}$.

To treat temporal dependence, we define the space $\mathscr{A}_{\mathcal{T}}$
of functions that are absolutely continuous [23, 21] on the interval
$\mathcal{T}$, having derivative belonging to $L^{2}(\mathcal{T})$, and the
space $\tilde{\mathscr{A}}_{\mathcal{T}}$ of absolutely continuous
$\mathcal{T}$–periodic functions with time derivatives belonging to
$L^{2}(\mathcal{T})$,

$\displaystyle\tilde{\mathscr{A}}_{\mathcal{T}}$
$\displaystyle=\\{\psi\in\mathscr{A}_{\mathcal{T}}\,|\,\psi(0)=\psi(T)\\},$
(13)

which is _not_ a Hilbert space but is instead an everywhere dense subset of
the Hilbert space $\mathscr{H}_{\mathcal{T}}$ [23]. To treat spatial
dependence, we now define the Sobolev space $\mathscr{H}^{1,2}_{\mathcal{V}}$
which is itself a Hilbert space [4, 8, 14],

$\displaystyle\mathscr{H}^{1,2}_{\mathcal{V}}=\big{\\{}\psi\in\mathscr{H}_{\mathcal{V}}\
|\ \|\mbox{\boldmath${\nabla}$}\psi\|_{\mathcal{V}}<\infty,\
\langle\psi\rangle_{\mathcal{V}}=0\big{\\}}.$ (14)

The condition $\langle\psi\rangle_{\mathcal{V}}=0$ in (14) is required to
eliminate non-zero constant $\psi$, which satisfies
$\|\mbox{\boldmath${\nabla}$}\psi\|_{\mathcal{V}}=0$. The
$\mathscr{H}^{1,2}_{\mathcal{V}}$–norm
$\|\mbox{\boldmath${\nabla}$}\cdot\|_{\mathcal{V}}$ is induced by the
$\mathscr{H}^{1,2}_{\mathcal{V}}$–inner-product:
$\|\mbox{\boldmath${\nabla}$}\psi\|_{\mathcal{V}}=\langle\mbox{\boldmath${\nabla}$}\psi\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\psi\rangle_{\mathcal{V}}^{1/2}$.

Finally, define the Hilbert space $\mathscr{H}$ and its everywhere dense
subset $\mathscr{F}$

$\displaystyle\mathscr{H}=\mathscr{H}_{\mathcal{T}}\otimes\mathscr{H}^{1,2}_{\mathcal{V}},\qquad\mathscr{F}=\tilde{\mathscr{A}}_{\mathcal{T}}\otimes\mathscr{H}^{1,2}_{\mathcal{V}}.$
(15)

Due to the presence of $\tilde{\mathscr{A}}_{\mathcal{T}}$ in the definition
of the function space $\mathscr{F}$, it is _not_ a complete Hilbert space, and
is instead an everywhere dense subset of the complete Hilbert space
$\mathscr{H}$. Recall that $\langle\cdot\rangle$ denotes space-time average
over $\mathcal{T}\times\mathcal{V}$ and
$\boldsymbol{\psi}\mbox{\boldmath${\cdot}$}\boldsymbol{\zeta}=\boldsymbol{\psi}^{T}\,\overline{\boldsymbol{\zeta}}$.
The sesquilinear $\mathscr{H}$–inner-product is given by
$\langle\psi,\varphi\rangle_{1,2}=\left\langle\mbox{\boldmath${\nabla}$}\psi\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\varphi\right\rangle$
with associated norm $\|\cdot\|_{1,2}$ given by
$\|\psi\|_{1,2}=\langle|\mbox{\boldmath${\nabla}$}\psi|^{2}\rangle^{1/2}$. We
emphasize that in the case of a time-dependent fluid velocity field, it is
necessary that $\psi\in\mathscr{H}$ satisfy
$\langle\psi\rangle_{\mathcal{V}}=0$, as required by the definition of
$\mathscr{H}^{1,2}_{\mathcal{V}}$ in (14). Otherwise,
$\|\cdot\|_{1,2}=|\mathcal{T}\times\mathcal{V}|^{-1}\int_{\mathcal{T}\times\mathcal{V}}\mathrm{d}t\,\mathrm{d}\boldsymbol{x}\,|\mbox{\boldmath${\nabla}$}\cdot|^{2}$
is not a norm, since a strictly positive function
$\psi(t,\boldsymbol{x})=\psi(t)$ on $\mathcal{T}\times\mathcal{V}$ satisfies
$\|\psi\|_{1,2}=0$, where $|\mathcal{T}\times\mathcal{V}|$ denotes Lebesgue
measure of the set $\mathcal{T}\times\mathcal{V}$. In the case of a time-
independent fluid velocity field
$\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x})$ we set
$\mathscr{H}\equiv\mathscr{F}\equiv\mathscr{H}^{1,2}_{\mathcal{V}}$, and
$\langle\cdot\rangle=\langle\cdot\rangle_{\mathcal{V}}$.

###  2.2 Integral representations for the effective diffusivity

In this section we summarize the results of [15], which provides Stieltjes
integral representations for both the symmetric $\mathfrak{S}^{*}$ and
antisymmetric $\mathfrak{A}^{*}$ parts of $\mathfrak{D}^{*}$. Since the
analysis in this section involves only the fast variables
$(\tau,\boldsymbol{\xi})$ defined in equation (5), for notational simplicity,
we will drop the subscripts $\xi$ shown in equation (8) and use $\partial_{t}$
to denote $\partial_{\tau}$.

Inserting the expression for $u_{j}$ on the right side of (8) into equation
(9) leads to the following functional representations for the components
$\mathfrak{S}^{*}_{jk}$ and $\mathfrak{A}^{*}_{jk}$, $j,k=1,\ldots,d$, of
$\mathfrak{S}^{*}$ and $\mathfrak{A}^{*}$ [18]

$\displaystyle\mathfrak{S}^{*}_{jk}=\varepsilon(\delta_{jk}+\langle\chi_{j},\chi_{k}\rangle_{1,2}),\quad\mathfrak{A}^{*}_{jk}=\langle
A\chi_{j},\chi_{k}\rangle_{1,2}\,,\quad
A=(-\Delta)^{-1}(\partial_{t}-\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$})\,.$
(16)

Here, $\langle
f,h\rangle_{1,2}=\langle\mbox{\boldmath${\nabla}$}f\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}h\rangle$
is a Sobolev-type _sesquilinear_ inner-product [14] and the operator
$(-\Delta)^{-1}$ is based on convolution with respect to the Green’s function
for the Laplacian $\Delta$ [22, 8]. Since the function $\chi_{j}$ is _real-
valued_ we have
$\langle\chi_{j},\chi_{k}\rangle_{1,2}=\langle\chi_{k},\chi_{j}\rangle_{1,2}$,
which implies that $\mathfrak{S}^{*}$ is a symmetric matrix. The function
$A\chi_{j}$ is also real-valued. The operator $A$ is skew-adjoint on the
Hilbert space $\mathscr{H}$ [15], which implies that
$\mathfrak{A}^{*}_{kj}=\langle
A\chi_{k},\chi_{j}\rangle_{1,2}=-\langle\chi_{k},A\chi_{j}\rangle_{1,2}=-\langle
A\chi_{j},\chi_{k}\rangle_{1,2}=-\mathfrak{A}^{*}_{jk}$ which, in turn,
implies that $\mathfrak{A}^{*}$ is an antisymmetric matrix, hence
$\mathfrak{A}^{*}_{kk}=\langle A\chi_{k},\chi_{k}\rangle_{1,2}=0$.

Applying the linear operator $(-\Delta)^{-1}$ to both sides of the cell
problem in equation (8) yields the following resolvent formula for $\chi_{j}$

$\displaystyle\chi_{j}=(\varepsilon+A)^{-1}g_{j},\qquad
g_{j}=(-\Delta)^{-1}u_{j}.$ (17)

From equations (16) and (17) we have the following functional formulas for
$\mathfrak{S}^{*}_{jk}$ and $\mathfrak{A}^{*}_{jk}$ involving the skew-adjoint
operator $A$

$\displaystyle\mathfrak{S}^{*}_{jk}$
$\displaystyle=\varepsilon\left(\delta_{jk}+\langle(\varepsilon+A)^{-1}g_{j},(\varepsilon+A)^{-1}g_{k}\rangle_{1,2}\right),$
(18) $\displaystyle\mathfrak{A}^{*}_{jk}$ $\displaystyle=\langle
A(\varepsilon+A)^{-1}g_{j},(\varepsilon+A)^{-1}g_{k}\rangle_{1,2}.$

Since $A$ is a skew-adjoint operator, it can be written as $A=\imath M$ where
$M$ is a symmetric operator [23]. In [15] it is shown that $M$ is _self-
adjoint_ on the Hilbert space $\mathscr{H}$.

The spectral theorem for self-adjoint operators states that there is a one-to-
one correspondence between the self-adjoint operator $M$ and a family of self-
adjoint projection operators $\\{Q(\lambda)\\}_{\lambda\in\Sigma}$ — the
resolution of the identity — that satisfies
$\lim_{\lambda\to\,\inf{\Sigma}}Q(\lambda)=0$ and
$\lim_{\lambda\to\,\sup{\Sigma}}Q(\lambda)=I$ [23]. Here, $\Sigma$ is the
_spectrum_ of the operator $M$, while $0$ and $I$ denote the null and identity
operators. Define the _complex valued_ function $\mu_{jk}(\lambda)=\langle
Q(\lambda)g_{j},g_{k}\rangle_{1,2}$, $j,k=1,\ldots,d$, where
$g_{j}=(-\Delta)^{-1}u_{j}$ is defined in (17). The real,
$\mbox{Re}\,\mu_{jk}(\lambda)$, and imaginary, $\mbox{Im}\,\mu_{jk}(\lambda)$,
parts of the function $\mu_{jk}(\lambda)$ are of bounded variation, and
therefore have Stieltjes measures $\mbox{Re}\,\mu_{jk}$ and
$\mbox{Im}\,\mu_{jk}$ associated with them [23]. The function
$\mu_{kk}(\lambda)$ is positive hence $\mu_{kk}$ is a positive measure, while
$\mbox{Re}\,\mu_{jk}$ and $\mbox{Im}\,\mu_{jk}$, $j\neq k$, are signed
measures. Given certain regularity conditions [15] on the components $u_{j}$
of the fluid velocity field $\boldsymbol{u}$, the functional formulas for
$\mathfrak{S}^{*}_{jk}$ and $\mathfrak{A}^{*}_{jk}$ in (18) have the following
Radon–Stieltjes integral representations, for all $0<\varepsilon<\infty$

$\displaystyle\mathfrak{S}^{*}_{jk}=\varepsilon\left(\delta_{jk}+\int_{-\infty}^{\infty}\frac{\mathrm{d}\mbox{Re}\,\,\mu_{jk}(\lambda)}{\varepsilon^{2}+\lambda^{2}}\right),\qquad\mathfrak{A}^{*}_{jk}=-\int_{-\infty}^{\infty}\frac{\lambda\,\mathrm{d}\mbox{Im}\,\,\mu_{jk}(\lambda)}{\varepsilon^{2}+\lambda^{2}}\,.$
(19)

The integration in (19) is over the spectrum $\Sigma\subseteq\mathbb{R}$ of
the self-adjoint operator $M=-\imath A$ [23, 20]. In the setting of a time-
independent flow, $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x})$, the
operator
$A=(-\Delta)^{-1}[\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}]$
and the self-adjoint operator $M=-\imath A$ is compact [4]. Therefore, the
spectrum $\Sigma$ is discrete outside a neighborhood of $\lambda=0$ with a
limit point at $\lambda=0$ [22]. In the setting of a time-dependent flow,
$\boldsymbol{u}=\boldsymbol{u}(t,\boldsymbol{x})$, $M=-\imath A$ is an
unbounded operator [15, 23]. Therefore, in general, the spectrum $\Sigma$ can
be an unbounded subset of $\mathbb{R}$, it can have discrete and continuous
components, and can even coincide with $\mathbb{R}$ itself [23].

##  3 Iterative moment method

In this section we provide an iterative method which may be used to calculate,
in principle, an arbitrary number of moments for spectral measures associated
with the effective diffusivity for spatially and space-time periodic fluid
velocity fields. The spectral theorem shows that the mass $\mu^{0}_{jk}$ and
the moments $\mu^{n}_{jk}$, $n=1,2,3,\ldots$, of the spectral measure
$\mu_{jk}$ are given by

$\displaystyle\mu^{0}_{jk}=\int_{-\infty}^{\infty}\mathrm{d}\mu_{jk}(\lambda)=\langle
g_{j},g_{k}\rangle_{1,2},\qquad\mu^{n}_{jk}=\int_{-\infty}^{\infty}\lambda^{n}\;\mathrm{d}\mu_{jk}(\lambda)=\langle
M^{n}g_{j},g_{k}\rangle_{1,2},\quad n=1,2,3,\ldots\;,$

where $M=-\imath A$ and the operator $M^{n}$ is defined through composition.
Integration by parts and properties of the operator $(-\Delta)^{-1}$ show that
the mass $\mu^{0}_{jk}$ of the measure $\mu_{jk}$ is given by [15]

$\displaystyle\mu^{0}_{jk}=\langle
g_{j},g_{k}\rangle_{1,2}=\langle\mbox{\boldmath${\nabla}$}(-\Delta)^{-1}u_{j}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}(-\Delta)^{-1}u_{k}\rangle=\langle(-\Delta)^{-1}u_{j},u_{k}\rangle_{2}=\langle
g_{j},u_{k}\rangle_{2},$ (20)

where we have denoted $\langle\cdot,\cdot\rangle_{2}$ the sesquilinear
$L^{2}(\mathcal{T}\times\mathcal{V})$ inner-product. Similarly, denoting the
_material derivative_
$D_{t}=\partial_{t}+\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}$
for time-dependent $\boldsymbol{u}$ and
$D_{t}=\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}$ for
time-independent $\boldsymbol{u}$, with $A=(-\Delta)^{-1}D_{t}$, integration
by parts and the anti-symmetry of $A$ yield [15]

$\displaystyle\mu^{1}_{jk}$ $\displaystyle=\langle
Mg_{j},g_{k}\rangle_{1,2}=-\imath\,\langle D_{t}g_{j},g_{k}\rangle_{2}\,,$
(21) $\displaystyle\mu^{2}_{jk}$ $\displaystyle=\langle
M^{2}g_{j},g_{k}\rangle_{1,2}=\langle Ag_{j},Ag_{k}\rangle_{1,2}=\langle
D_{t}g_{j},Ag_{k}\rangle_{2}\,.$

Higher moments are found in a similar way, for $n=1,2,3,\cdots$,

$\displaystyle\mu^{2n+1}_{jk}$ $\displaystyle=\langle
M^{2n+1}g_{j},g_{k}\rangle_{1,2}=-\imath\,\langle
A^{n+1}g_{j},A^{n}g_{k}\rangle_{1,2}=-\imath\,\langle
D_{t}A^{n}g_{j},A^{n}g_{k}\rangle_{2}\,,$ (22) $\displaystyle\mu^{2n+2}_{jk}$
$\displaystyle=\langle M^{2n+2}g_{j},g_{k}\rangle_{1,2}=\langle
A^{n+1}g_{j},A^{n+1}g_{k}\rangle_{1,2}=\langle
D_{t}A^{n}g_{j},A^{n+1}g_{k}\rangle_{2}\,.$

From equation (22) and the asymmetry of $A$, we have $\mu^{2n+1}_{jk}=0$ for
all $n=0,1,2,\ldots$ [15]. Moreover, $\mu^{2n}_{kj}=\mu^{2n}_{jk}$ as $g_{k}$,
$D_{t}g_{k}$, $A^{n}g_{k}$, and $D_{t}A^{n}g_{k}$ are real-valued functions
for all $n=1,2,\ldots$. In summary,

$\displaystyle\mu^{2n}_{kj}=\mu^{2n}_{jk}\,,\qquad\mu^{2n+1}_{jk}=0\,.$ (23)

A key property of the operators $(-\Delta)^{-1}$, $D_{t}$, and
$A=(-\Delta)^{-1}D_{t}$ is that they are _linear_. In order to take advantage
of this property in our calculation of the moments $\mu^{n}_{jk}$ of the
measure $\mu_{jk}$, we take $\mathcal{T}=[0,2\pi]$ and
$\mathcal{V}=[0,2\pi]^{2}$, and write the fluid velocity field
$\boldsymbol{u}$ in the complex, orthonormal Fourier basis for
$L^{2}(\mathcal{T}\times\mathcal{V})$ [9]

$\displaystyle\boldsymbol{u}=\sum_{\ell,\boldsymbol{k}}\boldsymbol{a}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}},\qquad\phi_{\ell,\boldsymbol{k}}(t,\boldsymbol{x})=\exp[\imath(\ell
t+\boldsymbol{k}\mbox{\boldmath${\cdot}$}\boldsymbol{x})]\,,\quad\boldsymbol{a}_{\ell,\boldsymbol{k}}=\langle\boldsymbol{u},\phi_{\ell,\boldsymbol{k}}\rangle_{2}\,,$
(24)

where $\ell\in\mathbb{Z}$, $\boldsymbol{k}\in\mathbb{Z}^{d}$,
$\boldsymbol{a}_{\ell,\boldsymbol{k}}=(a^{1}_{\ell,\boldsymbol{k}},\ldots,a^{d}_{\ell,\boldsymbol{k}})$,
and the average $\langle\boldsymbol{u},\phi_{\ell,\boldsymbol{k}}\rangle_{2}$
is understood to be component-wise, so
$(\boldsymbol{a}_{\ell,\boldsymbol{k}})_{i}=\langle
u_{i},\phi_{\ell,\boldsymbol{k}}\rangle_{2}$. The orthonormal basis vectors
satisfy

$\displaystyle\langle\phi_{\ell,\boldsymbol{k}},\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}\rangle_{2}=\langle\phi_{\ell-\ell^{\prime},\boldsymbol{k}-\boldsymbol{k}^{\prime}}\rangle=\delta_{\ell,\ell^{\prime}}\,\delta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}\,,$
(25)

where $\langle\cdot\rangle$ denotes averaging over
$\mathcal{T}\times\mathcal{V}$ and $\delta_{i,j}$ is the Kronecker delta, with
$\delta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}=\prod_{i}\delta_{k_{i},k_{i}^{\prime}}$.
For simplicity, we assume that there are a _finite_ number of terms in the
Fourier expansion of $\boldsymbol{u}$. In the setting of a time-independent
flow, $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x})$, the basis function
$\phi_{\ell,\boldsymbol{k}}$ is replaced by
$\phi_{\boldsymbol{k}}(\boldsymbol{x})=\exp(\imath\,\boldsymbol{k}\mbox{\boldmath${\cdot}$}\boldsymbol{x})$,
so that
$\boldsymbol{a}_{\ell,\boldsymbol{k}}=\boldsymbol{a}_{0,\boldsymbol{k}}=\boldsymbol{a}_{\boldsymbol{k}}$.
The basis functions $\\{\phi_{\ell,\boldsymbol{k}}\\}_{\ell,\boldsymbol{k}}$
are _eigenfunctions_ of the operators $\partial_{t}$, ${\nabla}$, and
$(-\Delta)^{-1}$, with [15]

$\displaystyle\partial_{t}\phi_{\ell,\boldsymbol{k}}=\imath\ell\,\phi_{\ell,\boldsymbol{k}}\,,\qquad\mbox{\boldmath${\nabla}$}\phi_{\ell,\boldsymbol{k}}=\imath\boldsymbol{k}\,\phi_{\ell,\boldsymbol{k}}\,,\qquad(-\Delta)^{-1}\phi_{\ell,\boldsymbol{k}}=|\boldsymbol{k}|^{-2}\phi_{\ell,\boldsymbol{k}}\,.$
(26)

It follows from equations (24) and (26), and the linearity of the operator
$(-\Delta)^{-1}$ that $g_{j}=(-\Delta)^{-1}u_{j}$ satisfies

$\displaystyle
g_{j}=\sum_{\ell,\boldsymbol{k}}|\boldsymbol{k}|^{-2}\,a^{\,j}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,.$
(27)

Moreover, from equation (26),
$D_{t}=\partial_{t}+\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}$,
and $A=(-\Delta)^{-1}D_{t}$, we have

$\displaystyle D_{t}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
$\displaystyle=(\partial_{t}+\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$})\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
(28)
$\displaystyle=\imath\ell^{\prime}\,\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}+[\boldsymbol{u}\,\mbox{\boldmath${\cdot}$}\,\imath\boldsymbol{k}^{\prime}]\,\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
$\displaystyle=\imath\ell^{\prime}\,\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}-\sum_{\ell,\boldsymbol{k}}\imath[\boldsymbol{a}_{\ell,\boldsymbol{k}}\,\mbox{\boldmath${\cdot}$}\,\boldsymbol{k}^{\prime}\,]\,\phi_{\ell^{\prime}+\ell,\boldsymbol{k}^{\prime}+\boldsymbol{k}}$
$\displaystyle A\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
$\displaystyle=\imath\ell^{\prime}|\boldsymbol{k}^{\prime}|^{-2}\,\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}-\sum_{\ell,\boldsymbol{k}}\imath[\boldsymbol{a}_{\ell,\boldsymbol{k}}\,\mbox{\boldmath${\cdot}$}\,\boldsymbol{k}^{\prime}\,]\,|\boldsymbol{k}^{\prime}+\boldsymbol{k}|^{-2}\,\phi_{\ell^{\prime}+\ell,\boldsymbol{k}^{\prime}+\boldsymbol{k}}\,,$

where the $-\imath$ comes from the sesquilinearity of the dot product.
Consequently, equations (27) and (28), and the linearity of the operators
$D_{t}$ and $A$ yield

$\displaystyle
D_{t}g_{j}=\sum_{\ell,\boldsymbol{k}}|\boldsymbol{k}|^{-2}a^{\,j}_{\ell,\boldsymbol{k}}\,D_{t}\phi_{\ell,\boldsymbol{k}}\,,\qquad
Ag_{j}=\sum_{\ell,\boldsymbol{k}}|\boldsymbol{k}|^{-2}a^{\,j}_{\ell,\boldsymbol{k}}\,A\phi_{\ell,\boldsymbol{k}}\,.$
(29)

The mass and moments of the spectral measure $\mu_{jk}$ can be expressed
explicitly in terms of just the Fourier coefficients
$\boldsymbol{a}_{\ell,\boldsymbol{k}}$ of $\boldsymbol{u}$, using the
orthonormality of the basis $\\{\phi_{\ell,\boldsymbol{k}}\\}$ and equations
(20)–(24), (27) and (29). For example equations (20), (24), and (27), and the
orthonormality of the $\phi_{\ell,\boldsymbol{k}}$ yield

$\displaystyle\mu^{0}_{jk}=\langle(-\Delta)^{-1}u_{j},u_{k}\rangle_{2}=\sum_{\ell,\boldsymbol{k}}|\boldsymbol{k}|^{-2}\,a^{\,j}_{\ell,\boldsymbol{k}}\,\bar{a}^{\,k}_{\ell,\boldsymbol{k}}\,,$
(30)

where $\bar{a}^{\,j}_{\ell,\boldsymbol{k}}$ is the complex conjugate of
$a^{\,j}_{\ell,\boldsymbol{k}}$. Similarly, the first $\mu^{1}_{jk}$ and
second $\mu^{2}_{jk}$ moments of the measure $\mu_{jk}$ can be expressed
explicitly in terms of just the Fourier coefficients
$\boldsymbol{a}_{\ell,\boldsymbol{k}}$ of $\boldsymbol{u}$. However, these
equations expressed in terms of the $\boldsymbol{a}_{\ell,\boldsymbol{k}}$ are
complicated and those for higher moments $\mu^{n}_{jk}$, $n=3,4,\ldots$,
expressed in terms of the $\boldsymbol{a}_{\ell,\boldsymbol{k}}$ become
unmanageable. However, the form of the formulas for $\mu^{n}_{jk}$ in
equations (21) and (22) show that an arbitrary number of the moments can be
found _iteratively_.

The key to doing so is determining how the operators $(-\Delta)^{-1}$,
$D_{t}=\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}$, and
$A=(-\Delta)^{-1}D_{t}$ map one iterate to the next. By linearity, it suffices
to understand how these operators transform a single Fourier mode
$p_{\ell,\boldsymbol{k}}\phi_{\ell,\boldsymbol{k}}$ to potentially multiple
modes
$\sum_{\ell^{\prime},\boldsymbol{k}^{\prime}}p^{\prime}_{\ell^{\prime},\boldsymbol{k}^{\prime}}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$.
By equations (24), (26) and (28), we have

$\displaystyle(-\Delta)^{-1}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
$\displaystyle=|\boldsymbol{k}^{\prime}|^{-2}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}\,,$
(31) $\displaystyle D_{t}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
$\displaystyle=\imath\ell^{\prime}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}-\sum_{\ell,\boldsymbol{k}}\imath[\boldsymbol{a}_{\ell,\boldsymbol{k}}\,\mbox{\boldmath${\cdot}$}\,\boldsymbol{k}^{\prime}\,]\,\phi_{\ell^{\prime}+\ell,\boldsymbol{k}^{\prime}+\boldsymbol{k}}\,,$
$\displaystyle A\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}$
$\displaystyle=\imath\ell^{\prime}|\boldsymbol{k}^{\prime}|^{-2}\phi_{\ell^{\prime},\boldsymbol{k}^{\prime}}-\sum_{\ell,\boldsymbol{k}}\imath[\boldsymbol{a}_{\ell,\boldsymbol{k}}\,\mbox{\boldmath${\cdot}$}\,\boldsymbol{k}^{\prime}\,]|\boldsymbol{k}^{\prime}+\boldsymbol{k}|^{-2}\,\phi_{\ell^{\prime}+\ell,\boldsymbol{k}^{\prime}+\boldsymbol{k}}\,,$

To make this explicit, we write

$\displaystyle
u_{j}=\sum_{\ell,\boldsymbol{k}}a^{\,j}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,,\qquad
g_{j}=\sum_{\ell,\boldsymbol{k}}b^{\,j}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,,\qquad
D_{t}g_{j}=\sum_{\ell,\boldsymbol{k}}c^{\,j}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,,\qquad
Ag_{j}=\sum_{\ell,\boldsymbol{k}}d^{\,j}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,.$
(32)

Therefore, equations (20), (21), and (32), and the orthonormality of the
$\phi_{\ell,\boldsymbol{k}}$ yield

$\displaystyle\mu^{0}_{jk}=\sum_{\ell,\boldsymbol{k}}b^{\,j}_{\ell,\boldsymbol{k}}\,\bar{a}^{\,k}_{\ell,\boldsymbol{k}}\,,\qquad\mu^{1}_{jk}=-\imath\sum_{\ell,\boldsymbol{k}}c^{\,j}_{\ell,\boldsymbol{k}}\,\bar{b}^{\,k}_{\ell,\boldsymbol{k}}\,,\qquad\mu^{2}_{jk}=\sum_{\ell,\boldsymbol{k}}c^{\,j}_{\ell,\boldsymbol{k}}\,\bar{d}^{\,k}_{\ell,\boldsymbol{k}}\,.$
(33)

In a similar way, for each $n=1,2,3,\ldots\,$, denote

$\displaystyle
D_{t}A^{n}g_{j}=\sum_{\ell,\boldsymbol{k}}c^{\,j\,,\,n}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,,\qquad
A^{n}g_{j}=\sum_{\ell,\boldsymbol{k}}d^{\,j\,,\,n}_{\ell,\boldsymbol{k}}\,\phi_{\ell,\boldsymbol{k}}\,,$
(34)

with $d^{\,j\,,\,1}_{\ell,\boldsymbol{k}}=d^{\,j}_{\ell,\boldsymbol{k}}$ and
$d^{\,j\,,\,n}_{\ell,\boldsymbol{k}}=|\boldsymbol{k}|^{-2}c^{\,j\,,\,n}_{\ell,\boldsymbol{k}}$.
Then, from equation (22) we have

$\displaystyle\mu^{2n+1}_{jk}=-\imath\sum_{\ell,\boldsymbol{k}}c^{\,j\,,\,n}_{\ell,\boldsymbol{k}}\,\bar{d}^{\,k\,,\,n}_{\ell,\boldsymbol{k}}\,,\qquad\mu^{2n+2}_{jk}=\sum_{\ell,\boldsymbol{k}}c^{\,j\,,\,n}_{\ell,\boldsymbol{k}}\,\bar{d}^{\,k\,,\,n+1}_{\ell,\boldsymbol{k}}\,.$
(35)

In the next section we will use the iterative mappings in (31) and formulas
for the spectral measure mass and moments in equations (33) and (35) to
calculate many moments for some steady and space-time periodic fluid flows.

##  4 Moment calculations

In this section we analytically and numerically calculate moments
$\mu^{2n}_{jk}$, for $n=0,1,2,\ldots$, and $j,k=1,\ldots,d$, for each of the
spectral measures $\mu_{jk}$. For shear flow we analytically calculate all of
the moments. These results show that each of the spectral measures is a
$\delta$-function centered at the spectral origin. For 2D steady BC-flow and a
related time-dependent 2D flow, we implement the iterative moment method
developed in Section 3 using Maple’s symbolic math toolbox and calculate
several moments analytically in closed form. We also implement the iterative
moment method in Matlab and compute many more moments to floating point
precision. In Section 5, we incorporate the positive moment values
$\mu^{2n}_{kk}$ in Padé approximant bounds for the diagonal components
$\mathfrak{D}^{*}_{kk}$ of the effective diffusivity matrix $\mathfrak{D}^{*}$
for these 2D fluid velocity fields.

###  4.1 Moments of shear flow

We now consider the special case of shear flow. The fluid velocity field for
shear flow is time-independent. In 2D shear flow in the $x$-direction has a
fluid velocity $\boldsymbol{u}=(0,\zeta(x))$ and $\boldsymbol{u}=(\xi(y),0)$
for flow in the $y$-direction, where in the current context $\xi$ and $\zeta$
are arbitrary mean-zero functions which are expressible by finite Fourier
series, as shown in equation in (24). Examples of 3D shear flow fluid velocity
fields are $\boldsymbol{u}=(\xi(y,z),0,0)$, $\boldsymbol{u}=(0,\zeta(x,z),0)$,
and $\boldsymbol{u}=(0,0,\gamma(x,y))$, where $\gamma$ is also an arbitrary
function of the same type as $\xi$ and $\zeta$.

The key property of these simple shear flows is that only one component of the
fluid velocity field is non-zero, the $i$th component say, and that component
is a function which is independent of $x_{i}$, the $i$th component of
$\boldsymbol{x}$. This property and equations (26) and (27) imply
$\mbox{\boldmath${\nabla}$}g_{j}=0$ for all $j=1,\ldots,d$, hence
$D_{t}g_{j}=\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}g_{j}=0$
for all $j=1,\ldots,d$ which, in turn, implies that
$Ag_{j}=(-\Delta)^{-1}D_{t}g_{j}=0$ for all $j=1,\ldots,d$ which, in turn,
implies that $A^{n}g_{j}=0$ and $D_{t}A^{n}g_{j}=0$ for all $j=1,\ldots,d$ and
$n=1,2,\ldots\,$. Therefore, by equations (20)–(22) we have

$\displaystyle\mu^{0}_{jk}=\langle
g_{j},u_{k}\rangle_{2}\,,\qquad\mu_{jk}^{n}=0\,,\quad j,k=1,\ldots,d\,,\quad
n=1,2,3,\ldots\;.$ (36)

Let’s first focus on the spectral measure $\mu_{kk}$ for some $k=1,\ldots,d$,
which is a _positive_ measure. It’s clear that the only positive Stieltjes
measure with all moments having value zero is a $\delta$-measure concentrated
at $\lambda=0$, $\delta_{0}(\mathrm{d}\lambda)$, where
$\delta_{a}(\mathrm{d}\lambda)$ is the $\delta$-measure concentrated at
$\lambda=a$. We therefore have the following result regarding the positive
measure $\mu_{kk}$ for simple shear flow,

$\displaystyle\mu_{kk}=\mu_{kk}^{0}\,\delta_{0}(\mathrm{d}\lambda)\,.$ (37)

We now show the signed measures $\mbox{Re}\,\mu_{jk}$ and
$\mbox{Im}\,\mu_{jk}$, $j\neq k$, also satisfy (37) in a weak sense. By the
Jordan decomposition theorem [9] there exist unique positive measures
$\mbox{Re}\,\mu_{jk}^{+}$ and $\mbox{Re}\,\mu_{jk}^{-}$ such that
$\mbox{Re}\,\mu_{jk}=\mbox{Re}\,\mu_{jk}^{+}-\mbox{Re}\,\mu_{jk}^{-}$ and
$\mbox{Re}\,\mu_{jk}^{+}\perp\mbox{Re}\,\mu_{jk}^{-}$, and similarly for the
signed measure $\mbox{Im}\,\mu_{jk}$. This and (36) imply the moments of these
measures satisfy

$\displaystyle[\mbox{Re}\,\mu_{jk}^{+}]^{n}=[\mbox{Re}\,\mu_{jk}^{-}]^{n}\,,\quad[\mbox{Im}\,\mu_{jk}^{+}]^{n}=[\mbox{Im}\,\mu_{jk}^{-}]^{n}\,,\quad
n=1,2,3,\ldots\,,$ (38)

where $[\mbox{Re}\,\mu_{jk}^{+}]^{n}$ and $[\mbox{Re}\,\mu_{jk}^{-}]^{n}$ are
the $n$th moments of the measures $\mbox{Re}\,\mu_{jk}^{+}$ and
$\mbox{Re}\,\mu_{jk}^{-}$, for example. Consequently, (36) implies for all
polynomials $P(\lambda)$ satisfying $P(0)=0$ that

$\displaystyle\int_{\Sigma}P(\lambda)\,\mathrm{d}\mbox{Re}\,\mu_{jk}(\lambda)=0\,,\qquad
j\neq k\,.$ (39)

Since $\mbox{Re}\,\mu_{jk}^{+}\perp\mbox{Re}\,\mu_{jk}^{-}$, it is clear that
this can only be the case if
$\mbox{Re}\,\mu_{jk}=\mbox{Re}\,\mu_{jk}^{0}\,\delta_{0}(\mathrm{d}\lambda)$.
A similar argument establishes that
$\mbox{Im}\,\mu_{jk}=\mbox{Im}\,\mu_{jk}^{0}\,\delta_{0}(\mathrm{d}\lambda)$.
Hence, we have that

$\displaystyle\mu_{jk}=\mu_{jk}^{0}\,\delta_{0}(\mathrm{d}\lambda)\,,\quad
j\neq k\,.$ (40)

This argument can be strengthened for the setting of a time-independent fluid
velocity field, $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{x})$. In this case,
the self-adjoint operator $M$ has bounded spectrum
$\Sigma\subseteq[-\|M\|,\|M\|]$ [20, 23]. The Stone Weierstrass theorem [9]
then extends the result in (39) to functions $\xi(\lambda)$ continuous on the
closed interval $[-\|M\|,\|M\|]$ satisfying $\xi(0)=0$, thus (40) holds in
this strengthened sense.

###  4.2 Moments for BC-flow

In this section we demonstrate how the iterative method developed in Section 3
can be used to compute an arbitrary number of moments for the spectral measure
$\mu_{jk}$ associated with $BC$-flow, which is given by the fluid velocity
field [15]

$\displaystyle\boldsymbol{u}=(C\cos{y},B\cos{x})\,.$ (41)

In particular, we implement the iterative method using Maple’s symbolic math
toolbox to calculate the mass $\mu_{jk}^{0}$ and the moments $\mu_{jk}^{2n}$,
$2n=2,4,\ldots,12$ analytically in closed form. We then discuss how higher
order moments can be numerically computed to floating point precision.

In terms of the Fourier basis $\\{\phi_{m,n}\\}$ for
$\mathscr{H}_{\mathcal{V}}$ [9], where $\phi_{m,n}(x,y)=\exp(\imath(mx+ny))$,
the components $u_{j}$, $j=1,2$, of the fluid velocity field
$\boldsymbol{u}=(C\cos{y},B\cos{x})$ are given by

$\displaystyle u_{1}$ $\displaystyle=(C/2)\,(\phi_{0,1}+\phi_{0,-1})\,,$ (42)
$\displaystyle u_{2}$ $\displaystyle=(B/2)\,(\phi_{1,0}+\phi_{-1,0})\,.$

From equations (26), (27), and (42) we have

$\displaystyle g_{j}=u_{j}\,,\quad j=1,2\,.$ (43)

The operator
$D_{t}=\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}=u_{1}\partial_{x}+u_{2}\partial_{y}$
is given by

$\displaystyle
D_{t}=(C/2)\,[\phi_{0,1}+\phi_{0,-1}]\,\partial_{x}+(B/2)\,[\phi_{1,0}+\phi_{-1,0}]\,\partial_{y}\,.$
(44)

Since $\partial_{x}u_{1}=0$ and $\partial_{y}u_{2}=0$, equations (26) and
(42)–(44) yield

$\displaystyle D_{t}\,g_{1}$ $\displaystyle=(C/2)(\imath
B/2)(\phi_{1,1}+\phi_{-1,1}-\phi_{1,-1}-\phi_{-1,-1}),$ (45) $\displaystyle
D_{t}\,g_{2}$ $\displaystyle=(B/2)(\imath
C/2)(\phi_{1,1}+\phi_{1,-1}-\phi_{-1,1}-\phi_{-1,-1})\,.$

Since $A=(-\Delta)^{-1}D_{t}$, it follows from equations (26) and (45) that

$\displaystyle Ag_{j}=D_{t}\,g_{j}/2\,,\quad j=1,2\,.$ (46)

$2n$ | $\mu_{kk}^{2n}$
---|---
0 | $\frac{1}{2}$
2 | $\frac{1}{8}$
4 | $\frac{3}{80}$
6 | $\frac{381}{32000}$
8 | $\frac{26277}{6800000}$
10 | $\frac{47519559}{37570000000}$
12 | $\frac{2960164002865793}{7127269448000000000}$
14 | $\frac{56807418712571064717219}{416027270403097600000000000}$
16 | $\frac{845725433928943189960402643663087}{18830209775901005048070400000000000000}$
18 | $\frac{2652281628393653311493590026036436288914383079}{179505850850574462175090974199721600000000000000000}$
20 | $\frac{404455666246342112121617203918794294909069461346892222329513233}{83202120549989484527334438746964410459680581766400000000000000000000}$
22 | $\frac{61678397622238580001722366830219450097176873936306735282205457266492841578250541}{38564719930926020344578565530438076789148157444429069162842585600000000000000000000000}$
24 | $\frac{24295131630327856551337800584398454447680582549232596740141494261473890932914297972532742346690937057}{46170636927742076288794326714045567306816555799878186559451173786860303947656387400000000000000000000000000}$
Table 1: Exact values of the spectral measure moments $\mu_{kk}^{2n}$ for BC-
flow with fluid velocity field in equation (41) with $B=C=1$, computed using
Maple’s symbolic math toolbox. Off-diagonal moments satisfy $\mu_{jk}^{2n}=0$,
$j\neq k$.

We now use equations (42)–(46) to compute the mass and the first two moments
of the spectral measure $\mu_{jk}$. From equations (32) and (43) we have
$b_{m,n}^{j}=a_{m,n}^{j}$, where $a_{0,1}^{1}=a_{0,-1}^{1}=C/2$ and
$a_{1,0}^{2}=a_{-1,0}^{2}=B/2$, and all other coefficients $a_{m,n}^{j}=0$.
Consequently, from equation (20), $\mu_{jk}^{0}=\langle
g_{j},u_{k}\rangle_{2}$, or the first formula in (33),
$\mu_{jk}^{0}=\sum_{m,n}b^{j}_{m,n}\bar{a}^{k}_{m,n}$, we have

$\displaystyle\mu_{11}^{0}=\frac{C^{2}}{2}\,,\qquad\mu_{22}^{0}=\frac{B^{2}}{2}\,,\qquad\mu_{12}^{0}=0\,.$
(47)

Similarly, from the formulas in equation (21) $\mu^{1}_{jk}=-\imath\,\langle
D_{t}g_{j},g_{k}\rangle_{2}$ and $\mu^{2}_{jk}=\langle
D_{t}g_{j},Ag_{k}\rangle_{2}$ or the last two formulas in equation (33),
$\mu^{1}_{jk}=-\imath\sum_{\ell,\boldsymbol{k}}c^{\,j}_{\ell,\boldsymbol{k}}\,\bar{b}^{\,k}_{\ell,\boldsymbol{k}}\,,$
and
$\mu^{2}_{jk}=\sum_{\ell,\boldsymbol{k}}c^{\,j}_{\ell,\boldsymbol{k}}\,\bar{d}^{\,k}_{\ell,\boldsymbol{k}}\,,$
we have $\mu^{1}_{jk}=0$ and

$\displaystyle\mu_{11}^{2}=\mu_{22}^{2}=\frac{B^{2}C^{2}}{8}\,,\qquad\mu_{12}^{2}=0\,.$
(48)

While the calculations of the mass and the first two moments are quite
transparent, calculations of higher order moments become complicated quickly.
Analytical closed form expressions can be obtained using Maple’s symbolic math
toolbox. These calculations indicate that $\mu_{12}^{2n}=0$ for all
$2n=0,2,4,6,\ldots$ and

$\displaystyle\mu^{4}_{11}=(1/320)B^{2}C^{4}+(11/320)B^{4}C^{2}\,,$ (49)
$\displaystyle\mu^{6}_{11}=\frac{3B^{2}C^{2}(101B^{4}+25B^{2}C^{2}+C^{4})}{32000}\,,$
$\displaystyle\mu^{8}_{11}=\frac{B^{2}C^{2}(567567B^{6}+233070B^{4}C^{2}+39610B^{2}C^{4}+617C^{6})}{217600000}\,,$
$\displaystyle\mu^{10}_{11}=\frac{C^{2}B^{2}(567567B^{6}+233070B^{4}C^{2}+39610B^{2}C^{4}+617C^{6})}{217600000}\,,$
$\displaystyle\mu^{12}_{11}=\frac{C^{2}B^{2}(3455217211B^{8}+1985186775B^{6}C^{2}+556637120B^{4}C^{4}+85049835B^{2}C^{6}+412611C^{8})}{4808960000000}\,,$

where $\mu^{2n}_{22}$ as a function of $B$ and $C$, $\mu^{2n}_{22}(B,C)$
satisfies $\mu^{2n}_{22}(B,C)=\mu^{2n}_{11}(C,B)$. We were able to calculate
the moments $\mu^{2n}_{ij}(B,C)$ for BC-flow in closed form up to $2n=26$ for
arbitrary $B,C\in(0,1]$, given the computational resources used. Some of these
are displayed in Table 1 for $B=C=1$.

Numerically computing higher order moments using floating point arithmetic can
be accomplished using the mapping in equation (31) as follows. The action of
the inverse Laplacian on a Fourier mode $\phi_{m,n}(x,y)=\exp(\imath(mx+ny))$
is given by,

$\displaystyle(-\Delta)^{-1}\phi_{m,n}=\phi_{m,n}/(m^{2}+n^{2})\,,\qquad
m^{2}+n^{2}\neq 0\,.$ (50)

Also, since $\phi_{i,j}\phi_{k,l}=\phi_{i+k,j+l}$, from equations (26) and
(44) we have

$\displaystyle D_{t}\phi_{m,n}=(\imath mC/2)\,\phi_{m,n+1}+(\imath
mC/2)\,\phi_{m,n-1}+(\imath nB/2)\,\phi_{m+1,n}+(\imath
nB/2)\,\phi_{m-1,n}\,,\quad m^{2}+n^{2}\neq 0\,,$ (51)

with $A\phi_{m,n}=(-\Delta)^{-1}D_{t}\phi_{m,n}$. Equations (33) and (35)
together with (50) and (51) define an iterative method that can be used to
compute Fourier coefficients for the functions in the iterative mapping chain:

$\displaystyle u_{j}\mapsto g_{j}\mapsto D_{t}\,g_{j}\mapsto Ag_{j}\mapsto
D_{t}Ag_{j}\mapsto A^{2}g_{j}\mapsto D_{t}A^{2}g_{j}\mapsto
A^{3}g_{j}\mapsto\ldots\,.$ (52)

This can be used to numerically compute an arbitrary number of moments
$\mu_{jk}^{n}$ for $BC$-flow, the maximum number limited only by numerical
accuracy and computational resources.

We have used this iterative method to calculate hundreds of spectral measure
moments $\mu_{11}^{2n}$ for BC-cell-flow with fluid velocity field in equation
(41) and $B=C=1$ using Matlab — shown in Fig. 1. Since the operator $A$ is
bounded for this 2D steady flow, the moments decrease exponentially for large
$n$, and Fig. 1(a) suggests the power law behavior $\mu_{11}^{2n}\sim
10^{-0.24(2n)}$ for $n\gg 1$. Fig. 1(b) shows the error of the computed values
of $\mu_{11}^{2n}$ using Matlab, relative to the exact values in Table 1, with
error values $\lesssim 10^{-15}$ indicating quite accurate numerical
computation of the $\mu_{11}^{2n}$.

Figure 1: Spectral measure moment computations for BC-cell-flow. (a) Numerical
computation of the moments $\mu_{11}^{2n}$ displayed in $\log$-10 scale for
BC-flow in (41) with $B=C=1$ using the iterative scheme shown in equations
(50)–(52). The linear fit to the data $\log_{10}(\mu_{11}^{2n})$ suggests the
power law behavior $\mu_{11}^{2n}\sim 10^{-0.24(2n)}$ for $n\gg 1$. (b) Error
of the floating point computations in (a) relative to the exact values, i.e.,
$|\mu^{2n}_{\text{float}}-\mu^{2n}_{\text{exact}}|/\mu^{2n}_{\text{exact}}$,
showing excellent agreement with error values $\lesssim 10^{-15}$ for
$2n=0,2,\ldots,26$.

###  4.3 Moments for a 2D time-periodic flow

In this section we demonstrate how the iterative method developed in Section 3
can be used to compute an arbitrary number of moments for the spectral measure
$\mu_{jk}$ associated with the space-time periodic flow [15, 5]

$\displaystyle\boldsymbol{u}=(C\cos{y},B\cos{x})+\theta\cos{t}\,(\sin{y},\sin{x}).$
(53)

In particular, we implement the iterative method using Maple’s symbolic math
toolbox to calculate the first few moments analytically in closed form. We
then demonstrate how higher order moments can be numerically computed to
floating point precision.

Denote $\boldsymbol{u}=\boldsymbol{u}_{BC}+\boldsymbol{u}_{\theta}$, where
$\boldsymbol{u}_{BC}=(C\cos{y},B\cos{x})$ is the fluid velocity field in
equation (42) and $\boldsymbol{u}_{\theta}=\theta\cos{t}\,(\sin{y},\sin{x})$.
The components $u_{j}$, $j=1,2$, of the fluid velocity field $\boldsymbol{u}$
in (53) expressed in terms of the Fourier basis $\\{\phi_{\ell,m,n}\\}$ for
$\mathscr{H}_{\mathcal{T}\mathcal{V}}$ [9], where
$\phi_{\ell,m,n}(x,y,t)=\exp(\imath(\ell t+mx+ny))$, are given by

$\displaystyle[\boldsymbol{u}_{BC}]_{1}$
$\displaystyle=(C/2)\,(\phi_{0,0,1}+\phi_{0,0,-1})\,,$ (54)
$\displaystyle[\boldsymbol{u}_{BC}]_{2}$
$\displaystyle=(B/2)\,(\phi_{0,1,0}+\phi_{0,-1,0})\,,$
$\displaystyle[\boldsymbol{u}_{\theta}]_{1}$
$\displaystyle=(\theta/4\imath)(\phi_{1,0,1}+\phi_{-1,0,1}-\phi_{1,0,-1}-\phi_{-1,0,-1})\,,$
$\displaystyle[\boldsymbol{u}_{\theta}]_{2}$
$\displaystyle=(\theta/4\imath)(\phi_{1,1,0}+\phi_{-1,1,0}-\phi_{1,-1,0}-\phi_{-1,-1,0})\,.$

Writing $\boldsymbol{g}=\boldsymbol{g}_{BC}+\boldsymbol{g}_{\theta}$ with
$\boldsymbol{g}=(-\Delta)^{-1}\boldsymbol{u}$, using the notation that
$(-\Delta)^{-1}$ operates component-wise on $\boldsymbol{u}$, we have from
(54)

$\displaystyle\boldsymbol{g}=\boldsymbol{u}_{BC}+\boldsymbol{u}_{\theta}\,.$
(55)

In this time-independent setting,
$D_{t}=\partial_{t}+\boldsymbol{u}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}$,
hence by linearity

$\displaystyle
D_{t}=\partial_{t}+\boldsymbol{u}_{BC}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}+\boldsymbol{u}_{bc}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\;.$
(56)

From equation (26), we have

$\displaystyle\partial_{t}\phi_{\ell,m,n}=\imath\ell\phi_{\ell,m,n}\,,\quad\ell\neq
0.$ (57)

We have the following analogue of equation (51), for $m^{2}+n^{2}\neq 0$,

$\displaystyle\boldsymbol{u}_{BC}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\phi_{\ell,m,n}=(\imath
mC/2)\,\phi_{\ell,m,n+1}+(\imath mC/2)\,\phi_{\ell,m,n-1}+(\imath
nB/2)\,\phi_{\ell,m+1,n}+(\imath nB/2)\,\phi_{\ell,m-1,n}\,.$ (58)

From equation (54) we have, for $m^{2}+n^{2}\neq 0$,

$\displaystyle\boldsymbol{u}_{\theta}\mbox{\boldmath${\cdot}$}\mbox{\boldmath${\nabla}$}\phi_{\ell,m,n}$
$\displaystyle=(\imath mc/4\imath)\,\phi_{\ell+1,m,n+1}+(\imath
mc/4\imath)\,\phi_{\ell-1,m,n+1}-(\imath mc/4\imath)\,\phi_{\ell+1,m,n-1}$
(59) $\displaystyle-(\imath mc/4\imath)\,\phi_{\ell-1,m,n-1}+(\imath
nb/4\imath)\,\phi_{\ell+1,m+1,n}+(\imath nb/4\imath)\,\phi_{\ell-1,m+1,n}$
$\displaystyle-(\imath nb/4\imath)\,\phi_{\ell+1,m-1,n}-(\imath
nb/4\imath)\,\phi_{\ell-1,m-1,n}\,.$

Similarly, we have the following analogue of equation (50)

$\displaystyle(-\Delta)^{-1}\phi_{\ell,m,n}=\phi_{\ell,m,n}/(m^{2}+n^{2})\qquad
m^{2}+n^{2}\neq 0.$ (60)

Equations (33) and (35) together with equations (58)–(60) enable all the
Fourier coefficients in the iterative sequence shown in equation (52) to be
obtained analytically in closed form.

Dynamic BC flow
---
$2n$ | $\mu_{kk}^{2n}$ | $\mu_{jk}^{2n}$
0 | 3/4 | 0
2 | 35/64 | 0
4 | 1069/1280 | 5/32
6 | 282223/163840 | 473/800
8 | 15778075321/3481600000 | 83952089/51200000
Table 2: Exact values of the spectral measure moments $\mu_{kk}^{2n}$ and
$\mu_{jk}^{2n}$, $j\neq k$, for the fluid velocity field in equation (53) with
$B=C=\theta=1$, calculated using Maple’s symbolic math toolbox.

Using Maple’s symbolic math toolbox, this yields the spectral measure masses

$\displaystyle\mu_{11}^{0}$
$\displaystyle=\frac{C^{2}}{2}+\frac{\theta^{2}}{4},\qquad\mu_{22}^{0}=\frac{B^{2}}{2}+\frac{\theta^{2}}{4},\qquad\mu_{12}^{0}=0.$
(61)

Recall that all odd moments are identically zero. The second, fourth, and
sixth moments are given by

$\displaystyle\mu_{11}^{2}$
$\displaystyle=\frac{3\theta^{4}}{64}+\frac{(B^{2}+C^{2}+4)\theta^{2}}{16}+\frac{B^{2}C^{2}}{8},$
(62) $\displaystyle\mu_{12}^{2}$ $\displaystyle=0,$
$\displaystyle\mu_{11}^{4}$
$\displaystyle=\frac{5\theta^{6}}{3072}+\frac{(323B^{2}+705C^{2}+1090)\theta^{4}}{92160}+\frac{(180+22B^{4}+(366C^{2}+245)B^{2}+8C^{4}+80C^{2})\theta^{2}}{11520}$
$\displaystyle+\frac{11B^{4}C^{2}}{320}+\frac{C^{4}B^{2}}{320},$
$\displaystyle\mu_{12}^{4}$
$\displaystyle=-\frac{13BC\theta^{4}}{5760}-\Big{(}\frac{18(B^{3}C+BC^{3})+35BC}{11520}\Big{)}\theta^{2},$
$\displaystyle\mu_{11}^{6}$
$\displaystyle=\frac{20701\theta^{8}}{49152000}+\Big{(}\frac{31583}{4608000}+\frac{4861B^{2}}{3686400}+\frac{3119C^{2}}{1228800}\Big{)}\theta^{6}+\Big{(}\frac{49}{2304}+\frac{163871B^{2}}{9216000}+\frac{4343B^{4}}{3072000}$
$\displaystyle+\frac{48329C^{2}}{3072000}+\frac{108793B^{2}C^{2}}{9216000}+\frac{161C^{4}}{204800}\Big{)}\theta^{4}+\Big{(}\frac{1}{64}+\frac{337B^{2}}{9216}+\frac{33191B^{4}}{2304000}+\frac{101B^{6}}{192000}+\frac{C^{2}}{576}$
$\displaystyle-\frac{2569B^{2}C^{2}}{2304000}+\frac{773B^{4}C^{2}}{48000}+\frac{11C^{4}}{16000}+\frac{137B^{2}C^{4}}{36000}+\frac{C^{6}}{48000}\Big{)}\theta^{2}+\frac{303B^{6}C^{2}}{32000}+\frac{3B^{4}C^{4}}{1280}+\frac{3B^{2}C^{6}}{32000},$
$\displaystyle\mu_{12}^{6}$
$\displaystyle=\frac{437BC\theta^{6}}{4608000}+\Big{(}\frac{6509BC}{1536000}+\frac{521B^{3}C}{576000}+\frac{521BC^{3}}{576000}\Big{)}\theta^{4}+\Big{(}\frac{65BC}{9216}-\frac{12841B^{3}C}{2304000}-\frac{9B^{5}C}{32000}-\frac{12841BC^{3}}{2304000}$
$\displaystyle+\frac{91B^{3}C^{3}}{28800}-\frac{9BC^{5}}{32000}\Big{)}\theta^{2}\,,$

where $\mu_{22}^{2n}$ as a function of $B$, $C$, and $\theta$,
$\mu_{22}^{2n}(B,C,\theta)$ satisfies
$\mu_{22}^{2n}(B,C,\theta)=\mu_{11}^{2n}(C,B,\theta)$. Exact rational values
for $2n=0,2,\ldots,8$ are displayed in Table 2 for $B=C=\theta=1$.

Figure 2: Spectral measure moment computations for a 2D time-periodic flow.
(a) Numerical computation of moments $\mu_{11}^{2n}$ displayed in $\log$-10
scale for the time-periodic flow in (53) with $B=C=\theta=1$ using the
iterative scheme shown in equations (57)–(60). The linear fit to the data
$\log_{10}(\mu_{11}^{2n})$ suggests the power law behavior $\mu_{11}^{2n}\sim
10^{1.07(2n)}$ for $n\gg 1$. (b) Error of the numerical moment computations in
(a) relative to the exact values, showing excellent agreement with error
values $\lesssim 10^{-14.8}$ for $2n=0,2,4,6,8$.

We have used this iterative method to calculate hundreds of spectral measure
moments $\mu_{11}^{2n}$ for the 2D time-periodic fluid velocity field in
equation (53) with $B=C=\theta=1$ using Matlab — shown in Fig. 2. Since the
operator $A$ is unbounded for this time-periodic flow, the moments increase
exponentially for large $n$, and Fig. 2(a) suggests the power law behavior
$\mu_{11}^{2n}\sim 10^{0.54(2n)}$ for $n\gg 1$. Fig. 1(b) shows the error of
the numerically computed values of $\mu_{11}^{2n}$ using Matlab relative to
the exact values in Table 2, with error values $\lesssim 10^{-14.8}$
indicating quite accurate numerical computation of the $\mu_{11}^{2n}$.

##  5 Bounds for effective diffusivity

In the advection dominated regime, where $\varepsilon\ll 1$, it is known [6]
that the asymptotic behavior of the effective diffusivity for BC-cell flow
discussed in Section 4.2 is given by
$\mathfrak{D}^{*}_{kk}\sim\varepsilon^{1/2}$ as $\varepsilon\to 0$. The
physical idea for this vanishing behavior of $\mathfrak{D}^{*}_{kk}$ as
$\varepsilon\to 0$ is as follows. For $\varepsilon=0$, passive tracers are
trapped in closed streamlines of the cell flow [16, 15] and without the
presence of molecular diffusivity $\varepsilon$ there can be no large scale,
macroscopic transport. However, for $\varepsilon>0$, molecular diffusion can
move tracers to neighboring streamlines and ultimately to neighboring period
cells, giving rise to macroscopic transport on large spatial scales, with
$\mathfrak{D}^{*}_{kk}\sim\varepsilon^{1/2}$ as $\varepsilon\to 0$. In [5] it
was conjectured that, due to the chaotic dynamics of the time-dependent fluid
velocity field $\boldsymbol{u}$ in (53), the effective diffusivity
$\mathfrak{D}^{*}_{kk}$ has the asymptotic behavior $\mathfrak{D}^{*}_{kk}\sim
1$ for $\varepsilon\ll 1$, which is a phenomenon known as residual
diffusivity.

Both of the asymptotic behaviors $\mathfrak{D}^{*}_{kk}\sim\varepsilon^{1/2}$
for BC- flow and $\mathfrak{D}^{*}_{kk}\sim 1$ for the dynamic flow in (53)
were numerically verified in [15, 16]. This was accomplished in [15] by
developing Fourier methods for computing the effective diffusivity
$\mathfrak{D}^{*}_{kk}$ by directly computing the spectral measure $\mu$. The
cell problem in equation (17) was represented as an infinite system of
equations involving the Fourier coefficients of the fluid velocity field
$\boldsymbol{u}$. Truncating the system of equations results in a Hermitian
matrix representation for the operator $M=-\imath A$ in equations (17) and
(18). The spectral measure of the matrix is given by a finite sum of weighted
Dirac-$\delta$ functions
$\mathrm{d}\mu(\lambda)=\sum_{i}m_{i}\,\delta_{\lambda_{i}}(\mathrm{d}\lambda)$,
for $i=0,1,2,\ldots$. Here, the spectral weights $m_{i}$ are given explicitly
in terms of the eigenvectors of the matrix and $\lambda_{i}$ are the
eigenvalues of the matrix.

In this section, we describe these and other phenomena in terms of Padé
approximants of the diagonal coefficients $\mathfrak{D}^{*}_{kk}$,
$k=1,\ldots,d$, of the effective diffusivity matrix which involve the moments
$\mu_{kk}^{n}$ of the _positive_ measure $\mu_{kk}$. For simplicity, let’s
focus on one diagonal coefficient and set $\mu_{kk}^{n}=\mu^{n}$. With this
notation, equation (19) can be written as
$\mathfrak{D}^{*}_{kk}=\varepsilon(1+\varepsilon^{-2}f(\varepsilon^{-2}))$,
where $f(z)$ is a _Stieltjes function_ [3] of the variable
$z=\varepsilon^{-2}$ involving the positive measure $\mu$.

Padé approximants $\\{[N-1/N]\\}_{N=1}^{\infty}$ and
$\\{[N/N]\\}_{N=1}^{\infty}$ form rigorous, converging lower and upper bounds
for the Stieltjes function $f(z)$ [3]

$\displaystyle[N-1/N](z)\leq f(z)\leq[N/N](z)\,,\qquad
f(z)=\int_{-\infty}^{\infty}\frac{\mathrm{d}\mu(\lambda)}{1+z\lambda^{2}}\,.$
(63)

Since all of the terms of the expression in
$\mathfrak{D}^{*}_{kk}=\varepsilon(1+\varepsilon^{-2}f(\varepsilon^{-2}))$ are
positive, we therefore have the following rigorous lower and upper bounds for
$\mathfrak{D}^{*}_{kk}$

$\displaystyle\varepsilon(1+\varepsilon^{-2}[N-1/N](\varepsilon^{-2}))\leq\mathfrak{D}^{*}_{kk}(\varepsilon)\leq\varepsilon(1+\varepsilon^{-2}[N/N](\varepsilon^{-2}))\,.$
(64)

The theory of Padé approximants for $f(z)$ follows by expanding
$1/(1+z\lambda^{2})$ in a geometric series and writing $f(z)$ as a _series of
Stieltjes_

$\displaystyle f(z)=\sum_{n=0}^{\infty}c_{n}\,z^{n}\,,\qquad
c_{n}=(-1)^{n}\,\mu^{2n},\qquad\mu^{n}=\int_{-\infty}^{\infty}\lambda^{n}\,\mathrm{d}\mu(\lambda)\,.$
(65)

The $L$, $N$ Padé approximant of $f(z)$ is given by

$\displaystyle[L/N](z)=\frac{P^{[L/N]}(z)}{Q^{[L/N]}(z)}\,,$ (66)

where $P^{[L/N]}(z)$ is a polynomial of degree of at most $L$ and
$Q^{[L/N]}(z)$ is a polynomial of degree of at most $N$. The _formal_ power
series for $f(z)$ in equation (65) and the condition
$f(z)-P^{[L/N]}(z)/Q^{[L/N]}(z)=O(z^{L+N+1})$ determines the coefficients of
$P^{[L/N]}(z)$ and $Q^{[L/N]}(z)$ and these polynomials can be written in
terms of determinants involving the $c_{n}$ [3].

###  5.1 Numerical implementation

Figure 3: Padé bound orders for effective diffusivity. Attempted and actual
numerator and denominator polynomial orders for $[N/N]$ upper and $[N-1/N]$
lower bounds for various values of $N$. These bounds for Stieltjes functions
associated with the effective diffusivity $\mathfrak{D}^{*}_{kk}$ are for (a)
BC-flow and (b) 2D time-periodic flow with fluid velocity fields in equations
(41) and (53). The actual numerator and denominator polynomial orders returned
for the attempted $[N/N]$ upper and $[N-1/N]$ lower bounds are
$[\alpha_{upper}/\beta_{upper}]$ and $[\alpha_{lower}/\beta_{lower}]$.

Given the moments $\mu^{2n}$ obtained in Section 4.3, it might seem that Padé
approximants for $\mathfrak{D}^{*}_{kk}$ of any order could be attainable.
However, even in the absence of rounding errors on a computer, the theoretical
treatment of Padé approximants is subject to the appearance of seemingly
spurious pole-zero pairs or “Froissart doublets” in arbitrary locations that
prevent pointwise convergence [10]. Such anomalies become common in the
presence of rounding errors or other forms of noise [10]. Hence, computing the
Padé approximants directly in terms of the determinants is essentially ill-
posed and becomes unreliable for large polynomial degree. Padé approximation
uses information about the series of Stieltjes $f(z)$ in (65) at a single
point $z=0$ in the complex plane to gain information about regions of the
complex plane away from $z=0$ — similar to analytic continuation — and becomes
more ill-posed further away from $z=0$. This is a fundamental limitation of
using Padé approximant bounds to study the effective diffusivity in the
advection dominated regime, where $\varepsilon\ll 1$, i.e., $z\gg 1$. Despite
this, we will show for the case of BC-flow that Padé approximants capture the
known asymptotic behavior $\mathfrak{D}^{*}_{kk}\sim\varepsilon^{1/2}$ for
$\varepsilon\gg 1$ [6].

Such issues associated with numerical computation of Padé approximants are
addressed in [10], where a Matlab function padeapprox was released, which is
freely available as part of Chebfun, which provides robust Padé approximants
via singular value decomposition (SVD) [10]. This numerical method does not
enable Padé approximants of an arbitrary order to be computed, but instead
truncates the Padé approximants to an order that ensures numerical stability,
hence numerically dependable results — incorporating any additional moments
typically results in the same truncated Padé order. This automatic truncation
is accomplished by introducing a tolerance tol that is used to zero out
singular values less than tol. For most purposes involving problems perturbed
just by rounding errors, tol $=10^{-14}$ is a reasonable value [10]. This
regularization can be circumvented by setting tol=0 [10]. The desired degrees
of the numerator/denominator polynomials, in our case $N/N$ and $N-1/N$, are
provided along with a vector of the coefficients $c_{n}$ and tol; Numerator
and denominator polynomials of actual degree $\alpha$ and $\beta$ are
returned, as shown in Figure 3.

Figure 4:  Scaled spectral measure moments. Unscaled and scaled spectral
measure moments for (a) BC-flow and (b) space-time periodic flow. (a) The
scaled moments $\mu^{2n}_{11}\gamma^{n}$ for BC-flow with $\gamma=3.0387$
reduces the variation of the moments from $\sim 7$ orders of magnitude to
$\sim 1$ order of magnitude. (b) The scaled moments $\mu^{2n}_{11}\gamma^{n}$
for space-time periodic flow with $\gamma=0.0342$ reduces the variation of the
moments from $\sim 14$ orders of magnitude to $\sim 10$ orders of magnitude.

The numerical computation of Padé approximants can be further stabilized by
scaling the $c_{n}$ in (65) by picking a $\gamma>0$ so that
$f(z)=\sum_{n}(c_{n}\gamma^{n})(z/\gamma)^{n}$ has scaled coefficients
$c_{n}\gamma^{n}$. The goal is to obtain scaled coefficients that have roughly
comparable orders of magnitude, neither decreasing nor increasing at a rapid
geometric rate, for which the algorithm padeapprox is most effective [10].
Figure 4a shows this rescaling is effective for BC-flow, with the variation in
the moments decreasing from $\sim 7$ orders of magnitude to $\sim 1$ order of
magnitude. Figure 4b shows the moment scaling for the time-periodic flow is
less effective but does decrease the variation in the moments from $\sim 14$
orders of magnitude to $\sim 10$ orders of magnitude. Once Padé approximants
are calculated using these scaled moments the original approximants are
recovered by simply using the scaled parameter $z\mapsto z/\gamma$.

###  5.2 Bounds for BC-flow

In this section we discuss our results for Padé approximant upper $[N/N]$ and
lower $[N-1/N]$ bounds for the effective diffusivity $\mathfrak{D}^{*}_{kk}$
for BC-flow with fluid velocity field given in equation (41). In particular,
we incorporated into the Matlab function padeapprox the moments $\mu^{2n}$,
$2n=0,\ldots,26$, for the spectral measure $\mu$ for BC-cell-flow with
$B=C=1$. Figures 5(a) and (b) display, in linear and log-log scale
respectively, progressively higher order Padé approximants $[N/N]$ and
$[N-1/N]$ for the effective diffusivity $\mathfrak{D}^{*}_{kk}$ as a function
of $\varepsilon$ for $N=1,\ldots,7$. Here, the maximum order $N=7$ corresponds
to the upper $[6/6]$ and lower $[6/7]$ Padé approximants due to the truncation
which occurs as illustrated in Figure 3(a). From
$\mathfrak{D}^{*}_{kk}=\varepsilon(1+\varepsilon^{-2}f(\varepsilon^{-2}))$
with $f\geq 0$ we see the global lower bound is
$\mathfrak{D}^{*}_{kk}\geq\varepsilon$, which is shown in Figures 5(a) and (b)
in black-dash line style. The global upper bound is attained by shear flow
[2], discussed in Section 4.1, where the spectral measure is given by a delta
measure at the spectral origin, $\mu=\mu^{0}\,\delta_{0}(\mathrm{d}\lambda)$,
yielding the global upper bound
$\mathfrak{D}^{*}_{kk}\leq\varepsilon(1+\mu^{0}/\varepsilon^{2})$, also shown
in Figures 5(a) and (b) in black-dash line style. The difference
$[N/N]-[N-1/N]$ of upper and lower bounds are displayed in Figure 5(c),
showing that the bounds are quite accurate for large $\varepsilon$, i.e.,
small $z$, and they become less accurate as $\varepsilon\to 0$, as anticipated
and discussed in Section 5.1.

Figure 5:  Effective diffusivity bounds for BC-flow. Padé approximant upper
$[N/N]$ and lower $[N-1/N]$ bounds for $\mathfrak{D}^{*}_{kk}(\varepsilon)$ as
a function of $\varepsilon$ for $N=1,\ldots,7$ in (a) linear scale and (b)
log-log scale. (c) Differences of upper and lower bounds $[N/N]-[N-1/N]$ for
various values of $N$ in log-log scale. (d) Linear polynomial fits to the
highest order upper and lower bounds with $N=7$ demonstrating the asymptotic
behavior $\mathfrak{D}^{*}_{kk}\sim\varepsilon^{1/2}$ as $\varepsilon\to 0$.

The bounds in Figures 5(a) and (b) get progressively tighter with increasing
$N$, as indicated by the bound differences in Figure 5(c). Figures 5(a), (b),
and (c) show that for $\varepsilon\gtrsim 1$ all the bounds but the global
upper and lower bounds lie virtually right on top of each other. The highest
order bounds shown in Figure 5 incorporate the moments $2n=0,\ldots,24$ and
$2n=0,\ldots,26$ for the upper and lower bounds, respectively. They capture
the behavior of $\mathfrak{D}^{*}_{kk}$ in the advection dominated regime
where $\varepsilon\ll 1$. Linear polynomial fits to the bounds $[6/6]$ and
$[6/7]$ in log-log scale capture the asymptotic behavior of the effective
diffusivity $\mathfrak{D}^{*}_{kk}\sim\varepsilon^{1/2}$ discussed in Section
5, with estimated critical exponents 0.49916 and 0.50166, respectively, as
shown in Figure 5(d). The bound differences shown in Figure 5(c) indicate that
the highest order bounds are accurate to within 0.01 for
$\varepsilon\gtrapprox 0.053$ and within 0.1 for $\varepsilon\gtrapprox 0.025$
with $\mathfrak{D}^{*}_{kk}(0.053)\approx 0.35$ and
$\mathfrak{D}^{*}_{kk}(0.025)\approx 0.24$ — demonstrating $\approx 6.6$ and
$\approx 9.6$ times the enhancement of the effective diffusivity above
$\varepsilon$ even for the small values of $\varepsilon=0.053$ and
$\varepsilon=0.025$, respectively.

When more than $2n=26$ moments are supplied to the function padeapprox, in
order to ensure numerical accuracy and stability, the function just truncates
the Padé order to $N=7$, as shown in Figure 3(a); supplying more moments just
leads to the same repeated maximum order $N=7$ of the Padé bounds. A tolerance
value tol$=10^{-14}$ was used for all bounds in Figure 5. Reducing this value
does lead to bounds of a higher order than $N=7$ when more moments than
$2n=26$ are supplied, at the expense of numerical accuracy of the Padé bounds.

###  5.3 Bounds for space-time periodic flow

Figure 6:  Effective diffusivity bounds for space-time periodic flow. Padé
approximant upper $[N/N]$ and lower $[N-1/N]$ bounds for $N=1,\ldots,8$ in
linear scale (a) with corresponding absolute differences (b), and upper and
lower bounds in log-log scale (c) in the variable $\varepsilon$. (d) Linear
polynomial fits to the highest order upper and lower bounds with $N=8$.

In this section we discuss our results for Padé approximant upper $[N/N]$ and
lower $[N-1/N]$ bounds for the effective diffusivity $\mathfrak{D}^{*}_{kk}$
for the space-time periodic flow with fluid velocity field given in equation
(53) with $B=C=\theta=1$. In particular, we incorporated into the Matlab
function padeapprox the moments $\mu^{2n}$, $2n=0,\ldots,32$, for the
associated spectral measure $\mu$. Similar to Figure 5, Figures 6(a) and (b)
display progressively higher order Padé approximants $[N/N]$ and $[N-1/N]$ for
the effective diffusivity $\mathfrak{D}^{*}_{kk}$ as a function of
$\varepsilon$, with global lower and upper bounds displayed in black-dash line
style.

The bounds in Figures 6(a) and (b) get progressively tighter with increasing
$N$. For $\varepsilon\gtrsim 1$ all bounds but the global upper and lower
bounds lie virtually right on top of each other, as demonstrated
quantitatively by the bound differences in Figure 6(c). The highest order
bounds shown in Figure 6 correspond to $N=8$ which incorporates the moments
$2n=0,\ldots,32$. Figure 6(c) indicates the highest order bounds are accurate
within a difference of 0.01 for $\varepsilon\gtrapprox 0.59$ and within a
difference of 0.1 for $\varepsilon\gtrapprox 0.36$ with
$\mathfrak{D}^{*}_{kk}(0.59)\approx 1.22$ and
$\mathfrak{D}^{*}_{kk}(0.36)\approx 1.06$, indicating $\approx 2.07$ and
$\approx 2.94$ times the enhancement of the effective diffusivity above
$\varepsilon$ when $\varepsilon=0.59$ and $\varepsilon=0.36$, respectively.

Comparing this to the value of the effective diffusivity for BC-flow for
$\varepsilon=0.59$ and $\varepsilon=0.36$ shows the addition of the time
dependent term $\boldsymbol{u}_{\theta}=\theta\cos{t}\,(\sin{y},\sin{x})$ with
$\theta=1$ leads to a $\approx 5\%$ and $\approx 16\%$ increase in the
effective diffusivity, respectively — hence an even greater enhancement above
the molecular diffusivity value. In [5, 15] numerical results indicate the
addition of the time dependent term gives rise to an additional
$\mathfrak{D}^{*}_{kk}$ enhancement of $\approx 10^{1.5}$ times that for BC-
flow when $\varepsilon\approx 10^{-4}$, with $\mathfrak{D}^{*}_{kk}$
flattening out into a $\mathfrak{D}^{*}_{kk}\sim 1$ behavior for
$\varepsilon\approx 0.1$ with $\mathfrak{D}^{*}_{kk}(10^{-4})\approx
10^{-1/2}$, compared to $\mathfrak{D}^{*}_{kk}(10^{-4})\approx 10^{-2}$ for
BC-flow.

##  6 Conclusions

The effective diffusivity $\mathfrak{D}^{*}$ associated with both spatially
and space-time periodic fluid flows has a Stieltjes function representation
involving a spectral measure $\mu$ of a self-adjoint operator and the Péclet
number Pe of the flow. We considered a non-dimensionalization that _separates_
the geometry and dynamics of the fluid velocity field — encoded in $\mu$ —
from the strength of the fluid velocity field in Pe, which leads to the
following simple relationship between the molecular diffusivity $\varepsilon$
and Pe: $\mbox{{Pe}}=1/\varepsilon$. The theory of Padé approximants for
Stieltjes functions provides rigorous nested upper and lower bounds for
$\mathfrak{D}^{*}$ that incorporate the moments of $\mu$. The parameter
separation property makes $\mu$ and the bounds independent of the geometry and
dynamics of the fluid velocity field. As more moments are incorporated the
bounds get tighter and can converge to the true value of $\mathfrak{D}^{*}$
for certain values of Pe near $\mbox{{Pe}}=0$. Fundamental theoretical and
numerical limitations of calculating Padé approximants far away from
$\mbox{{Pe}}=0$ limit the accuracy of bounds for $\mathfrak{D}^{*}$ in the
advection dominated regime, where $\mbox{{Pe}}\gg 1$, i.e., $\varepsilon\ll
1$. These limitations also limit the number of bounds that can be stably
computed.

An iterative method was introduced that enables an arbitrary number of moments
of $\mu$ to be calculated. This method was implemented in Maple to calculate
several moments of both a steady and space-time periodic flow, in closed form
as a function of flow parameters. This method was also implemented in Matlab
to compute hundreds of moments to floating point precision, having excellent
agreement with the exact values. These moments were incorporated into an
existing numerical algorithm padeapprox [10] which computes Padé approximants
in a robust, stable way.

Despite the fundamental limitations of calculating Padé approximants far away
from $\mbox{{Pe}}=0$, for the steady cell flow, BC-cell-flow, the known [6, 7]
asymptotic behavior $\mathfrak{D}^{*}\sim\varepsilon^{1/2}$, for
$\varepsilon\ll 1$, in the advection dominated regime was accurately captured
by high order upper and lower bounds for $\mathfrak{D}^{*}$, for values of
$\varepsilon$ down to $\varepsilon\approx 0.053$. Bound differences for these
highest order bounds indicate $\mathfrak{D}^{*}$ was accurately computed
within $0.01$ for $\varepsilon\gtrapprox 0.053$ and within $0.1$ for
$\varepsilon\gtrapprox 0.025$.

For the space-time periodic flow, which is BC-cell-flow with an additional
space-time periodic term, the operator underlying the spectral measure $\mu$
is _unbounded_ which makes the computation of Padé approximant bounds for
$\mathfrak{D}^{*}$ less stable than the bounds for BC-cell-flow. This
restricted our analysis of $\mathfrak{D}^{*}$ in the advection dominated
regime. Bound differences for the highest order bounds indicate the bounds are
accurate to within $0.01$ for $\varepsilon\gtrapprox 0.59$ and within $0.1$
for $\varepsilon\gtrapprox 0.36$. Comparing the values of $\mathfrak{D}^{*}$
for BC-cell-flow and the space-time periodic flow for $\varepsilon=0.59$ and
$\varepsilon=0.36$ shows the addition of the time dependent term to BC-cell-
flow leads to a $\approx 5\%$ and $\approx 16\%$ increase in
$\mathfrak{D}^{*}$, respectively.

The bounds for $\mathfrak{D}^{*}$ for the space-time periodic flow are not
accurate for small enough $\varepsilon$ to capture the _residual diffusivity_
phenomenon suggested by direct numerical computations of $\mathfrak{D}^{*}$.
In particular, in [5, 15] numerical results indicate the addition of the time
dependent term to BC-cell-flow gives rise to an additional $\mathfrak{D}^{*}$
enhancement of $\approx 10^{1.5}\approx 31.62$ times that for BC-cell-flow
when $\varepsilon\approx 10^{-4}$, with $\mathfrak{D}^{*}$ flattening out into
a $\mathfrak{D}^{*}\sim 1$ behavior when $\varepsilon\lessapprox 0.1$ with
$\mathfrak{D}^{*}_{kk}(10^{-4})\approx 10^{-1/2}$, compared to
$\mathfrak{D}^{*}_{kk}(10^{-4})\approx 10^{-2}$ for BC-flow.

The iterative moment method developed in this manuscript enables analytic
calculation of the moments of $\mu$ in closed form for a large class of fluid
velocity fields $\boldsymbol{u}$, namely spatially and space-time periodic
$\boldsymbol{u}$ with finite Fourier series, which are dense in the associated
spaces of square integrable periodic functions [9]. This and the theory of
Padé approximants for Stieltjes functions provide a remarkable advance in
estimating the effective diffusivity, even in the advection dominated regime
where $\varepsilon$ is small for some $\boldsymbol{u}$.

##  7 Acknowledgements

We gratefully acknowledge support from the Applied and Computational Analysis
Program and the Arctic and Global Prediction Program at the US Office of Naval
Research through grants N00014-18-1-2552, N00014-21-1-2909, N00014-13-10291,
N00014-15-1-2455, and N00014-18-1-2041). We are also grateful for support from
the Division of Mathematical Sciences and the Division of Polar Programs at
the US National Science Foundation (NSF) through grants DMS-2309520,
DMS-0940249, DMS-2136198, DMS-2111117, DMS-2206171, DMS-1715680, and
DMS-1413454. Finally, we would like to thank the NSF Math Climate Research
Network (MCRN), and especially Chris Jones, for supporting this work.

## References

  * [1] M. Avellaneda and A. Majda.  Stieltjes Integral Representation and Effective Diffusivity Bounds for Turbulent Transport.  Phys. Rev. Lett., 62:753–755, 1989.
  * [2] M. Avellaneda and A. Majda.  An Integral Representation and Bounds on the Effective Diffusivity in Passive Advection by Laminar and Turbulent Flows.  Comm. Math. Phys., 138:339–391, 1991.
  * [3] G. A. Baker and P. R. Graves-Morris.  Padé Approximants.  Encyclopedia of Mathematics and its Applications. Cambridge University Press, 1996.
  * [4] R. Bhattacharya.  Multiscale Diffusion Processes with Periodic Coefficients and an Application to Solute Transport in Porous Media.  Ann. Appl. Probab., 9(4):951–1020, 1999.
  * [5] L. Biferale, A. Crisanti, M. Vergassola, and A. Vulpiani.  Eddy diffusivities in scalar transport.  Phys. of Fluids, 7:2725–2734, 1995.
  * [6] A. Fannjiang and G. Papanicolaou.  Convection Enhanced Diffusion for Periodic Flows.  SIAM J. Appl. Math., 54(2):333–408, 1994.
  * [7] A. Fannjiang and G. Papanicolaou.  Convection-enhanced diffusion for random flows.  J. Stat. Phys., 88(5-6):1033–1076, 1997.
  * [8] G. B. Folland.  Introduction to Partial Differential Equations.  Princeton University Press, Princeton, NJ, 1995.
  * [9] G. B. Folland.  Real Analysis: Modern Techniques and Their Applications.  Wiley–Interscience, New York, NY, 1999.
  * [10] P. Gonnet, S. Güttel, and L. N. Trefethen.  Robust Padé approximation via SVD.  SIAM Review, 55(1):101–117, 2013.
  * [11] A. Majda and P. R. Kramer.  Simplified Models for Turbulent Diffusion: Theory, Numerical Modelling, and Physical Phenomena.  Physics reports. North-Holland, 1999.
  * [12] D. McLaughlin, G. Papanicolaou, and O. Pironneau.  Convection of Microstructure and Related Problems.  SIAM J. Appl. Math., 45:780–797, 1985.
  * [13] R. M. McLaughlin and M. G. Forest.  An anelastic, scale-separated model for mixing, with application to atmospheric transport phenomena.  Phys. of Fluids, 11(4):880–892, 1999.
  * [14] R. C. McOwen.  Partial Differential Equations: Methods and Applications.  Prentice Hall PTR, 2003.
  * [15] N. B. Murphy, E. Cherkaev, J. Xin, J. Zhu, and K. M. Golden.  Spectral analysis and computation of effective diffusivities in space-time periodic incompressible flows.  Ann. Math. Sci. Appl., 2(1):3–66, 2017.
  * [16] N. B. Murphy, E. Cherkaev, J. Zhu, J. Xin, and K. M. Golden.  Spectral analysis and computation for homogenization of advection diffusion processes in steady flows.  Journal of Mathematical Physics, 61(1):013102, 2020.
  * [17] A. Novikov, G. Papanicolaou, and L. Ryzhik.  Boundary layers for cellular flows at high Péclet numbers.  Comm. Pure Appl. Math., 58(7):867–922, 2005.
  * [18] G. A. Pavliotis.  Homogenization theory for advection-diffusion equations with mean flow.  PhD thesis, Rensselaer Polytechnic Institute Troy, New York, 2002.
  * [19] G. A. Pavliotis and A. Stuart.  Multiscale Methods: Averaging and Homogenization.  Texts in Applied Mathematics. Springer New York, 2008.
  * [20] M. C. Reed and B. Simon.  Functional Analysis.  Academic Press, San Diego CA, 1980.
  * [21] H. L. Royden.  Real Analysis.  Prentice-Hall Of India Pvt. Limited, third edition, 1988.
  * [22] I. Stakgold.  Boundary Value Problems of Mathematical Physics.  Classics in Applied Mathematics. SIAM, 2000.  2-volume set.
  * [23] M. H. Stone.  Linear Transformations in Hilbert Space.  American Mathematical Society, Providence, RI, 1964.
  * [24] G. I. Taylor.  Eddy motion in the atmosphere.  Philos. Trans. Royal Soc. A, 215(523–537):1–26, 1915.
  * [25] G. I. Taylor.  Diffusion by continuous movements.  Proceedings of the London Mathematical Society. Third Series, 2:196–211, 1921.