Advective diffusion plays a key role in the transport of salt, heat, buoys, and markers in geophysical flows, in the dispersion of pollutants and trace gases in the atmosphere, and even in the dynamics of sea ice floes influenced by winds and ocean currents. The long time, large scale behavior of such systems is equivalent to an enhanced diffusion process with an effective diffusivity matrix $D*$. Three decades ago, a Stieltjes integral representation for the homogenized matrix, involving a spectral measure of a self-adjoint operator, was developed. However, analytical calculations of $D*$ have been obtained for only a few simple flows, and numerical computations of the effective behavior based on this spectral representation have apparently not been attempted. We overcome these limitations by providing a mathematical foundation for the computation of Stieltjes integral representations of $D*$. We explore two different approaches and for each approach we derive new Stieltjes integral representations and rigorous bounds for the symmetric and antisymmetric parts of $D*$, involving the molecular diffusivity and a spectral measure *μ* of a self-adjoint operator that depends on the characteristics of a randomly perturbed periodic flow. In discrete formulations of each approach, we express *μ* explicitly in terms of standard (or generalized) eigenvalues and eigenvectors of Hermitian matrices. We develop and implement an efficient numerical algorithm that combines beneficial numerical attributes of each approach. We use this method to compute the effective behavior for model flows and relate spectral characteristics to flow geometry and transport properties.

## I. INTRODUCTION

The advective enhancement of diffusive transport of passive scalars by complex fluid flows plays a key role in many important processes in the global climate system^{63} and Earth’s ecosystems.^{15} Advection by geophysical fluids intensifies the dispersion and large scale transport of heat,^{41} pollutants,^{55} and nutrients^{15,24} diffusing in their environment. Advective processes also enhance the large scale transport of plankton,^{24} which is an important component of the food web that sustains life in the polar oceans. The transport of sea ice floes in the polar oceans is diffusive in nature over short time scales as floes bump against each other, yet can be influenced and enhanced by eddy fluxes, storms, and prevailing winds in the atmosphere as well as ocean currents.^{31,51,52,63,64} Thermal transport *through* sea ice, coupling the temperature field in the upper ocean to the lower atmosphere, can be enhanced due to the presence of convective fluid flow in the porous brine microstructure of sea ice.^{21,30,32,62}

It was discovered in the early 1900s^{60} that complex fluid flows transport passive scalars in much the same way as molecular diffusion. The mathematical description of this phenomenon^{61} demonstrated that the long time, large scale behavior of a diffusing particle or tracer being advected by an *incompressible* fluid velocity field is equivalent to an enhanced diffusion process with an effective diffusivity matrix $D*$. Describing the enhancement of the effective transport properties by fluid advection is a challenging problem with theoretical and practical importance in many fields of science and engineering, ranging from turbulent combustion^{2,12,50,59,65,66} to mass, heat, and salt transport in geophysical flows.^{41}

A broad range of mathematical techniques have been developed that reduce the analysis of complex fluid flows, with rapidly varying structures in space and time, to solving averaged or *homogenized* equations that do not have rapidly varying data, but involve effective parameters.^{11,17,18,33–35,48,66} Homogenization of the advection diffusion equation for passive scalar transport by stochastically stationary, time-independent, *mean-zero* fluid velocity fields was treated in Ref. 35. This reduced the analysis of advective diffusion to solving a diffusion equation involving a homogenized temperature field and a constant effective diffusivity matrix $D*$.

An important consequence of this analysis is that the effective diffusivity $D*$ is given in terms of a solution of the so-called “cell problem,”^{35} which can be written as a steady state diffusion equation involving a skew-symmetric random matrix H.^{3,4,17,18} By adapting the analytic continuation method of homogenization theory for composite materials,^{8,22,39} it was shown^{3,4} that the cell problem could be written in resolvent form involving a self-adjoint random operator acting on the Hilbert space of *curl-free vector fields*. This, in turn, led to a Stieltjes integral representation for the *symmetric* part of $D*$, involving the Péclet number *Pe* of the flow and a *spectral measure* of the operator. A key feature of the integral representation for $D*$ is that parameter information in *Pe* is *separated* from the geometry of the fluid velocity field. The velocity field geometry is encoded in the spectral measure through its moments. This parameter separation has led^{3,4,17,18} to rigorous forward bounds for the diagonal components of $D*$.

The mathematical framework developed in Ref. 35 was adapted^{33,48} to the case of a periodic, time-dependent fluid velocity field with *nonzero mean*. It was shown^{48} that, depending on the strength of the fluctuations relative to the mean flow, the effective diffusivity matrix $D*$ can be constant or a function of both space and time. When $D*$ is constant, only its symmetric part appears in the homogenized equation as an enhancement in the diffusivity. However, when $D*$ is a function of space and time, its antisymmetric part also plays a key role in the homogenized equation. Based on Refs. 9 and 10, the cell problem associated with a *time-independent* flow was transformed^{48} into a resolvent formula involving a self-adjoint operator, acting on a Sobolev space^{19,37} of spatially periodic *scalar fields*, which is also a Hilbert space. This, in turn, led to a discrete Stieltjes integral representation for the *antisymmetric* part of $D*$, involving the Péclet number of the flow and a spectral measure of the operator.

Such methods have been extended to steady flows where particles diffuse according to linear collisions,^{49} solute transport in porous media,^{9,10} anelastic (weakly compressible) flows,^{36} as well as to the setting of a time-dependent fluid velocity field.^{5,33,44} All these approaches yield integral representations of the symmetric and, when appropriate, the antisymmetric part of $D*$. Variational formulations of the effective parameter problem for $D*$ are given in Refs. 4, 17, and 18.

We now discuss the main contributions of this work. In Refs. 3 and 4, a Stieltjes integral representation and the associated rigorous bounds were developed for the diagonal components of $D*$ (which are components of the *symmetric* part of $D*$); the off-diagonal components were ignored. In Ref. 48, using a different method, a Stieltjes integral representation was developed for the off-diagonal components of the *antisymmetric* part of $D*$ (the diagonal components are zero); no bounds were developed using the spectral representation. Here, we adapt and extend the mathematical frameworks developed in Refs. 3, 4, and 48 to the case of a time-independent randomly perturbed periodic flow. We obtain *new* Stieltjes integral representations and rigorous bounds for both the symmetric and antisymmetric parts of $D*$—for both its diagonal and off-diagonal components. These integral representations achieve a separation of the geometry of the fluid velocity—through the spectral measures of self-adjoint operators—from the relevant material property of the fluid, its molecular diffusivity.

For each approach, we provide a mathematical foundation for the computation of Stieltjes integral representations for $D*$, by developing discrete, matrix formulations of the effective parameter problem. The spectral measure in each of these two approaches is given explicitly in terms of the eigenvalues and eigenvectors of either a standard or generalized Hermitian eigenvalue problem. In general, the numerical implementation of the generalized eigenvalue problem is more computationally expensive than a standard eigenvalue problem.^{47} However, the approach involving the standard problem involves a matrix that is larger by a factor of the dimension *d* compared to the matrix of the other approach involving the generalized problem. The numerical eigenvalue decomposition of a matrix is quite expensive, and the high resolution of the discretized domain results in very large matrices. Therefore, it is important to develop an efficient way to numerically compute the spectral measure. We provide a detailed matrix analysis that demonstrates both approaches can be formulated in terms of a common *standard* eigenvalue problem involving a matrix with the *smaller size* encountered in the generalized eigenvalue problem, thus combining beneficial numerical attributes of both approaches.

In the continuum setting, the self-adjoint operators used in Refs. 3, 4, and 48 involve the inverse of the negative Laplacian (−Δ)^{−1}, which is given in terms of convolution with the Green’s function for Δ = ∇^{2}.^{56} In the discrete setting, the negative matrix Laplacian is full-rank, hence invertible when Dirichlet boundary conditions are considered, for example. However, the matrix Laplacian is rank-deficient, hence noninvertible when periodic boundary conditions are considered. The matrix analysis discussed above in the full-rank setting reveals useful structure with minimal effort. However, the matrix analysis of the rank-deficient setting is quite a bit more involved, but necessary, since here we investigate advection enhanced diffusion by periodic flows. This analysis demonstrates that the two approaches considered yield equivalent spectral representations for the effective diffusivity matrix $D*$.

We utilize this unified mathematical framework to compute the effective diffusivity matrix $D*$ for some model periodic and randomly perturbed periodic flows and describe the behavior of the enhancement of diffusive transport in terms of the behavior of the spectral measure. There are several approaches to computing $D*$, including numerical solutions of the underlying cell problem partial differential equation,^{48} Monte Carlo methods,^{13} and a method accurate for large Péclet numbers.^{23} Our work here was not motivated by a goal of finding a faster or more accurate method of computing $D*$, although our approach does provide an alternative way of computing $D*$ and is quite robust.

Instead, for randomly perturbed periodic flows, our spectral method for computing $D*$ is set apart from more traditional methods in that it was developed not only to study homogenized behavior but to investigate *spectral statistics*. Indeed, the effective behavior of the system, as encapsulated by $D*$, is closely connected to the statistical behavior of the spectral measure which, in turn, is determined by the statistics of the random eigenvalues and eigenvectors. Consequently, the spectral method enables the homogenization of advection diffusion processes to be viewed through the lens of random matrix theory. In the theory of two-phase composite materials, this approach has led to a new understanding of critical behavior of transport in high contrast composite materials as a percolation threshold is approached. The focus on calculating spectral measures in Stieltjes representations for composite materials through computation of eigenvalues and eigenvectors of random matrices led to the realization that critical behavior of classical transport at a percolation threshold can be viewed as a type of Anderson transition.^{45} In the analysis enabled by the spectral approach for composites, we investigated the appearance, for example, of localization phenomena, mobility edges, and universal Wigner-Dyson statistics of the Gaussian orthogonal ensemble for the eigenvalue spacing distribution. The results in this manuscript lay the mathematical and computational groundwork for such investigations in the context of advection diffusion processes, where the geometry of the fluid velocity field plays the role of the microstructural properties of the composite. Our results in this direction are somewhat outside the scope of this paper and will be published elsewhere, although we consider here some spectral characteristics and their relations to flow geometry and transport properties.

As another example of how this work enables applications to geophysics, we consider advection diffusion within sea ice. The enhancement of the effective thermal conductivity of sea ice due to the presence of convective fluid flow has long been known from an observational perspective.^{21,32,62} The authors are unaware, though, of any predictive, theoretical works on this enhancement. In Ref. 30, the effective thermal conductivity of sea ice in the presence of *bulk* fluid convection is investigated, by applying the results developed here. Using Stieltjes integral representations, a series of bounds on the effective thermal conductivity were obtained for model flows using Pad$e\u0301$ approximants and the analytic continuation method, in terms of the P$e\u0301$clet number.

Motivated by the theoretical findings in this work, in Ref. 44, we generalized the results given here to the setting of a time-dependent fluid velocity field ** u** =

**(**

*u**t*,

**). Furthermore, we used different abstract methods of functional analysis to generalize the equivalence results summarized in Theorem 5, Corollary 6, Lemma 9, and Theorem 10 of this work to the continuum, steady, and dynamic settings.**

*x*In order to streamline the presentation leading to the numerical results, we have placed in Appendixes D and E the development of integral representations for effective diffusivities using the approach introduced in Refs. 3 and 4—for both the continuum and discrete settings. The matrix analysis of the rank deficient setting and the equivalence results discussed above have also been placed in Appendixes F and G. Comments on the notation used throughout this manuscript are given in Appendix B.

## II. HOMOGENIZATION OF THE ADVECTION DIFFUSION EQUATION

The dispersion of a passive scalar with density *ϕ* diffusing in a fluid with molecular diffusivity *ɛ* and being advected by a mean-zero incompressible velocity field ** u** satisfying

*∇*

*⋅***= 0 and ⟨**

*u***⟩ = 0 is described by the advection diffusion equation**

*u*with initial density *ϕ*_{0}(** x**) given. Here, ⟨

**⟩ denotes volume averaging over the period cell $V$,**

*⋅**ϕ*

_{t}denotes partial differentiation of

*ϕ*with respect to time

*t*, Δ =

*∇*

*⋅***= ∇**

*∇*^{2}is the Laplacian,

*ɛ*> 0, and

*d*is the system dimension, and we denote by 0 the null element in all linear spaces in question. Moreover,

*ξ*

*⋅***=**

*ζ*

*ξ*^{†}

**and † is the operation of complex-conjugate-transpose, with**

*ζ*

*ξ*

*⋅***= |**

*ξ***|**

*ξ*^{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*. We assume for now that the time-independent fluid velocity field

**is spatially periodic on the region $V\u2282Rd$. Later, we will discuss the case of an array of randomly perturbed, periodic flows.**

*u*The long time, large scale dispersion of passive scalars can be described^{61} by an effective diffusivity matrix $D*$. An explicit representation for $D*$ can be found using methods of homogenization theory.^{7,46} These methods have demonstrated that the averaged or *homogenized* behavior of the advection diffusion equation in (1) is determined by a diffusion equation involving an averaged scalar density $\varphi \xaf$ and a (constant) effective diffusivity matrix $D*$,^{35}

The components $Djk*$, *j*, *k* = 1, …, *d*, of $D*$ are given by^{35}

The function *χ*_{j} in (3) satisfies a cell problem which is a steady state advection diffusion equation with a forcing term involving *u*_{j}, the *j*th component of the fluid velocity field ** u**,

^{18,35}

Equations (2)–(4) follow from the assumption that the length scale associated with spatial variations of the initial density *ϕ*_{0} is much larger than the length scale of spatial variations associated with the fluid velocity field *u*^{18,35} (separation of scales). This information is incorporated into Eq. (1) by introducing a small dimensionless parameter *δ* ≪ 1 and writing^{35}

Anticipating that *ϕ* will have diffusive dynamics as *t* → *∞*, space and time are rescaled by ** x** ↦

**/**

*x**δ*and

*t*↦

*t*/

*δ*

^{2}. As

*δ*→ 0, the associated solution

*ϕ*

^{δ}(

*t*,

**) =**

*x**ϕ*(

*t*/

*δ*

^{2},

**/**

*x**δ*) of Eq. (1) (in the rescaled variables) converges to $\varphi \xaf(t,x)$ which satisfies Eq. (2). The convergence is in an

*L*

^{2}sense that depends on the technical assumptions made about the fluid velocity field

**.**

*u*^{4,17,18,33,35,48}

We emphasize that the cell problem in (4) involves only the fast variable ** x**/

*δ*.

^{35}Other space-time scalings have also been considered, which have led to space-time dependent $D*$

^{48}and even anomalous diffusive dynamics.

^{33}Homogenization theorems for space-time dependent fluid velocity fields are treated in Refs. 11, 33, and 48.

In our analysis of the effective diffusivity matrix $D*$, it is beneficial to use nondimensional parameters. We therefore assume that Eq. (1) has been nondimensionalized as follows: Let *ℓ* and *τ* be the typical length and time scales associated with the problem of interest. Mapping to the nondimensional variables *t* ↦ *t*/*τ* and ** x** ↦

**/**

*x**ℓ*, one finds that

*ϕ*satisfies the advection diffusion equation in (1) with a nondimensional molecular diffusivity and fluid velocity field,

This demonstrates that by nondimensionalizing Eq. (1), the fluid velocity field ** u** is, in turn, divided by a quantity with dimensions of velocity and the molecular diffusivity is divided by a quantity with dimensions of velocity multiplied by the spatial length. It is convenient to choose the rescaled

**and**

*u**ɛ*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

*ɛ*from the

*geometry*of the flow; this leads to mathematically and physically meaningful properties of rigorous bounds for $D*$ which follow from the analytic structure of Stieltjes integral representations for $D*$

^{4,6}—discussed in Sec. III.

We accomplish both of these goals as follows: Define the dimensional fluid velocity field by ** u** =

*u*

_{0}

**, where the parameter**

*v**u*

_{0}has dimensions of velocity and represents the “

*flow strength*” of

**which is independent of the geometry of the flow; the flow geometry is encapsulated in the nondimensional vector field**

*u***. With these definitions, we choose reference scales**

*v**τ*and

*ℓ*in Eq. (6) to satisfy

*u*

_{0}=

*ℓ*/

*τ*so that

**↦**

*u***and**

*v**ɛ*↦

*ɛ*/

*u*

_{0}

*ℓ*. For example, in Sec. V, we compute $D*$ for

*BC*-flow

^{11}having dimensional fluid velocity field

**=**

*u**u*

_{0}(

*C*cos

*y*,

*B*cos

*x*), where the flow strength

*u*

_{0}∈ (0,

*∞*) is independent of the nondimensional parameters

*B*,

*C*∈ [0, 1] which determine the streamline geometry of

**.**

*u*An example of a nondimensional, parameter that compares the rate of scalar advection to the rate of diffusion is the Péclet number. We define it by the ratio *Pe* = *ℓu*_{0}/*ɛ*, although other definitions have been used.^{33,36} The advection and diffusion dominated regimes are characterized by *Pe* ≫ 1 and *Pe* ≪ 1, respectively. Therefore, our choice of the rescaled *ɛ* satisfies *Pe* = 1/*ɛ*.

The *parameter separation* between *Pe* and the geometry of the flow is important for rigorous upper and lower Padé approximant bounds for $D*$.^{4,43} Padé approximants of $D*$ are given in terms of ratios of polynomials^{6} *P*(*z*)/*Q*(*z*), where *z* = *Pe*^{2}, 0 < *z* < *∞*, and the coefficients of these polynomials depend on the moments of a *spectral measure* that, in turn, depend on the fluid velocity field ** u**.

^{4,43}For example, when

**is given by**

*u**BC*-flow, the moments of the measure depend

^{43}on the parameters

*B*and

*C*. Our numerical investigations have shown if the nondimensionalization of Eq. (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

^{43}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 $D*$, which is not physically or mathematically consistent with the known behavior of $D*$.

^{11,17,33,48}This demonstrates the importance of

*parameter separation*between

*z*and the flow geometry for Padé approximant bounds for $D*$.

This way of nondimensionalizing Eq. (1) is also convenient in the case of a time-dependent fluid velocity field,^{44} where the parameter *u*_{0} again represents the flow strength and the vector field ** v** encapsulates the

*geometric and dynamical*properties of the flow. For example, the space-time periodic flow with velocity field

**=**

*u**u*

_{0}((

*C*cos

*y*,

*B*cos

*x*) + cos

*t*(

*γ*sin

*y*,

*β*sin

*x*)) has dynamical behavior exhibiting Lagrangian chaos.

^{11,44}Here, the flow strength

*u*

_{0}∈ (0,

*∞*) is independent of the parameters

*B*,

*C*,

*γ*,

*β*∈ [0, 1] which determine the geometric and dynamical properties of

**. This choice of nondimensionalization gives a clearer interpretation of the advection and diffusion dominated regimes in terms of**

*u**Pe*= 1/

*ɛ*than that given in Ref. 44. A detailed discussion of various nondimensionalizations of Eq. (1) is given in Refs. 33 and 36.

## III. HILBERT SPACE AND INTEGRAL REPRESENTATIONS

In this section, we adapt and extend a method^{9,10,48} which provides Stieltjes integral representations for the effective diffusivity matrix $D*$. We do so by providing functional formulas for the symmetric $S*$ and antisymmetric $A*$ parts of $D*$, involving the scalar field *χ*_{j} in Eq. (4). We also provide a Sobolev space formulation of the effective parameter problem^{48} which yields a resolvent formula for *χ*_{j}, involving a self-adjoint operator that depends only on the fluid velocity field ** u**. This and the spectral theorem

^{53,58}yield Stieltjes integral representations for $S*$ and $A*$ involving a spectral measure of the operator.

Consider the Hilbert space $H$,

where *ν* is the Lebesgue measure on $Rd$, restricted to $V$, and the *σ*-algebra associated with the underlying measure space is generated by the Lebesgue measurable subsets of $Rd$. The Hilbert space $H$ is equipped with a sesquilinear inner product ⟨** ⋅**,

**⟩ defined by $\u27e8f,h\u27e9=\u27e8\u2009f\xaf\u2009h\u2009\u27e9$, with $\u27e8h,f\u27e9=\u27e8f,h\u27e9\xaf$ for $\u2009f,h\u2208H$, which induces a norm ∥**

*⋅***∥ defined by ∥**

*⋅**f*∥ = ⟨

*f*,

*f*⟩

^{1/2}and $f\u2208H$ implies ∥

*f*∥ <

*∞*.

One could also consider a random fluid flow filling all of $Rd$, with a velocity field ** u** determined by the probability space (Ω,

*P*) with

*σ*-algebra generated by the sets {

**(**

*u***) ∈**

*x**B*}, where $x\u2208Rd$ and

*B*is a Borel subset of $Rd$.

^{4}Here, Ω is the set of all geometric realizations of

**, which is indexed by the parameter**

*u**ω*∈ Ω representing one particular geometric realization, and

*P*is the associated probability measure. The underlying Hilbert space in this case can be taken to be $H=L2(\Omega ,P)$, i.e., the space of all

*P*-measurable complex-valued scalar functions

*ξ*satisfying $\Vert \xi \Vert =\u27e8|\xi |2\u27e91/2<\u221e$, where ⟨

**⟩ denotes ensemble averaging and the underlying sesquilinear inner-product is defined by $\u27e8\xi ,\zeta \u27e9=\u27e8\u2009\xi \xaf\u2009\zeta \u27e9$. In this case, one could consider a random fluid flow with a velocity field**

*⋅***that is stationary**

*u*^{35}or ergodic,

^{3,4}for example, with regularity conditions at infinity, i.e., as |

**| →**

*x**∞*. In these cases, one works with an infinite medium directly, which presents substantial computational difficulties.

A more computationally tractable random system is given by an *n* × *n* array of randomly perturbed periodic flows.^{18} In this case, the *σ*-algebra associated with the underlying probability space is generated by the Lebesgue measurable subsets of $Rd$. Here, the Hilbert space $H$ is given by Eq. (7) and the averaged quantities depend on the realization of the random medium because ⟨⋅⟩ is given by volume averaging over the period cell $V$.^{18} The effective diffusivity matrix $D*$ is obtained by taking an infinite volume limit, $D*=limn\u2192\u221eDn*\u2009$, of the finite volume effective diffusivity matrix $Dn*$ and evoking an ergodic theorem.^{18,22} Numerically, by the law of large numbers,^{16} it is natural to spatially average each statistical trial and then ensemble average over all of the sampled random realizations. This is the approach we take here.

In any case, once the Hilbert space $H$ is established, with associated average ⟨⋅⟩, inner-product ⟨⋅, ⋅⟩, and norm ∥⋅∥, the spectral theory presented in the remainder of this section progresses independent of the underlying details, as it lays on an axiomatic foundation.^{58} For the sake of numerical tractability, we will assume that the Hilbert space $H$ is given by Eq. (7). The fluid velocity field ** u** can be assumed to represent a periodic or randomly perturbed periodic flow. Now, consider the Sobolev space $H1,2$, which itself is a Hilbert space,

^{9,10}

where ∥⋅∥_{1,2} is the norm induced by the underlying sesquilinear inner-product ⟨⋅,⋅⟩_{1,2} defined by ⟨*f*,*h*⟩_{1,2} = ⟨*∇**f* ⋅ *∇**h*⟩, with $\u27e8h,f\u27e91,2=\u27e8f,h\u27e9\xaf1,2\u2009$ (recall ** ξ** ⋅

**=**

*ζ*

*ξ*^{†}

**includes complex conjugation).**

*ζ*Recall the definition of the components $Djk*=\epsilon \delta jk+\u27e8uj\chi k\u27e9$, *j*, *k* = 1, …, *d*, of the effective diffusivity matrix $D*$ in (3). Rewrite the functional ⟨*u*_{j}*χ*_{k}⟩ as^{48}

where we have integrated by parts and used the periodicity of the functions *u*_{j} and *χ*_{k}. Here, Δ^{−1} is based on convolution with the Green’s function for the Laplacian Δ.^{56} The symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$ are defined by

Substituting into Eq. (9) the expression − *u*_{j} = ** u** ⋅

*∇**χ*

_{j}+

*ɛ*Δ

*χ*

_{j}for −

*u*

_{j}in (4) yields the following functional formulas for the components $Sjk*$ and $Ajk*$,

*j*,

*k*= 1, …,

*d*, of $S*$ and $A*$,

Due to the incompressibility of the fluid velocity field, ** ∇** ⋅

**= 0, the operator**

*u**A*is antisymmetric on $H1,2$, i.e., ⟨

*Aξ*,

*ζ*⟩

_{1,2}= −⟨

*ξ*,

*Aζ*⟩

_{1,2}for all $\xi ,\zeta \u2009\u2208H1,2$ [see Eq. (C1)]. Since the scalar fields

*χ*

_{j}and

*Aχ*

_{j}are

*real-valued*, we have that $\u27e8\chi j,\chi k\u27e91,2=\u27e8\chi k,\chi j\u27e91,2$ and $Akj*=\u27e8A\chi k,\chi j\u27e91,2=\u2212\u27e8\chi k,A\chi j\u27e91,2=\u2212\u27e8A\chi j,\chi k\u27e91,2=\u2212Ajk*$, confirming that $S*$ is symmetric and $A*$ is antisymmetric, with $Akk*=0$.

Applying the operator Δ^{−1} to both sides of the cell problem *ɛ*Δ*χ*_{j} + ** u** ⋅

*∇**χ*

_{j}= −

*u*

_{j}in Eq. (4) yields the following resolvent formula for

*χ*

_{j}involving the operator

*A*in (11):

Substituting the resolvent formula for *χ*_{j} in (12) into the bilinear functionals in Eq. (11) yields

Since $V$ is a bounded domain, the linear operator Δ^{−1} is bounded on $H$ with bounded operator norm ∥Δ^{−1}∥ < *∞*.^{56} When |** u**| is uniformly bounded on the period cell $V$, $supx\u2208V|u(x)|<\u221e$, the linear operator

*A*is bounded on $H1,2$, with [see Eqs. (C2) and (C3)]

where ∥*A*∥_{1,2} denotes the operator norm of *A* induced by the norm ∥⋅∥_{1,2} in Eq. (8). All of the fluid velocity fields that we consider in our numerical computations of $D*$ in Sec. V are uniformly bounded. More generally, for $uk\u2208H$, *k* = 1, …, *d*, the operator *A* is compact on $H1,2$,^{9,10,48} hence bounded.^{56} Consequently, *M* = −*ıA*, where $\u0131=\u22121$, is a bounded linear operator on $H1,2$ with ∥*M*∥_{1,2} = ∥*A*∥_{1,2} < *∞*. Moreover, the skew-symmetry of *A* and the sesquilinearity of the $H1,2$-inner-product imply that *M* is also symmetric, ⟨*Mf*, *h*⟩_{1,2} = ⟨*f*, *Mh*⟩_{1,2}. The bounded linear symmetric operator *M* with domain $H1,2$ is *self-adjoint*.^{53,58}

The spectrum Σ of the self-adjoint operator *M* is real-valued, with the spectral radius equal to its operator norm,^{53} i.e.,

The spectral theorem for bounded linear self-adjoint operators in Hilbert space states that there is a one-to-one correspondence between the operator *M* and a family of self-adjoint projection operators {Q(*λ*)}_{λ∈Σ}—the resolution of the identity—satisfying^{58}

Furthermore, for all $\xi ,\zeta \u2009\u2208H1,2$, the *complex-valued* function of the spectral variable *λ* defined by *μ*_{ξζ}(*λ*) = ⟨Q(*λ*)*ξ*,*ζ*⟩_{1,2} has real and imaginary parts that are of bounded variation.^{58} Therefore, there is a *complex* Stieltjes measure *μ*_{ξζ} associated with *μ*_{ξζ}(*λ*).^{20,57,58}

The spectral theorem also states, for all complex-valued functions *f* ∈ *L*^{2}(*μ*_{ξξ}) and *h* ∈ *L*^{2}(*μ*_{ζζ}), there exist linear operators denoted by *f*(*M*) and *h*(*M*) which are defined in terms of the bilinear functional ⟨*f*(*M*)*ξ*, *h*(*M*) *ζ*⟩_{1,2}.^{58} In particular, this functional has the following Radon–Stieltjes integral representation involving the Stieltjes measure *μ*_{ξζ}, for all $\xi ,\zeta \u2009\u2208H1,2$:

where the integration is over the spectrum Σ of *M*^{53,58} and $f\xaf$ denotes the complex conjugation of the scalar function *f*. Since the Stieltjes measure *μ*_{ξζ} has the property^{58} $\u222babd\mu \xi \zeta (\lambda )=\mu \xi \zeta (b)\u2212\mu \xi \zeta (a)$, Eq. (16) implies that the mass $\mu \xi \zeta 0$ of the measure *μ*_{ξζ} is given by

which is bounded in the sense that $|\mu \xi \zeta 0|\u2264\Vert \xi \Vert 1,2\u2009\Vert \zeta \Vert 1,2<\u221e$ for all $\xi ,\zeta \u2009\u2208H1,2$. Due to the sesquilinearity of the inner-product and the fact that the projection operator Q(*λ*) is self-adjoint on $H1,2$, the complex-valued function *μ*_{ζξ}(*λ*) satisfies $\mu \zeta \xi (\lambda )=\mu \xaf\xi \zeta (\lambda )$ and *μ*_{ξξ}(*λ*) = ∥Q(*λ*)*ξ*∥^{2}, so *μ*_{ξξ} is a *positive* measure. The formulas in Eqs. (16)–(18) will be used several times throughout this manuscript to clarify and streamline the development of various Stieltjes integral representations for the effective diffusivity matrix $D*$. They will be used both in the continuum setting and the discrete setting, which leads to an efficient numerical algorithm for our computation of spectral representations for $D*$ for various model fluid flows in Sec. V.

We are now ready to present the main results of this section. For notational simplicity, denote the complex-valued function *μ*_{jk}(*λ*) = ⟨Q(*λ*)*g*_{j}, *g*_{k}⟩, instead of $\mu gjgk(\lambda )$, where *g*_{j} = (−Δ)^{−1}*u*_{j} is defined in (12). Denote the real and imaginary parts of the function *μ*_{jk}(*λ*) by Re *μ*_{jk}(*λ*) and Im *μ*_{jk}(*λ*), respectively.

*Let*

*Q(λ)*

*denote the resolution of the identity corresponding to the self-adjoint operator*

*M*

*. Then, the components*$Sjk*$

*and*$Ajk*$

*,*

*j, k = 1, …, d*

*, of the symmetric*$S*$

*and antisymmetric*$A*$

*parts of the effective diffusivity matrix*$D*$

*have the following Radon-Stieltjes integral representations:*

*For*

*j ≠ k*

*,*

*Re*

*μ*

_{jk}

*and*

*Im*

*μ*

_{jk}

*are signed measures.*

^{20}*For*

*j = k*

*,*

*Re*

*μ*

_{kk}

*= μ*

_{kk}

*is a positive measure*

^{20}*, which demonstrates that the effective transport of the scalar density*

*ϕ*

*is enhanced by the presence of an incompressible fluid velocity field*

*u**, above the bare diffusive value*

*ɛ*,

*The mass*$\mu jk0$

*of the measure*

*μ*

_{jk}

*is real-valued and satisfies*

We first provide integral representations for the bilinear functional formulas in Eq. (13) for $Sjk*$ and $Ajk*$. Since we have already established the operator *M* is self-adjoint, we only need to identify the appropriate functions *f*(*λ*) and *h*(*λ*) as well as the Hilbert space members $\xi ,\zeta \u2009\u2208H1,2$ in (17). Using *A* = *ıM*, comparison of the functionals in Eqs. (13) and (17) identifies *f*(*λ*) = *h*(*λ*) = (*ɛ* + *ıλ*)^{−1} for the first formula in (13), while *f*(*λ*) = *ıλ*(*ɛ* + *ıλ*)^{−1} and *h*(*λ*) = (*ɛ* + *ıλ*)^{−1} for the second formula, with *ξ* = *g*_{j} and *ζ* = *g*_{k} for both of these formulas.

*f*,

*h*∈

*L*

^{2}(

*μ*

_{kk}) for all

*k*= 1, …,

*d*and 0 <

*ɛ*<

*∞*. From Eq. (18), the mass $\mu jk0$ of the measure

*μ*

_{jk}is given by $\mu jk0=\u27e8gj,gk\u27e91,2=\u27e8\u2207\Delta \u22121uj\u22c5\u2207\Delta \u22121uk\u27e9$. Integration by parts and the Cauchy-Schwartz inequality

^{20}then yield Eq. (21). In particular, $|\mu kk0|\u2264\Vert \Delta \u22121\Vert \Vert uk\Vert 2<\u221e$, as Δ

^{−1}has bounded operator norm ∥Δ

^{−1}∥ on $L2(V)$.

^{56}Consequently, since $0<(\epsilon 2+\lambda 2)\u22121\u22641/\epsilon 2<\u221e$ and $0<\lambda 2(\epsilon 2+\lambda 2)\u22121<1$ for all 0 <

*ɛ*<

*∞*, it is clear that

*f*,

*h*∈

*L*

^{2}(

*μ*

_{kk}). Consequently, the spectral theorem in (17) implies that the functional formulas for $Sjk*$ and $Ajk*$ in Eq. (13) have the following Stieltjes integral representations:

*complex measure*

*μ*

_{jk}.

*signed measures*Re

*μ*

_{jk}and Im

*μ*

_{jk}. Since

*ɛ*,

*χ*

_{k},

*u*

_{k}, and

*g*

_{k}in (13) are

*real-valued*, we have from (12) the following symmetry conditions:

^{58}of the Stieltjes integrals in Eq. (22) with respect to the function

*μ*

_{jk}(

*λ*), and the two identities $Re\u2009\mu jk(\lambda )=(\mu jk(\lambda )+\mu \xafjk(\lambda ))/2$ and $Im\u2009\u2009\mu jk(\lambda )=(\mu jk(\lambda )\u2212\mu \xafjk(\lambda ))/(2\u0131)$ yield Eq. (19). This concludes our Proof of Theorem 1.

*Let*Σ

*denote the spectrum of the operator*

*M*

*, denote*

*λ*

_{+}= sup

_{λ∈Σ}

*|λ|*

*, and assume*

*0 < ɛ < ∞*

*. Then, the diagonal components*$Skk*=Dkk*$

*of the effective diffusivity matrix satisfy the following, rigorous upper*

^{4}*and lower bounds:*

*For*

*j ≠ k*

*, let*$Re\u2009\mu jk=Re\u2009\mu jk+\u2212Re\u2009\mu jk\u2212$

*and*$Im\u2009\u2009\mu jk=Im\u2009\mu jk+\u2212Im \u2009\mu jk\u2212$

*denote the Jordan decomposition of the signed measures*

*Re*

*μ*

_{jk}

*and*

*Im*

*μ*

_{jk}

*, respectively, and denote the total variation of these measures by*$|Re\u2009\mu |jk=Re\u2009\mu jk++Re\u2009\mu jk\u2212$

*and*$|Im \u2009\mu |jk=Im\u2009\u2009\mu jk++Im \u2009\mu jk\u2212$

*, respectively. Finally, denote the masses of the measures*$Re\u2009\mu jk+$

*,*$Re\u2009\mu jk\u2212$

*, and*

*|Im*

*μ|*

_{jk}

*by*$[Re\u2009\mu jk+]0$

*,*$[Re\u2009\mu jk\u2212]0$

*, and*$|Im\u2009\mu |jk0$

*, respectively. Then, the off-diagonal components*$Sjk*$

*and*$Ajk*$

*,*

*j ≠ k*

*, satisfy the following, rigorous upper and lower bounds:*

Assume that 0 < *ɛ* < *∞*. From Eq. (15), the spectrum Σ of the *compact*^{9,10,33} operator *M* = −*ıA* is a bounded subset of $R$. The spectrum is discrete^{56} away from the spectral origin *λ* = 0 and comes in ± pairs^{48} with an accumulation point at *λ* = 0.^{56} Denote *λ*_{+} = sup_{λ∈Σ}|*λ*| and note that inf_{λ∈Σ}*λ*^{2} = 0. The inequalities $1/(\epsilon 2+\lambda +2)\u22641/(\epsilon 2+\lambda 2)\u22641/\epsilon 2$ hold for all *λ* ∈ Σ. Consequently, since *μ*_{kk} is a positive measure with finite mass $\mu kk0$, the inequalities in (24) hold.^{20} It may happen that $\mu kk0=0$, hence $Skk*=\epsilon $, e.g., shear flow orthogonal to the *k*th direction.^{4,17}

*j*≠

*k*, Re

*μ*

_{jk}is a signed measure. By the Jordan decomposition of Re

*μ*

_{jk}, there are unique,

*positive*measures $Re\u2009\mu jk+$ and $Re\u2009\mu jk\u2212$ such that $Re\u2009\mu jk=Re\u2009\mu jk+\u2212Re\u2009\mu jk\u2212\u2009$. Moreover, associated with the signed measure, Re

*μ*

_{jk}is its

*total variation*|Re

*μ*|

_{jk},

^{20}

*μ*|

_{jk}is also bounded. Since we have

^{20}that $|Sjk*|\u2264\u222bd|Re\u2009\mu |jk(\lambda )/(\epsilon 2+\lambda 2)$, the upper bound in Eq. (24) with $\mu kk0$ replaced by $|Re\u2009\mu |jk0$ holds for the positive quantity $|Sjk*|$. These bounds for $Sjk*$ can be improved upon by separately considering the positive and negative contributions of the integral representation for $Sjk*$, yielding Eq. (25). In a similar way, we obtain the bounds for $Ajk*$ in Eq. (26). This concludes our Proof of Theorem 2.

We conclude this section by noting that bounds on $Skk*$ can also be obtained using variational methods^{4,17,18} as well as Padé approximants^{4,6} of Stieltjes functions.

## IV. DISCRETE SETTING: SOBOLEV SPACE OF SCALAR FIELDS

For our numerical computations of the effective diffusivity $D*$ in Sec. V, we consider a discrete approximation of the cell problem in Eq. (4) written as (*ɛ* + *ıM*)*χ*_{j} = *g*_{j}. Here, *M* = −*ıA*, *A* = Δ^{−1}[** u** ⋅

**], and**

*∇**g*

_{j}= −Δ

^{−1}

*u*

_{j}, as defined in Eqs. (11) and (12). In this section, we manipulate these formulas in order to develop a numerical algorithm which enables numerical computations of $D*$ by directly computing a discrete representation of the spectral measure

*μ*

_{jk}in Eq. (19), in terms of the eigenvalues and eigenvectors of a

*generalized*eigenvalue problem. This is not a trivial extension of the spectral theory for the continuum setting discussed in Sec. III, as the matrix representation of the operator (−Δ)

^{−1}[

**⋅**

*u***] is**

*∇**not*Hermitian. The special structure of the generalized eigenvalue problem itself makes these matrix operators Hermitian only with respect to a discrete Sobolev-like inner-product analogous to the inner-product ⟨⋅,⋅⟩

_{1,2}introduced after Eq. (8). Moreover, the eigenvector orthogonality is only with respect to this inner-product.

Towards this goal, we begin by noting that since ** u**(

**) is incompressible,**

*x***⋅**

*∇***= 0, there is a real (nondimensional)**

*u**antisymmetric*matrix H(

**)**

*x*^{3,4}such that

where H^{T} denotes transposition of the matrix H. Due to the antisymmetry of the matrix H in Eq. (28) and the symmetry of the Hessian operator ** ∇∇** when acting on a sufficiently smooth space of functions, we have H

*: ∇∇**φ*= 0 for all such smooth functions

*φ*, where

**denotes matrix contraction. Consequently,**

*:***⋅ [H**

*∇*

*∇**φ*] = [

**⋅ H] ⋅**

*∇*

*∇**φ*+ H

*: ∇∇**φ*= [

**⋅ H] ⋅**

*∇*

*∇**φ*, or

as operators acting on such functions. With (29) we can write the operator *M* = *ı*(−Δ)^{−1}[** u** ⋅

**] as**

*∇**M*=

*ı*(−Δ)

^{−1}

**⋅ [H**

*∇***].**

*∇*We now discuss our discrete formulation of the effective parameter problem, which leads to Stieltjes integral representations for the symmetric $S*$ and antisymmetric $A*$ components of the effective diffusivity matrix $D*$, involving a discrete spectral measure. For notational brevity, assume that the period cell $V$ is square. In the discrete setting,^{42} $V$ is represented by a square grid of size *L*, which is bijectively mapped to a vector with *L*^{d} components. The functions *u*_{j}(** x**) and

*χ*

_{j}(

**) are mapped to vectors**

*x*

*u*_{j}and

*χ*_{j}with

*L*

^{d}components, respectively, and the matrix H(

**) is mapped to a square banded antisymmetric matrix of size**

*x**N*=

*L*

^{d}

*d*(see Sec. V A for details). For simplicity, we will not make a notational distinction for H between the discrete and continuum settings, as the context will be clear.

The differential operator ** ∇** is represented by a finite difference matrix ∇,

^{14,42}where $\u2207T=(\u22071T,\u2026,\u2207dT)$ and ∇

_{j},

*j*= 1, …,

*d*, are also finite difference matrices. Moreover, the divergence operator

**⋅ is given by −∇**

*∇*^{T}and the matrix representation of the negative Laplacian −Δ is given by ∇

^{T}∇. Consequently, we may write the discrete, matrix representation M of the self-adjoint operator

*M*=

*ı*(−Δ)

^{−1}

**⋅ [H**

*∇***] as $M=(\u2207T\u2207)\u22121[\u2212\u0131\u2207TH\u2207]$. This composition of the Hermitian matrix −**

*∇**ı*∇

^{T}H∇ and the real-symmetric matrix $(\u2207T\u2207)\u22121$ is neither real-symmetric nor Hermitian symmetric. From Eq. (C1), we see that the symmetry properties of the integro-differential operator

*M*depend intimately on the inner-product ⟨

*f*,

*h*⟩

_{1,2}= ⟨

*∇**f*⋅

*∇**h*⟩ of the underlying Sobolev space $H1,2$ defined in Eq. (8). Therefore, we anticipate that the properties of this inner-product must be incorporated into the discrete formulation of integral representations for $D*$.

We now provide a matrix formulation of the effective parameter problem introduced in Sec. III, which involves a *generalized eigenvalue problem* that has the Sobolev-type inner-product as a central feature. In particular, this formulation retains the key properties of the weak form of the eigenvalue problem $\u27e8M\phi n,\phi n\u27e91,2=\lambda n$. Namely, the operator *M* is self-adjoint *with respect to the inner-product* ⟨⋅, ⋅⟩_{1,2}, its eigenfunctions $\phi n\u2208H1,2$ are orthonormal $\u27e8\phi n,\phi m\u27e91,2=\delta nm$, *n*, *m* = 1, 2, 3, …, with respect to the inner-product ⟨⋅, ⋅⟩_{1,2}, and the spectrum Σ of *M* is real valued, $\Sigma \u2282R$. Toward this goal, consider the eigenvalue problem *Mφ*_{n} = *λ*_{n}*φ*_{n} in the following “strong” form:

Using a discrete version of Eq. (30), we will establish the integral representations in (19) for discrete versions of the functionals $Sjk*=\epsilon (\delta jk+\u27e8\chi j,\chi k\u27e91,2)$ and $Ajk*=\u27e8\u0131M\chi j,\chi k\u27e91,2$ in (11), involving a discrete spectral measure.

By our discussion in this section, the matrix representation of the eigenvalue problem in (30) is

The first formula in Eq. (31) is a *generalized eigenvalue problem*^{47} associated with the pencil B − *λ*C, where B and C are Hermitian and real-symmetric matrices, respectively, of size *K* = *L*^{d}. The *λ*_{n} and *z*_{n}, *n* = 1, …, *K*, are known as generalized eigenvalues and eigenvectors, respectively. The matrix C = ∇^{T}∇ is clearly positive semidefinite. In this section, we will assume that C is positive definite, hence invertible. We will discuss the case where C is positive semidefinite in Appendix G, which is necessary for the setting where ∇ incorporates *periodic* boundary conditions—needed for our study of advection enhanced diffusion by a spatially periodic fluid velocity field ** u**.

Since B and C are Hermitian and real-symmetric, respectively, and C is positive definite, the generalized eigenvalue problem in (31) has properties which are similar to the properties of the standard symmetric eigenvalue problem.^{47} In particular, the generalized eigenvalues *λ*_{n} are all real, the generalized eigenvectors *z*_{n} form a basis for $CK$, and the *z*_{n} are orthonormal in the following sense: $zn\u2009\u2020Czm=\delta nm$, *n*, *m* = 1, …, *K*.^{47} Since C = ∇^{T}∇ is real-valued, this is equivalent to the Sobolev-type orthogonality condition,

In other words, the generalized eigenvectors *z*_{n} are orthonormal with respect to the “discrete inner-product” ⟨⋅, ⋅⟩_{1,2} defined by ⟨** ξ**,

**⟩**

*ζ*_{1,2}= ⟨∇

**⋅ ∇**

*ξ***⟩, for $\xi ,\zeta \u2208CK$ and ⟨⋅⟩ denotes ensemble averaging. Denoting by Z the matrix with columns consisting of the generalized eigenvectors**

*ζ*

*z*_{n}, Eq. (32) is seen to be equivalent to [∇Z]

^{†}[∇Z] = I, or Z

^{†}CZ = I. A key feature of the generalized eigenvalue problem is that the matrix Z simultaneously diagonalizes B and C. Specifically, if Λ is the diagonal matrix whose elements on the main diagonal are the generalized eigenvalues

*λ*

_{n}, then

^{47}

We now derive the discrete version of Eqs. (16)–(18) comprising the key results of the spectral theorem, for the generalized eigenvalue problem setting. These derivations will streamline and clarify our results for the current section, showing how the derived series representations of the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$ can be written as the Stieltjes integrals in Eq. (19), involving a discrete spectral measure. This derivation will also clarify and streamline our results in Appendixes E, F, and G, which lead to the numerical algorithm used in Sec. V for spectral computations of $D*$ for periodic fluid flows. The numerical algorithm used in Sec. V is analogous to the algorithm that we develop in this section, which is elegant for the full-rank setting—revealing a great deal of structure with minimal effort. The matrix analysis of the rank-deficient setting developed in Appendix G is quite a bit more involved. Developing the full-rank setting in this section first makes the results of the rank-deficient setting more transparent and the final results more anticipated.

Since the *z*_{n}, *n* = 1, …, *K*, form a basis for $CK$ and satisfy the orthogonality relation in (32), for all $\xi \u2208CK$, we have $\xi =\u2211n(\u2207zn\u22c5\u2207\xi )zn=\u2211n(zn[\u2207zn]\u2009\u2020\u2207)\xi $, which implies

where the matrix operators Q_{n}, *n* = 1, …, *K*, are self-adjoint with respect to the *discrete* inner-product ⟨⋅, ⋅⟩_{1,2} defined above, i.e., $\u27e8Qn\xi ,\zeta \u27e91,2=\u27e8\xi ,Qn\zeta \u27e91,2$ for all $\xi ,\zeta \u2208CK$. From B*z*_{n} = *λ*_{n}C*z*_{n} and Eq. (34), we have that BQ_{n} = *λ*_{n}CQ_{n} which is equivalent to C^{−1}BQ_{n} = *λ*_{n}Q_{n}, since the matrix C is invertible. This formula and (34) then imply that the matrix C^{−1}B has the spectral decomposition C^{−1}B = *∑*_{n}*λ*_{n}Q_{n}. By the mutual orthogonality of the projection matrices Q_{n} and by induction, we have that $[C\u22121B]m=\u2211n\lambda nmQn$ for all $m\u2208N$. This, in turn, implies that *f*(C^{−1}B) = *∑*_{n}*f*(*λ*_{n})Q_{n} for any polynomial $f:R\u21a6C$.

From the mutual orthogonality of the projection matrices Q_{n} and their symmetry with respect to the discrete inner-product ⟨⋅, ⋅⟩_{1,2} it follows that, for all $\xi ,\zeta \u2208CK$ and all complex-valued polynomials *f*(*λ*) and *h*(*λ*), the bilinear functional $\u27e8f(C\u22121B)\xi ,h(C\u22121B)\zeta \u27e91,2$ has the integral representation in (17), with *M* substituted by C^{−1}B and other appropriate notational changes. Moreover, the complex-valued function *μ*_{ξζ}(*λ*) = ⟨Q(*λ*)*ξ*,*ζ*⟩_{1,2} in (17) is now given by *μ*_{ξζ}(*λ*) = ⟨Q(*λ*)** ξ**,

**⟩**

*ζ*_{1,2}, where the associated matrix representation Q(

*λ*) of the projection operator Q(

*λ*) and the discrete spectral measure d

*μ*

_{ξζ}(

*λ*) are given by

Here, *θ*(*λ*) is the Heaviside function, satisfying *θ*(*λ*) = 0 for *λ* < 0 and *θ*(*λ*) = 1 for *λ* ≥ 0, and $\delta \lambda n(d\lambda )=d\theta (\lambda \u2212\lambda n)$ is the *δ*-measure centered at *λ*_{n}. From Eq. (34) and the well-known properties of *θ*(*λ*), we have that Q(*λ*) satisfying Eq. (16). From Eq. (34) and the well-known properties of the *δ*-measure, the mass $\mu \xi \zeta 0$ of the spectral measure in (35) satisfies Eq. (18) with appropriate notational changes. We are now ready to present the main result of this section.

*Consider the generalized eigenvalue problem*

*B*

*z*_{n}

*= λ*

_{n}

*C*

*z*_{n}

*in*(

*31*).

*Let*

*Z*

*be the matrix with columns consisting of the eigenvectors*

*z*_{n}

*and*Λ

*be the diagonal matrix with eigenvalues*

*λ*

_{n}

*on the diagonal, which satisfy Eq.*(

*33*).

*The discrete, matrix representations of the bilinear functional formulas for*$Sjk*$

*and*$Ajk*$

*in Eq.*(

*11*)

*are given by*

*Also, the discrete representation of the resolvent formula for*

*χ*_{j}

*in Eq.*(

*12*)

*is given by*

*The discrete representations of the bilinear functional formulas for*$Sjk*$

*and*$Ajk*$

*in*(

*13*)

*are given by*

*Consequently, from the discrete analog of Eq. (*

*23*

*) and the formulas for*$Sjk*$

*and*$Ajk*$

*in*(

*38*),

*we have the following series representations:*

*Finally, defining*$gj=(\u2207T\u2207)\u22121uj$

*and recalling the projection matrix*

*Q*

_{n}

*in*(

*34*),

*we have*

*It follows from Eqs.*(39)

*and*(40)

*that the series representations for*$Sjk*$

*and*$Ajk*$

*in*(

*39*)

*have the Stieltjes integral representations in Eq. (*

*19*

*), involving the discrete spectral measure*

*μ*

_{jk}

*in (*

*35*

*) with*$\xi =gj$

*and*$\zeta =gk\u2009$

*.*

*A*= Δ

^{−1}

**⋅ [H**

*∇***] and Eq. (31), the matrix representation of the functional formulas $Sjk*=\epsilon (\delta jk+\u27e8\chi j,\chi k\u27e91,2)$ and $Ajk*=\u27e8A\chi j,\chi k\u27e91,2$ in Eq. (11) are given by (36). Moreover, the matrix representation of the cell problem**

*∇**ɛ*Δ

*χ*

_{j}+

**⋅ [H**

*∇***]**

*∇**χ*

_{j}= −

*u*

_{j}in (4) is given by

*z*_{n}form a basis for $CK$. Consequently, by Eq. (33), we have that B = Z

^{−†}ΛZ

^{−1}, C = Z

^{−†}Z

^{−1}. It now follows from Eq. (41) that Z

^{−†}(

*ɛ*I +

*ı*Λ)Z

^{−1}

*χ*_{j}=

*u*_{j}, which implies the resolvent formula for

*χ*_{j}in Eq. (37).

*χ*_{j}in (37) into Eq. (36) yields Eq. (38), where we have used that [∇Z]

^{†}= [∇Z]

^{−1}. The quadratic form Z

^{†}

*u*_{j}⋅ Z

^{†}

*u*_{k}arising in (38) can be written in terms of the projection matrices Q

_{n}defined in (34) as follows:

In Sec. V, we use a generalization of Eq. (39) to compute Stieltjes integral representations for the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$ for some model periodic fluid velocity fields. This generalization, discussed in Theorem 10, is a direct analog of Eq. (39) and holds for the setting where the matrix gradient ∇ with periodic boundary conditions is rank-deficient, so the negative matrix Laplacian ∇^{T}∇ is *noninvertible*.

We conclude this section with a discussion that helps reduce the amount of memory required to store the eigenvalues and spectral weights comprising the discrete spectral measure, which is useful when computing a large number of statistical realizations associated with randomly perturbed periodic flows. From the formula for *μ*_{ξζ}(*λ*) above Eq. (35), the fact that the matrix ∇ and vectors $gk$ are real-valued and the two identities $Re\u2009\mu jk(\lambda )=(\mu jk(\lambda )+\mu \xafjk(\lambda ))/2$ and $Im\u2009\mu jk(\lambda )=(\mu jk(\lambda )\u2212\mu \xafjk(\lambda ))/(2\u0131)\u2009$, we have

with $[\u2207(Qn+Q\xafn)gj\u22c5\u2207gk\u2009]=2Re\u2009[\u2207Qn\u2009gj\u22c5\u2207gk]$ and $[\u2207(Qn\u2212Q\xafn)gj\u22c5\u2207gk\u2009]=2\u0131Im\u2009[\u2207Qn\u2009gj\u22c5\u2207gk]$.

Consider the generalized eigenvalue problem in Eq. (31) written as [*ı*B]*z*_{n} = *ıλ*_{n}C*z*_{n}, with B = −*ı*∇^{T}H∇ and C = ∇^{T}∇. Since the matrices *ı*B and C are real-valued, we have $[\u0131B]z\xafn=\u2212\u0131\lambda nC\u2009z\xafn$. Consequently, the (generalized) eigenvectors *z*_{n} come in complex conjugate pairs and the *λ*_{n} come in ± pairs. Therefore, if the size *K* of these matrices is even, then we may renumber the index set $IK$ as $IK={\u2212K/2,\u2026,\u22121,1,\u2026,K/2}$ such that *λ*_{−n} = −*λ*_{n} and $z\u2212n=zn\xaf$. If *K* is odd, then *λ*_{0} = 0 is also an eigenvalue with a *real-valued* eigenvector *z*_{0}. Consequently, the representations of the measures Re*μ*_{jk} and Im*μ*_{jk}, following from the functions in Eq. (44), can be simplified^{48} to depend only on the restricted index set {*n* ≥ 0: *λ*_{n} ≤ *λ*}. This is clear from Eqs. (19) and (44), since for *n* ≥ 0 we have $\lambda \u2212n2=(\u2212\lambda n)2=\lambda n2$ and $z\u2212n=zn\xaf$, thus $Q\u2212n=Q\xafn$. Consequently,

with $\lambda 0Im[Q0gj\u22c5gk]\u22610\u2009$. For numerical computations of statistical properties of advection enhanced diffusion by randomly perturbed periodic flows, a useful consequence of this is only *half* of the eigenvalues and spectral weights need to be held in memory, as the other half are redundant, which saves memory consumption.

## V. SPECTRAL COMPUTATIONS OF THE EFFECTIVE DIFFUSIVITY MATRIX

In Sec. IV, we developed a mathematical framework that provides discrete Stieltjes integral representations for the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$. These discrete integral representations can be written as the series shown in Eq. (39) or the integrals shown in (19), involving the discrete spectral measure in (35) with the appropriate notational changes described in Theorem 3. This framework assumes that the matrix gradient ∇ has Dirichlet boundary conditions, for example, so the matrix is *full-rank* and the negative matrix Laplacian ∇^{T}∇ is invertible. However, for our studies of advection enhanced diffusion by spatially periodic fluid velocity fields, we need to use a matrix gradient ∇ with periodic boundary conditions and this matrix is *rank-deficient*. In order to streamline our presentation leading to the numerical results of this section, in Appendix G, we extend the discrete mathematical framework developed in Sec. IV to the rank-deficient setting. This analysis, in the Proof of Theorem 10, shows that the discrete Stieltjes integral representations for $S*$ and $A*$ are given by direct analogs of the formulas in Eq. (39).

In this section, we use these direct analogs of the formulas in Eq. (39) to compute $D*$ for various model periodic flows and randomly perturbed periodic flows. As a benchmark test, we compute the spectral measure and $D*$ for a shear flow, which are known in closed form,^{4} and a cell flow with closed streamlines, for which it is known^{17,18} that $D*\u223c\epsilon 1/2$ for *ɛ* ≪ 1. Our numerical results are in good agreement with the theoretical results. We also consider a fluid velocity field that has a free parameter. As the parameter varies, the flow transitions from cell flow with closed streamlines to shear flow in the diagonal *x* − *y* direction. This gives rise to transitional behavior in the spectral measure, which governs transitional behavior in the effective diffusivity matrix $D*$. For the sake of brevity, we will focus our attention on the *ɛ*-behavior of the components $Sjk*$, *j*, *k* = 1, …, *d*, of $S*$. Also, for computational simplicity, we have restricted our computations to dimension *d* = 2.

### A. Numerical methods

By Eq. (28), the time-independent fluid velocity field ** u** =

**(**

*u***) is given in terms of an antisymmetric matrix H = H(**

*x***),**

*x***=**

*u***⋅ H. For dimension**

*∇**d*= 2, the matrix H is determined by a stream function Ψ = Ψ(

**),**

*x*yielding ** u** = [−

*∂*

_{y}Ψ,

*∂*

_{x}Ψ], where

*∂*

_{x}denotes partial differentiation in the variable

*x*, for example. In this section, we consider two flows with free parameters which transition from cell flow to shear flow as parameters vary. In particular, we consider BC-flow

^{11}and “cat’s eye” flow,

^{17}which are defined by the following stream functions Ψ

_{BC}and Ψ

_{CE}, respectively:

where we have denoted ** x** = (

*x*,

*y*). The corresponding fluid velocity fields are

The flow geometry of these fluid velocity fields transition from the shear flow to cell flow structure as the system parameters *α*, *B*, *C* ∈ [0, 1] vary.

The streamlines of a 2D flow are level sets of the stream function Ψ, which define a family of curves that are instantaneously tangent to the fluid velocity field ** u**, since

**= [−**

*u**∂*

_{y}Ψ,

*∂*

_{x}Ψ] implies that

**⋅ $\u2207$Ψ = 0. In Fig. 1, we display streamlines of the flows in Eq. (47) for various parameter values. The streamlines for cat’s eye flow are symmetric about the line**

*u**y*=

*x*, which follows from the symmetry of the stream function Ψ

_{CE}(

*x*,

*y*) = Ψ

_{CE}(

*y*,

*x*). The stream function for BC-flow has a more complicated symmetry Ψ(

*x*,

*y*;

*B*,

*C*) = −Ψ(

*y*,

*x*;

*C*,

*B*). This symmetry indicates that if the values of

*B*and

*C*are interchanged,

*B*↔

*C*, then the original flow is recovered from a 90° rotation (

*x*→

*y*,

*y*→ −

*x*) followed by a reflection about the

*x*-axis (

*y*→ −

*y*). Consequently, flows elongated in the

*y*-direction become flows elongated in the

*x*-direction under the interchange

*B*↔

*C*. Consistently, our numerical computations of

*μ*

_{jk}and $Sjk*$,

*j*,

*k*= 1, 2, exhibit these symmetries. For

*BC*-flow, these symmetries allow us to restrict our attention to the behavior of

*μ*

_{jk}and $Sjk*$ as only one parameter varies. Here, we discuss our results for cat’s eye flow for various (deterministic) values of the parameter

*α*as well as

*α*uniformly distributed on the interval [0,

*p*] for various values of

*p*. For the sake of brevity, we consider

*BC*-shear flow in the

*x*and

*y*directions only, which are obtained for parameter values (

*B*,

*C*) = (0, 1) and (

*B*,

*C*) = (1, 0), respectively, and do not display our results for the transitional behavior from one extreme to the other. The spectral measure and $D*$ were computed for

*BC*-cell flow in Ref. 44.

In Sec. III, we gave an overview of the effective parameter problem for the setting of randomly perturbed $V$-periodic flows and introduced the Hilbert space $H$ in Eq. (7). Numerically, it is natural to set $H$ to be the space of randomly perturbed $V$-periodic functions, $H=L2(\nu \xd7P)$, where *P* is the probability measure associated with the random variable *α*.^{29} In this case, by Fubini’s theorem,^{20} ⟨⋅⟩ can be considered to denote spatial averaging followed by statistical averaging.

We now discuss in more detail our discrete, matrix formulation of the effective parameter problem. To illustrate how to generalize these ideas to dimension *d* larger than *d* = 2, we will maintain the aspects of our general notation. In this discrete setting, the spatial region $V=[0,2\pi ]d$, for example, is replaced by a square lattice $VLd$ of size *L* containing *L*^{d} equally spaced points in $V$. As discussed in Sec. IV, the differential operators ** ∇** and

**⋅ are replaced by finite difference, matrix operators ∇ and −∇**

*∇*^{T}, respectively, with suitable boundary conditions. Periodic boundary conditions will be assumed throughout this section. Since these matrices operate on vectors, the

*d*-dimensional lattice $VLd$ must be bijectively mapped to a one dimensional lattice $VN$ of size

*N*=

*L*

^{d}

*d*. The specific structure of $VN$ and ∇ depend on the bijective mapping Θ chosen. In our computations discussed in this section, we used the mapping Θ described in Ref. 42.

The spatially dependent *d*-dimensional vector field ** v**(

**) = (**

*x**v*

_{1}(

**), …,**

*x**v*

_{d}(

**)), say, is bijectively mapped by Θ to a discretized**

*x**constant*vector (

*v*_{1}, …,

*v*_{d}) with

*N*elements, where the vectors

*v*_{i},

*i*= 1, …,

*d*, each have

*L*

^{d}elements and contain all of the spatial information about the

*v*

_{i}(

**) for $x\u2208VLd$. Similarly, the**

*x**d*-dimensional standard basis vector

*e*_{1}= (1, 0, …, 0) is mapped to the

*N*-dimensional vector (

**1**,

**0**, …,

**0**), where

**1**and

**0**are vectors of ones and zeros with

*L*

^{d}elements, and similarly for the

*e*_{j}for

*j*= 2, …,

*d*. Therefore, the vectors $e^j$,

*j*= 1, …,

*d*, satisfying

serve as “lattice standard basis vectors.” With this convention, the division by *L*^{d/2} in (49) takes care of the uniform *L*-scaling in discrete approximations of spatial integrals; instead of the (2*π*)^{d} normalized Lebesgue measure d** x**, we have Δ

**= (2**

*x**π*/

*L*)

^{d}/(2

*π*)

^{d}when $V=[0,2\pi ]d$ and the spatial average $\u27e8\xi (x)\u2009\zeta (x)\u27e9V$ becomes

**⋅**

*ξ***/**

*ζ**L*

^{d}. In a similar way, the 2 × 2 matrix H(

**) in (46) becomes a**

*x**N*×

*N*antisymmetric banded matrix, where the stream function Ψ(

**) is represented by a diagonal**

*x**L*

^{d}×

*L*

^{d}matrix and the zero element 0 is represented by a matrix of zeros. As in previous sections, for notational simplicity, we will not make a notational distinction for the matrix H between the continuum and discrete settings as the context will be clear.

In Theorem 10 of Appendix G, we extend our results developed in Theorem 1 of Sec. III to the setting where the matrix gradient ∇ has periodic boundary conditions and is rank-deficient. This is accomplished by considering the singular value decomposition (SVD) ∇ = UΣV^{T}. Here, ∇ is of size *N* × *K*, say, where *K* = *L*^{d} and *N* = *Kd*. Also, Σ = diag(*σ*_{1}, …, *σ*_{K}), where 0 ≤ *σ*_{1} ≤ ⋯ ≤ *σ*_{K}. The matrices U and V are of size *N* × *K* and *K* × *K*, respectively, and satisfy U^{T}U = I and V^{T}V = VV^{T} = I, where I is the identity matrix of size *K*.^{14} When ∇ is of full rank, then the singular values satisfy *σ*_{j} > 0 for all *j* = 1, …, *K* and the matrix Laplacian ∇^{T}∇ = VΣ^{2}V^{T} is invertible. When ∇ is rank-deficient, then there are *K*_{1} nonzero singular values and *K*_{0} = *K* − *K*_{1} zero singular values, say. For example, when *d* = 2, the nullity of ∇ is 1 so *K*_{1} = *K* − 1. In this case, we write U = [U_{0} U_{1}], Σ = diag(O, Σ_{1}), and V = [V_{0} V_{1}], where O is a matrix of zeros so that $\u2207=U1\Sigma 1V1T$ and the negative matrix Laplacian $\u2207T\u2207=V1\Sigma 12V1T$ is noninvertible, since $V1V1T\u2260I\u2009$.

The associated matrix analysis in Appendix G demonstrates that the spectral measure *μ*_{jk} underlying the Stieltjes integral representation of $Sjk*$ is given by

where $\lambda n1$, *n* = 1, …, *K*_{1}, are the eigenvalues of the matrix $\u2212\u0131U1THU1$, while various equivalent representations of the spectral weights *m*_{jk}(*n*), *j*, *k* = 1, …, *d*, are given in Eq. (G7) of Theorem 10. For notational simplicity, in this section, we denote Re*μ*_{jk} by *μ*_{jk}. In our computations of *μ*_{jk}, we used

which follows from Eqs. (43), (44), (G6), and (G7). Here, $rn1$, *n* = 1, …, *K*_{1}, are the *complex* eigenvectors of the matrix $\u2212\u0131U1THU1$. Consequently, *m*_{kk}(*n*) ≥ 0, so *μ*_{kk} is a positive measure and *μ*_{jk}, *j* ≠ *k*, is a signed measure, where *m*_{jk}(*n*) in (51) can take positive or negative values.

To reveal the structure of *μ*_{12} and $S12*$ in our numerical computations discussed in Sec. V B, we denote the spectral weights *m*_{jk}(*n*) associated with the Jordan decomposition $\mu jk=\mu jk+\u2212\mu jk\u2212$ in (27) by $mjk+(n)$ and $mjk\u2212(n)$, where $mjk\xb1(n)\u22650$. Also, we define the functions $[S12*]+$ and $[S12*]\u2212$,

for each 0 < *ɛ* < *∞*, so that $S12*=[S12*]+\u2212[S12*]\u2212$, $[S12*]\xb1(\epsilon )=S12*(\epsilon ;\mu 12\xb1)$, and $[S12*]\xb1\u22650$.

In the case of a nonrandom fluid velocity field ** u**, we used

*L*= 200 so that

*K*

_{1}= 39 999. The eigenvalues $\lambda n1$ and eigenvectors $rn1$ of the nonrandom Hermitian matrix $\u2212\u0131U1THU1$ were computed using the Matlab function

*eig()*. In this case, the averaging ⟨⋅⟩ in (50) is interpreted as spatial averaging over the period cell $V$. In the setting of a randomly perturbed flow, we used

*L*= 100 so that

*K*

_{1}= 9999. In this case, the averaging ⟨⋅⟩ in (50) is interpreted as spatial averaging followed by ensemble averaging over ∼10

^{3}statistical trials.

The numerical accuracy of the eigenvalue problem is determined by the *eigenvalue condition numbers* $K(\lambda n1)$, *n* = 1, …, *K*_{1}, which are the reciprocals of the cosines of the angles between the left and right eigenvectors. Large eigenvalue condition numbers of a Hermitian matrix implies that it is near a matrix with multiple eigenvalues, while eigenvalue condition numbers close to 1 imply that the eigenvalue problem is well-conditioned. The eigenvalue problem for the matrix $\u2212\u0131U1THU1$ is extremely well conditioned with $maxn|1\u2212K(\lambda n1)|\u223c10\u221214$ typical, computed using Matlab’s function *condeig()*.

To the best of our knowledge, Matlab does not provide a function that describes the accuracy of the computed SVD of the matrix ∇ = UΣV^{T}. In order to better understand the numerical accuracy in the entries of the matrix U, which is central to our computational method, we performed the following tests. For the case of Dirichlet boundary conditions, the matrix ∇ is full-rank; hence, the matrix Laplacian ∇^{T}∇ is invertible. We computed the matrix $\Gamma =\u2207(\u2207T\u2207)\u22121\u2207T$ directly using Matlab’s *mldivide* function, i.e., Γ = ∇((∇^{T}∇)\∇^{T}), and also using the SVD of the matrix ∇, with Γ = UU^{T}. We then computed the componentwise maximum difference $maxlm|[\u2207((\u2207T\u2207)\\u2207T)\u2212UUT]lm|$. When *L* = 100 and *L* = 200, this difference is ∼10^{−15}, which gives an idea of the accuracy of the SVD of ∇ for the rank-deficient, periodic setting. The matrix Γ is used extensively throughout Appendix E. In all of our computations, Matlab’s sparse architecture was employed wherever possible to reduce roundoff errors.

### B. Numerical results

Before we discuss our numerical results in this section, it is helpful to first describe the relationship between the spectral measure *μ*_{jk} and the *ɛ*-behavior of $Sjk*$. This relationship is easiest to understand when *j* = *k*, in terms of the enhancement in scalar transport above the bare diffusive value *ɛ*, as shown in Eq. (20). Consider *μ*_{kk} and $Skk*$, for some *k* = 1, …, *d*, and write the formula in (20) as

where $Ekk(\epsilon )\u22650$ denotes the *enhancement* above *ɛ*. By Eq. (15), the spectrum Σ of the self-adjoint operator *M* is a bounded subset of $R$. Consequently, in the *diffusion dominated regime* where *ɛ* ≫ |*λ*| for all *λ* ∈ Σ, we have^{20} $Skk*\u223c\epsilon +\mu kk0/\epsilon \u2009$. The enhancement $Ekk(\epsilon )\u223c\mu kk0/\epsilon $ is only nominal in this regime where *ɛ* ≫ 1 and is independent of the distribution of measure mass—dependent only the total mass. However, in the *advection dominated regime* where *ɛ* ≪ 1, if the spectral measure *μ*_{kk} has significant mass near the spectral origin *λ* = 0, e.g., if the spectral weights *m*_{kk}(*n*) in (50) have values significantly greater than zero for |*λ*_{n}| ≪ 1, then the integrand associated with $Ekk(\epsilon )$ can introduce singular behavior that competes with the small *ɛ* prefactor in front of the integral, giving rise to a significant enhancement $Ekk(\epsilon )$ for 0 < *ɛ* ≪ 1. Although, if the mass of the measure is zero or extremely small in a neighborhood of *λ* = 0, then the enhancement $Ekk(\epsilon )$ can be less pronounced, since the small *ɛ* prefactor dominates the (lack of) singular behavior in the integrand near *ɛ* = 0. This illustrates that in the advection dominated regime, where 0 < *ɛ* ≪ 1, the *ɛ*-behavior of the enhancement $Ekk(\epsilon )$ depends strongly on the details of the spectral measure *μ*_{kk} near *λ* = 0.

We emphasize that, due to roundoff errors and finite resolution (*L* < *∞*) effects, our numerical approximations of *μ*_{jk} and the *ɛ*-behavior of $Sjk$ breakdown for extremely small values of *ɛ*. For example, in the discrete setting, it is highly unlikely that *λ*_{n} = 0 is (exactly) a numerical eigenvalue solution of Eq. (30), even though in the continuum setting *λ* = 0 *is an accumulation point*^{28,56} of the discrete spectra for the compact^{9,10} operator *M*. Therefore, our numerical simulations can probe the *ɛ*-behavior of $Sjk$ for moderately small values of *ɛ*, but the approximation ultimately breaks down in the limit *ɛ* → 0. However, for moderately small values of *ɛ*, our description above regarding the relationship between *μ*_{jk} and the *ɛ*-behavior of $Sjk$ is still valid—illustrating that the details of the spectral measure *μ*_{jk} near the spectral origin *λ* = 0 strongly influence the *ɛ*-behavior of $Sjk$ when *ɛ* ≪ 1.

These concepts are illustrated in our computations of *μ*_{jk} and $Sjk$ for “cat’s eye” flow displayed in Sec. V B 2. As the free parameter *α* increases from 0 to 1, the flow transitions from cell flow with closed streamlines to shear flow in the *y* = *x* direction, as shown in Fig. 1. Our computations of *μ*_{jk} display a transitional behavior near *λ* = 0 which gives rise to a pronounced change in the *ɛ*-behavior of $Sjk$ near *ɛ* = 0, as well as a significant enhancement in $Skk$ above the bare diffusive value *ɛ*.

As a benchmark result, we demonstrate in Sec. V B 1 that our computations accurately capture the known behavior of *μ*_{kk} and $Dkk*$ for shear flow in the *k*th direction,^{4} where $\mu kk=\mu kk0\u2009\delta 0(d\lambda )$ and $Dkk*=\epsilon +\mu kk0/\epsilon $. As another benchmark result, we demonstrate in Sec. V B 2 that our computations accurately capture the known^{17,18} asymptotic behavior $Dkk*\u223c\epsilon a$ with critical exponent *a* = 1/2, for *ɛ* ≪ 1. In particular, the numerical methods developed in this manuscript compute *a* ≈ 0.54, with an 8% error relative to the true value. For the sake of comparison, our Fourier approach to computing the spectral measure *μ*_{kk} discussed in Ref. 44 computes *a* ≈ 0.52, with a 4% relative error, and our implementation of the linear systems approach to computing $Dkk*$ discussed in Ref. 48 computes *a* ≈ 0.49 with a 2% relative error. This marked increase in the accuracy of the linear systems approach is largely due to the ability to handle larger matrix sizes which, in turn, is due to the *O*(*N*^{2}) numerical complexity of the method compared to the *O*(*N*^{3}) numerical complexity of the spectral measure method.

#### 1. *BC*-shear flow

In the continuum setting, it is known^{4} for shear flow in the *x*-direction that the measure *μ*_{11} is given by a *δ*-measure concentrated at the spectral origin, while *μ*_{22} ≡ 0, and similarly for shear flow in the *y*-direction. As a baseline result, we computed the spectral measures and effective diffusivities for BC-shear-flow in both the *x* and *y* directions, which are obtained for parameter values (*B*, *C*) = (0, 1) and (*B*, *C*) = (1, 0), respectively. Our computations for the components *μ*_{jk}, *j*, *k* = 1, 2, of the spectral measure for BC-shear-flow displayed in Fig. 2 are in good agreement with the theoretical prediction in Ref. 4.

Figure 2 displays the streamlines for *BC*-shear-flow in (a) the *x*-direction and (b) the *y*-direction. In Fig. 2(c), the components $Sjk*$, *j*, *k* = 1, 2, of the effective diffusivity matrix are displayed for *BC*-shear-flow in the *x*-direction. The analogous result for *BC*-shear-flow in the *y*-direction is visually identical to Fig. 2(c) under the mapping $S11*\u2194S22*$, i.e., under the exchange of the colors black ↔ blue. The components *μ*_{jk}, *j*, *k* = 1, 2, of the spectral measure are displayed for *BC*-shear-flow in (d) the *x*-direction and (e) the *y*-direction.

We focus our discussion on the results for *BC*-shear-flow in the *x*-direction, as the discussion regarding *BC*-shear-flow in the *y*-direction is analogous. For all *n* = 1, …, *K*_{1}, the spectral weights *m*_{22}(*n*) in Fig. 2(d) satisfy *m*_{22}(*n*) ≲ 10^{−29}. With the effects of finite resolution (*L* < *∞*) and roundoff error associated with a machine epsilon of ∼10^{−16}, these spectral weights can be considered “numerically zero.” The spectral weights $m12\xb1(n)$ satisfy $m12\xb1(n)\u227210\u221228$ for *λ*_{n} away from the spectral origin *λ* = 0 with a peak near *λ* = 0 having magnitudes $m12\xb1(n)\u227210\u221216$. The spectral weights for the *x*-direction satisfy *m*_{11}(*n*) ≲ 10^{−28} away from *λ* = 0, while the weights near *λ* = 0 satisfy 10^{−9} ≲ *m*_{11} ≲ 10^{−1}, resembling a *δ*-measure with virtually all of its mass concentrated near *λ* = 0. This is consistent with theoretical predictions.^{4} Due to the antisymmetry of the *real-valued* matrix $U1THU1$, its complex eigenvectors and purely imaginary eigenvalues come in complex conjugate pairs.^{25} Consequently, the eigenvalues of the Hermitian matrix $\u2212\u0131U1THU1$ come in ± pairs with identical spectral weights, resulting in the symmetry about *λ* = 0 displayed by the spectral measures in Fig. 2.

Due to the high concentration of measure mass in *μ*_{11} very near the spectral origin, our computation of $S11*$ displayed in Fig. 2(c) behaves like it is being governed by a delta function concentrated at the origin. In particular, Fig. 2(c) shows that the computed *ɛ*-behavior of $S11*$, displayed in black color with solid line-style, lays right on top of the graph of its upper bound $\epsilon \u2009[1+\mu 110/\epsilon 2]$ given in (24), with $\mu 110\u22484.975\xd710\u22121$, displayed in black color and dashed-dotted line-style. (We had to increase the linewidth of the upper bound to be able to distinguish between the two black lines.) Due to the extremely small magnitudes of the spectral weights *m*_{22} and $m12\xb1$, with measure masses $\mu 220\u22485.33\xd710\u221229$, $[\mu 120]+\u22481.03\xd710\u221216$, and $[\mu 120]\u2212\u22483.34\xd710\u221215$, both the upper and lower bounds for $S22*$ and $S12*$ in Eqs. (24) and (25) are very close to *ɛ* and 0, respectively; the graph of $S22*$ is virtually right on top of the lower bound *ɛ* in cyan color and solid line-style, and the magnitudes of $[S12*]+$ and $[S12*]\u2212$ are so small they do not even appear. Since the support of the spectral measure is contained in the interval [−1, 1], the components of the effective diffusivity approach their bare diffusive value *ɛ δ*_{jk} for large *ɛ*, as discussed above.

In Ref. 44, we developed Fourier methods for the computation of the spectral measure *μ*_{jk} for *BC*-cell flow, with *B* = *C* = 1. In particular, the eigenvalue problem *Mφ*_{n} = *λ*_{n}*φ*_{n} associated with the operator *M* = −*ı*Δ^{−1}[** u** ⋅

**] was transformed into an infinite algebraic system of equations defining a discrete, generalized eigenvalue problem. The Fourier coefficients of the eigenfunctions**

*∇**φ*

_{n},

*n*= 1, 2, 3, …, for the continuum setting comprise the components of the generalized eigenvectors in the discrete setting. Since we already treated

*BC*-cell flow in Ref. 44, and for the sake of brevity, we now turn our attention to a discussion regarding our numerical results for “cat’s eye flow” displayed in Fig. 3.

#### 2. Cat’s eye flow

Since the streamlines for cat’s eye flow in Fig. 1 are symmetric about the line *y* = *x* for all *α* ∈ [0, 1], as discussed above, we anticipate that *μ*_{11} = *μ*_{22}. Our computations of the components *μ*_{jk}, *j*, *k* = 1, 2, of the spectral measure shown in Figs. 3 and 5 display this symmetry. A closer look at these figures reveals a deeper symmetry, namely, that *μ*_{11} = *μ*_{22} = |*μ*_{12}|, where $|\mu 12|=\mu 12++\mu 12\u2212$ is the total variation of the signed measure *μ*_{12} introduced in Eq. (27), i.e., superimposing the panels for $m12+$ and $m12\u2212$ in Figs. 3 and 5 yields the figure panels for *m*_{11} and *m*_{22}. We have also numerically verified the behavior *μ*_{11} = *μ*_{22} = |*μ*_{12}|.

Since the operator *A* = Δ^{−1}[** ∇** ⋅

**] is compact,**

*u*^{9,10}its spectrum is discrete with an accumulation point at the spectral origin

*λ*= 0.

^{56}This accumulation point behavior of the measures

*μ*

_{jk},

*j*,

*k*= 1, 2, can be seen in all of the panels of Fig. 3. When the parameter

*α*= 0, the streamlines of cat’s eye flow are closed cell structures, as shown in Fig. 1, so that large scale transport occurs only when

*ɛ*> 0.

^{17}In this case, the computed “accumulation point” has eigenvalues

*λ*

_{n}with very small magnitude 10

^{−19}≲ |

*λ*

_{n}| ≲ 10

^{−14}. (It is clear that a finite number of eigenvalues do not constitute an accumulation point, but we will use this terminology to identify the concentration of eigenvalues near

*λ*= 0 shown in Figs. 3–5 and to set ideas.) However, the spectral measure weights

*m*

_{kk}(

*n*) and $mjk\xb1(n)$ of the accumulation point have even smaller magnitudes with 10

^{−35}≲

*m*

_{kk}(

*n*) ≲ 10

^{−29}and similarly for $mjk\xb1(n)$, as shown in Fig. 3. Consequently, this component of the spectral measure does not contribute significantly to the enhancement $Ejk(\epsilon )=\epsilon \u2211n\u2009mjk(n)/(\epsilon 2+\lambda n2)$ of $Sjk*(\epsilon )=\epsilon \delta jk+Ejk(\epsilon )$. Plotting the panels of Fig. 3 with the log scale

*x*-axis reveals that there is a

*gap*in the spectral measure with no spectra in the interval 10

^{−14}≲ |

*λ*

_{n}| ≲ 10

^{−7}, as shown in Fig. 4. The other component of the spectral measure has spectra in the interval 10

^{−7}≲ |

*λ*

_{n}| ≲ 10

^{0}and weights 10

^{−37}≲

*m*

_{kk}(

*n*) ≲ 10

^{−1}. However, the part of this component of the spectral measure with spectral weights having more significant magnitudes 10

^{−4}≲

*m*

_{kk}(

*n*) ≲ 10

^{−1}are associated with eigenvalues with magnitude |

*λ*

_{n}| ≳ 10

^{−1}and consequently also do not contribute significantly to the enhancement $Ejk$.

When 0 < *α* ≪ 1, open channels connect neighboring cells and large scale transport takes place both in thin boundary layers and within the channels.^{17} This is reflected in the spectral measure $d\mu jk(\lambda )=\u2211n\u27e8mjk(n)\u2009\delta \lambda n1(d\lambda )\u27e9$ in Eq. (50) by a dramatic increase in the magnitude of the spectral weights *m*_{kk}(*n*) and $m12\xb1(n)$ associated with the accumulation point—by more than *14 orders of magnitude* with 10^{−19} ≲ *m*_{kk}(*n*) ≲ 10^{−15}—corresponding to a change of only 10^{−6} in the magnitude of *α*, as shown in Figs. 3 and 4. The associated change in the spectral weights away from *λ* = 0 is nominal. However, Fig. 4 reveals that as *α* increases from 0 to 10^{−6} a localized portion of the accumulation point begins to migrate away from *λ* = 0 to the other component of the spectral measure with spectra in the interval 10^{−7} ≲ |*λ*_{n}| ≲ 10^{0}—decreasing the mass of the accumulation point for *μ*_{kk}. Moreover, all of the positive masses $m12+(n)$ migrate away from the accumulation point so the accumulation point of *μ*_{12} becomes comprised with purely negative weights. From this discussion, it is clear that the increased magnitudes of the masses comprising the accumulation point provide an increased contribution to the enhancement $Ejk(\epsilon )$ for *ɛ* ∼ 10^{−14}, for example, but only a moderate increase, and the increase in $Ejk(\epsilon )$ for 10^{−7} ≲ *ɛ* ≲ 10^{0} is also not significant.

As the value of *α* increases to *α* = 1, the magnitudes of the masses comprising the accumulation point increase until they reach maximum values with 10^{−11} ≲ *m*_{kk}(*n*) ≲ 10^{−5}, associated with eigenvalues with magnitudes 10^{−19} ≲ |*λ*_{n}| ≲ 10^{−14}, as shown in Fig. 4. This marked increase in the magnitudes of the spectral weights provide a significant contribution to the enhancement $Ejk(\epsilon )$ for *ɛ* ∼ 10^{−14}, for example. Moreover, as *α* increases in the range (10^{−1}, 10^{0}), a significant transitional behavior arises in the other component of the spectral measure away from *λ* = 0, as shown in Figs. 4 and 5. In particular, a bulge of measure weights with significant magnitude forms for eigenvalues in the range 10^{−4} ≲ |*λ*_{n}| ≲ 10^{−1}, with a marked increase in weight magnitude from 10^{−7} ≲ *m*_{kk}(*n*) ≲ 10^{−4} to 10^{−7} ≲ *m*_{kk}(*n*) ≲ 10^{−2}. This provides a significant contribution to the enhancement $Ejk(\epsilon )$ even for 10^{−4} ≲ *ɛ* ≲ 10^{−1}.

The behavior of $Sjk*$ that we deduced from the behavior of the spectral measure *μ*_{jk} is consistent with our computations of the components $Sjk*$, *j*, *k* = 1, 2, of the matrix $S*$, which are displayed in Fig. 6. Since the support of the spectral measure *μ*_{jk} is contained in the interval [−1, 1], the components $Sjk*$ of the effective diffusivity approach their bare diffusive value *ɛ δ*_{jk} as *ɛ* surpasses *ɛ* = 1, as discussed in the beginning of Sec. V B.

For cat’s eye cell-flow, when *α* = 0 the log-log plot of $Skk*$ displays a linear trend for 10^{−3} ≲ *ɛ* ≲ 10^{−1}, capturing the known^{17,18} power law behavior $Skk*\u223c\epsilon a$ as *ɛ* → 0. The polynomial fit *P*(*ɛ*) to this line is given by *P*(*ɛ*) = 0.54*ɛ* + 0.08. This calculation of the critical exponent *a* = 0.54 has an 8% error relative to the value of the theoretical result *a* = 1/2.^{17,18} The presence of spectral weights with magnitudes 10^{−7} ≲ *m*_{kk}(*n*) ≲ 10^{−4} and associated eigenvalues 10^{−4} ≲|*λ*_{n}|≲ 10^{−1} gives rise to a moderate enhancement in $Skk*$. This enhancement increases from a fraction of an order of magnitude to one and a half orders of magnitude above its bare diffusive value *ɛ*, as *ɛ* decreases from 10^{−1}–10^{−3}, as shown in Fig. 3 for *α* = 0.

We deduced from Fig. 3 that for *α* ∈ (0, 10^{−2}) the behavior of the spectral measure gives rise to only a moderate enhancement $Ejk(\epsilon )=\epsilon \u2211n\u2009mjk(n)/(\epsilon 2+\lambda n2)$ of the effective diffusivity $Sjk*(\epsilon )=\epsilon \delta jk+Ejk(\epsilon )$ both for the advection dominated regime where *ɛ* ∼ 10^{−14}, for example, and for the transitional regime 10^{−3} ≲ *ɛ* ≲ 10^{−1}. This is consistent with the behavior of $Sjk*$ shown in Fig. 6 for *α* = 0 and *α* = 0.05. We further deduced from Figs. 4 and 5 that for *α* ∈ (10^{−2}, 1) the behavior of the spectral measure gives rise to a significant enhancement for both *ɛ* ∼ 10^{−14} and 10^{−3} ≲ *ɛ* ≲ 10^{−1} where there is a marked increase in spectral weight magnitude from 10^{−7} ≲ *m*_{kk}(*n*) ≲ 10^{−4} to 10^{−7} ≲ *m*_{kk}(*n*) ≲ 10^{−2} for eigenvalues satisfying 10^{−4} ≲ |*λ*_{n}| ≲ 10^{−1}. This is consistent with the panels of Fig. 6 corresponding to 0.1 ≤ *α* ≤ 1, with $Sjk*$ enhanced many orders of magnitude above its bare diffusive value *δ*_{jk} *ɛ*, and $Skk*$ as well as $S12*$ closely following their upper bounds for *ɛ* ≳ 10^{−0.5} when *α* = 1.

We conclude this section with a description of various symmetries arising in our numerical computations and their consequences. We discussed above that our computations of *μ*_{jk}, *j* = 1, 2, display the symmetry *μ*_{11} = *μ*_{22} = |*μ*_{12}|. This gives rise to the symmetry $S11*(\epsilon )=S22*(\epsilon )=\epsilon +E12(\epsilon ;\mu 12+)+E12(\epsilon ;\mu 12\u2212)$ between the components of the effective diffusivity, where we have denoted $Ejk(\epsilon ;\mu jk)=\epsilon \u222bd\mu jk(\lambda )/(\epsilon 2+\lambda 2)$, e.g., $S11*(\epsilon )=\epsilon +E(\epsilon ;\mu 11)$. The symmetry $S11*=S22*$ can be clearly seen in our computations of $Sjk*$, *j*, *k* = 1, 2, displayed in Fig. 6; the two curves lay right on top of one another, as do their upper and lower bounds as $\mu 110=\mu 220$. We have also numerically explored the empirical relationship $S11*\u2248\epsilon +[S12*]++[S12*]\u2212$ by plotting $S11*$ and $\epsilon +[S12*]++[S12*]\u2212$ on one graph. For most values of *α* and *ɛ* considered, the three curves lay virtually on top of each other (not shown), and when there is a deviation of $\epsilon +[S12*]++[S12*]\u2212$ from $S11*$, it is slight. This property also leads to the inequalities $S11*\u2265\epsilon +[S12*]+$ and $S11*\u2265\epsilon +[S12*]\u2212$, with $S11*=S22*$, which is consistent with the behavior of $Sjk*$ shown in Fig. 6.

#### 3. Randomly perturbed cat’s eye flow

We now discuss our computations of the components *μ*_{jk} and $Sjk*$, *j*, *k* = 1, 2, of the spectral measure and effective diffusivity matrix, respectively, for randomly perturbed cat’s eye flow with *α* uniformly distributed on the interval [0, *p*]. For each statistical trial of a sample space Ω_{0}, with |Ω_{0}| ∼ 10^{3} and *L* = 100, we computed *every* eigenvalue $\lambda n1$ and eigenvector $rn1$, *n* = 1, …, *K*_{1}, of the matrix $\u2212\u0131U1THU1$ to form the spectral measure *μ*_{jk} in Eq. (50). In order to visually determine the behavior of the function $\mu jk(\lambda )=\u27e8Q(\lambda )e^j,e^k\u27e9$ underlying the spectral measure *μ*_{jk}, we plot a histogram representation of *μ*_{jk}(*λ*) called the *spectral function*, which we will also denote by *μ*_{jk}(*λ*). We now describe how we computed this graphical representation of the measure *μ*_{jk}. First, the spectral interval *I* ⊇ Σ was divided into *V* subintervals *I*_{v}, *v* = 1, …, *V*, of equal length. In our computations of the spectral functions, we typically used *V* ∼ 10^{2}. Second, for fixed *v*, we identified all of the eigenvalues that satisfy $\lambda n1(\omega )\u2208Iv$, for *n* = 1, …, *K*_{1} and *ω* ∈ Ω_{0}. The assigned value of *μ*_{jk}(*λ*) at the midpoint *λ* of the interval *I*_{v} is the sum of the spectral weights *m*_{jk}(*ω*) associated with all such $\lambda n1(\omega )\u2208Iv$, normalized by |Ω_{0}|. In this way, the area underneath the curve is the measure mass $\mu jk0$.

Consistent with the symmetries of the randomly perturbed flow, our computations of the spectral function satisfy *μ*_{11}(*λ*) = *μ*_{22}(*λ*); hence, the ensemble averaged components $Sjk*$ of the effective diffusivity also satisfy $S11*=S22*$, as shown in Fig. 7. Similar to our computations for nonrandom *α*, when *p* = 0.1, the measure mass of *μ*_{jk}, *j*, *k* = 1, 2, near the spectral origin *λ* = 0 is quite small and, on average, the region with a significant magnitude increases in breadth as *p* increases, with the formation of a high concentration of measure mass at the spectral origin *λ* = 0 as *p* → 1—due to the incorporation of statistical realizations with near shear-flow characteristics associated with *α* ≈ 1. This average increase in the breadth of the region with significant mass and the formation of the high concentration of measure mass at *λ* = 0 give rise to a substantial enhancement of the components $Sjk*$ of the effective diffusivity above the bare molecular diffusivity values *ɛ δ*_{jk}. The sign changes in *μ*_{12}(*λ*) give rise to sign changes in $S12*=[S12*]+\u2212[S12*]\u2212$. In the log-log plot of $S12*$, a negative singularity forms in $S12*$ at the location of sign changes.

## VI. CONCLUSIONS

We adapted and extended two methods previously introduced in Refs. 3, 4, 9, 10, and 48 to provide new Stieltjes integral representations for the symmetric $S*$ and antisymmetric $A*$ parts of effective diffusivity matrix $D*$ in (19)—for all components of these homogenized matrices. Each integral representation involves the nondimensionalized molecular diffusivity *ɛ* and a spectral measure of a self-adjoint operator acting on an appropriate Hilbert space. We utilized these integral representations to derive rigorous bounds for the off-diagonal components of the matrices $S*$ and $A*$. We also proved that the spectral measures of both methods are identical, establishing that the two approaches yield equivalent spectral representations for $D*$.

We developed discrete formulations of these two mathematical frameworks involving matrix representations of the self-adjoint operators and developed a standard and a *nonstandard* spectral theorem in terms of a standard and *generalized* eigenvalue problem, respectively. This matrix analysis provided the Stieltjes integral representations for $S*$ and $A*$ in (19), involving discrete spectral measures given explicitly in terms of the eigenvalues and eigenvectors of the matrices. We developed these discrete frameworks for both of the settings where the matrix gradient ∇ has Dirichlet boundary conditions, for example, and is therefore *full-rank*, and where ∇ has periodic boundary conditions and is therefore *rank-deficient*. In our studies of advection enhanced diffusion by a *periodic* fluid velocity field ** u** here, it is necessary to use a matrix gradient ∇ with periodic boundary conditions. We also proved, in both the full-rank and rank-deficient settings, that the spectral measures of both methods are identical, establishing that the two approaches yield equivalent spectral representations for the effective diffusivity matrix $D*$. More specifically, this matrix analysis demonstrates that both approaches can be formulated in terms of a common

*standard*eigenvalue problem involving a matrix with the

*smaller size*encountered in the generalized eigenvalue problem, thus combining the beneficial numerical attributes of both approaches.

We employed these discrete formulations to compute the components $Sjk*$, *j*, *k* = 1, …, *d*, of the matrix $S*$ for some model 2D (*d* = 2) periodic flows and randomly perturbed periodic flows, by directly computing the associated discrete spectral measure Re *μ*_{jk}. As a baseline result, we computed $Sjk*$ and Re *μ*_{jk} for *BC*-shear-flow, for which the spectral measure is known.^{4} Our numerical results are in good agreement with the theoretical result. We also computed $Sjk*$ and Re *μ*_{jk} for the “cat’s eye” flow, for both the nonrandom and randomly perturbed settings, as a function of a free parameter. As the parameter varies, the flow transitions from cell-flow to shear-flow in the diagonal direction *y* = *x*. For cat’s eye cell-flow, our computations capture the known^{17,18} power law behavior $Skk*\u223c\epsilon 1/2$ for *ɛ* ≪ 1. The spectral measure Re *μ*_{jk} and $Sjk*$ have transitional behavior as the parameter varies. This reveals how the details of Re *μ*_{kk} near the spectral origin govern the enhancement of $Skk*$ above its bare diffusive value *ɛ* in the advection dominated regime where *ɛ* ≪ 1 and similarly for $S12*$. Consistent with the symmetries of the flow, our computations indicate that Re *μ*_{11} = Re *μ*_{22}. Our computations of Re *μ*_{12} for cat’s eye flow also suggest a deeper symmetry, namely, |Re *μ*_{12}| = Re *μ*_{11} = Re *μ*_{22}, where |Re *μ*_{12}| is the total variation of the signed measure Re *μ*_{12}. Our computations of $Sjk*$ are consistent with these symmetries and rigorous bounds derived in Theorem 2. In order to streamline the presentation of the main theoretical and numerical results in the body of the manuscript, we have placed more detailed developmental material in Appendixes D-- G.

In almost 30 years since the initial formulation^{3,4} of Stieltjes integral representations for the effective diffusivity matrix $D*$, analytical calculations of $D*$ have been obtained for only a few simple flows, such as shear flow. Our results help further advance the applicability of this spectral measure approach by providing a mathematical foundation for computation of spectral representations of $D*$. For randomly perturbed periodic flows, the spectral method differs from more traditional methods in that it enables statistical investigation of the random eigenvalues and eigenvectors, thus connecting homogenization of advection diffusion processes to random matrix theory.^{45} The results in this manuscript lay the groundwork for such investigations.

## ACKNOWLEDGMENTS

We gratefully acknowledge support from the Applied and Computational Analysis Program and the Arctic and Global Prediction Program at the U.S. Office of Naval Research through Grant Nos. N00014-12-10861, N00014-13-10291, and N00014-18-1-2552. We are also grateful for support from the Division of Mathematical Sciences at the U.S. National Science Foundation (NSF) through Grant Nos. DMS-0940249, DMS-1413454, DMS-1715680, DMS-1211179, and DMS-1522383. Finally, we would like to thank the NSF Math Climate Research Network (MCRN) for their support of this work.

### APPENDIX A: APPENDIX OVERVIEW

In order to streamline the presentation of the main theoretical and numerical results in the body of the manuscript, we have placed more detailed developmental material in a series of appendixes here. We now give an overview of the topics covered in this appendix. In Appendix B, we comment on the notation used throughout this manuscript. In Appendix C, we derive important properties of the linear operator *A* = Δ^{−1}[** u** ⋅

**] defined in Eq. (11).**

*∇*In Theorem 1 of Sec. III, we adapted and extended a method^{9,10,48} involving a self-adjoint operator *M* acting on a Sobolev space of scalar-valued functions which provides the Stieltjes integral representations for $S*$ and $A*$ in Eq. (19). In Appendix D, we adapt and extend a different method^{3,4} involving a self-adjoint operator **M** acting on the Hilbert space of curl-free vector-valued functions which also leads to the Stieltjes integral representations for $S*$ and $A*$ in Eq. (19). These results are summarized in Corollary 4 of Theorem 1.

In Theorem 5 of Appendix D, we prove that the spectral measures arising in Theorem 1 and Corollary 4 are identical, establishing that the two approaches yield equivalent spectral representations of $D*$. This is accomplished by proving that the masses and all the moments of the two spectral measures are equal and citing the Hausdorff moment problem for measures with bounded support.^{1,58} In Corollary 6 of Theorem 5, we utilize a one-to-one isometric correspondence^{44} between the Sobolev space arising in Theorem 1 and the Hilbert space arising in Corollary 4, to extend the results of Theorem 5 to every spectral measure associated with the two self-adjoint operators *M* and **M**. This corollary also proves that the masses and all the moments of the two spectral measures are equal in the generalized setting of a space-time periodic flow, possibly associated with chaotic dynamics and unbounded spectrum.^{44} In the setting of unbounded spectrum, the Hamburger or Stieltjes moment problems are instead relevant, and more conditions must be met beyond the equality of the masses and moments to ensure that the measures are identical, such as Carleman’s criterion.^{1} In Ref. 44, an alternate method was used to determine that the spectral measures arising in the two different approaches are indeed identical in the time-dependent setting.

In Appendix E, we develop a discrete formulation of the mathematical framework given in Appendix D, involving a Hermitian matrix representation for the self-adjoint operator **M**. We also briefly review the standard spectral theorem for Hermitian matrices. This provides the Stieltjes integral representations for $S*$ and $A*$ in (19) involving a discrete spectral measure, given explicitly in terms of the eigenvalues and eigenvectors of the matrix. This discrete framework holds for the setting where the matrix gradient has Dirichlet boundary conditions, for example, and is therefore *full-rank*. These results are analogs of those in Theorem 3 and are summarized in Corollary 7. In Theorem 8, we develop a *projection method* that is used to generalize the results in Theorem 3 and Corollary 7 to the setting where the matrix gradient has periodic boundary conditions, for example, and is therefore *rank-deficient*. The results of Theorem 8 are also used to prove that the discrete spectral representations arising in Theorem 3 and Corollary 7 are equivalent in this rank-deficient setting.

Specifically, in Lemma 9 of Appendix F, we use the properties of the singular value decomposition of the matrix gradient ∇ to reveal symmetries between the two discrete approaches formulated in Sec. IV and Appendix E, establishing that the two approaches yield equivalent spectral representations of the effective diffusivity matrix $D*$ when ∇ is full-rank. In particular, we establish in the Proof of Lemma 9 that the eigenvalues and generalized eigenvalues underlying the spectral measures for each method are in fact eigenvalues of a Hermitian matrix arising in both methods. Moreover, the eigenvectors *w*_{n} and generalized eigenvectors *z*_{n} of the two methods are related by *w*_{n} = ∇*z*_{n}, which leads to the equivalence of the discrete spectral measures of the two approaches.

In Theorem 10 of Appendix G, we generalize Lemma 9 to the rank-deficient setting. This generalizes both the numerical algorithms developed in Sec. IV and Appendix E to the setting of periodic boundary conditions and combines beneficial numerical attributes of each algorithm. In Sec. V, this common method is used to compute $S*$ for model flows and relate spectral characteristics to flow geometry and transport properties.

### APPENDIX B: NOTATION

We now briefly discuss the notation used throughout the manuscript. Operators are denoted by capital letters, while functions comprising the domains of these operators are denoted by lowercase letters. (Capital letters are also used to denote the size of matrices.) Furthermore, the vector-valued functions are denoted by lowercase bold font, e.g., ** ξ**, while the scalar-valued functions are denoted by lowercase nonbold font, e.g.,

*ξ*. Matrices denoted by capital Latin letters are in math-serif font, e.g., H, B, C, Z, Q, A, U, V, W, etc. Integro-differential operators mapping scalar-valued functions to scalar-valued functions are in standard math-italic font, e.g.,

*A*and

*M*. Integro-differential operators mapping vector-valued functions to vector-valued functions are in math-boldface font, e.g.,

**A**and

**M**. We use a similar notation for differential operators, e.g.,

*∇**ξ*and −Δ

*ξ*and their discrete, matrix counterparts ∇

**and ∇**

*ξ*^{T}∇

**, respectively. All homogenized matrices are in Gothic font and have a superscript asterisk, e.g., $D*$, $S*$, and $A*$.**

*ξ*### APPENDIX C: PROPERTIES OF THE LINEAR OPERATOR *A*

In this section, we derive various properties of the linear operator *A* = Δ^{−1}[** u** ⋅

**] defined in Eq. (11). In particular, we demonstrate that**

*∇**A*is antisymmetric on the Hilbert space $H1,2$ defined in (8). Moreover, we show that

*A*is bounded on $H1,2$ and we provide an upper bound for ∥

*A*∥

_{1,2}when

**is uniformly bounded on the period cell $V$.**

*u*We first show that the incompressibility condition ** ∇** ⋅

**= 0 implies that the operator**

*u**A*is antisymmetric on $H1,2$,

^{9,10}i.e., ⟨

*Af*,

*h*⟩

_{1,2}= −⟨

*f*,

*Ah*⟩

_{1,2}. On the Hilbert space $H$ defined in Eq. (7), the linear operator Δ

^{−1}satisfies ⟨ΔΔ

^{−1}

*f*,

*h*⟩ = ⟨

*f*,

*h*⟩ in a distributional sense, for all $f,h\u2208H$.

^{19,37,44}Consequently, integration by parts and

**⋅**

*∇***= 0 yields**

*u*^{9,10,44,48}

for all $f,h\u2009\u2208H1,2$ and real-valued incompressible ** u** (see Ref. 44 for more details).

Now, we derive the bound for ∥*A*∥_{1,2} given in Eq. (14). From the Cauchy-Schwartz inequality |⟨*f*, *h*⟩| ≤ ∥*f*∥ ∥*h*∥, we have

We now provide an upper bound for ∥** u** ⋅

*∇**f*∥ when the components

*u*

_{k},

*k*= 1, …,

*d*, of the fluid velocity field

**are uniformly bounded on the period cell $V$. By the Cauchy-Schwartz inequality, |**

*u***⋅**

*ξ***| ≤ |**

*ζ***| |**

*ξ***|, we have**

*ζ*The result in Eq. (14) is now clear.

### APPENDIX D: CURL-FREE FIELDS AND THE EFFECTIVE DIFFUSIVITY MATRIX

In this section, we adapt and extend an alternative method^{3,4} to the method discussed in Sec. III which leads to the integral representations for the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$ in Eq. (19). More specifically, in Appendix D 1, we provide functional formulas for $S*$ and $A*$ that are analogous to the formulas in (11). We review a Hilbert space formulation of this effective parameter problem^{3,4,17,18,33} in Appendix D 2 which leads to a resolvent formula for *∇**χ*_{j} that is analogous to the resolvent formula for *χ*_{j} in (12), involving a self-adjoint operator. We use this result and the spectral theorem^{53,58} to provide the Stieltjes integral representations for $S*$ and $A*$ in (19), involving a spectral measure of the operator. Finally, we prove that the spectral measure corresponding to the Stieltjes integral representation for $D*$ in Eq. (19) of Theorem 1 is identical to the spectral measure corresponding the Stieltjes integral representation for $D*$ developed in this section. This establishes that the two different approaches yield equivalent spectral representations of $D*$.

#### 1. Functional formulas for the effective diffusivity matrix

In this section, we derive functional formulas for the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$ that are analogous to the formulas in Eq. (11). Using Eq. (29) and the representation of the fluid velocity field ** u** in (28), the advection diffusion equation in (1) can be written as a diffusion equation,

^{17,18}

Here, *e*_{k} is a standard basis vector, *k* = 1, …, *d*, and D(** x**) =

*ɛ*I + H(

**) can be viewed as a local diffusivity matrix with coefficients**

*x*where *δ*_{jk} is the Kronecker delta and I is the identity operator on $Rd$.

Substituting into Eq. (3) the expression for *u*_{j} in (4) and using Eq. (28), ** u** =

**⋅ H, show that the components $Sjk*$ and $Ajk*$ of $S*$ and $A*$ can be written in terms of the following functional formulas involving the**

*∇**real-valued*vector field

*∇**χ*

_{k}:

The functional formulas in (D4) are analogous to the functional formulas in Eq. (11). The symmetry $Sjk*=Skj*\u2009$ of the matrix $S*$ follows from (D4) and the fact that the vector field *∇**χ*_{k} is real-valued so that ⟨*∇**χ*_{j} ⋅ *∇**χ*_{k}⟩ = ⟨*∇**χ*_{k} ⋅ *∇**χ*_{j}⟩. Moreover, the positivity condition ⟨*∇**χ*_{k} ⋅ *∇**χ*_{k}⟩ = ⟨|*∇**χ*_{k}|^{2}⟩≥ 0 demonstrates that the effective transport of the scalar density *ϕ* is always *enhanced* by the presence of an incompressible velocity field, i.e., $Dkk*=Skk*\u2265\epsilon $. The equality $Dkk*=Skk*$ follows from the skew-symmetry of the matrix $A*$, $Akj*=\u2212Ajk*$, hence $Akk*=0$. The skew-symmetry of $A*$ follows from the skew-symmetry of the *real-valued* matrix H, $Ajk*=\u27e8H\u2207\chi j\u22c5\u2207\chi k\u27e9=\u2212\u27e8\u2207\chi j\u22c5H\u2207\chi k\u27e9=\u2212\u27e8H\u2207\chi k\u22c5\u2207\chi j\u27e9=\u2212Akj*\u2009$.

#### 2. The analytic continuation method and integral representations of $D*$

In this section, we begin by noting that the cell problem in Eq. (D2) is equivalent^{3,4,17,18} to the quasi-static limit of Maxwell’s equations,^{22,27,40} which describe the transport properties of an electromagnetic wave in a composite material,

Here, *E*_{k} = *∇**χ*_{k} + *e*_{k} plays the role of the local electric field, *J*_{k} = D*E*_{k} plays the role of the local current density, and D = *ɛ*I + H plays the role of the local conductivity matrix of the medium. Since H is skew-symmetric, the intensity-flux relation *J*_{k} = D*E*_{k} is not the usual Fourier law, but instead resembles that of a Hall medium.^{17,18,26,40}

In Refs. 3 and 4, the analytic continuation method for representing transport in composites was adapted to provide a Stieltjes integral representation for the symmetric part of the effective diffusivity matrix $D*$, involving a spectral measure of a self-adjoint operator. This method provides Stieltjes integral representations for the bulk transport coefficients of composite media,^{22} such as the effective electrical conductivity. This method is based on the spectral theorem of Hilbert space theory and a resolvent formula for, say, the electric field, involving a self-adjoint operator^{22} or matrix^{42} which depends only on the composite geometry. In this section, we adapt the method developed in Refs. 3 and 4 to provide Stieltjes integral representations for *both* the symmetric and antisymmetric parts of the effective diffusivity matrix $D*$, which encodes the flow geometry of the fluid velocity field in a spectral measure of a self-adjoint operator.

Following the discussion leading to Eq. (8), we consider a fluid velocity field that is spatially periodic on a region $V\u2282Rd$ and define the Hilbert space $H$,^{18,33}

which is analogous to the Hilbert space $H$ defined in Eq. (7), where $H$ is defined over vector-fields instead of scalar-fields as in $H$. The Hilbert space $H$ is equipped with a sesquilinear inner-product ⟨⋅, ⋅⟩ defined by ⟨** ξ**,

**⟩ = ⟨**

*ζ***⋅**

*ξ***⟩, with**

*ζ***⋅**

*ξ***=**

*ζ*

*ξ*^{†}

**, † is the operation of complex-conjugate-transpose, and $\u27e8\zeta ,\xi \u27e9=\u27e8\xi ,\zeta \u27e9\xaf$ for $\u2009\xi ,\zeta \u2208H$. This inner-product induces a norm ∥⋅∥ defined by ∥**

*ζ***∥ = ⟨**

*ξ***,**

*ξ***⟩**

*ξ*^{1/2}and $\xi \u2208H$ implies that ∥

**∥ <**

*ξ**∞*. Now consider the associated Hilbert space $H\xd7$ of curl-free fields,

^{4,17,18,22,40,44}

The curl-free vector field *∇**χ*_{k} in the cell problem in (D2) is mean-zero ⟨*∇**χ*_{k}⟩ = 0. When the matrix D in Eq. (D1) is bounded in the operator norm ∥⋅∥ induced by the $H$-inner-product,^{20} ∥D∥ < *∞*, then there exists unique $\u2207\chi k\u2208H\xd7$ satisfying Eq. (D2).^{22,46} We assume that

which together imply ∥D∥ < *∞*.

The linear operator ** Γ** =

**(Δ**

*∇*^{−1})

**⋅ is a projection onto the Hilbert space $H\xd7$ in the sense that $\Gamma :H\u21a6H\xd7$ and**

*∇***=**

*Γξ***(weakly) for all $\xi \u2208H\xd7$, in particular,**

*ξ*

*Γ ∇**χ*

_{k}=

*∇**χ*

_{k}.

^{17,44}It is based on convolution with respect to the Green’s function for the Laplacian Δ = ∇

^{2}.

^{38,56}Applying the integro-differential operator

**Δ**

*∇*^{−1}to the cell problem in Eq. (D2) yields

**[(**

*Γ**ɛ*I + H)(

*∇**χ*

_{k}+

*e*_{k})] = 0. Since

*Γe*_{k}= 0 and

*Γ∇**χ*

_{k}=

*∇**χ*

_{k}, this formula is equivalent to (

*ɛ*I +

**H**

*Γ***)**

*Γ*

*∇**χ*

_{k}= −

**H**

*Γ*

*e*_{k}, which yields the following resolvent formula for

*∇**χ*

_{k}:

which is analogous to Eq. (12). Since ** Γ** is a projection operator onto $H\xd7\u2282H$, it is bounded by unity in operator norm on $H$, ∥

**∥ ≤ 1.**

*Γ*^{20,54}Integration by parts and the symmetry of the operator Δ

^{−1}

^{56}(or the projective nature of

**itself) shows that**

*Γ***is also a**

*Γ**symmetric*operator, i.e., ⟨

**⋅**

*Γξ***⟩ = ⟨**

*ζ***⋅**

*ξ***⟩ for all $\xi ,\zeta \u2208H$.**

*Γζ*^{44}These two properties show that

**with domain $H$ is a**

*Γ**self-adjoint*operator.

^{53}Since

**is self-adjoint and**

*Γ*

*Γ∇**χ*

_{k}=

*∇**χ*

_{k}, we may write $Ajk*$ in Eq. (D4) as $Ajk*=\u27e8A\u2009\u2207\chi j\u22c5\u2207\chi k\u27e9$. Consequently, substituting the resolvent formula for

*∇**χ*

_{k}in Eq. (D9) into the functional formulas in Eq. (D4) yields

which is a direct analog of Eq. (13).

Since ** Γ** is self-adjoint on $H$, the antisymmetry of the matrix H implies that

**A**=

**H**

*Γ***is an**

*Γ**antisymmetric*operator on $H$, i.e., ⟨

**A**⋅

*ξ***⟩ = −⟨**

*ζ***⋅**

*ξ***A**⟩. We emphasize that the operator

*ζ***A**depends only on the fluid velocity field via Eq. (28). By Eq. (D8), the operator

**A**is bounded on $H$ with ∥

**A**∥ ≤ ∥H∥ <

*∞*. This, the skew-symmetry of

**A**, and the sesquilinearity of the $H$-inner-product imply that

**M**= −

*ı*

**A**, where $\u0131=\u22121$, is a bounded symmetric operator, hence self-adjoint on $H$ with ∥

**M**∥ = ∥

**A**∥ <

*∞*. The spectrum Σ of the self-adjoint operator

**M**is real-valued with the spectral radius equal to its operator norm,

^{53}thus

We are now ready to present the main results of this section. We start with the following corollary of Theorem 1.

*Let*

*Q**(λ)*

*denote the resolution of the identity corresponding to the self-adjoint operator*

*M**. Then, the components*$Sjk*$

*and*$Ajk*$

*,*

*j, k = 1, …, d*

*, of the symmetric*$S*$

*and antisymmetric*$A*$

*parts of the effective diffusivity matrix*$D*$

*have the Stieltjes-Radon integral representations in (*

*19*

*) with*

*μ*

_{jk}

*replaced by the spectral measure*$\mu \u0303jk$

*of*

*M**in the*$(gj,gk)$

*state. The bounds in Theorem 2 hold for these integral representations for*$Sjk*$

*and*$Ajk*$

*. The mass*$\mu \u0303jk0$

*of the measure*$\mu \u0303jk$

*is real-valued and satisfies*

*λ*defined by $\mu \u0303jk(\lambda )=\u27e8Q(\lambda )gj,gk\u27e9$. Here, $gk=\u2212\Gamma Hek$ is defined in (D9) and {

**Q**(

*λ*)}

_{λ∈Σ}is the family of self-adjoint projection operators that is in one-to-one correspondence with the bounded linear self-adjoint operator

**M**on the Hilbert space $H\xd7$.

^{53,58}From Eq. (18) and the fact that

**is a self-adjoint projection operator on $H\xd7$, the mass $\mu \u0303jk0$ of the measure $\mu \u0303jk$ is real-valued and satisfies**

*Γ*It is worth making the following observation. In the current setting, the formulas for $Sjk*$ and $Ajk*$ in (19) are computed with respect to the standard basis {*e*_{k}}, through the definition of $\mu \u0303jk(\lambda )=\u27e8Q(\lambda )gj,gk\u27e9$ with $gk=\u2212\Gamma Hek$. We now show that, given $Sjk*$ and $Ajk*$, *j*, *k* = 1, …, *d*, the effective diffusivity matrix can be computed relative to any directions. This is due to the bilinearity of the inner-product underlying the definition of $\mu \u0303jk(\lambda )$. More specifically, if $\xi ,\zeta \u2208Rd$ are arbitrary directions of interest, then $\u27e8Q(\lambda )\Gamma H\xi ,\Gamma H\zeta \u27e9=\u2211jkajbk\u27e8Q(\lambda )gj,gk\u27e9$, where the constants *a*_{j} and *b*_{k}, *j*, *k* = 1, …, *d*, are the coordinates of the vectors ** ξ** and

**with respect to the standard basis. This immediately leads to integral representations for the effective diffusivity matrix relative to any desired directions. This observation had useful consequences in Refs. 17 and 18.**

*ζ**The spectral measure*

*μ*

_{jk}

*in Theorem 1 is identical to the spectral measure*$\mu \u0303jk$

*in Corollary 4,*

Before we prove Theorem 5, we give a brief outline of the proof and display the key formulas that lead to the result in (D14). The vector field $gk=\u2212\Gamma Hek$ defined in Eq. (D9) and the scalar field *g*_{k} = (−Δ)^{−1}*u*_{k} defined in (12) are (weakly) related by^{44}

The operator **A** = ** Γ**H

**defined in Eq. (D9) and the operator**

*Γ**A*= Δ

^{−1}[

**⋅**

*u***] defined in (12), hence**

*∇***M**= −

*ı*

**A**and

*M*= −

*ıA*, are (weakly) related by

Consequently, the masses and moments of the spectral measure *μ*_{jk}, *j*, *k* = 1, …, *d*, in Theorem 1 and the spectral measure $\mu \u0303jk$ in Corollary 4 are identical,

Since the self-adjoint operators **M** and *M* are bounded, hence the spectra of these operators are bounded subsets of $R$ and the Hausdorff moment problem is determinate,^{1,58} and this establishes Eq. (D14).

**=**

*u***(**

*u***,**

*x**t*) is periodic in both space and time

*t*. We will see that the proof of Eq. (D17) essentially depends only on Eqs. (D15) and (D16). Consequently, en route, we will also establish that Eq. (D17) holds for this more general, time-dependent context. Here, we will only sketch the key ideas that were developed at length in Ref. 44. Equation (28),

**=**

*u***⋅ H, and the definitions**

*∇***=**

*Γ***(Δ**

*∇*^{−1})

**⋅, $gk=\u2212\Gamma Hek$, and**

*∇**g*

_{k}= (−Δ)

^{−1}

*u*

_{k}, together yield

^{44}Eq. (D15),

*μ*

_{jk}and $\mu \u0303jk$ are equal,

*Γ∇**ξ*=

*∇**ξ*(weakly),

^{44}Eqs. (28) and (29),

**⋅ [H**

*∇***] = [**

*∇***⋅ H] ⋅**

*∇***, imply that for all $\xi \u2009\u2208H1,2$,**

*∇*^{44}which establishes Eq. (D16). It follows from (D15) and (D16) that

We conclude this section with the following corollary of Theorem 5.

*For each*$\u2009\xi \u2009\u2208H1,2$,

*we have*$\u2207\xi \u2208H\xd7\u2009$

*, and conversely, for each*$\xi \u2208H\xd7$,

*there exists unique*$\xi \u2009\u2208H1,2$

*such that*

*ξ**=*

*∇**ξ*

*. Consequently, by Theorem 5, for every*$\xi ,\zeta \u2009\u2208H1,2$

*and*$\xi ,\zeta \u2208H\xd7\u2009$

*related by*

*ξ**=*

*∇**ξ*

*and*

*ζ**=*

*∇**ζ*,

*we have*

*This establishes that the spectral measures*

*μ*

_{ξζ}

*and*$\mu \u0303\xi \zeta $

*are identical,*

*Moreover,*Eq. (

*D22*)

*also holds for the class of space-time periodic fluid velocity fields*

*u**=*

*u**(*

*x**, t)*

*described in Ref.*

*44*

*.*

The Hilbert spaces $H1,2$ and $H\xd7$ are in one-to-one isometric correspondence.^{44} More specifically, for every $\xi \u2009\u2208H1,2$, we have $\u2207\xi \u2208H\xd7$ satisfying ∥*ξ*∥_{1,2} = ∥*∇**ξ*∥ and ∥*Aξ*∥_{1,2} = ∥**A ∇**

*ξ*∥. Conversely, for each $\xi \u2208H\xd7$, there exists unique

^{44}$\xi \u2009\u2208H1,2$ (up to equivalence class) such that

**=**

*ξ*

*∇**ξ*, ∥

**∥ = ∥**

*ξ**ξ*∥

_{1,2}, and ∥

**A**∥ = ∥

*ξ**Aξ*∥

_{1,2}. By this and Theorem 5, we have Eqs. (D22) and (D23).

The proof of Eq. (D17) depends only on Eqs. (D15) and (D16), which also hold^{44} for the class of space-time periodic fluid velocity fields ** u** =

**(**

*u***,**

*x**t*) described in Ref. 44. The argument in the previous paragraph above also holds for this time-dependent setting.

^{44}Therefore, Eq. (D22) holds for the time-dependent setting as well. This concludes our Proof of Corollary 6.

In order to establish Eq. (D23) for the time-dependent setting associated with operators having unbounded spectra,^{44} one must first establish Carleman’s criterion,^{1} for example, associated with the Stieltjes or Hamburger moment problems. This requires an involved analysis that is outside of the scope of the current work. Our results in this direction will be published elsewhere.

### APPENDIX E: DISCRETE SETTING: HILBERT SPACE OF CURL-FREE FIELDS

Here, we provide a matrix formulation for the effective parameter problem discussed in Appendix D 2, which provides discrete versions of the Stieltjes integral representations for $D*$ shown in Eq. (19). These integrals involve a discrete spectral measure *μ*_{jk} analogous to the discrete measure in Eq. (35). More specifically, we use a discretized version of the cell problem in Eqs. (D2) and (D9), written as $(\epsilon I+A)\u2207\chi k=gk$, to express the discrete spectral measure *μ*_{jk} explicitly in terms of the eigenvalues and eigenvectors of a matrix representation A of the operator **A** = ** Γ**H

**. We then develop a**

*Γ**projection method*which additionally shows that the discrete measure

*μ*

_{jk}is actually determined by the eigenvalues and eigenvectors of a matrix that is much smaller than the matrix A. This projection method is used in Appendix G to establish the equivalence of the two discrete approaches given in this appendix and Sec. IV, for the setting where the matrix gradient ∇ is rank-deficient. This equivalence proof establishes, en route, that the common discrete spectral measure can be computed by a method that combines the computational benefits of both approaches.

From the discussion in Sec. IV, the discrete representation of the projection operator ** Γ** =

**(Δ**

*∇*^{−1})

**⋅ is given by the symmetric projection matrix $\Gamma =\u2207(\u2207T\u2207)\u22121\u2207T$ satisfying Γ**

*∇*^{2}= Γ and Γ∇ = ∇, where $(\u2207T\u2207)\u22121$ is now interpreted as a matrix inversion. We assume here that the matrix ∇ is of full-rank so that $(\u2207T\u2207)\u22121$ exists. The rank-deficient case, where the matrix ∇

^{T}∇ is singular, is examined in Appendix G. In this way, the integro-differential operator

**A**=

**H**

*Γ***is represented by an antisymmetric matrix A = ΓHΓ satisfying A**

*Γ*^{T}= −A. In a similar way, the vectors $gk=\u2212\Gamma Hek$,

*k*= 1, …,

*d*, are redefined for this matrix setting. For simplicity, we will not make a notational distinction between the continuum and discrete cases for the vectors $gk$ and

*e*_{k}as well as the matrix H, as the context will be clear.

The spectrum of the antisymmetric matrix A of size *N*, say, is composed of eigenvalues *υ*_{n}, *n* = 1, …, *N*, with the corresponding eigenvectors *w*_{n} satisfying A*w*_{n} = *υ*_{n}*w*_{n}. Since A is skew-symmetric, the eigenvalues *υ*_{n} are purely imaginary,^{25} *υ*_{n} = *ıλ*_{n} with $\lambda n\u2208R$. Therefore, the matrix M = −*ı*A is Hermitian (M^{†} = M) and it has the same eigenvectors *w*_{n} as the matrix A and real eigenvalues given by *λ*_{n} = Im *υ*_{n}. The eigenvectors *w*_{n}, *n* = 1, …, *N*, of the Hermitian matrix M form an orthonormal basis for $CN$,^{25,28} i.e., $wn\u2009\u2020wm=\delta nm$ and for every $\xi \u2208CN$ we have $\xi =\u2211n(wn\u2020\xi )wn=\u2211nwnwn\u2020\xi \u2009$. Consequently, defining $Qn=wnwn\u2020$, *n* = 1, …, *N*, to be the mutually orthogonal Hermitian projection matrices onto the eigenspaces spanned by the *w*_{n}, we have the following analog of Eq. (34):

With an argument similar to the one in the paragraph following Eq. (34), one can show that for all $\xi ,\zeta \u2208CN$ and complex-valued polynomials *f*(*λ*) and *h*(*λ*), the bilinear functional ⟨*f*(M)** ξ** ⋅

*h*(M)

**⟩ has the integral representation in Eq. (17), with**

*ζ**M*substituted by M and the scalars

*ξ*and

*ζ*replaced by vectors

**and**

*ξ***. Moreover, the complex-valued function**

*ζ**μ*

_{ξζ}(

*λ*) = ⟨Q(

*λ*)

*ξ*,

*ζ*⟩ in Eq. (17) is now given by $\mu \xi \zeta (\lambda )=\u27e8Q(\lambda )\xi \u22c5\zeta \u27e9$, where the associated matrix representation $Q(\lambda )$ of the projection operator Q(

*λ*) and the discrete spectral measure d

*μ*

_{ξζ}(

*λ*) are given by the following analog of Eq. (35):

We are now ready to present the main results of this section, given in the following corollary of Theorem 3.

*Consider the standard eigenvalue problem*

*M*

*w*_{n}

*= λ*

_{n}

*w*_{n}

*associated with the Hermitian matrix*

*M*

*= −ı*

*A*

*. Let*

*W*

*be the matrix with columns consisting of the eigenvectors*

*w*_{n}

*and*Λ

*be the diagonal matrix with eigenvalues*

*λ*

_{n}

*on its diagonal. The discrete, matrix representations of the bilinear functional formulas for*$Sjk*$

*and*$Ajk*$

*in*

*Eq. (D4)*

*with*

*H*

*replaced by*

*A*

*, as discussed right before*Eq. (

*D10*)

*, are given by*

*which is analogous to Eq. (*

*36*

*). Also, the discrete representation of the resolvent formula for*

*∇*

*χ*

_{j}

*in Eq. (*

*D9*

*) is given by*

*which is analogous to Eq. (*

*37*

*). The discrete representations of the bilinear functional formulas for*$Sjk*$

*and*$Ajk*$

*in (*

*D10*

*) are given by*

*which are analogous to those in Eq. (*

*38*

*). Consequently, the formulas for*$Sjk*$

*and*$Ajk*$

*in (*

*E5*

*) have the following series representations:*

*which is analogous to Eq. (*

*39*

*). Finally, recalling the projection matrix*$Qn=wnwn\u2020$

*in (*

*E1*

*), we have*

*which is analogous to Eq. (*

*40*

*). It follows from Eq. (*

*E7*

*) that the series representations for*$Sjk*$

*and*$Ajk*$

*in (*

*E6*

*) have the Stieltjes integral representations in Eq. (*

*19*

*), involving the discrete spectral measure*

*μ*

_{jk}

*in (*

*E2*

*) with*$\xi =gj=\u2212\Gamma Hej$

*and*$\zeta =gk=\u2212\Gamma Hek\u2009$

*.*

*χ*_{j}⋅ ∇

*χ*_{k}⟩ = ⟨A ∇

*χ*_{j}⋅ ∇

*χ*_{k}⟩. Consider the spectral decomposition M = WΛW

^{†}of the Hermitian matrix M = −

*ı*A involving the real eigenvalues

*λ*

_{n}comprising the main diagonal of the diagonal matrix Λ and the eigenvectors

*w*_{n}comprising the columns of the unitary matrix W satisfying W

^{†}W = WW

^{†}= I.

^{25}Equation (E4) follows from this spectral decomposition of A =

*ı*M, the discrete version of the resolvent formula in Eq. (D9) written as $\u2207\chi j=(\epsilon I+A)\u22121gj$, the fact that W is unitary satisfying W

^{†}= W

^{−1}, and the elementary properties of matrix inversion.

^{25}Substituting the resolvent formula for ∇

*χ*_{j}in (E4) into Eq. (E3) and using W

^{†}= W

^{−1}yield Eq. (E5). The quadratic form $W\u2020gj\u22c5W\u2020gk$ arising in (E5) can be written in terms of the projection matrices $Qn=wnwn\u2020$ defined in (E1) as follows:

*μ*

_{jk}(

*λ*) in Eq. (E2) with $\xi =gj=\u2212\Gamma Hej$ and $\zeta =gk=\u2212\Gamma Hek$. This concludes our Proof of Corollary 7.

We conclude this section with the development of a *projection method* that is used in the Proof of Theorem 10. This method shows that the presence of the projection matrix Γ in the operator ΓHΓ and $gk=\u2212\Gamma Hek$ projects out contributions of the null space of Γ in the series representation for $D*$ in (E6). This is reminiscent of problems encountered for many elliptic equations on periodic domains, where the analysis requires a deficient dimension to be projected out.^{7} Theorem 10 establishes that the discrete framework developed in Sec. IV and this section yield equivalent Stieltjes integral representations for the effective diffusivity matrix $D*$, involving a discrete spectral measure, and also generalizes this result to the setting where the matrix ∇ is rank-deficient. The Proof of Theorem 10 demonstrates that the common discrete spectral measure can be computed by a method that combines the computational benefits of both approaches. This, in turn, is used in our numerical computations of Stieltjes integral representations of $D*$ in Sec. V.

*The real-symmetric projection matrix*

*Γ*

*of size*

*N*

*has eigenvalues satisfying*

*γ*

_{n}

*= 0, 1*

*, hence the spectral decomposition*

*Here,*$0N0$

*and*$1N1$

*are vectors of zeros and ones with*

*N*

_{0}

*and*

*N*

_{1}

*components, respectively, where*

*N = N*

_{0}

*+ N*

_{1}

*. Moreover, the columns comprising the*

*N × N*

_{0}

*matrix*

*P*

_{0}

*and the*

*N × N*

_{1}

*matrix*

*P*

_{1}

*are orthonormal eigenvectors that span the null space and range of*

*Γ*

*, respectively. Consequently, the matrix*

*Γ*

*can be written as*

*Consider the spectral composition of the antisymmetric matrix*$P1THP1$

*of size*

*N*

_{1}

*,*

*where*

*R*

_{11}

*is a unitary matrix,*$R11\u2020R11=R11R11\u2020=I11$

*,*

*Λ*

_{11}

*is a real-valued diagonal matrix, and*

*I*

_{11}

*is the identity matrix of size*

*N*

_{1}

*× N*

_{1}

*. Consequently, the matrix*

*A*

*= ıΓ*

*H*

*Γ*

*has the spectral decomposition*

*Equations (*

*E10*

*) and (*

*E12*

*), and the mutual orthogonality of the matrices*

*P*

_{0}

*and*

*P*

_{1}

*yield*

*where*

*O*

_{0}

*is a matrix of zeros of size*

*N × N*

_{0}

*. Consequently, Eq. (*

*E8*

*) can be written as*

*Moreover, the spectral weights in Eq. (*

*E8*

*) associated with*

*γ*

_{n}

*= 0*

*are identically zero*$[Qngj\u22c5gk]=0\u2009$

*.*

Using the spectral decomposition Γ = PGP^{T} in (E9), we write the matrix A = ΓHΓ as $A=P[G(PTHP)G]PT=P\u2009diag(O00\u2009P1THP1)\u2009PT$, where O_{00} is a matrix of zeros of size *N*_{0} × *N*_{0}. Due to the skew-symmetry of H, the *N*_{1} × *N*_{1} matrix $P1THP1$ is also skew-symmetric. Consequently, it has the spectral decomposition in Eq. (E11); hence, $A=P\u2009diag(O00\u2009\u0131R11\Lambda 11R11\u2020)\u2009PT$. Writing this in block matrix form^{42} and multiplying yield Eq. (E12). Equations (E10) and (E12), and the mutual orthogonality of the matrices P_{0} and P_{1} yield Eq. (E13). This concludes our Proof of Theorem 8.

### APPENDIX F: DISCRETE EQUIVALENCE OF THE EFFECTIVE PARAMETER PROBLEMS

In Sec. III, we provided Stieltjes integral representations for the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$, associated with an incompressible fluid velocity field. A discrete version of this mathematical framework was formulated in Sec. IV. An alternate approach to the effective parameter problem was formulated in Appendix D, and its discrete version was formulated in Appendix E. In this section, we demonstrate that the discrete versions of these effective parameter problems yield equivalent spectral representations of $D*$ when the matrix ∇ is of full-rank, as in the case of Dirichlet boundary conditions, so that the matrix Laplacian is invertible. In Appendix G, this result is extended to the setting where ∇ is rank-deficient, as in the case of periodic boundary conditions.

Let ∇ = UΣV^{T} be the singular value decomposition (SVD) of the matrix ∇ of size *N* × *K*, where *K* = *L*^{d} and *N* = *Kd*. Here, Σ = diag(*σ*_{1}, …, *σ*_{K}), where 0 ≤ *σ*_{1} ≤ ⋯ ≤ *σ*_{K}, and the matrices U and V are of size *N* × *K* and *K* × *K*, respectively, and satisfy^{14}

where I is the *K* × *K* identity matrix. The columns of U are called left singular vectors, the columns of V are called right singular vectors, and the *σ*_{i} are called singular values.

It follows from ∇ = UΣV^{T} and Eq. (F1) that the spectral decomposition of the negative matrix Laplacian ∇^{T}∇ is given by ∇^{T}∇ = VΣ^{2}V^{T}.^{14} We assume that ∇ is of full-rank so that *σ*_{i} > 0 for all *i* = 1, …, *K*. This implies that Σ^{−1} exists so that the matrix Laplacian is invertible. In this case, it follows from ∇ = UΣV^{T} and Eq. (F1) that the projection matrix $\Gamma =\u2207(\u2207T\u2207)\u22121\u2207T$ is given by

which is an *N* × *N* symmetric projection matrix satisfying Γ^{2} = Γ and Γ∇ = ∇. A key property of the SVD of the *full-rank* matrix ∇ is that its range is spanned by the columns of U;^{14} hence, Γ = UU^{T} projects subspaces of $RN$ onto the range of ∇.

From Eqs. (F1) and (F2), we can write the eigenvalue problem −*ı*ΓHΓ*w*_{n} = *λ*_{n}*w*_{n} discussed in Appendix E as

Now consider the generalized eigenvalue problem −*ı*∇^{T}H∇*z*_{n} = *α*_{n}∇^{T}∇*z*_{n} discussed in Sec. IV and recall that ∇ = UΣV^{T} and ∇^{T}∇ = VΣ^{2}V^{T}. Since Σ is invertible, by Eq. (F1), we can write this generalized eigenvalue problem as the following standard eigenvalue problem:

Comparing the formulas in Eqs. (F3) and (F4) indicates that spectrum associated with each of these eigenvalue problems is identical, *α*_{n} = *λ*_{n}, and that the eigenvectors are related by U^{T}*w*_{n} = ΣV^{T}*z*_{n}. Since Γ is a projection matrix, Γ^{2} = Γ, the eigenvalue problem ΓHΓ*w*_{n} = *ıλ*_{n}*w*_{n} can be written as ΓHΓ[Γ*w*_{n}] = *ıλ*_{n}[Γ*w*_{n}], which implies that Γ*w*_{n} = *w*_{n}. Consequently, applying the matrix U to both sides of the formula U^{T}*w*_{n} = ΣV^{T}*z*_{n} and recalling that Γ = UU^{T} and ∇ = UΣV^{T}, we have

In the following lemma, we make precise the correspondence between the standard eigenvalue problem −*ı*ΓHΓ*w*_{n} = *λ*_{n}*w*_{n} and the generalized eigenvalue problem −*ı*∇^{T}H∇*z*_{n} = *α*_{n}∇^{T}∇*z*_{n}, as well as the associated spectral measures in Eqs. (E2) and (35), respectively.

*Consider the standard eigenvalue problem and the generalized eigenvalue problem given, respectively, in the following equations:*

*Let*

*∇ =*

*U*

*Σ*

*V*

^{T}

*be the SVD of the matrix*

*∇*

*, which we assume to be of full-rank. Then, Eq. (*

*F6*

*) implies and is implied by Eq. (*

*F7*

*), with*

*w*_{n}

*and*

*z*_{n}

*related as in Eq. (*

*F5*

*). This implies that the spectrum associated with each of these eigenvalue problems is identical. Moreover, the spectral weights in Eqs. (*E8

*) and (*40 )

*are identical; specifically,*

*This, in turn, implies that the associated spectral measures in*Eqs. (

*E2*)

*and (*

*35*

*) are identical.*

^{T}, ∇

^{T}∇ = VΣ

^{2}V

^{T}, and Γ = UU

^{T}, where Σ is invertible, and the matrices V and U satisfy Eq. (F1). First, consider Eq. (F6) written as in Eq. (F3), [−

*ı*U

^{T}HU][U

^{T}

*w*_{n}] =

*λ*

_{n}[U

^{T}

*w*_{n}]. Since the matrix Σ is invertible and V

^{T}V = I, we can rewrite Eq. (F3) as

^{T}with

*z*_{n}= VΣ

^{−1}U

^{T}

*w*_{n}. This formula for

*z*_{n}, Eq. (F1), and the formula Γ

*w*_{n}=

*w*_{n}above Eq. (F5) imply that

*w*_{n}= UΣV

^{T}

*z*_{n}= ∇

*z*_{n}, which is the formula in (F5). Now consider Eq. (F7) written as in (F4), [−

*ı*U

^{T}HU][ΣV

^{T}

*z*_{n}] =

*λ*

_{n}[ΣV

^{T}

*z*_{n}]. Since U

^{T}U = I, we can rewrite Eq. (F4) as

^{T}with

*w*_{n}= UΣV

^{T}

*z*_{n}= ∇

*z*_{n}.

**=**

*u***⋅ H in (28), for the continuum setting, we have that**

*∇**u*

_{j}= (

**⋅ H) ⋅**

*∇*

*e*_{j}=

**⋅ (H**

*∇*

*e*_{j}). Since ∇ = UΣV

^{T}, the discrete version of this formula is given by

^{T}and $(\u2207T\u2207)\u22121=V\Sigma \u22122VT$, we have $\u2207(\u2207T\u2207)\u22121=U\Sigma \u22121VT$. Consequently, Γ = UU

^{T}, and Eqs. (F1) and (F11) yield $\u2207(\u2207T\u2207)\u22121uj=\u2212\Gamma Hej$. Equation (F8) now follows from the formula

*w*_{n}= ∇

*z*_{n}in (F5),

^{25}The equivalence of Eqs. (F8) and (F12) now follows from Eqs. (E1), (34), and (40). This concludes our Proof of Lemma 9.

We conclude this section with a discussion regarding numerical computations of the effective diffusivity matrix $D*$. The approach discussed in this section and the projection method discussed in Theorem 8 demonstrate that a hybrid of these two approaches leads to the most efficient algorithm for numerical computations of spectral representations for $D*$—combining the computational advantages of both the methods discussed in Appendix E and Sec. IV. More specifically, in the full rank setting, the spectral measure underlying the discrete integral representation for $D*$ was calculated in Appendix E in terms of the *standard* eigenvalue problem −*ı*ΓHΓ*w*_{m} = *λ*_{n}*w*_{n}, where the matrix −*ı*ΓHΓ is of size *N* × *N*. In Sec. IV, $D*$ was calculated in terms of the *generalized* eigenvalue problem −*ı*∇^{T}H∇*z*_{n} = *λ*_{n}∇^{T}∇*z*_{n}, involving the *K* × *K* matrices −*ı*∇^{T}H∇ and ∇^{T}∇. Since $\u2207T=(\u22071T,\u2026,\u2207dT)$ is of size *K* × *N*, we have that *K* = *N*/*d* < *N*. However, the generalized eigenvalue problem is more computationally intensive than the standard eigenvalue problem.^{47} For the case of randomly perturbed flows, many statistical iterations are necessary to compute $D*$ and the efficiency of the numerical algorithm is key. Neither of the methods discussed in Appendix E and Sec. IV are optimal.

The projection method developed in Theorem 8 demonstrates that, by first computing the standard eigenvalue decomposition of the nonrandom matrix Γ, the spectral statistics of the eigenvalue problem −*ı*ΓHΓ*w*_{n} = *λ*_{n}*w*_{n} can then be obtained by repeatedly computing the standard eigenvalue decomposition of smaller matrices. They are of size *K* × *K* by Eqs. (E10)–(E12) and (F2). We emphasize that in the setting of full-rank ∇, the parameter *K* in this section and *N*_{1} in Theorem 8 both denote the rank of the matrix Γ, i.e., *K* = *N*_{1}. Note, computing the matrix $\Gamma =\u2207(\u2207T\u2207)\u22121\u2207T$ itself involves the cost of numerically solving *N* linear systems of size *K* × *K*. Alternatively, the Proof of Lemma 9 illustrates that by first computing the SVD of the matrix gradient, ∇ = UΣV^{T}, the spectral statistics of the generalized eigenvalue problem − *ı*∇^{T}H∇*z*_{n} = *λ*_{n}∇^{T}∇*z*_{n} can then be obtained by repeatedly computing the *standard* eigenvalue decomposition of the matrix −*ı*U^{T}HU which is of size *K* × *K*. When *N* is large, these equivalent methods are more numerically efficient than the other approaches discussed in Appendix E and Sec. IV.

In Appendix G, we generalize Lemma 9 to the case where ∇ is rank-deficient with rank *K*_{1} satisfying *K*_{1} < *K*. Our analysis demonstrates that the two formulations discussed in Appendix E and Sec. IV yield equivalent spectral representations of the effective diffusivity matrix $D*$ in this rank-deficient setting. Moreover, en route, a more efficient numerical algorithm for computations of $D*$ is revealed. More specifically, we demonstrate that by first computing the SVD of the matrix gradient, $D*$ can be computed via a *standard* eigenvalue problem for matrices of size *K*_{1} × *K*_{1}. Consequently, the rank deficiency of the problem actually increases the numerical efficiency of computations.

### APPENDIX G: RANK DEFICIENCY AND A UNIFYING STANDARD EIGENVALUE PROBLEM

In Sec. IV and Appendix E, we provided two discrete, matrix formulations of the effective parameter problem for advection enhanced diffusion, which led to discrete Stieltjes integral representations for the effective diffusivity matrix $D*$ involving spectral measures of Hermitian matrices. These two formulations assume that the *N* × *K* matrix ∇ is of full-rank *K* so that the negative matrix Laplacian ∇^{T}∇ is invertible. Lemma 9 of Appendix F shows that, given this condition, the two formulations yield equivalent spectral representations of $D*$. This analysis also demonstrates that a hybrid of the two formulations leads to a numerical algorithm for computing $D*$ that is more efficient than the numerical algorithms for either approach. In this section, we generalize Lemma 9 to the case where ∇ is rank-deficient, with rank *K*_{1} satisfying *K*_{1} < *K*. We demonstrate that the two formulations are equivalent in this rank-deficient setting and we also derive an efficient hybrid numerical algorithm for computations of $D*$. This framework is used in Sec. V to compute $D*$ for periodic flows, for which the matrix ∇ with periodic boundary conditions is rank-deficient.

Toward this goal, let ∇ = UΣV^{T} be the SVD of the matrix gradient ∇ of size *N* × *K*, introduced in Sec. IV and Appendix F. We assume that ∇ is rank deficient so that ∇^{T}∇ = VΣ^{2}V^{T} is singular, with *K*_{1} nonzero eigenvalues and *K*_{0} = *K* − *K*_{1} zero eigenvalues, and write

Here, we denote O_{ab}, *a*, *b* = 0, 1, to be matrices of zeros of size *K*_{a} × *K*_{b}, U_{a} is of size *N* × *K*_{a}, V_{a} is *K* × *K*_{a}, and Σ_{1} is a *K*_{1} × *K*_{1} diagonal, *invertible* matrix. By Eq. (F1), the matrices U_{1} and V_{1} satisfy $U1TU1=I1$ and $V1TV1=I1\u2009$, where I_{1} is the *K*_{1} × *K*_{1} identity matrix, but $V1V1T\u2260I$. Due to the blocks of zeros in the matrix Σ in Eq. (G1), we can write the matrix gradient as $\u2207=U1\Sigma 1V1T$ and the negative matrix Laplacian as $\u2207T\u2207=V1\Sigma 12V1T$. An important property of the SVD of the matrix ∇ is that its null space is spanned by the columns of V_{0} and its range is spanned by the columns of U_{1}.^{14} We emphasize that in the setting where ∇ is full-rank, we have U_{1} = U, Σ_{1} = Σ, and V_{1} = V satisfying $V1V1T=V1TV1=I$.

Consider the cell problem in Eq. (4) written, via (28) and [** ∇** ⋅ H] ⋅

*∇**φ*=

**⋅ [H**

*∇*

*∇**φ*], as

**⋅ H**

*∇*

*∇**χ*

_{j}+

*ɛ*Δ

*χ*

_{j}= −

*u*

_{j}. Discretizing this formula yields (see the discussion following Eq. (29) for details regarding the discretization of these differential operators, etc.)

where *u*_{j} is the discrete, vector representation of the *j*th component of the fluid velocity field *u*_{j} and similarly for *χ*_{j}. Substituting the formula for *u*_{j} in (G2) into the discrete version $Djk*=\epsilon \delta jk+\u27e8uj\u22c5\chi k\u27e9$ of Eq. (3) yields

where, as before, $Skj*=Sjk*$ and $Akj*=\u2212Ajk*$. We are now ready to state the key result of this appendix.

*Let the matrix gradient*

*∇*

*be rank-deficient and let*$\u2207=U1\Sigma 1V1T$

*be its SVD. Also, let*$U1THU1=\u0131R1\Lambda 1R1\u2020$

*be the spectral decomposition of the antisymmetric matrix*$U1THU1$

*. Consider the discrete formulation of the effective parameter problem developed in*4

*Sec. IV*

*. We have the following generalization of Eq. (*

*34*

*):*

*where the matrices*$Qn1$

*,*

*n = 1, …, K*

_{1}

*, are self-adjoint with respect to the discrete inner-product*

*⟨*⋅

*,*⋅

*⟩*

_{1,2}

*defined by*

*⟨*

*ξ**,*

*ζ**⟩*

_{1,2}

*= ⟨∇*

**⋅**

*ξ**∇*

*ζ**⟩*

*, i.e.,*$\u27e8Qn1\xi ,\zeta \u27e91,2=\u27e8\xi ,Qn1\zeta \u27e91,2\u2009$

*for*$\xi ,\zeta \u2208CK1$

*. Moreover, the generalization of the resolvent formula*

*χ*_{j}

*=*

*Z*

*(ɛ*

*I*

*+ ıΛ)*

^{−1}

*Z*

^{†}

*u*_{j}

*in Eq. (*

*37*

*) is given by*

*Now consider the discrete formulation of the effective parameter problem developed in* *Appendix E**. Let* $\Gamma 1H\Gamma 1=\u0131W1\Lambda \u0303W1\u2020$ *be the spectral decomposition of the antisymmetric matrix* *Γ*_{1}*H**Γ*_{1}*, where* $\Gamma 1=U1U1T$*,* $\Lambda \u0303$ *is a diagonal real-valued matrix with* $\lambda \u0303n1$*,* *n = 1, …, N**, on its diagonal, and* *W*_{1} *is a unitary matrix with columns* $wn1\u2009$*. Then, Eqs. (**E1**) and (**E2**) hold with* $Qn$ *and* *λ*_{n} *replaced by* $Qn1=[wn1][wn1]\u2020$ *and* $\lambda \u0303n1\u2009$*. Also, the resolvent formula in Eq. (**E4** holds with* *W* *replaced by* *W*_{1} *and* Λ *replaced by* $\Lambda \u0303\u2009$*.*

*These two discrete formulations of the effective parameter problem are related as follows: The diagonal eigenvalue matrices*$\Lambda \u0303$

*and*Λ

_{1}

*are related by*$\Lambda \u0303=diag(O00,\Lambda 1)$

*. The eigenvector matrices*

*W*

_{1}

*and*

*Z*

_{1}

*are related by the following generalization of Eq. (*

*F5*

*) [also see (*

*E12*

*)]:*

*where the columns of*

*P*

_{0}

*are eigenvectors of*

*Γ*

_{1}

*which span its null space. Moreover, the following formulas generalize Eqs. (*

*E14*

*), (*

*F8*

*), and (*

*F12*

*):*

*In this rank-deficient setting, these two approaches both yield discrete Stieltjes integral representations for the functional formulas in* (*G3*) *for the symmetric* $S*$ *and antisymmetric* $A*$ *parts of the effective diffusivity matrix* $D*$*. The two representations are equivalent by the relations discussed above and are given by* Eq. (*39)**, with* *z*_{n} *and* *λ*_{n} *replaced by* $zn1=V1\Sigma 1\u22121rn1$ *and* $\lambda n1$*,* *n = 1, …, K*_{1}*, where* $(\lambda n1,rn1)$ *are eigenpairs of matrix* $\u2212\u0131U1THU1$*. The results discussed here also hold for the setting where* *∇* *is of full-rank and therefore generalize the discrete mathematical frameworks developed in Sec. IV* and Appendixes E and F*.*

We first work with Eq. (G2) directly and develop a mathematical framework which parallels the framework of Sec. IV for the rank-deficient setting. We then transform Eq. (G2) into a discrete analog of Eq. (D9) written as $\epsilon I+\Gamma H\Gamma \u2207\chi j=\u2212\Gamma Hej$, with a suitable generalization of the formula for the matrix Γ in (F2), and develop a mathematical framework which parallels the framework of Appendix E. We then generalize Lemma 9 of Appendix F, establishing the equivalence of these two formulations for the rank-deficient setting and, en route, derive a hybrid numerical algorithm for computing spectral representations of $D*$ which is more efficient than both of the other numerical algorithms.

Since Σ = diag(O_{00}, Σ_{1}), we have $\u2207=U1\Sigma 1V1T$ and $\u2207T\u2207=V1\Sigma 12V1T$, so the cell problem in (G2) can be written as $[V1\Sigma 1][U1THU1][\Sigma 1V1T]\chi j+\epsilon V1\Sigma 12V1T\chi j=uj$. The *K*_{1} × *K*_{1} antisymmetric matrix $U1THU1$ has the spectral decomposition $U1THU1=\u0131R1\Lambda 1R1\u2020$, where Λ_{1} is a diagonal real-valued matrix and R_{1} is a unitary matrix, $R1\u2020R1=R1R1\u2020=I1\u2009$. Equation (G5) follows from these formulas, which is an analog of Eq. (37). The formula $\u2207=U1\Sigma 1V1T$ and Eq. (G5), in turn, imply that $\u2207\chi j=U1R1(\epsilon I1+\u0131\Lambda 1)\u22121[V1\Sigma 1\u22121R1]\u2020uj$.

*χ*_{j}into the formula for $Sjk*$ in Eq. (G3) and using $U1TU1=I1$ and $R1\u2020R1=I1$ yield the functional formula for $Sjk*$ in (38) with Λ replaced by Λ

_{1}and Z replaced by $Z1=V1\Sigma 1\u22121R1$. We also have

_{1}and Z replaced by $Z1=V1\Sigma 1\u22121R1$. The quadratic form $Z1\u2020uj\u22c5Z1\u2020uk$ arising in these functional formulas yields Eq. (42) with

*z*_{n}replaced by $zn1=V1\Sigma 1\u22121rn1$ and other appropriate notational changes, where $rn1$,

*n*= 1, …,

*K*

_{1}, are the orthonormal eigenvectors of the matrix $U1THU1$ which comprise the columns of R

_{1}. The formula $zn1=V1\Sigma 1\u22121rn1$ can be written as $\u2207zn1=U1rn1$. The orthogonality properties of R

_{1}and U

_{1}then imply that the vectors $zn1$ satisfy the Sobolev-type orthogonality condition in Eq. (32), $\u2207zn1\u22c5\u2207zm1=\delta nm$. Moreover, since $\u2211nrn1[rn1]\u2020=I1$, we also have Eq. (G4). Consequently, the functional representations of $Sjk*$ and $Ajk*$ in Eq. (G3) have the discrete integral representations in Eq. (39) with

*z*_{n}and

*λ*

_{n}replaced by $zn1$ and $\lambda n1$, and $(\lambda n1,rn1)$ are eigenpairs of Hermitian matrix $\u2212\u0131U1THU1$ of size

*K*

_{1}. These summations have the same properties as the summations discussed in the Proof of Theorem 3 and Eqs. (44) and (45).

We now argue that the mathematical framework developed so far in this proof generalizes the full-rank case in Sec. IV. Indeed, in the full-rank setting, the matrix V_{1} = V is orthogonal, Σ_{1} = Σ is invertible, and R_{1} = R is orthogonal so that the matrix Z_{1} = Z defined in (G5) is given by Z = VΣ^{−1}R and is invertible with Z^{−1} = R^{†}ΣV^{T}. Equation (33) is satisfied with Z^{−†} = VΣR and Λ = Λ_{1} which, in turn, establishes the current rank-deficient setting generalizes the full-rank setting summarized by Eqs. (30)–(42).

We now generalize the mathematical framework developed in Appendix E to the case that the matrix ∇ is rank deficient. Using $\u2207=U1\Sigma 1V1T$, $\u2207T\u2207=V1\Sigma 12V1T$, $U1TU1=V1TV1=I1$, and the invertibility of the matrix Σ_{1}, the cell problem in (G2) can be written as $U1[U1THU1][U1TU1][\Sigma 1V1T]\chi j+\epsilon U1\Sigma 1V1T\chi j=U1\Sigma 1\u22121V1Tuj$. Substituting *u*_{j} = −∇^{T}H*e*_{j} into this expression yields $\epsilon I+\Gamma 1H\Gamma 1\u2207\chi j=gj1$, where $\Gamma 1=U1U1T$ and $gj1=\u2212\Gamma 1Hej$ which is analogous to Eq. (D9). As in Appendix F, the matrix $\Gamma 1=U1U1T$ projects subspaces of $RN$ onto the *range* of ∇. The antisymmetric matrix Γ_{1}HΓ_{1} has the spectral decomposition $\Gamma 1H\Gamma 1=\u0131W1\Lambda \u0303W1\u2020$, where $\Lambda \u0303$ is a diagonal real-valued matrix and W_{1} is a unitary matrix $W1\u2020W1=W1W1\u2020=I$, which yields a generalization of (E4). From Γ_{1}∇ = ∇, we have ⟨∇^{T}H∇*χ*_{j} ⋅ *χ*_{k}⟩ = ⟨Γ_{1}HΓ_{1}∇*χ*_{j} ⋅ ∇*χ*_{k}⟩. Exactly as in Appendix E, these formulas lead to generalizations of Eqs. (E1)–(E14), with appropriate notational changes for the rank-deficient setting. In the case that the matrix ∇ is of full-rank, we have U_{1} = U hence W_{1} = W, which establishes that the mathematical framework discussed in this paragraph generalizes the mathematical framework in Appendix E for the full-rank setting.

We now establish that these two different approaches provide equivalent spectral representations for the effective diffusivity matrix $D*$ for the rank-deficient setting. In Eq. (E9) of Theorem 8, we wrote the eigenvalue decomposition Γ = PGP^{T}, where P is an orthogonal matrix and G is a diagonal matrix. Moreover, we wrote P = [P_{0} P_{1}], where the columns of the matrices P_{0} and P_{1} are orthonormal eigenvectors that span the null space and range of Γ, respectively. Since the eigenvalues *γ*_{n} associated with the eigenvectors in the matrix P_{1} satisfy *γ*_{n} = 1, any linear combination of the corresponding eigenvectors is also an eigenvector of Γ with eigenvalue *γ*_{n} = 1. Therefore, since the orthonormal columns of the matrix U_{1} span the range of $\Gamma 1=U1U1T$, we may take P_{1} = U_{1} so that P = [P_{0} U_{1}]. Consequently, we can rewrite Eq. (E11) as $U1THU1=\u0131R11\Lambda 11R11\u2020$. Identifying the notations of Appendix E with this section, we have R_{11} = R_{1} and Λ_{11} = Λ_{1}. This equation and Eq. (E12) establish that $\Lambda \u0303=diag(O00,\Lambda 1)$.

We now establish that the spectral weights associated with each approach are identical. Using the notation of this section, Eq. (E12) yields W_{1} = [P_{0} U_{1}R_{1}]. Consequently, from $\u2207=U1\Sigma 1V1T$, $U1TU1=V1TV1=I1$, and the formula $Z1=V1\Sigma 1\u22121R1$ in (G5), we have that ∇Z_{1} = U_{1}R_{1}, implying Eq. (G6), which is a generalization of Eq. (F5). Equation (G7) now follows from $\Gamma 1T=\Gamma 1$, Γ_{1}P_{0} = O, Γ_{1}U_{1} = U_{1}, and the formula *u*_{j} = −∇^{T}H*e*_{j} in (F11). This establishes that the spectral weights associated with both of the approaches discussed in Sec. IV and Appendix E are identical which, in turn, establishes that both approaches provide equivalent spectral representations of the effective diffusivity matrix $D*$. This concludes our Proof of Theorem 10.

In Sec. V, spectral representations of the symmetric $S*$ and antisymmetric $A*$ parts of the effective diffusivity matrix $D*$ are computed for various periodic flows, and spectral characteristics are related to flow geometry and transport properties. We accomplish this by using Eq. (39) with *z*_{n} and *λ*_{n} replaced by $zn1=V1\Sigma 1\u22121rn1$ and $\lambda n1$, respectively, where $(\lambda n1,rn1)$ are eigenpairs of matrix $\u2212\u0131U1THU1$. We emphasize that this matrix is of size *K*_{1} which is more than a factor of *d* smaller than the matrix −*ı*Γ_{1}HΓ_{1}. Moreover, we have established that the discrete spectral measure at the heart of the integral representation for $D*$ is given in terms of the *standard* eigenvalues and eigenvectors of $\u2212\u0131U1THU1$, which is a less costly numerical computation than the computation associated with the *generalized* eigenvalue problem^{47} in (31). Consequently, in the process of establishing the equivalence of the effective parameter problems discussed in Sec. IV and Appendix E for the rank-deficient setting, we have also developed a hybrid numerical algorithm that is more efficient at computing the spectral representation of $D*$ than both of the other approaches.