Spectral analysis and computation for homogenization of advection diffusion processes in steady flows

Murphy, N. B.; Cherkaev, E.; Zhu, J.; Xin, J.; Golden, K. M.

doi:10.1063/1.5127457

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⁶³ and Earth’s ecosystems.¹⁵ Advection by geophysical fluids intensifies the dispersion and large scale transport of heat,⁴¹ pollutants,⁵⁵ and nutrients^15,24 diffusing in their environment. Advective processes also enhance the large scale transport of plankton,²⁴ 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⁶⁰ that complex fluid flows transport passive scalars in much the same way as molecular diffusion. The mathematical description of this phenomenon⁶¹ 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.⁴¹

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,”³⁵ 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⁴⁸ 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⁴⁸ 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,⁴⁹ solute transport in porous media,^9,10 anelastic (weakly compressible) flows,³⁶ 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.⁴⁷ 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 (−Δ)⁻¹, which is given in terms of convolution with the Green’s function for Δ = ∇².⁵⁶ 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,⁴⁸ Monte Carlo methods,¹³ and a method accurate for large Péclet numbers.²³ 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.⁴⁵ 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 $\overset{́}{e}$ approximants and the analytic continuation method, in terms of the P $\overset{́}{e}$ 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, x). 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.

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 ∇ ⋅ u = 0 and ⟨u⟩ = 0 is described by the advection diffusion equation

ϕ_{t} (t, x) = u (x) \cdot \nabla ϕ (t, x) + ε Δ ϕ (t, x), ϕ (0, x) = ϕ_{0} (x), x \in R^{d}, t > 0,

(1)

with initial density ϕ₀(x) given. Here, ⟨⋅⟩ denotes volume averaging over the period cell $V$ ⁠, ϕ_t denotes partial differentiation of ϕ with respect to time t, Δ = ∇ ⋅ ∇ = ∇² 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 ξ ⋅ ξ = |ξ|². 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 u is spatially periodic on the region $V \subset R^{d}$ ⁠. Later, we will discuss the case of an array of randomly perturbed, periodic flows.

The long time, large scale dispersion of passive scalars can be described⁶¹ 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 $\bar{ϕ}$ and a (constant) effective diffusivity matrix $D^{*}$ ⁠,³⁵

{\bar{ϕ}}_{t} (t, x) = \nabla \cdot [D^{*} \nabla \bar{ϕ} (t, x)], \bar{ϕ} (0, x) = ϕ_{0} (x) .

(2)

The components $D_{j k}^{*}$ ⁠, j, k = 1, …, d, of $D^{*}$ are given by³⁵

D_{j k}^{*} = ε δ_{j k} + ⟨ u_{j} χ_{k} ⟩ .

(3)

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 jth component of the fluid velocity field u,^18,35

u \cdot \nabla χ_{j} + ε Δ χ_{j} = - u_{j}, ⟨ \nabla χ_{j} ⟩ = 0 .

(4)

Equations (2)–(4) follow from the assumption that the length scale associated with spatial variations of the initial density ϕ₀ 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³⁵

ϕ (0, x) = ϕ_{0} (δ x) .

(5)

Anticipating that ϕ will have diffusive dynamics as t → ∞, space and time are rescaled by x ↦ x/δ and t ↦ t/δ². As δ → 0, the associated solution ϕ^δ(t, x) = ϕ(t/δ², x/δ) of Eq. (1) (in the rescaled variables) converges to $\bar{ϕ} (t, x)$ which satisfies Eq. (2). The convergence is in an L² 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/δ.³⁵ Other space-time scalings have also been considered, which have led to space-time dependent $D^{*}$ ⁴⁸ and even anomalous diffusive dynamics.³³ 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,

ε \mapsto τ ε / ℓ^{2}, u \mapsto τ u / ℓ .

(6)

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 u and ɛ 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₀v, where the parameter u₀ has dimensions of velocity and represents the “flow strength” of u which is independent of the geometry of the flow; the flow geometry is encapsulated in the nondimensional vector field v. With these definitions, we choose reference scales τ and ℓ in Eq. (6) to satisfy u₀ = ℓ/τ so that u ↦ v and ɛ ↦ ɛ/u₀ ℓ. For example, in Sec. V, we compute $D^{*}$ for BC-flow¹¹ having dimensional fluid velocity field u = u₀ (C cos y, B cos x), where the flow strength u₀ ∈ (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₀/ɛ, 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⁶ P(z)/Q(z), where z = Pe², 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 u is given by BC-flow, the moments of the measure depend⁴³ on the parameters B and C. Our numerical investigations have shown if the 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⁴³ 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,⁴⁴ where the parameter u₀ 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₀((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, ∞) is independent of the parameters B, C, γ, β ∈ [0, 1] which determine the geometric and dynamical properties of u. This choice of nondimensionalization gives a clearer interpretation of the advection and diffusion dominated regimes in terms of 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⁴⁸ 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$ ⁠,

H = {f \in L^{2} (V, ν) : f (x) is periodic in V and ⟨ f ⟩ = 0},

(7)

where ν is the Lebesgue measure on $R^{d}$ ⁠, restricted to $V$ ⁠, and the σ-algebra associated with the underlying measure space is generated by the Lebesgue measurable subsets of $R^{d}$ ⁠. The Hilbert space $H$ is equipped with a sesquilinear inner product ⟨⋅, ⋅⟩ defined by $⟨ f, h ⟩ = ⟨ \bar{f} h ⟩$ ⁠, with $⟨ h, f ⟩ = \bar{⟨ f, h ⟩}$ for $f, h \in H$ ⁠, which induces a norm ∥⋅∥ defined by ∥f∥ = ⟨f,f⟩^1/2 and $f \in H$ implies ∥f∥ < ∞.

One could also consider a random fluid flow filling all of $R^{d}$ ⁠, with a velocity field u determined by the probability space (Ω, P) with σ-algebra generated by the sets {u(x) ∈ B}, where $x \in R^{d}$ and B is a Borel subset of $R^{d}$ ⁠.⁴ Here, Ω is the set of all geometric realizations of u, which is indexed by the parameter ω ∈ Ω 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 = L^{2} (Ω, P)$ ⁠, i.e., the space of all P-measurable complex-valued scalar functions ξ satisfying $‖ ξ ‖ = {⟨ | ξ |^{2} ⟩}^{1 / 2} < \infty$ ⁠, where ⟨⋅⟩ denotes ensemble averaging and the underlying sesquilinear inner-product is defined by $⟨ ξ, ζ ⟩ = ⟨ \bar{ξ} ζ ⟩$ ⁠. In this case, one could consider a random fluid flow with a velocity field u that is stationary³⁵ 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.¹⁸ In this case, the σ-algebra associated with the underlying probability space is generated by the Lebesgue measurable subsets of $R^{d}$ ⁠. 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$ ⁠.¹⁸ The effective diffusivity matrix $D^{*}$ is obtained by taking an infinite volume limit, $D^{*} = \lim_{n \to \infty} D_{n}^{*}$ ⁠, of the finite volume effective diffusivity matrix $D_{n}^{*}$ and evoking an ergodic theorem.^18,22 Numerically, by the law of large numbers,¹⁶ 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.⁵⁸ 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 $H^{1, 2}$ ⁠, which itself is a Hilbert space,^9,10

H^{1, 2} = {f \in H : {‖ f ‖}_{1, 2} < \infty, ⟨ f ⟩ = 0}, {‖ f ‖}_{1, 2} = {⟨ | \nabla f |^{2} ⟩}^{1 / 2},

(8)

where ∥⋅∥_1,2 is the norm induced by the underlying sesquilinear inner-product ⟨⋅,⋅⟩_1,2 defined by ⟨f,h⟩_1,2 = ⟨∇f ⋅ ∇h⟩, with ${⟨ h, f ⟩}_{1, 2} = {\bar{⟨ f, h ⟩}}_{1, 2}$ (recall ξ ⋅ ζ = ξ^†ζ includes complex conjugation).

Recall the definition of the components $D_{j k}^{*} = ε δ_{j k} + ⟨ u_{j} χ_{k} ⟩$ ⁠, j, k = 1, …, d, of the effective diffusivity matrix $D^{*}$ in (3). Rewrite the functional ⟨u_jχ_k⟩ as⁴⁸

⟨ u_{j} χ_{k} ⟩ = ⟨ [Δ Δ^{- 1} u_{j}] χ_{k} ⟩ = - ⟨ \nabla Δ^{- 1} u_{j} \cdot \nabla χ_{k} ⟩ = - {⟨ Δ^{- 1} u_{j}, χ_{k} ⟩}_{1, 2},

(9)

where we have integrated by parts and used the periodicity of the functions u_j and χ_k. Here, Δ⁻¹ is based on convolution with the Green’s function for the Laplacian Δ.⁵⁶ The symmetric $S^{*}$ and antisymmetric $A^{*}$ parts of the effective diffusivity matrix $D^{*}$ are defined by

D^{*} = S^{*} + A^{*}, S^{*} = \frac{1}{2} (D^{*} + {[D^{*}]}^{T}), A^{*} = \frac{1}{2} (D^{*} - {[D^{*}]}^{T}) .

(10)

Substituting into Eq. (9) the expression − u_j = u ⋅ ∇χ_j + ɛΔχ_j for −u_j in (4) yields the following functional formulas for the components $S_{j k}^{*}$ and $A_{j k}^{*}$ ⁠, j, k = 1, …, d, of $S^{*}$ and $A^{*}$ ⁠,

S_{j k}^{*} = ε (δ_{j k} + {⟨ χ_{j}, χ_{k} ⟩}_{1, 2}), A_{j k}^{*} = {⟨ A χ_{j}, χ_{k} ⟩}_{1, 2}, A = Δ^{- 1} [u \cdot \nabla] .

(11)

Due to the incompressibility of the fluid velocity field, ∇ ⋅ u = 0, the operator A is antisymmetric on $H^{1, 2}$ ⁠, i.e., ⟨Aξ, ζ⟩_1,2 = −⟨ξ, Aζ⟩_1,2 for all $ξ, ζ {\in H}^{1, 2}$ [see Eq. (C1)]. Since the scalar fields χ_j and Aχ_j are real-valued, we have that ${⟨ χ_{j}, χ_{k} ⟩}_{1, 2} = {⟨ χ_{k}, χ_{j} ⟩}_{1, 2}$ and $A_{k j}^{*} = {⟨ A χ_{k}, χ_{j} ⟩}_{1, 2} = - {⟨ χ_{k}, A χ_{j} ⟩}_{1, 2} = - {⟨ A χ_{j}, χ_{k} ⟩}_{1, 2} = - A_{j k}^{*}$ ⁠, confirming that $S^{*}$ is symmetric and $A^{*}$ is antisymmetric, with $A_{k k}^{*} = 0$ ⁠.

Applying the operator Δ⁻¹ 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):

χ_{j} = {(ε + A)}^{- 1} g_{j}, g_{j} = {(- Δ)}^{- 1} u_{j} .

(12)

Substituting the resolvent formula for χ_j in (12) into the bilinear functionals in Eq. (11) yields

\begin{align} S_{j k}^{*} & = ε (δ_{j k} + {⟨ {(ε + A)}^{- 1} g_{j}, {(ε + A)}^{- 1} g_{k} ⟩}_{1, 2}), \\ A_{j k}^{*} & = {⟨ A {(ε + A)}^{- 1} g_{j}, {(ε + A)}^{- 1} g_{k} ⟩}_{1, 2} . \end{align}

(13)

Since $V$ is a bounded domain, the linear operator Δ⁻¹ is bounded on $H$ with bounded operator norm ∥Δ⁻¹∥ < ∞.⁵⁶ When |u| is uniformly bounded on the period cell $V$ ⁠, $\sup_{x \in V} | u (x) | < \infty$ ⁠, the linear operator A is bounded on $H^{1, 2}$ ⁠, with [see Eqs. (C2) and (C3)]

{‖ A ‖}_{1, 2} \leq [) ‖ Δ^{- 1} ‖ \sup_{x \in V} | u (x) |^{2} {(]}^{1 / 2} < \infty,

(14)

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 $u_{k} \in H$ ⁠, k = 1, …, d, the operator A is compact on $H^{1, 2}$ ⁠,^9,10,48 hence bounded.⁵⁶ Consequently, M = −ıA, where $ı = \sqrt{- 1}$ ⁠, is a bounded linear operator on $H^{1, 2}$ with ∥M∥_1,2 = ∥A∥_1,2 < ∞. Moreover, the skew-symmetry of A and the sesquilinearity of the $H^{1, 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 $H^{1, 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,⁵³ i.e.,

Σ \subseteq [- {‖ A ‖}_{1, 2}, {‖ A ‖}_{1, 2}] .

(15)

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⁵⁸

\lim_{λ \to \inf Σ} Q (λ) = 0, \lim_{λ \to \sup Σ} Q (λ) = I .

(16)

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

The spectral theorem also states, for all complex-valued functions f ∈ L²(μ_ξξ) and h ∈ L²(μ_ζζ), 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.⁵⁸ In particular, this functional has the following Radon–Stieltjes integral representation involving the Stieltjes measure μ_ξζ, for all $ξ, ζ {\in H}^{1, 2}$ ⁠:

{⟨ f (M) ξ, h (M) ζ ⟩}_{1, 2} = \int_{- \infty}^{\infty} \bar{f} (λ) h (λ) d μ_{ξ ζ} (λ), μ_{ξ ζ} (λ) = {⟨ Q (λ) ξ, ζ ⟩}_{1, 2},

(17)

where the integration is over the spectrum Σ of M^53,58 and $\bar{f}$ denotes the complex conjugation of the scalar function f. Since the Stieltjes measure μ_ξζ has the property⁵⁸ $\int_{a}^{b} d μ_{ξ ζ} (λ) = μ_{ξ ζ} (b) - μ_{ξ ζ} (a)$ ⁠, Eq. (16) implies that the mass $μ_{ξ ζ}^{0}$ of the measure μ_ξζ is given by

μ_{ξ ζ}^{0} = \int_{- \infty}^{\infty} d μ_{ξ ζ} (λ) = \int_{- \infty}^{\infty} d {⟨ Q (λ) ξ, ζ ⟩}_{1, 2} = {⟨ ξ, ζ ⟩}_{1, 2},

(18)

which is bounded in the sense that $| μ_{ξ ζ}^{0} | \leq {‖ ξ ‖}_{1, 2} {‖ ζ ‖}_{1, 2} < \infty$ for all $ξ, ζ {\in H}^{1, 2}$ ⁠. Due to the sesquilinearity of the inner-product and the fact that the projection operator Q(λ) is self-adjoint on $H^{1, 2}$ ⁠, the complex-valued function μ_ζξ(λ) satisfies $μ_{ζ ξ} (λ) = {\bar{μ}}_{ξ ζ} (λ)$ and μ_ξξ(λ) = ∥Q(λ)ξ∥², 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 $μ_{g_{j} g_{k}} (λ)$ ⁠, where g_j = (−Δ)⁻¹u_j is defined in (12). Denote the real and imaginary parts of the function μ_jk(λ) by Re μ_jk(λ) and Im μ_jk(λ), respectively.

Theorem 1.

Let Q(λ) denote the resolution of the identity corresponding to the self-adjoint operator M. Then, the components

S_{j k}^{*}

and

A_{j k}^{*}

, 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:

S_{j k}^{*} = ε (δ_{j k} + \int_{- \infty}^{\infty} \frac{d R e μ_{j k} (λ)}{ε^{2} + λ^{2}}), A_{j k}^{*} = \int_{- \infty}^{\infty} \frac{λ d I m μ_{j k} (λ)}{ε^{2} + λ^{2}} .

(19)

For j ≠ k, Re μ_jk and Im μ_jk are signed measures.²⁰ For j = k, Re μ_kk = μ_kk is a positive measure²⁰ , 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 ɛ,

S_{k k}^{*} = ε (1 + \int_{- \infty}^{\infty} \frac{d μ_{k k} (λ)}{ε^{2} + λ^{2}}) \geq ε .

(20)

The mass

μ_{j k}^{0}

of the measure μ_jk is real-valued and satisfies

μ_{j k}^{0} = ⟨ [{(- Δ)}^{- 1} u_{j}] u_{k} ⟩, | μ_{j k}^{0} | \leq ‖ Δ^{- 1} u_{j} ‖ ‖ u_{k} ‖ < \infty .

(21)

Proof of Theorem 1.

We first provide integral representations for the bilinear functional formulas in Eq. (13) for $S_{j k}^{*}$ and $A_{j k}^{*}$ ⁠. 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 $ξ, ζ {\in H}^{1, 2}$ in (17). Using A = ıM, comparison of the functionals in Eqs. (13) and (17) identifies f(λ) = h(λ) = (ɛ + ıλ)⁻¹ for the first formula in (13), while f(λ) = ıλ(ɛ + ıλ)⁻¹ and h(λ) = (ɛ + ıλ)⁻¹ for the second formula, with ξ = g_j and ζ = g_k for both of these formulas.

Now we just need to verify f, h ∈ L²(μ_kk) for all k = 1, …, d and 0 < ɛ < ∞. From Eq. (), the mass

μ_{j k}^{0}

of the measure μ_jk is given by

μ_{j k}^{0} = {⟨ g_{j}, g_{k} ⟩}_{1, 2} = ⟨ \nabla Δ^{- 1} u_{j} \cdot \nabla Δ^{- 1} u_{k} ⟩

⁠. Integration by parts and the Cauchy-Schwartz inequality²⁰ then yield Eq. (). In particular,

| μ_{k k}^{0} | \leq ‖ Δ^{- 1} ‖ {‖ u_{k} ‖}^{2} < \infty

⁠, as Δ⁻¹ has bounded operator norm ∥Δ⁻¹∥ on

L^{2} (V)

⁠.⁵⁶ Consequently, since

0 < {(ε^{2} + λ^{2})}^{- 1} \leq 1 / ε^{2} < \infty

and

0 < λ^{2} {(ε^{2} + λ^{2})}^{- 1} < 1

for all 0 < ɛ < ∞, it is clear that f, h ∈ L²(μ_kk). Consequently, the spectral theorem in () implies that the functional formulas for

S_{j k}^{*}

and

A_{j k}^{*}

in Eq. () have the following Stieltjes integral representations:

S_{j k}^{*} = ε (δ_{j k} + \int_{- \infty}^{\infty} \frac{d μ_{j k} (λ)}{ε^{2} + λ^{2}}), A_{j k}^{*} = - \int_{- \infty}^{\infty} \frac{ı λ d μ_{j k} (λ)}{ε^{2} + λ^{2}},

(22)

which involve the complex measure μ_jk.

We now show how the integrals for

S_{j k}^{*}

and

A_{j k}^{*}

in () can be represented in terms of the signed measures Re μ_jk and Im μ_jk. Since ɛ, χ_k, u_k, and g_k in () are real-valued, we have from () the following symmetry conditions:

\begin{align} {⟨ {(ε I + A)}^{- 1} g_{j}, {(ε I + A)}^{- 1} g_{k} ⟩}_{1, 2} & = {⟨ {(ε I + A)}^{- 1} g_{k}, {(ε I + A)}^{- 1} g_{j} ⟩}_{1, 2}, \\ {⟨ A {(ε I + A)}^{- 1} g_{j}, {(ε I + A)}^{- 1} g_{k} ⟩}_{1, 2} & = {⟨ {(ε I + A)}^{- 1} g_{k}, A {(ε I + A)}^{- 1} g_{j} ⟩}_{1, 2} . \end{align}

(23)

These symmetries, the sesquilinearity of the

H

-inner-product, the linearity⁵⁸ of the Stieltjes integrals in Eq. () with respect to the function μ_jk(λ), and the two identities

R e μ_{j k} (λ) = (μ_{j k} (λ) + {\bar{μ}}_{j k} (λ)) / 2

and

Im μ_{j k} (λ) = (μ_{j k} (λ) - {\bar{μ}}_{j k} (λ)) / (2 ı)

yield Eq. (). This concludes our Proof of Theorem 1.

Theorem 2.

Let Σ denote the spectrum of the operator M, denote λ₊ = sup_λ∈Σ|λ|, and assume 0 < ɛ < ∞. Then, the diagonal components

S_{k k}^{*} = D_{k k}^{*}

of the effective diffusivity matrix satisfy the following, rigorous upper⁴ and lower bounds:

ε [1 + μ_{k k}^{0} / (ε^{2} + λ_{+}^{2})] \leq S_{k k}^{*} \leq ε [1 + μ_{k k}^{0} / ε^{2}] .

(24)

For j ≠ k, let

R e μ_{j k} = R e μ_{j k}^{+} - R e μ_{j k}^{-}

and

Im μ_{j k} = Im μ_{j k}^{+} - Im μ_{j k}^{-}

denote the Jordan decomposition of the signed measures Re μ_jk and Im μ_jk, respectively, and denote the total variation of these measures by

| R e μ |_{j k} = R e μ_{j k}^{+} + R e μ_{j k}^{-}

and

| Im μ |_{j k} = Im μ_{j k}^{+} + Im μ_{j k}^{-}

, respectively. Finally, denote the masses of the measures

R e μ_{j k}^{+}

,

R e μ_{j k}^{-}

, and |Im μ|_jk by

{[R e μ_{j k}^{+}]}^{0}

,

{[R e μ_{j k}^{-}]}^{0}

, and

| Im μ |_{j k}^{0}

, respectively. Then, the off-diagonal components

S_{j k}^{*}

and

A_{j k}^{*}

, j ≠ k, satisfy the following, rigorous upper and lower bounds:

ε \frac{{[R e μ_{j k}^{+}]}^{0}}{ε^{2} + λ_{+}^{2}} - \frac{{[R e μ_{j k}^{-}]}^{0}}{ε} \leq S_{j k}^{*} \leq \frac{{[R e μ_{j k}^{+}]}^{0}}{ε} - ε \frac{{[R e μ_{j k}^{-}]}^{0}}{ε^{2} + λ_{+}^{2}}, j \neq k,

(25)

- \frac{λ_{+} | Im μ |_{j k}^{0}}{ε^{2}} \leq A_{j k}^{*} \leq \frac{λ_{+} | Im μ |_{j k}^{0}}{ε^{2}}, j \neq k .

(26)

Proof of Theorem 2.

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⁵⁶ away from the spectral origin λ = 0 and comes in ± pairs⁴⁸ with an accumulation point at λ = 0.⁵⁶ Denote λ₊ = sup_λ∈Σ|λ| and note that inf_λ∈Σλ² = 0. The inequalities $1 / (ε^{2} + λ_{+}^{2}) \leq 1 / (ε^{2} + λ^{2}) \leq 1 / ε^{2}$ hold for all λ ∈ Σ. Consequently, since μ_kk is a positive measure with finite mass $μ_{k k}^{0}$ ⁠, the inequalities in (24) hold.²⁰ It may happen that $μ_{k k}^{0} = 0$ ⁠, hence $S_{k k}^{*} = ε$ ⁠, e.g., shear flow orthogonal to the kth direction.^4,17

When j ≠ k, Re μ_jk is a signed measure. By the Jordan decomposition of Re μ_jk, there are unique, positive measures

R e μ_{j k}^{+}

and

R e μ_{j k}^{-}

such that

R e μ_{j k} = R e μ_{j k}^{+} - R e μ_{j k}^{-}

⁠. Moreover, associated with the signed measure, Re μ_jk is its total variation |Re μ|_jk,²⁰

R e μ_{j k} = R e μ_{j k}^{+} - R e μ_{j k}^{-}, | R e μ |_{j k} = R e μ_{j k}^{+} + R e μ_{j k}^{-} .

(27)

From Eq. (), the measures

Re μ_{j k}^{+}

and

R e μ_{j k}^{-}

have bounded mass, which we denote

{[Re μ_{j k}^{+}]}^{0}

and

{[Re μ_{j k}^{-}]}^{0}

⁠, respectively; thus, the mass

| R e μ |_{j k}^{0}

of the measure |Re μ|_jk is also bounded. Since we have²⁰ that

| S_{j k}^{*} | \leq \int d | Re μ |_{j k} (λ) / (ε^{2} + λ^{2})

⁠, the upper bound in Eq. () with

μ_{k k}^{0}

replaced by

| Re μ |_{j k}^{0}

holds for the positive quantity

| S_{j k}^{*} |

⁠. These bounds for

S_{j k}^{*}

can be improved upon by separately considering the positive and negative contributions of the integral representation for

S_{j k}^{*}

⁠, yielding Eq. (). In a similar way, we obtain the bounds for

A_{j k}^{*}

in Eq. (). This concludes our Proof of Theorem 2.

We conclude this section by noting that bounds on $S_{k k}^{*}$ 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 = Δ⁻¹[u ⋅ ∇], and g_j = −Δ⁻¹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 (−Δ)⁻¹[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(x) is incompressible, ∇ ⋅ u = 0, there is a real (nondimensional) antisymmetric matrix H(x)^3,4 such that

u = \nabla \cdot H, H^{T} = - H,

(28)

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

\nabla \cdot [H \nabla] = [\nabla \cdot H] \cdot \nabla,

(29)

as operators acting on such functions. With (29) we can write the operator M = ı(−Δ)⁻¹[u ⋅ ∇] as M = ı(−Δ)⁻¹∇ ⋅ [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,⁴² $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(x) are mapped to vectors u_j and χ_j with L^d components, respectively, and the matrix H(x) is mapped to a square banded antisymmetric matrix of size N = L^dd (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 $\nabla^{T} = (\nabla_{1}^{T}, \dots, \nabla_{d}^{T})$ 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 = ı(−Δ)⁻¹∇ ⋅ [H∇] as $M = {(\nabla^{T} \nabla)}^{- 1} [- ı \nabla^{T} H \nabla]$ ⁠. This composition of the Hermitian matrix −ı∇^TH∇ and the real-symmetric matrix ${(\nabla^{T} \nabla)}^{- 1}$ 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 $H^{1, 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 ${⟨ M φ_{n}, φ_{n} ⟩}_{1, 2} = λ_{n}$ ⁠. Namely, the operator M is self-adjoint with respect to the inner-product ⟨⋅, ⋅⟩_1,2, its eigenfunctions $φ_{n} {\in H}^{1, 2}$ are orthonormal ${⟨ φ_{n}, φ_{m} ⟩}_{1, 2} = δ_{n m}$ ⁠, n, m = 1, 2, 3, …, with respect to the inner-product ⟨⋅, ⋅⟩_1,2, and the spectrum Σ of M is real valued, $Σ \subset R$ ⁠. Toward this goal, consider the eigenvalue problem Mφ_n = λ_nφ_n in the following “strong” form:

ı \nabla \cdot [H \nabla φ_{n}] = λ_{n} (- Δ) φ_{n} .

(30)

Using a discrete version of Eq. (30), we will establish the integral representations in (19) for discrete versions of the functionals $S_{j k}^{*} = ε (δ_{j k} + {⟨ χ_{j}, χ_{k} ⟩}_{1, 2})$ and $A_{j k}^{*} = {⟨ ı M χ_{j}, χ_{k} ⟩}_{1, 2}$ in (11), involving a discrete spectral measure.

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

B z_{n} = λ_{n} C z_{n}, B = - ı \nabla^{T} H \nabla, C = \nabla^{T} \nabla .

(31)

The first formula in Eq. (31) is a generalized eigenvalue problem⁴⁷ 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.⁴⁷ In particular, the generalized eigenvalues λ_n are all real, the generalized eigenvectors z_n form a basis for $C^{K}$ ⁠, and the z_n are orthonormal in the following sense: $z_{n}^{†} C z_{m} = δ_{n m}$ ⁠, n, m = 1, …, K.⁴⁷ Since C = ∇^T∇ is real-valued, this is equivalent to the Sobolev-type orthogonality condition,

\nabla z_{n} \cdot \nabla z_{m} = δ_{n m} .

(32)

In other words, the generalized eigenvectors z_n are orthonormal with respect to the “discrete inner-product” ⟨⋅, ⋅⟩_1,2 defined by ⟨ξ,ζ⟩_1,2 = ⟨∇ξ ⋅ ∇ζ⟩, for $ξ, ζ \in C^{K}$ 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⁴⁷

Z^{†} B Z = Λ, Z^{†} C Z = I .

(33)

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 $C^{K}$ and satisfy the orthogonality relation in (32), for all $ξ \in C^{K}$ ⁠, we have $ξ = \sum_{n} (\nabla z_{n} \cdot \nabla ξ) z_{n} = \sum_{n} (z_{n} {[\nabla z_{n}]}^{†} \nabla) ξ$ ⁠, which implies

\sum_{n = 1}^{K} Q_{n} = I, Q_{n} = z_{n} {[\nabla z_{n}]}^{†} \nabla, Q_{l} Q_{m} = Q_{l} δ_{l m},

(34)

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., ${⟨ Q_{n} ξ, ζ ⟩}_{1, 2} = {⟨ ξ, Q_{n} ζ ⟩}_{1, 2}$ for all $ξ, ζ \in C^{K}$ ⁠. From Bz_n = λ_nCz_n and Eq. (34), we have that BQ_n = λ_nCQ_n which is equivalent to C⁻¹BQ_n = λ_nQ_n, since the matrix C is invertible. This formula and (34) then imply that the matrix C⁻¹B has the spectral decomposition C⁻¹B = ∑_nλ_nQ_n. By the mutual orthogonality of the projection matrices Q_n and by induction, we have that ${[C^{- 1} B]}^{m} = \sum_{n} λ_{n}^{m} Q_{n}$ for all $m \in N$ ⁠. This, in turn, implies that f(C⁻¹B) = ∑_nf(λ_n)Q_n for any polynomial $f : R \mapsto C$ ⁠.

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 $ξ, ζ \in C^{K}$ and all complex-valued polynomials f(λ) and h(λ), the bilinear functional ${⟨ f (C^{- 1} B) ξ, h (C^{- 1} B) ζ ⟩}_{1, 2}$ has the integral representation in (17), with M substituted by C⁻¹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

Q (λ) = \sum_{n : λ_{n} \leq λ} θ (λ - λ_{n}) Q_{n}, d μ_{ξ ζ} (λ) = \sum_{n : λ_{n} \leq λ} ⟨ δ_{λ_{n}} (d λ) [\nabla Q_{n} ξ \cdot \nabla ζ] ⟩ .

(35)

Here, θ(λ) is the Heaviside function, satisfying θ(λ) = 0 for λ < 0 and θ(λ) = 1 for λ ≥ 0, and $δ_{λ_{n}} (d λ) = d θ (λ - λ_{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 $μ_{ξ ζ}^{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.

Theorem 3.

Consider the generalized eigenvalue problem Bz_n = λ_nCz_n in (). 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. (). The discrete, matrix representations of the bilinear functional formulas for

S_{j k}^{*}

and

A_{j k}^{*}

in Eq. () are given by

S_{j k}^{*} = ε (δ_{j k} + ⟨ \nabla χ_{j} \cdot \nabla χ_{k} ⟩), A_{j k}^{*} = ⟨ \nabla C^{- 1} [ı B] χ_{j} \cdot \nabla χ_{k} ⟩ .

(36)

Also, the discrete representation of the resolvent formula for χ_j in Eq. () is given by

χ_{j} = Z {(ε I + ı Λ)}^{- 1} Z^{†} u_{j} .

(37)

The discrete representations of the bilinear functional formulas for

S_{j k}^{*}

and

A_{j k}^{*}

in () are given by

\begin{align} S_{j k}^{*} & = ε (δ_{j k} + ⟨ {(ε I + ı Λ)}^{- 1} Z^{†} u_{j} \cdot {(ε I + ı Λ)}^{- 1} Z^{†} u_{k} ⟩), \\ A_{j k}^{*} & = ⟨ ı Λ {(ε I + ı Λ)}^{- 1} Z^{†} u_{j} \cdot {(ε I + ı Λ)}^{- 1} Z^{†} u_{k} ⟩ . \end{align}

(38)

Consequently, from the discrete analog of Eq. () and the formulas for

S_{j k}^{*}

and

A_{j k}^{*}

in (), we have the following series representations:

S_{j k}^{*} / ε - δ_{j k} = \sum_{n = 1}^{K} \frac{R e [\bar{(z_{n}^{†} u_{j})} (z_{n}^{†} u_{k})]}{ε^{2} + λ_{n}^{2}}, A_{j k}^{*} = \sum_{n = 1}^{K} \frac{λ_{n} Im [\bar{(z_{n}^{†} u_{j})} (z_{n}^{†} u_{k})]}{ε^{2} + λ_{n}^{2}} .

(39)

Finally, defining

g_{j} = {(\nabla^{T} \nabla)}^{- 1} u_{j}

and recalling the projection matrix Q_n in (), we have

\bar{(z_{n}^{†} u_{j})} (z_{n}^{†} u_{k}) = \nabla Q_{n} g_{j} \cdot \nabla g_{k} .

(40)

It follows from Eqs. () and () that the series representations for

S_{j k}^{*}

and

A_{j k}^{*}

in () have the Stieltjes integral representations in Eq. (), involving the discrete spectral measure μ_jk in () with

ξ = g_{j}

and

ζ = g_{k}

.

Proof of Theorem 3.

From the matrix representation

A = {(\nabla^{T} \nabla)}^{- 1} (\nabla^{T} H \nabla)

of the operator A = Δ⁻¹∇ ⋅ [H∇] and Eq. (), the matrix representation of the functional formulas

S_{j k}^{*} = ε (δ_{j k} + {⟨ χ_{j}, χ_{k} ⟩}_{1, 2})

and

A_{j k}^{*} = {⟨ A χ_{j}, χ_{k} ⟩}_{1, 2}

in Eq. () are given by (). Moreover, the matrix representation of the cell problem ɛΔχ_j + ∇ ⋅ [H∇]χ_j = −u_j in () is given by

(ε C + ı B) χ_{j} = u_{j} .

(41)

The matrix Z is invertible, as the generalized eigenvectors z_n form a basis for

C^{K}

⁠. Consequently, by Eq. (), we have that B = Z^−†ΛZ⁻¹, C = Z^−†Z⁻¹. It now follows from Eq. () that Z^−†(ɛI + ıΛ)Z⁻¹χ_j = u_j, which implies the resolvent formula for χ_j in Eq. ().

Substituting the resolvent formula for χ_j in () into Eq. () yields Eq. (), where we have used that [∇Z]^† = [∇Z]⁻¹. The quadratic form Z^†u_j ⋅ Z^†u_k arising in () can be written in terms of the projection matrices Q_n defined in () as follows:

Z^{†} u_{j} \cdot Z^{†} u_{k} = \sum_{n = 1}^{K} \bar{(z_{n}^{†} u_{j})} (z_{n}^{†} u_{k}) = \sum_{n = 1}^{K} z_{n} z_{n}^{†} u_{j} \cdot u_{k} .

(42)

We now demonstrate that

z_{n} z_{n}^{†} u_{j} \cdot u_{k} = \nabla Q_{n} g_{j} \cdot \nabla g_{k}

⁠, where

g_{j} = {(\nabla^{T} \nabla)}^{- 1} u_{j}

⁠,

z_{n} z_{n}^{†} u_{j} \cdot u_{k} = z_{n} z_{n}^{†} [\nabla^{T} \nabla] g_{j} \cdot [\nabla^{T} \nabla] g_{k} = [\nabla z_{n}] {[\nabla z_{n}]}^{†} \nabla g_{j} \cdot \nabla g_{k} = \nabla Q_{n} g_{j} \cdot \nabla g_{k},

(43)

which establishes Eq. (). This concludes our Proof of Theorem 3.

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 $g_{k}$ are real-valued and the two identities $R e μ_{j k} (λ) = (μ_{j k} (λ) + {\bar{μ}}_{j k} (λ)) / 2$ and $Im μ_{j k} (λ) = (μ_{j k} (λ) - {\bar{μ}}_{j k} (λ)) / (2 ı)$ ⁠, we have

\begin{align} Re μ_{j k} (λ) & = \frac{1}{2} \sum_{n : λ_{n} \leq λ} ⟨ θ (λ - λ_{n}) [\nabla (Q_{n} + {\bar{Q}}_{n}) g_{j} \cdot \nabla g_{k}] ⟩, \\ Im μ_{j k} (λ) & = \frac{1}{2 ı} \sum_{n : λ_{n} \leq λ} ⟨ θ (λ - λ_{n}) [\nabla (Q_{n} - {\bar{Q}}_{n}) g_{j} \cdot \nabla g_{k}] ⟩, \end{align}

(44)

with $[\nabla (Q_{n} + {\bar{Q}}_{n}) g_{j} \cdot \nabla g_{k}] = 2 Re [\nabla Q_{n} g_{j} \cdot \nabla g_{k}]$ and $[\nabla (Q_{n} - {\bar{Q}}_{n}) g_{j} \cdot \nabla g_{k}] = 2 ı Im [\nabla Q_{n} g_{j} \cdot \nabla g_{k}]$ ⁠.

Consider the generalized eigenvalue problem in Eq. (31) written as [ıB]z_n = ıλ_nCz_n, with B = −ı∇^TH∇ and C = ∇^T∇. Since the matrices ıB and C are real-valued, we have $[ı B] {\bar{z}}_{n} = - ı λ_{n} C {\bar{z}}_{n}$ ⁠. 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 $I_{K}$ as $I_{K} = {- K / 2, \dots, - 1, 1, \dots, K / 2}$ such that λ_−n = −λ_n and $z_{- n} = \bar{z_{n}}$ ⁠. If K is odd, then λ₀ = 0 is also an eigenvalue with a real-valued eigenvector z₀. Consequently, the representations of the measures Reμ_jk and Imμ_jk, following from the functions in Eq. (44), can be simplified⁴⁸ 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 $λ_{- n}^{2} = {(- λ_{n})}^{2} = λ_{n}^{2}$ and $z_{- n} = \bar{z_{n}}$ ⁠, thus $Q_{- n} = {\bar{Q}}_{n}$ ⁠. Consequently,

\begin{align} Re [\nabla Q_{n} g_{j} \cdot \nabla g_{k}] + R e [\nabla Q_{- n} g_{j} \cdot \nabla g_{k}] = 2 R e [\nabla Q_{n} g_{j} \cdot \nabla g_{k}], \\ λ_{n} Im [\nabla Q_{n} g_{j} \cdot \nabla g_{k}] + λ_{- n} Im [\nabla Q_{- n} g_{j} \cdot \nabla g_{k}] = 2 λ_{n} Im [\nabla Q_{n} g_{j} \cdot \nabla g_{k}], \end{align}

(45)

with $λ_{0} Im [Q_{0} g_{j} \cdot g_{k}] \equiv 0$ ⁠. 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,⁴ and a cell flow with closed streamlines, for which it is known^17,18 that $D^{*} \sim ε^{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 $S_{j k}^{*}$ ⁠, 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(x) is given in terms of an antisymmetric matrix H = H(x), u = ∇ ⋅ H. For dimension d = 2, the matrix H is determined by a stream function Ψ = Ψ(x),

H = [\begin{matrix} 0 & Ψ \\ - Ψ & 0 \end{matrix}],

(46)

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¹¹ and “cat’s eye” flow,¹⁷ which are defined by the following stream functions Ψ_BC and Ψ_CE, respectively:

Ψ_{B C} (x) = B \sin x - C \sin y, Ψ_{C E} (x) = \sin x \sin y + α \cos x \cos y,

(47)

where we have denoted x = (x, y). The corresponding fluid velocity fields are

\begin{align} u_{B C} (x) & = (C \cos y, B \cos x), \\ u_{C E} (x) & = (- \sin x \cos y + α \cos x \sin y, \cos x \sin y - α \sin x \cos y) . \end{align}

(48)

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 u ⋅ $\nabla$ Ψ = 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 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 $S_{j k}^{*}$ ⁠, j, k = 1, 2, exhibit these symmetries. For BC-flow, these symmetries allow us to restrict our attention to the behavior of μ_jk and $S_{j k}^{*}$ 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.

FIG. 1.

View large Download slide

Transitions in the flow structure. The streamlines for BC-flow and cat’s eye flow are displayed for various parameter values, transitioning from the shear to cell flow structure. From left to right and top to bottom, the parameter values associated with the BC-flow are B = 1 fixed and C = 0, 0.01, 0.1, 0.3, 0.5, and 1, while those for cat’s eye flow are α = 0, 0.1, 0.3, 0.5, 0.7, and 1.

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 = L^{2} (ν \times P)$ ⁠, where P is the probability measure associated with the random variable α.²⁹ In this case, by Fubini’s theorem,²⁰ ⟨⋅⟩ 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 π]}^{d}$ ⁠, for example, is replaced by a square lattice $V_{L}^{d}$ 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 $V_{L}^{d}$ must be bijectively mapped to a one dimensional lattice $V_{N}$ of size N = L^dd. The specific structure of $V_{N}$ 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₁(x), …, v_d(x)), say, is bijectively mapped by Θ to a discretized constant vector (v₁, …, 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(x) for $x {\in V}_{L}^{d}$ ⁠. Similarly, the d-dimensional standard basis vector e₁ = (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 ${\hat{e}}_{j}$ ⁠, j = 1, …, d, satisfying

{\hat{e}}_{j} = Θ (e_{j}) / L^{d / 2}, {\hat{e}}_{j} \cdot {\hat{e}}_{k} = δ_{j k},

(49)

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 dx, we have Δx = (2π/L)^d/(2π)^d when $V = {[0, 2 π]}^{d}$ and the spatial average ${⟨ ξ (x) ζ (x) ⟩}_{V}$ becomes ξ ⋅ ζ/L^d. In a similar way, the 2 × 2 matrix H(x) in (46) becomes a N × N antisymmetric banded matrix, where the stream function Ψ(x) is represented by a diagonal 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(σ₁, …, σ_K), where 0 ≤ σ₁ ≤ ⋯ ≤ σ_K. The matrices U and V are of size N × K and K × K, respectively, and satisfy U^TU = I and V^TV = VV^T = I, where I is the identity matrix of size K.¹⁴ When ∇ is of full rank, then the singular values satisfy σ_j > 0 for all j = 1, …, K and the matrix Laplacian ∇^T∇ = VΣ²V^T is invertible. When ∇ is rank-deficient, then there are K₁ nonzero singular values and K₀ = K − K₁ zero singular values, say. For example, when d = 2, the nullity of ∇ is 1 so K₁ = K − 1. In this case, we write U = [U₀ U₁], Σ = diag(O, Σ₁), and V = [V₀ V₁], where O is a matrix of zeros so that $\nabla = U_{1} Σ_{1} V_{1}^{T}$ and the negative matrix Laplacian $\nabla^{T} \nabla = V_{1} Σ_{1}^{2} V_{1}^{T}$ is noninvertible, since $V_{1} V_{1}^{T} \neq I$ ⁠.

The associated matrix analysis in Appendix G demonstrates that the spectral measure μ_jk underlying the Stieltjes integral representation of $S_{j k}^{*}$ is given by

d μ_{j k} (λ) = \sum_{n : λ_{n}^{1} \leq λ} ⟨ m_{j k} (n) δ_{λ_{n}^{1}} (d λ) ⟩,

(50)

where $λ_{n}^{1}$ ⁠, n = 1, …, K₁, are the eigenvalues of the matrix $- ı U_{1}^{T} H U_{1}$ ⁠, 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

m_{j k} (n) = Re [\bar{({[r_{n}^{1}]}^{†} U_{1}^{T} H {\hat{e}}_{j})} ({[r_{n}^{1}]}^{†} U_{1}^{T} H {\hat{e}}_{k})], n = 1, \dots, K_{1},

(51)

which follows from Eqs. (43), (44), (G6), and (G7). Here, $r_{n}^{1}$ ⁠, n = 1, …, K₁, are the complex eigenvectors of the matrix $- ı U_{1}^{T} H U_{1}$ ⁠. 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 μ₁₂ and $S_{12}^{*}$ in our numerical computations discussed in Sec. V B, we denote the spectral weights m_jk(n) associated with the Jordan decomposition $μ_{j k} = μ_{j k}^{+} - μ_{j k}^{-}$ in (27) by $m_{j k}^{+} (n)$ and $m_{j k}^{-} (n)$ ⁠, where $m_{j k}^{\pm} (n) \geq 0$ ⁠. Also, we define the functions ${[S_{12}^{*}]}^{+}$ and ${[S_{12}^{*}]}^{-}$ ⁠,

{[S_{12}^{*}]}^{+} (ε) = \max {S_{12}^{*} (ε), 0}, {[S_{12}^{*}]}^{-} (ε) = \max {- S_{12}^{*} (ε), 0},

(52)

for each 0 < ɛ < ∞, so that $S_{12}^{*} = {[S_{12}^{*}]}^{+} - {[S_{12}^{*}]}^{-}$ ⁠, ${[S_{12}^{*}]}^{\pm} (ε) = S_{12}^{*} (ε; μ_{12}^{\pm})$ ⁠, and ${[S_{12}^{*}]}^{\pm} \geq 0$ ⁠.

In the case of a nonrandom fluid velocity field u, we used L = 200 so that K₁ = 39 999. The eigenvalues $λ_{n}^{1}$ and eigenvectors $r_{n}^{1}$ of the nonrandom Hermitian matrix $- ı U_{1}^{T} H U_{1}$ 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₁ = 9999. In this case, the averaging ⟨⋅⟩ in (50) is interpreted as spatial averaging followed by ensemble averaging over ∼10³ statistical trials.

The numerical accuracy of the eigenvalue problem is determined by the eigenvalue condition numbers $K (λ_{n}^{1})$ ⁠, n = 1, …, K₁, 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 $- ı U_{1}^{T} H U_{1}$ is extremely well conditioned with $\max_{n} | 1 - K (λ_{n}^{1}) | \sim 1 0^{- 14}$ 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 $Γ = \nabla {(\nabla^{T} \nabla)}^{- 1} \nabla^{T}$ 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 $\max_{l m} | {[\nabla ((\nabla^{T} \nabla) \ \nabla^{T}) - U U^{T}]}_{l m} |$ ⁠. When L = 100 and L = 200, this difference is ∼10⁻¹⁵, 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 $S_{j k}^{*}$ ⁠. 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 $S_{k k}^{*}$ ⁠, for some k = 1, …, d, and write the formula in (20) as

S_{k k}^{*} (ε) = ε + E_{k k} (ε), E_{k k} (ε) = ε \int_{- \infty}^{\infty} \frac{d μ_{k k} (λ)}{ε^{2} + λ^{2}},

(53)

where $E_{k k} (ε) \geq 0$ 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²⁰ $S_{k k}^{*} \sim ε + μ_{k k}^{0} / ε$ ⁠. The enhancement $E_{k k} (ε) \sim μ_{k k}^{0} / ε$ 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 $E_{k k} (ε)$ can introduce singular behavior that competes with the small ɛ prefactor in front of the integral, giving rise to a significant enhancement $E_{k k} (ε)$ for 0 < ɛ ≪ 1. Although, if the mass of the measure is zero or extremely small in a neighborhood of λ = 0, then the enhancement $E_{k k} (ε)$ 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 $E_{k k} (ε)$ 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 $S_{j k}$ 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 $S_{j k}$ 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 $S_{j k}$ is still valid—illustrating that the details of the spectral measure μ_jk near the spectral origin λ = 0 strongly influence the ɛ-behavior of $S_{j k}$ when ɛ ≪ 1.

These concepts are illustrated in our computations of μ_jk and $S_{j k}$ 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 $S_{j k}$ near ɛ = 0, as well as a significant enhancement in $S_{k k}$ 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 $D_{k k}^{*}$ for shear flow in the kth direction,⁴ where $μ_{k k} = μ_{k k}^{0} δ_{0} (d λ)$ and $D_{k k}^{*} = ε + μ_{k k}^{0} / ε$ ⁠. As another benchmark result, we demonstrate in Sec. V B 2 that our computations accurately capture the known^17,18 asymptotic behavior $D_{k k}^{*} \sim ε^{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 $D_{k k}^{*}$ 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²) numerical complexity of the method compared to the O(N³) numerical complexity of the spectral measure method.

1. BC-shear flow

In the continuum setting, it is known⁴ for shear flow in the x-direction that the measure μ₁₁ is given by a δ-measure concentrated at the spectral origin, while μ₂₂ ≡ 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.

FIG. 2.

View large Download slide

Shear flow baseline result. The streamlines of BC-shear-flow in (a) the x-direction and (b) the y-direction. (c) The ɛ behavior of $S_{j k}^{*}$ ⁠, j, k = 1, 2, for shear flow in the x-direction. The weights m_jk of the spectral measure Re μ_jk for shear flow in (d) the x-direction and (e) the y-direction as a function of λ. Consistent with theoretical predictions, the measure associated with the direction of the flow resembles a delta measure centered at the origin, while the other two components have spectral weights m_jk with very small magnitudes.

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 $S_{j k}^{*}$ ⁠, 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 $S_{11}^{*} \leftrightarrow S_{22}^{*}$ ⁠, 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₁, the spectral weights m₂₂(n) in Fig. 2(d) satisfy m₂₂(n) ≲ 10⁻²⁹. With the effects of finite resolution (L < ∞) and roundoff error associated with a machine epsilon of ∼10⁻¹⁶, these spectral weights can be considered “numerically zero.” The spectral weights $m_{12}^{\pm} (n)$ satisfy $m_{12}^{\pm} (n) ≲ 1 0^{- 28}$ for λ_n away from the spectral origin λ = 0 with a peak near λ = 0 having magnitudes $m_{12}^{\pm} (n) ≲ 1 0^{- 16}$ ⁠. The spectral weights for the x-direction satisfy m₁₁(n) ≲ 10⁻²⁸ away from λ = 0, while the weights near λ = 0 satisfy 10⁻⁹ ≲ m₁₁ ≲ 10⁻¹, resembling a δ-measure with virtually all of its mass concentrated near λ = 0. This is consistent with theoretical predictions.⁴ Due to the antisymmetry of the real-valued matrix $U_{1}^{T} H U_{1}$ ⁠, its complex eigenvectors and purely imaginary eigenvalues come in complex conjugate pairs.²⁵ Consequently, the eigenvalues of the Hermitian matrix $- ı U_{1}^{T} H U_{1}$ come in ± pairs with identical spectral weights, resulting in the symmetry about λ = 0 displayed by the spectral measures in Fig. 2.

FIG. 3.

View large Download slide

Transition away from cat’s eye cell flow. The spectral weights m_jk for the components Re μ_jk, j, k = 1, 2, of the spectral measure are displayed with increasing values of the free parameter α from left to right. As the parameter α increases, the streamlines of the flow transition away from cell structures to open channels. This is reflected in the measure by a dramatic increase in the magnitude of the spectral weights m_jk associated with the “accumulation point” of the measure at λ = 0, while the other weights change only slightly.

Due to the high concentration of measure mass in μ₁₁ very near the spectral origin, our computation of $S_{11}^{*}$ 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 $S_{11}^{*}$ ⁠, displayed in black color with solid line-style, lays right on top of the graph of its upper bound $ε [1 + μ_{11}^{0} / ε^{2}]$ given in (24), with $μ_{11}^{0} \approx 4.975 \times 1 0^{- 1}$ ⁠, 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₂₂ and $m_{12}^{\pm}$ ⁠, with measure masses $μ_{22}^{0} \approx 5.33 \times 1 0^{- 29}$ ⁠, ${[μ_{12}^{0}]}^{+} \approx 1.03 \times 1 0^{- 16}$ ⁠, and ${[μ_{12}^{0}]}^{-} \approx 3.34 \times 1 0^{- 15}$ ⁠, both the upper and lower bounds for $S_{22}^{*}$ and $S_{12}^{*}$ in Eqs. (24) and (25) are very close to ɛ and 0, respectively; the graph of $S_{22}^{*}$ is virtually right on top of the lower bound ɛ in cyan color and solid line-style, and the magnitudes of ${[S_{12}^{*}]}^{+}$ and ${[S_{12}^{*}]}^{-}$ 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 = −ıΔ⁻¹[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 μ₁₁ = μ₂₂. 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 μ₁₁ = μ₂₂ = |μ₁₂|, where $| μ_{12} | = μ_{12}^{+} + μ_{12}^{-}$ is the total variation of the signed measure μ₁₂ introduced in Eq. (27), i.e., superimposing the panels for $m_{12}^{+}$ and $m_{12}^{-}$ in Figs. 3 and 5 yields the figure panels for m₁₁ and m₂₂. We have also numerically verified the behavior μ₁₁ = μ₂₂ = |μ₁₂|.

Since the operator A = Δ⁻¹[∇ ⋅ u] is compact,^9,10 its spectrum is discrete with an accumulation point at the spectral origin λ = 0.⁵⁶ 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.¹⁷ In this case, the computed “accumulation point” has eigenvalues λ_n with very small magnitude 10⁻¹⁹ ≲ |λ_n| ≲ 10⁻¹⁴. (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 $m_{j k}^{\pm} (n)$ of the accumulation point have even smaller magnitudes with 10⁻³⁵ ≲ m_kk(n) ≲ 10⁻²⁹ and similarly for $m_{j k}^{\pm} (n)$ ⁠, as shown in Fig. 3. Consequently, this component of the spectral measure does not contribute significantly to the enhancement $E_{j k} (ε) = ε \sum_{n} m_{j k} (n) / (ε^{2} + λ_{n}^{2})$ of $S_{j k}^{*} (ε) = ε δ_{j k} + E_{j k} (ε)$ ⁠. 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⁻¹⁴ ≲ |λ_n| ≲ 10⁻⁷, as shown in Fig. 4. The other component of the spectral measure has spectra in the interval 10⁻⁷ ≲ |λ_n| ≲ 10⁰ and weights 10⁻³⁷ ≲ m_kk(n) ≲ 10⁻¹. However, the part of this component of the spectral measure with spectral weights having more significant magnitudes 10⁻⁴ ≲ m_kk(n) ≲ 10⁻¹ are associated with eigenvalues with magnitude |λ_n| ≳ 10⁻¹ and consequently also do not contribute significantly to the enhancement $E_{j k}$ ⁠.

FIG. 4.

View large Download slide

Migration of positive measure mass from the computed “accumulation point” near the spectral origin. As the free parameter α of cat’s eye flow increases from zero, the magnitude of the spectral weights $m_{12}^{\pm}$ comprising the accumulation point increases dramatically. Moreover, the spectra associated with the positive weights $m_{12}^{+}$ migrate away from the spectral origin until the accumulation point is composed only of negative valued mass. The corresponding behavior for m_kk, k = 1, 2, is determined by the relation μ_kk = |μ₁₂| between the components of the spectral measure $d μ_{j k} (λ) = \sum_{n} ⟨ m_{j k} (n) δ_{λ_{n}^{1}} (d λ) ⟩$ ⁠.

FIG. 5.

View large Download slide

Transition toward cat’s eye shear flow. The spectral weights m_jk for the components Re μ_jk, j, k = 1, 2, of the spectral measure are displayed with increasing values of the free parameter α from left to right. As the value of the parameter α increases, the streamlines become more elongated in the x-y diagonal direction, becoming shear flow when α = 1. This is reflected in the spectral measure by an increase in the breadth of the spectral region with significant measure mass.

When 0 < α ≪ 1, open channels connect neighboring cells and large scale transport takes place both in thin boundary layers and within the channels.¹⁷ This is reflected in the spectral measure $d μ_{j k} (λ) = \sum_{n} ⟨ m_{j k} (n) δ_{λ_{n}^{1}} (d λ) ⟩$ in Eq. (50) by a dramatic increase in the magnitude of the spectral weights m_kk(n) and $m_{12}^{\pm} (n)$ associated with the accumulation point—by more than 14 orders of magnitude with 10⁻¹⁹ ≲ m_kk(n) ≲ 10⁻¹⁵—corresponding to a change of only 10⁻⁶ 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⁻⁶ 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⁻⁷ ≲ |λ_n| ≲ 10⁰—decreasing the mass of the accumulation point for μ_kk. Moreover, all of the positive masses $m_{12}^{+} (n)$ migrate away from the accumulation point so the accumulation point of μ₁₂ 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 $E_{j k} (ε)$ for ɛ ∼ 10⁻¹⁴, for example, but only a moderate increase, and the increase in $E_{j k} (ε)$ for 10⁻⁷ ≲ ɛ ≲ 10⁰ 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⁻¹¹ ≲ m_kk(n) ≲ 10⁻⁵, associated with eigenvalues with magnitudes 10⁻¹⁹ ≲ |λ_n| ≲ 10⁻¹⁴, as shown in Fig. 4. This marked increase in the magnitudes of the spectral weights provide a significant contribution to the enhancement $E_{j k} (ε)$ for ɛ ∼ 10⁻¹⁴, for example. Moreover, as α increases in the range (10⁻¹, 10⁰), 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⁻⁴ ≲ |λ_n| ≲ 10⁻¹, with a marked increase in weight magnitude from 10⁻⁷ ≲ m_kk(n) ≲ 10⁻⁴ to 10⁻⁷ ≲ m_kk(n) ≲ 10⁻². This provides a significant contribution to the enhancement $E_{j k} (ε)$ even for 10⁻⁴ ≲ ɛ ≲ 10⁻¹.

The behavior of $S_{j k}^{*}$ that we deduced from the behavior of the spectral measure μ_jk is consistent with our computations of the components $S_{j k}^{*}$ ⁠, 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 $S_{j k}^{*}$ of the effective diffusivity approach their bare diffusive value ɛ δ_jk as ɛ surpasses ɛ = 1, as discussed in the beginning of Sec. V B.

FIG. 6.

View large Download slide

Transitional behavior of the effective diffusivity from cat’s eye cell flow to shear flow. The behavior of the components $S_{j k}^{*}$ ⁠, j, k = 1, 2, of the effective diffusivity as a function of the nondimensionalized molecular diffusivity ɛ and increasing values of the free parameter α from left to right and top to bottom. The upper and lower bounds corresponding to $S_{j k}^{*}$ are in dashed-dotted and dashed line-style, respectively, and are the same line colors as the line colors for $S_{k k}^{*}$ ⁠, k = 1, 2, and red for $S_{12}^{*}$ ⁠. The trivial lower bound ɛ for $S_{k k}^{*}$ is in cyan color and solid line-style. When the lower bound for $S_{12}^{*}$ becomes negative, it tends to −∞ in this log-scale. For α = 0, a polynomial fit to $S_{k k}^{*}$ is also displayed in magenta color and dashed line-style. As the value of the parameter α increases and the flow transitions from the cell to shear structure, there is a substantial enhancement in the effective diffusivity for small values of ɛ.

For cat’s eye cell-flow, when α = 0 the log-log plot of $S_{k k}^{*}$ displays a linear trend for 10⁻³ ≲ ɛ ≲ 10⁻¹, capturing the known^17,18 power law behavior $S_{k k}^{*} \sim ε^{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⁻⁷ ≲ m_kk(n) ≲ 10⁻⁴ and associated eigenvalues 10⁻⁴ ≲|λ_n|≲ 10⁻¹ gives rise to a moderate enhancement in $S_{k k}^{*}$ ⁠. 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⁻¹–10⁻³, as shown in Fig. 3 for α = 0.

We deduced from Fig. 3 that for α ∈ (0, 10⁻²) the behavior of the spectral measure gives rise to only a moderate enhancement $E_{j k} (ε) = ε \sum_{n} m_{j k} (n) / (ε^{2} + λ_{n}^{2})$ of the effective diffusivity $S_{j k}^{*} (ε) = ε δ_{j k} + E_{j k} (ε)$ both for the advection dominated regime where ɛ ∼ 10⁻¹⁴, for example, and for the transitional regime 10⁻³ ≲ ɛ ≲ 10⁻¹. This is consistent with the behavior of $S_{j k}^{*}$ shown in Fig. 6 for α = 0 and α = 0.05. We further deduced from Figs. 4 and 5 that for α ∈ (10⁻², 1) the behavior of the spectral measure gives rise to a significant enhancement for both ɛ ∼ 10⁻¹⁴ and 10⁻³ ≲ ɛ ≲ 10⁻¹ where there is a marked increase in spectral weight magnitude from 10⁻⁷ ≲ m_kk(n) ≲ 10⁻⁴ to 10⁻⁷ ≲ m_kk(n) ≲ 10⁻² for eigenvalues satisfying 10⁻⁴ ≲ |λ_n| ≲ 10⁻¹. This is consistent with the panels of Fig. 6 corresponding to 0.1 ≤ α ≤ 1, with $S_{j k}^{*}$ enhanced many orders of magnitude above its bare diffusive value δ_jk ɛ, and $S_{k k}^{*}$ as well as $S_{12}^{*}$ closely following their upper bounds for ɛ ≳ 10^−0.5 when α = 1.

FIG. 7.

View large Download slide

Spectral functions and effective diffusivities for randomly perturbed cat’s eye flow. The random parameter α is uniformly distributed on the interval [0, p]. The spectral functions μ_jk(λ) are displayed with the corresponding effective diffusivities $S_{j k}^{*}$ directly below for various values of p, increasing from left to right. As p increases and the streamlines of the flow become more elongated in the x-y direction, on average, the region about the spectral origin λ = 0 with substantial measure mass increases in breadth and magnitude. This gives rise to a substantial enhancement in the components $S_{j k}^{*}$ of the effective diffusivity for larger values of the nondimensionalized molecular diffusivity ɛ. The color scheme of the panels for $S_{j k}^{*}$ is the same as that in Fig. 6.

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 μ₁₁ = μ₂₂ = |μ₁₂|. This gives rise to the symmetry $S_{11}^{*} (ε) = S_{22}^{*} (ε) = ε {+ E}_{12} (ε; μ_{12}^{+}) {+ E}_{12} (ε; μ_{12}^{-})$ between the components of the effective diffusivity, where we have denoted $E_{j k} (ε; μ_{j k}) = ε \int d μ_{j k} (λ) / (ε^{2} + λ^{2})$ ⁠, e.g., $S_{11}^{*} (ε) = ε + E (ε; μ_{11})$ ⁠. The symmetry $S_{11}^{*} = S_{22}^{*}$ can be clearly seen in our computations of $S_{j k}^{*}$ ⁠, 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 $μ_{11}^{0} = μ_{22}^{0}$ ⁠. We have also numerically explored the empirical relationship $S_{11}^{*} \approx ε + {[S_{12}^{*}]}^{+} + {[S_{12}^{*}]}^{-}$ by plotting $S_{11}^{*}$ and $ε + {[S_{12}^{*}]}^{+} + {[S_{12}^{*}]}^{-}$ 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 $ε + {[S_{12}^{*}]}^{+} + {[S_{12}^{*}]}^{-}$ from $S_{11}^{*}$ ⁠, it is slight. This property also leads to the inequalities $S_{11}^{*} \geq ε + {[S_{12}^{*}]}^{+}$ and $S_{11}^{*} \geq ε + {[S_{12}^{*}]}^{-}$ ⁠, with $S_{11}^{*} = S_{22}^{*}$ ⁠, which is consistent with the behavior of $S_{j k}^{*}$ shown in Fig. 6.

3. Randomly perturbed cat’s eye flow

We now discuss our computations of the components μ_jk and $S_{j k}^{*}$ ⁠, 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 Ω₀, with |Ω₀| ∼ 10³ and L = 100, we computed every eigenvalue $λ_{n}^{1}$ and eigenvector $r_{n}^{1}$ ⁠, n = 1, …, K₁, of the matrix $- ı U_{1}^{T} H U_{1}$ to form the spectral measure μ_jk in Eq. (50). In order to visually determine the behavior of the function $μ_{j k} (λ) = ⟨ Q (λ) {\hat{e}}_{j}, {\hat{e}}_{k} ⟩$ 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². Second, for fixed v, we identified all of the eigenvalues that satisfy $λ_{n}^{1} (ω) \in I_{v}$ ⁠, for n = 1, …, K₁ and ω ∈ Ω₀. 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 $λ_{n}^{1} (ω) \in I_{v}$ ⁠, normalized by |Ω₀|. In this way, the area underneath the curve is the measure mass $μ_{j k}^{0}$ ⁠.

Consistent with the symmetries of the randomly perturbed flow, our computations of the spectral function satisfy μ₁₁(λ) = μ₂₂(λ); hence, the ensemble averaged components $S_{j k}^{*}$ of the effective diffusivity also satisfy $S_{11}^{*} = S_{22}^{*}$ ⁠, 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 $S_{j k}^{*}$ of the effective diffusivity above the bare molecular diffusivity values ɛ δ_jk. The sign changes in μ₁₂(λ) give rise to sign changes in $S_{12}^{*} = {[S_{12}^{*}]}^{+} - {[S_{12}^{*}]}^{-}$ ⁠. In the log-log plot of $S_{12}^{*}$ ⁠, a negative singularity forms in $S_{12}^{*}$ 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 $S_{j k}^{*}$ ⁠, 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 $S_{j k}^{*}$ and Re μ_jk for BC-shear-flow, for which the spectral measure is known.⁴ Our numerical results are in good agreement with the theoretical result. We also computed $S_{j k}^{*}$ 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 $S_{k k}^{*} \sim ε^{1 / 2}$ for ɛ ≪ 1. The spectral measure Re μ_jk and $S_{j k}^{*}$ have transitional behavior as the parameter varies. This reveals how the details of Re μ_kk near the spectral origin govern the enhancement of $S_{k k}^{*}$ above its bare diffusive value ɛ in the advection dominated regime where ɛ ≪ 1 and similarly for $S_{12}^{*}$ ⁠. Consistent with the symmetries of the flow, our computations indicate that Re μ₁₁ = Re μ₂₂. Our computations of Re μ₁₂ for cat’s eye flow also suggest a deeper symmetry, namely, |Re μ₁₂| = Re μ₁₁ = Re μ₂₂, where |Re μ₁₂| is the total variation of the signed measure Re μ₁₂. Our computations of $S_{j k}^{*}$ 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.⁴⁵ 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 = Δ⁻¹[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⁴⁴ 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.⁴⁴ 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.¹ 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 = Δ⁻¹[u ⋅ ∇] defined in Eq. (11). In particular, we demonstrate that A is antisymmetric on the Hilbert space $H^{1, 2}$ defined in (8). Moreover, we show that A is bounded on $H^{1, 2}$ and we provide an upper bound for ∥A∥_1,2 when u is uniformly bounded on the period cell $V$ ⁠.

We first show that the incompressibility condition ∇ ⋅ u = 0 implies that the operator A is antisymmetric on $H^{1, 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 Δ⁻¹ satisfies ⟨ΔΔ⁻¹f, h⟩ = ⟨f, h⟩ in a distributional sense, for all $f, h \in H$ ⁠.^19,37,44 Consequently, integration by parts and ∇ ⋅ u = 0 yields^9,10,44,48

\begin{align} {⟨ A f, h ⟩}_{1, 2} & = ⟨ [\nabla (Δ^{- 1}) (u \cdot \nabla) f] \cdot \nabla h ⟩ \\ = - ⟨ [(u \cdot \nabla) f], h ⟩ \\ = - ⟨ [\nabla \cdot (u f)], h ⟩ \\ = ⟨ f, [(u \cdot \nabla) h] ⟩ \\ = ⟨ f, [Δ (Δ^{- 1}) (u \cdot \nabla) h] ⟩ \\ = - ⟨ \nabla f \cdot [\nabla (Δ^{- 1}) (u \cdot \nabla) h] ⟩ \\ = - {⟨ f, A h ⟩}_{1, 2}, \end{align}

(C1)

for all $f, h {\in H}^{1, 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

\begin{align} {‖ A f ‖}_{1, 2}^{2} & = | ⟨ \nabla [Δ^{- 1} (u \cdot \nabla f)] \cdot \nabla [Δ^{- 1} (u \cdot \nabla f)] ⟩ | \\ = | - ⟨ [Δ^{- 1} (u \cdot \nabla f)], (u \cdot \nabla f) ⟩ | \\ \leq ‖ Δ^{- 1} (u \cdot \nabla f) ‖ ‖ u \cdot \nabla f ‖ \\ \leq ‖ Δ^{- 1} ‖ {‖ u \cdot \nabla f ‖}^{2} . \end{align}

(C2)

We now provide an upper bound for ∥u ⋅ ∇f∥ when the components u_k, k = 1, …, d, of the fluid velocity field u are uniformly bounded on the period cell $V$ ⁠. By the Cauchy-Schwartz inequality, |ξ ⋅ ζ| ≤ |ξ| |ζ|, we have

{‖ u \cdot \nabla f ‖}^{2} = ⟨ | u \cdot \nabla f |^{2} ⟩ \leq ⟨ | u |^{2} | \nabla f |^{2} ⟩ \leq \sup_{x \in V} | u (x) |^{2} {‖ f ‖}_{1, 2}^{2} .

(C3)

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

ϕ_{t} = \nabla \cdot [D \nabla ϕ], ϕ (0, x) = ϕ_{0} (x), D = ε I + H .

(D1)

Moreover, the cell problem in (4) can be written as a steady-state diffusion equation,^17,18

\nabla \cdot [D (\nabla χ_{k} + e_{k})] = 0, ⟨ \nabla χ_{k} ⟩ = 0, k = 1, \dots, d .

(D2)

Here, e_k is a standard basis vector, k = 1, …, d, and D(x) = ɛI + H(x) can be viewed as a local diffusivity matrix with coefficients

D_{j k} = ε δ_{j k} + H_{j k}, j, k = 1, \dots, d,

(D3)

where δ_jk is the Kronecker delta and I is the identity operator on $R^{d}$ ⁠.

Substituting into Eq. (3) the expression for u_j in (4) and using Eq. (28), u = ∇ ⋅ H, show that the components $S_{j k}^{*}$ and $A_{j k}^{*}$ of $S^{*}$ and $A^{*}$ can be written in terms of the following functional formulas involving the real-valued vector field ∇χ_k:

S_{j k}^{*} = ε (δ_{j k} + ⟨ \nabla χ_{j} \cdot \nabla χ_{k} ⟩), A_{j k}^{*} = ⟨ H \nabla χ_{j} \cdot \nabla χ_{k} ⟩ .

(D4)

The functional formulas in (D4) are analogous to the functional formulas in Eq. (11). The symmetry $S_{j k}^{*} = S_{k j}^{*}$ 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|²⟩≥ 0 demonstrates that the effective transport of the scalar density ϕ is always enhanced by the presence of an incompressible velocity field, i.e., $D_{k k}^{*} = S_{k k}^{*} \geq ε$ ⁠. The equality $D_{k k}^{*} = S_{k k}^{*}$ follows from the skew-symmetry of the matrix $A^{*}$ ⁠, $A_{k j}^{*} = - A_{j k}^{*}$ ⁠, hence $A_{k k}^{*} = 0$ ⁠. The skew-symmetry of $A^{*}$ follows from the skew-symmetry of the real-valued matrix H, $A_{j k}^{*} = ⟨ H \nabla χ_{j} \cdot \nabla χ_{k} ⟩ = - ⟨ \nabla χ_{j} \cdot H \nabla χ_{k} ⟩ = - ⟨ H \nabla χ_{k} \cdot \nabla χ_{j} ⟩ = - A_{k j}^{*}$ ⁠.

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,

\nabla \times E_{k} = 0, \nabla \cdot J_{k} = 0, J_{k} = D E_{k}, ⟨ E_{k} ⟩ = e_{k}, D = ε I + H .

(D5)

Here, E_k = ∇χ_k + e_k plays the role of the local electric field, J_k = DE_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 = DE_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,²² 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²² or matrix⁴² 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 \subset R^{d}$ and define the Hilbert space $H$ ⁠,^18,33

H = {ξ \in \otimes_{n = 1}^{d} L^{2} (V, ν) : ξ (x) is periodic in V and ⟨ ξ ⟩ = 0},

(D6)

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 $⟨ ζ, ξ ⟩ = \bar{⟨ ξ, ζ ⟩}$ for $ξ, ζ \in H$ ⁠. This inner-product induces a norm ∥⋅∥ defined by ∥ξ∥ = ⟨ξ,ξ⟩^1/2 and $ξ \in H$ implies that ∥ξ∥ < ∞. Now consider the associated Hilbert space $H_{\times}$ of curl-free fields,^{4,17,18,22,40,44}

H_{\times} = \{ξ \in H : \nabla \times ξ = 0 weakly and ⟨ ξ ⟩ = 0\} .

(D7)

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,²⁰ ∥D∥ < ∞, then there exists unique $\nabla χ_{k} {\in H}_{\times}$ satisfying Eq. (D2).^22,46 We assume that

0 < ε < \infty, ‖ H ‖ < \infty,

(D8)

which together imply ∥D∥ < ∞.

The linear operator Γ = ∇(Δ⁻¹)∇ ⋅ is a projection onto the Hilbert space $H_{\times}$ in the sense that $Γ : H {\mapsto H}_{\times}$ and Γξ = ξ (weakly) for all $ξ {\in H}_{\times}$ ⁠, in particular, Γ ∇χ_k = ∇χ_k.^17,44 It is based on convolution with respect to the Green’s function for the Laplacian Δ = ∇².^38,56 Applying the integro-differential operator ∇Δ⁻¹ 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 = −ΓHe_k, which yields the following resolvent formula for ∇χ_k:

\nabla χ_{k} = {(ε I + A)}^{- 1} g_{k}, A = Γ H Γ, g_{k} = - Γ H e_{k},

(D9)

which is analogous to Eq. (12). Since Γ is a projection operator onto $H_{\times} \subset H$ ⁠, it is bounded by unity in operator norm on $H$ ⁠, ∥Γ∥ ≤ 1.^20,54 Integration by parts and the symmetry of the operator Δ⁻¹⁵⁶ (or the projective nature of Γ itself) shows that Γ is also a symmetric operator, i.e., ⟨Γξ ⋅ ζ⟩ = ⟨ξ ⋅ Γζ⟩ for all $ξ, ζ \in H$ ⁠.⁴⁴ These two properties show that Γ with domain $H$ is a self-adjoint operator.⁵³ Since Γ is self-adjoint and Γ∇χ_k = ∇χ_k, we may write $A_{j k}^{*}$ in Eq. (D4) as $A_{j k}^{*} = ⟨ A \nabla χ_{j} \cdot \nabla χ_{k} ⟩$ ⁠. Consequently, substituting the resolvent formula for ∇χ_k in Eq. (D9) into the functional formulas in Eq. (D4) yields

\begin{align} S_{j k}^{*} & = ε (δ_{j k} + ⟨ {(ε I + A)}^{- 1} g_{j} \cdot {(ε I + A)}^{- 1} g_{k} ⟩), \\ A_{j k}^{*} & = ⟨ A {(ε I + A)}^{- 1} g_{j} \cdot {(ε I + A)}^{- 1} g_{k} ⟩, \end{align}

(D10)

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 $ı = \sqrt{- 1}$ ⁠, 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,⁵³ thus

Σ \subseteq [- ‖ H ‖, ‖ H ‖] .

(D11)

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

Corollary 4.

Let Q(λ) denote the resolution of the identity corresponding to the self-adjoint operator M. Then, the components

S_{j k}^{*}

and

A_{j k}^{*}

, 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 () with μ_jk replaced by the spectral measure

{\tilde{μ}}_{j k}

of M in the

(g_{j}, g_{k})

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

S_{j k}^{*}

and

A_{j k}^{*}

. The mass

{\tilde{μ}}_{j k}^{0}

of the measure

{\tilde{μ}}_{j k}

is real-valued and satisfies

{\tilde{μ}}_{j k}^{0} = ⟨ H^{T} Γ H e_{j} \cdot e_{k} ⟩, | μ_{j k}^{0} | \leq {‖ H ‖}^{2} < \infty .

(D12)

Proof of Corollary 4.

Exactly as in the Proof of Theorem 1, Eq. () and a direct analog of Eq. () lead to the Stieltjes integral representations in (), involving a Stieltjes measure

{\tilde{μ}}_{j k}

associated with the function of the spectral variable λ defined by

{\tilde{μ}}_{j k} (λ) = ⟨ Q (λ) g_{j}, g_{k} ⟩

⁠. Here,

g_{k} = - Γ H e_{k}

is defined in () 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_{\times}

⁠.^53,58 From Eq. () and the fact that Γ is a self-adjoint projection operator on

H_{\times}

⁠, the mass

{\tilde{μ}}_{j k}^{0}

of the measure

{\tilde{μ}}_{j k}

is real-valued and satisfies

{\tilde{μ}}_{j k}^{0} = ⟨ g_{j}, g_{k} ⟩ = ⟨ Γ H e_{j} \cdot Γ H e_{k} ⟩ = ⟨ H^{T} Γ H e_{j} \cdot e_{k} ⟩, | {\tilde{μ}}_{j k}^{0} | \leq {‖ H ‖}^{2} < \infty .

(D13)

By analogy, the bounds in Theorem 2 also hold for these integral representations for

S_{j k}^{*}

and

A_{j k}^{*}

⁠. This concludes our Proof of Corollary 4.

It is worth making the following observation. In the current setting, the formulas for $S_{j k}^{*}$ and $A_{j k}^{*}$ in (19) are computed with respect to the standard basis {e_k}, through the definition of ${\tilde{μ}}_{j k} (λ) = ⟨ Q (λ) g_{j}, g_{k} ⟩$ with $g_{k} = - Γ H e_{k}$ ⁠. We now show that, given $S_{j k}^{*}$ and $A_{j k}^{*}$ ⁠, 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 ${\tilde{μ}}_{j k} (λ)$ ⁠. More specifically, if $ξ, ζ \in R^{d}$ are arbitrary directions of interest, then $⟨ Q (λ) Γ H ξ, Γ H ζ ⟩ = \sum_{j k} a_{j} b_{k} ⟨ Q (λ) g_{j}, g_{k} ⟩$ ⁠, 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.

Theorem 5.

The spectral measure μ_jk in Theorem 1 is identical to the spectral measure

{\tilde{μ}}_{j k}

in Corollary 4,

μ_{j k} \equiv {\tilde{μ}}_{j k} .

(D14)

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 $g_{k} = - Γ H e_{k}$ defined in Eq. (D9) and the scalar field g_k = (−Δ)⁻¹u_k defined in (12) are (weakly) related by⁴⁴

g_{k} = \nabla g_{k} .

(D15)

The operator A = ΓHΓ defined in Eq. (D9) and the operator A = Δ⁻¹[u ⋅ ∇] defined in (12), hence M = −ıA and M = −ıA, are (weakly) related by

A \nabla = \nabla A, M \nabla = \nabla M .

(D16)

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

μ_{j k}^{n} = {\tilde{μ}}_{j k}^{n}, n = 0, 1, 2, \dots .

(D17)

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).

Proof of Theorem 5.

In Ref. 44, Eqs. () and () were established in a more general context than our considerations here, where u = u(x, t) is periodic in both space and time t. We will see that the proof of Eq. () essentially depends only on Eqs. () and (). Consequently, en route, we will also establish that Eq. () 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 (), u = ∇ ⋅ H, and the definitions Γ = ∇(Δ⁻¹)∇⋅,

g_{k} = - Γ H e_{k}

⁠, and g_k = (−Δ)⁻¹u_k, together yield⁴⁴ Eq. (),

g_{j} = - Γ H e_{j} = \nabla {(- Δ)}^{- 1} u_{j} = \nabla g_{j} .

(D18)

Consequently, by Eq. (), the masses

μ_{j k}^{0}

and

{\tilde{μ}}_{j k}^{0}

of the spectral measures μ_jk and

{\tilde{μ}}_{j k}

are equal,

{\tilde{μ}}_{j k}^{0} = ⟨ g_{j} \cdot g_{k} ⟩ = ⟨ \nabla g_{j} \cdot \nabla g_{k} ⟩ = {⟨ g_{j}, g_{k} ⟩}_{1, 2} = μ_{j k}^{0} .

(D19)

Since Γ∇ξ = ∇ξ (weakly),⁴⁴ Eqs. () and (), ∇ ⋅ [H∇] = [∇ ⋅ H] ⋅ ∇, imply that for all

ξ {\in H}^{1, 2}

⁠,

A \nabla ξ = Γ H Γ \nabla ξ = Γ H \nabla ξ = \nabla [(Δ^{- 1}) u \cdot \nabla] ξ = \nabla A ξ,

(D20)

in a weak sense,⁴⁴ which establishes Eq. (). It follows from () and () that

\begin{align} {\tilde{μ}}_{j k}^{1} & = ⟨ M g_{j} \cdot g_{k} ⟩ = ⟨ M \nabla g_{j} \cdot \nabla g_{k} ⟩ = ⟨ \nabla M g_{j} \cdot \nabla g_{k} ⟩ = {⟨ M g_{j}, g_{k} ⟩}_{1, 2} = μ_{j k}^{1}, \\ {\tilde{μ}}_{j k}^{2} & = ⟨ M^{2} g_{j} \cdot g_{k} ⟩ = ⟨ M \nabla M g_{j} \cdot \nabla g_{k} ⟩ = ⟨ \nabla M^{2} g_{j} \cdot \nabla g_{k} ⟩ = {⟨ M^{2} g_{j}, g_{k} ⟩}_{1, 2} = μ_{j k}^{2} . \end{align}

(D21)

An inductive argument establishes Eq. (). This concludes our Proof of Theorem 5.

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

Corollary 6.

For each

ξ {\in H}^{1, 2}

⁠, we have

\nabla ξ {\in H}_{\times}

, and conversely, for each

ξ {\in H}_{\times}

⁠, there exists unique

ξ {\in H}^{1, 2}

such that ξ = ∇ξ. Consequently, by Theorem 5, for every

ξ, ζ {\in H}^{1, 2}

and

ξ, ζ {\in H}_{\times}

related by ξ=∇ξ and ζ=∇ζ, we have

μ_{ξ ζ}^{n} = {\tilde{μ}}_{ξ ζ}^{n}, n = 0, 1, 2, \dots .

(D22)

This establishes that the spectral measures μ_ξζ and

{\tilde{μ}}_{ξ ζ}

are identical,

μ_{ξ ζ} \equiv {\tilde{μ}}_{ξ ζ} .

(D23)

Moreover, Eq. () also holds for the class of space-time periodic fluid velocity fields u=u(x, t) described in Ref. 44 .

Proof of Corollary 6.

The Hilbert spaces $H^{1, 2}$ and $H_{\times}$ are in one-to-one isometric correspondence.⁴⁴ More specifically, for every $ξ {\in H}^{1, 2}$ ⁠, we have $\nabla ξ {\in H}_{\times}$ satisfying ∥ξ∥_1,2 = ∥∇ξ∥ and ∥Aξ∥_1,2 = ∥A∇ξ∥. Conversely, for each $ξ {\in H}_{\times}$ ⁠, there exists unique⁴⁴ $ξ {\in H}^{1, 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⁴⁴ 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.⁴⁴ 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,⁴⁴ one must first establish Carleman’s criterion,¹ 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 $(ε I + A) \nabla χ_{k} = g_{k}$ ⁠, 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 Γ = ∇(Δ⁻¹)∇⋅ is given by the symmetric projection matrix $Γ = \nabla {(\nabla^{T} \nabla)}^{- 1} \nabla^{T}$ satisfying Γ² = Γ and Γ∇ = ∇, where ${(\nabla^{T} \nabla)}^{- 1}$ is now interpreted as a matrix inversion. We assume here that the matrix ∇ is of full-rank so that ${(\nabla^{T} \nabla)}^{- 1}$ 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 $g_{k} = - Γ H e_{k}$ ⁠, 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 $g_{k}$ 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 Aw_n = υ_nw_n. Since A is skew-symmetric, the eigenvalues υ_n are purely imaginary,²⁵ υ_n = ıλ_n with $λ_{n} \in R$ ⁠. 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 $C^{N}$ ⁠,^25,28 i.e., $w_{n}^{†} w_{m} = δ_{n m}$ and for every $ξ \in C^{N}$ we have $ξ = \sum_{n} (w_{n}^{†} ξ) w_{n} = (\sum_{n} w_{n} w_{n}^{†}) ξ$ ⁠. Consequently, defining $Q_{n} = w_{n} w_{n}^{†}$ ⁠, 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):

\sum_{n = 1}^{N} Q_{n} = I, Q_{n} = w_{n} w_{n}^{†}, Q_{l} Q_{m} = Q_{l} δ_{l m} .

(E1)

With an argument similar to the one in the paragraph following Eq. (34), one can show that for all $ξ, ζ \in C^{N}$ 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 $μ_{ξ ζ} (λ) = ⟨ Q (λ) ξ \cdot ζ ⟩$ ⁠, where the associated matrix representation $Q (λ)$ of the projection operator Q(λ) and the discrete spectral measure dμ_ξζ(λ) are given by the following analog of Eq. (35):

Q (λ) = \sum_{n : λ_{n} \leq λ} θ (λ - λ_{n}) Q_{n}, d μ_{ξ ζ} (λ) = \sum_{n : λ_{n} \leq λ} ⟨ δ_{λ_{n}} (d λ) [Q_{n} ξ \cdot ζ] ⟩ .

(E2)

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

Corollary 7.

Consider the standard eigenvalue problem Mw_n = λ_nw_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

S_{j k}^{*}

and

A_{j k}^{*}

in Eq. (D4) with H replaced by A, as discussed right before Eq. (), are given by

S_{j k}^{*} = ε (δ_{j k} + ⟨ \nabla χ_{j} \cdot \nabla χ_{k} ⟩), A_{j k}^{*} = ⟨ A \nabla χ_{j} \cdot \nabla χ_{k} ⟩,

(E3)

which is analogous to Eq. (). Also, the discrete representation of the resolvent formula for ∇χ_j in Eq. () is given by

\nabla χ_{j} = W {(ε I + ı Λ)}^{- 1} W^{†} g_{j}, g_{j} = - Γ H e_{j},

(E4)

which is analogous to Eq. (). The discrete representations of the bilinear functional formulas for

S_{j k}^{*}

and

A_{j k}^{*}

in () are given by

\begin{align} S_{j k}^{*} & = ε (δ_{j k} + ⟨ {(ε I + ı Λ)}^{- 1} W^{†} g_{j} \cdot {(ε I + ı Λ)}^{- 1} W^{†} g_{k} ⟩), \\ A_{j k}^{*} & = ⟨ ı Λ {(ε I + ı Λ)}^{- 1} W^{†} g_{j} \cdot {(ε I + ı Λ)}^{- 1} W^{†} g_{k} ⟩, \end{align}

(E5)

which are analogous to those in Eq. (). Consequently, the formulas for

S_{j k}^{*}

and

A_{j k}^{*}

in () have the following series representations:

S_{j k}^{*} / ε - δ_{j k} = \sum_{n = 1}^{N} \frac{R e [\bar{(w_{n}^{†} g_{j})} (w_{n}^{†} g_{k})]}{ε^{2} + λ_{n}^{2}}, A_{j k}^{*} = \sum_{n = 1}^{N} \frac{λ_{n} Im [\bar{(w_{n}^{†} g_{j})} (w_{n}^{†} g_{k})]}{ε^{2} + λ_{n}^{2}},

(E6)

which is analogous to Eq. (). Finally, recalling the projection matrix

Q_{n} = w_{n} w_{n}^{†}

in (), we have

\bar{(w_{n}^{†} g_{j})} (w_{n}^{†} g_{k}) = Q_{n} g_{j} \cdot g_{k},

(E7)

which is analogous to Eq. (). It follows from Eq. () that the series representations for

S_{j k}^{*}

and

A_{j k}^{*}

in () have the Stieltjes integral representations in Eq. (), involving the discrete spectral measure μ_jk in () with

ξ = g_{j} = - Γ H e_{j}

and

ζ = g_{k} = - Γ H e_{k}

.

Proof of Corollary 7.

Equation () follows from the discrete version of () and Γ∇ = ∇ so ⟨H ∇χ_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.²⁵ Equation () follows from this spectral decomposition of A = ıM, the discrete version of the resolvent formula in Eq. () written as

\nabla χ_{j} = {(ε I + A)}^{- 1} g_{j}

⁠, the fact that W is unitary satisfying W^† = W⁻¹, and the elementary properties of matrix inversion.²⁵ Substituting the resolvent formula for ∇χ_j in () into Eq. () and using W^† = W⁻¹ yield Eq. (). The quadratic form

W^{†} g_{j} \cdot W^{†} g_{k}

arising in () can be written in terms of the projection matrices

Q_{n} = w_{n} w_{n}^{†}

defined in () as follows:

W^{†} g_{j} \cdot W^{†} g_{k} = \sum_{n = 1}^{N} \bar{(w_{n}^{†} g_{j})} (w_{n}^{†} g_{k}) = \sum_{n = 1}^{N} Q_{n} g_{j} \cdot g_{k} .

(E8)

Exactly as in Theorem 3, we have Eq. (), which can be written in terms of the integral representations in Eq. () involving the discrete spectral measure dμ_jk(λ) in Eq. () with

ξ = g_{j} = - Γ H e_{j}

and

ζ = g_{k} = - Γ H e_{k}

⁠. 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 $g_{k} = - Γ H e_{k}$ 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.⁷ 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.

Theorem 8

(Projection method). The real-symmetric projection matrix Γ of size N has eigenvalues satisfying γ_n = 0, 1, hence the spectral decomposition

Γ = P G P^{T}, G = diag (0_{N_{0}}, 1_{N_{1}}), P = [P_{0} P_{1}] .

(E9)

Here,

0_{N_{0}}

and

1_{N_{1}}

are vectors of zeros and ones with N₀ and N₁ components, respectively, where N = N₀ + N₁. Moreover, the columns comprising the N × N₀ matrix P₀ and the N × N₁ matrix P₁ are orthonormal eigenvectors that span the null space and range of Γ, respectively. Consequently, the matrix Γ can be written as

Γ = P_{1} P_{1}^{T} .

(E10)

Consider the spectral composition of the antisymmetric matrix

P_{1}^{T} H P_{1}

of size N₁,

P_{1}^{T} H P_{1} = ı R_{11} Λ_{11} R_{11}^{†},

(E11)

where R₁₁ is a unitary matrix,

R_{11}^{†} R_{11} = R_{11} R_{11}^{†} = I_{11}

, Λ₁₁ is a real-valued diagonal matrix, and I₁₁ is the identity matrix of size N₁ × N₁. Consequently, the matrix A = ıΓHΓ has the spectral decomposition

A = ı W Λ W^{†}, W = [P_{0} P_{1} R_{11}], Λ = diag (O_{00} Λ_{11}) .

(E12)

Equations () and (), and the mutual orthogonality of the matrices P₀ and P₁ yield

W^{†} g_{k} = W^{†} Γ H e_{k} = {[O_{0} P_{1} R_{11}]}^{†} H e_{k},

(E13)

where O₀ is a matrix of zeros of size N × N₀ . Consequently, Eq. () can be written as

W^{†} g_{j} \cdot W^{†} g_{k} = {(P_{1} R_{11})}^{†} H e_{k} \cdot {(P_{1} R_{11})}^{†} H e_{k} .

(E14)

Moreover, the spectral weights in Eq. () associated with γ_n = 0 are identically zero

[Q_{n} g_{j} \cdot g_{k}] = 0

.

Proof of Theorem 8.

Using the spectral decomposition Γ = PGP^T in (E9), we write the matrix A = ΓHΓ as $A = P [G (P^{T} H P) G] P^{T} = P diag (O_{00} P_{1}^{T} H P_{1}) P^{T}$ ⁠, where O₀₀ is a matrix of zeros of size N₀ × N₀. Due to the skew-symmetry of H, the N₁ × N₁ matrix $P_{1}^{T} H P_{1}$ is also skew-symmetric. Consequently, it has the spectral decomposition in Eq. (E11); hence, $A = P diag (O_{00} ı R_{11} Λ_{11} R_{11}^{†}) P^{T}$ ⁠. Writing this in block matrix form⁴² and multiplying yield Eq. (E12). Equations (E10) and (E12), and the mutual orthogonality of the matrices P₀ and P₁ 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(σ₁, …, σ_K), where 0 ≤ σ₁ ≤ ⋯ ≤ σ_K, and the matrices U and V are of size N × K and K × K, respectively, and satisfy¹⁴

U^{T} U = I, V^{T} V = V V^{T} = I,

(F1)

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Σ²V^T.¹⁴ We assume that ∇ is of full-rank so that σ_i > 0 for all i = 1, …, K. This implies that Σ⁻¹ exists so that the matrix Laplacian is invertible. In this case, it follows from ∇ = UΣV^T and Eq. (F1) that the projection matrix $Γ = \nabla {(\nabla^{T} \nabla)}^{- 1} \nabla^{T}$ is given by

Γ = U U^{T},

(F2)

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

From Eqs. (F1) and (F2), we can write the eigenvalue problem −ıΓHΓw_n = λ_nw_n discussed in Appendix E as

[- ı U^{T} H U] [U^{T} w_{n}] = λ_{n} [U^{T} w_{n}] .

(F3)

Now consider the generalized eigenvalue problem −ı∇^TH∇z_n = α_n∇^T∇z_n discussed in Sec. IV and recall that ∇ = UΣV^T and ∇^T∇ = VΣ²V^T. Since Σ is invertible, by Eq. (F1), we can write this generalized eigenvalue problem as the following standard eigenvalue problem:

[- ı U^{T} H U] [Σ V^{T} z_{n}] = α_{n} [Σ V^{T} z_{n}] .

(F4)

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^Tw_n = ΣV^Tz_n. Since Γ is a projection matrix, Γ² = Γ, the eigenvalue problem ΓHΓw_n = ıλ_nw_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^Tw_n = ΣV^Tz_n and recalling that Γ = UU^T and ∇ = UΣV^T, we have

w_{n} = \nabla z_{n} .

(F5)

In the following lemma, we make precise the correspondence between the standard eigenvalue problem −ıΓHΓw_n = λ_nw_n and the generalized eigenvalue problem −ı∇^TH∇z_n = α_n∇^T∇z_n, as well as the associated spectral measures in Eqs. (E2) and (35), respectively.

Lemma 9.

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

- ı Γ H Γ w_{n} = λ_{n} w_{n},

(F6)

- ı \nabla^{T} H \nabla z_{n} = λ_{n} \nabla^{T} \nabla z_{n} .

(F7)

Let ∇ = UΣV^T be the SVD of the matrix ∇, which we assume to be of full-rank. Then, Eq. () implies and is implied by Eq. (), with w_n and z_n related as in Eq. (). This implies that the spectrum associated with each of these eigenvalue problems is identical. Moreover, the spectral weights in Eqs. () and ( ) are identical; specifically,

Q_{n} Γ H e_{j} \cdot Γ H e_{k} = \nabla Q_{n} {[\nabla^{T} \nabla]}^{- 1} u_{j} \cdot \nabla {[\nabla^{T} \nabla]}^{- 1} u_{k} .

(F8)

This, in turn, implies that the associated spectral measures in Eqs. () and () are identical.

Proof of Lemma 9.

Recall that ∇ = UΣV^T, ∇^T∇ = VΣ²V^T, and Γ = UU^T, where Σ is invertible, and the matrices V and U satisfy Eq. (). First, consider Eq. () written as in Eq. (), [−ıU^THU][U^Tw_n] = λ_n[U^Tw_n]. Since the matrix Σ is invertible and V^TV = I, we can rewrite Eq. () as

V Σ [- ı U^{T} H U] (Σ V^{T}) (V Σ^{- 1}) [U^{T} w_{n}] = λ_{n} (V Σ^{2} V^{T}) (V Σ^{- 1}) [U^{T} w_{n}],

(F9)

which is precisely Eq. () written in terms of ∇ = UΣV^T with z_n = VΣ⁻¹U^Tw_n. This formula for z_n, Eq. (), and the formula Γw_n = w_n above Eq. () imply that w_n = UΣV^Tz_n = ∇z_n, which is the formula in (). Now consider Eq. () written as in (), [−ıU^THU][ΣV^Tz_n] = λ_n[ΣV^Tz_n]. Since U^TU = I, we can rewrite Eq. () as

U [- ı U^{T} H U] (U^{T} U) [Σ V^{T} z_{n}] = λ_{n} U [Σ V^{T} z_{n}],

(F10)

which is precisely () written in terms of Γ = UU^T with w_n = UΣV^Tz_n = ∇z_n.

We now establish Eq. (). From the formula u = ∇ ⋅ H in (), for the continuum setting, we have that u_j = (∇ ⋅ H) ⋅ e_j = ∇ ⋅ (He_j). Since ∇ = UΣV^T, the discrete version of this formula is given by

u_{j} = - \nabla^{T} H e_{j} = - V Σ U^{T} H e_{j} .

(F11)

From ∇ = UΣV^T and

{(\nabla^{T} \nabla)}^{- 1} = V Σ^{- 2} V^{T}

⁠, we have

\nabla {(\nabla^{T} \nabla)}^{- 1} = U Σ^{- 1} V^{T}

⁠. Consequently, Γ = UU^T, and Eqs. () and () yield

\nabla {(\nabla^{T} \nabla)}^{- 1} u_{j} = - Γ H e_{j}

⁠. Equation () now follows from the formula w_n = ∇z_n in (),

\begin{align} w_{n} w_{n}^{†} \nabla {(\nabla^{T} \nabla)}^{- 1} u_{j} \cdot \nabla {(\nabla^{T} \nabla)}^{- 1} u_{k} & = [{(\nabla^{T} \nabla)}^{- 1} \nabla^{T} w_{n}] {[{(\nabla^{T} \nabla)}^{- 1} \nabla^{T} w_{n}]}^{†} u_{j} \cdot u_{k} \\ = [{(\nabla^{T} \nabla)}^{- 1} \nabla^{T} \nabla z_{n}] {[{(\nabla^{T} \nabla)}^{- 1} \nabla^{T} \nabla z_{n}]}^{†} u_{j} \cdot u_{k} \\ = z_{n} z_{n}^{†} u_{j} \cdot u_{k}, \end{align}

(F12)

where we have used that the inverse of a symmetric matrix is also symmetric.²⁵ The equivalence of Eqs. () and () now follows from Eqs. (), (), and (). 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 = λ_nw_n, where the matrix −ıΓHΓ is of size N × N. In Sec. IV, $D^{*}$ was calculated in terms of the generalized eigenvalue problem −ı∇^TH∇z_n = λ_n∇^T∇z_n, involving the K × K matrices −ı∇^TH∇ and ∇^T∇. Since $\nabla^{T} = (\nabla_{1}^{T}, \dots, \nabla_{d}^{T})$ 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.⁴⁷ 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 = λ_nw_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₁ in Theorem 8 both denote the rank of the matrix Γ, i.e., K = N₁. Note, computing the matrix $Γ = \nabla {(\nabla^{T} \nabla)}^{- 1} \nabla^{T}$ 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 − ı∇^TH∇z_n = λ_n∇^T∇z_n can then be obtained by repeatedly computing the standard eigenvalue decomposition of the matrix −ıU^THU 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₁ satisfying K₁ < 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₁ × K₁. 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₁ satisfying K₁ < 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Σ²V^T is singular, with K₁ nonzero eigenvalues and K₀ = K − K₁ zero eigenvalues, and write

U = [U_{0} U_{1}], Σ = diag (O_{00}, Σ_{1}), V = [V_{0} V_{1}] .

(G1)

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 Σ₁ is a K₁ × K₁ diagonal, invertible matrix. By Eq. (F1), the matrices U₁ and V₁ satisfy $U_{1}^{T} U_{1} = I_{1}$ and $V_{1}^{T} V_{1} = I_{1}$ ⁠, where I₁ is the K₁ × K₁ identity matrix, but $V_{1} V_{1}^{T} \neq I$ ⁠. Due to the blocks of zeros in the matrix Σ in Eq. (G1), we can write the matrix gradient as $\nabla = U_{1} Σ_{1} V_{1}^{T}$ and the negative matrix Laplacian as $\nabla^{T} \nabla = V_{1} Σ_{1}^{2} V_{1}^{T}$ ⁠. An important property of the SVD of the matrix ∇ is that its null space is spanned by the columns of V₀ and its range is spanned by the columns of U₁.¹⁴ We emphasize that in the setting where ∇ is full-rank, we have U₁ = U, Σ₁ = Σ, and V₁ = V satisfying $V_{1} V_{1}^{T} = V_{1}^{T} V_{1} = 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.)

\nabla^{T} H \nabla χ_{j} + ε \nabla^{T} \nabla χ_{j} = u_{j},

(G2)

where u_j is the discrete, vector representation of the jth component of the fluid velocity field u_j and similarly for χ_j. Substituting the formula for u_j in (G2) into the discrete version $D_{j k}^{*} = ε δ_{j k} + ⟨ u_{j} \cdot χ_{k} ⟩$ of Eq. (3) yields

D_{j k}^{*} = S_{j k}^{*} + A_{j k}^{*}, S_{j k}^{*} = ε (δ_{j k} + ⟨ \nabla χ_{j} \cdot \nabla χ_{k} ⟩), A_{j k}^{*} = ⟨ \nabla^{T} H \nabla χ_{j} \cdot χ_{k} ⟩,

(G3)

where, as before, $S_{k j}^{*} = S_{j k}^{*}$ and $A_{k j}^{*} = - A_{j k}^{*}$ ⁠. We are now ready to state the key result of this appendix.

Theorem 10.

Let the matrix gradient ∇ be rank-deficient and let

\nabla = U_{1} Σ_{1} V_{1}^{T}

be its SVD. Also, let

U_{1}^{T} H U_{1} = ı R_{1} Λ_{1} R_{1}^{†}

be the spectral decomposition of the antisymmetric matrix

U_{1}^{T} H U_{1}

. Consider the discrete formulation of the effective parameter problem developed in 4Sec. IV. We have the following generalization of Eq. ():

\sum_{n = 1}^{K_{1}} Q_{n}^{1} = V_{1} V_{1}^{T}, Q_{n}^{1} = z_{n}^{1} {[\nabla z_{n}^{1}]}^{†} \nabla, Q_{l}^{1} Q_{m}^{1} = Q_{l}^{1} δ_{l m},

(G4)

where the matrices

Q_{n}^{1}

, n = 1, …, K₁, are self-adjoint with respect to the discrete inner-product ⟨⋅, ⋅⟩_1,2 defined by ⟨ξ,ζ⟩_1,2 = ⟨∇ξ ⋅ ∇ζ⟩, i.e.,

{⟨ Q_{n}^{1} ξ, ζ ⟩}_{1, 2} = {⟨ ξ, Q_{n}^{1} ζ ⟩}_{1, 2}

for

ξ, ζ \in C^{K_{1}}

. Moreover, the generalization of the resolvent formula χ_j = Z(ɛI + ıΛ)⁻¹Z^†u_j in Eq. () is given by

V_{1}^{T} χ_{j} = V_{1}^{T} Z_{1} {(ε I_{1} + ı Λ_{1})}^{- 1} Z_{1}^{†} u_{j}, Z_{1} = V_{1} Σ_{1}^{- 1} R_{1} .

(G5)

Now consider the discrete formulation of the effective parameter problem developed in Appendix E. Let $Γ_{1} H Γ_{1} = ı W_{1} \tilde{Λ} W_{1}^{†}$ be the spectral decomposition of the antisymmetric matrix Γ₁HΓ₁, where $Γ_{1} = U_{1} U_{1}^{T}$ , $\tilde{Λ}$ is a diagonal real-valued matrix with ${\tilde{λ}}_{n}^{1}$ , n = 1, …, N, on its diagonal, and W₁ is a unitary matrix with columns $w_{n}^{1}$ . Then, Eqs. (E1) and (E2) hold with $Q_{n}$ and λ_n replaced by $Q_{n}^{1} = [w_{n}^{1}] {[w_{n}^{1}]}^{†}$ and ${\tilde{λ}}_{n}^{1}$ . Also, the resolvent formula in Eq. (E4 holds with W replaced by W₁ and Λ replaced by $\tilde{Λ}$ .

These two discrete formulations of the effective parameter problem are related as follows: The diagonal eigenvalue matrices

\tilde{Λ}

and Λ₁ are related by

\tilde{Λ} = diag (O_{00}, Λ_{1})

. The eigenvector matrices W₁ and Z₁ are related by the following generalization of Eq. () [also see ()]:

W_{1} = [P_{0} \nabla Z_{1}], \nabla Z_{1} = U_{1} R_{1},

(G6)

where the columns of P₀ are eigenvectors of Γ₁ which span its null space. Moreover, the following formulas generalize Eqs. (), (), and ():

W_{1}^{†} Γ_{1} H e_{j} \cdot W_{1}^{†} Γ_{1} H e_{k} = {[\nabla Z_{1}]}^{†} H e_{j} \cdot {[\nabla Z_{1}]}^{†} H e_{k} = Z_{1}^{†} u_{j} \cdot Z_{1}^{†} u_{k} .

(G7)

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 $z_{n}^{1} = V_{1} Σ_{1}^{- 1} r_{n}^{1}$ and $λ_{n}^{1}$ , n = 1, …, K₁, where $(λ_{n}^{1}, r_{n}^{1})$ are eigenpairs of matrix $- ı U_{1}^{T} H U_{1}$ . 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.

Proof of Theorem 10.

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 $() ε I + Γ H Γ () \nabla χ_{j} = - Γ H e_{j}$ ⁠, 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₀₀, Σ₁), we have $\nabla = U_{1} Σ_{1} V_{1}^{T}$ and $\nabla^{T} \nabla = V_{1} Σ_{1}^{2} V_{1}^{T}$ ⁠, so the cell problem in (G2) can be written as $[V_{1} Σ_{1}] [U_{1}^{T} H U_{1}] [Σ_{1} V_{1}^{T}] χ_{j} + ε V_{1} Σ_{1}^{2} V_{1}^{T} χ_{j} = u_{j}$ ⁠. The K₁ × K₁ antisymmetric matrix $U_{1}^{T} H U_{1}$ has the spectral decomposition $U_{1}^{T} H U_{1} = ı R_{1} Λ_{1} R_{1}^{†}$ ⁠, where Λ₁ is a diagonal real-valued matrix and R₁ is a unitary matrix, $R_{1}^{†} R_{1} = R_{1} R_{1}^{†} = I_{1}$ ⁠. Equation (G5) follows from these formulas, which is an analog of Eq. (37). The formula $\nabla = U_{1} Σ_{1} V_{1}^{T}$ and Eq. (G5), in turn, imply that $\nabla χ_{j} = U_{1} R_{1} {(ε I_{1} + ı Λ_{1})}^{- 1} {[V_{1} Σ_{1}^{- 1} R_{1}]}^{†} u_{j}$ ⁠.

Substituting this formula for ∇χ_j into the formula for

S_{j k}^{*}

in Eq. () and using

U_{1}^{T} U_{1} = I_{1}

and

R_{1}^{†} R_{1} = I_{1}

yield the functional formula for

S_{j k}^{*}

in () with Λ replaced by Λ₁ and Z replaced by

Z_{1} = V_{1} Σ_{1}^{- 1} R_{1}

⁠. We also have

⟨ \nabla^{T} H \nabla χ_{j} \cdot χ_{k} ⟩ = ⟨ ı [V_{1} Σ_{1} R_{1}] Λ_{1} [R_{1}^{†} Σ_{1} V_{1}^{T}] χ_{j} \cdot χ_{k} ⟩ = ⟨ ı [Σ_{1} R_{1} Λ_{1} R_{1}^{†} Σ_{1}] V_{1}^{T} χ_{j} \cdot V_{1}^{T} χ_{k} ⟩ .

(G8)

This formula,

R_{1}^{†} R_{1} = I_{1}

⁠,

Σ_{1}^{T} = Σ_{1}

⁠, and Eqs. () and () yield the functional formula for

A_{j k}^{*}

in () with Λ replaced by Λ₁ and Z replaced by

Z_{1} = V_{1} Σ_{1}^{- 1} R_{1}

⁠. The quadratic form

Z_{1}^{†} u_{j} \cdot Z_{1}^{†} u_{k}

arising in these functional formulas yields Eq. () with z_n replaced by

z_{n}^{1} = V_{1} Σ_{1}^{- 1} r_{n}^{1}

and other appropriate notational changes, where

r_{n}^{1}

⁠, n = 1, …, K₁, are the orthonormal eigenvectors of the matrix

U_{1}^{T} H U_{1}

which comprise the columns of R₁. The formula

z_{n}^{1} = V_{1} Σ_{1}^{- 1} r_{n}^{1}

can be written as

\nabla z_{n}^{1} = U_{1} r_{n}^{1}

⁠. The orthogonality properties of R₁ and U₁ then imply that the vectors

z_{n}^{1}

satisfy the Sobolev-type orthogonality condition in Eq. (),

\nabla z_{n}^{1} \cdot \nabla z_{m}^{1} = δ_{n m}

⁠. Moreover, since

\sum_{n} r_{n}^{1} {[r_{n}^{1}]}^{†} = I_{1}

⁠, we also have Eq. (). Consequently, the functional representations of

S_{j k}^{*}

and

A_{j k}^{*}

in Eq. () have the discrete integral representations in Eq. () with z_n and λ_n replaced by

z_{n}^{1}

and

λ_{n}^{1}

⁠, and

(λ_{n}^{1}, r_{n}^{1})

are eigenpairs of Hermitian matrix

- ı U_{1}^{T} H U_{1}

of size K₁. These summations have the same properties as the summations discussed in the Proof of Theorem 3 and Eqs. () and ().

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₁ = V is orthogonal, Σ₁ = Σ is invertible, and R₁ = R is orthogonal so that the matrix Z₁ = Z defined in (G5) is given by Z = VΣ⁻¹R and is invertible with Z⁻¹ = R^†ΣV^T. Equation (33) is satisfied with Z^−† = VΣR and Λ = Λ₁ 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 $\nabla = U_{1} Σ_{1} V_{1}^{T}$ ⁠, $\nabla^{T} \nabla = V_{1} Σ_{1}^{2} V_{1}^{T}$ ⁠, $U_{1}^{T} U_{1} = V_{1}^{T} V_{1} = I_{1}$ ⁠, and the invertibility of the matrix Σ₁, the cell problem in (G2) can be written as $U_{1} [U_{1}^{T} H U_{1}] [U_{1}^{T} U_{1}] [Σ_{1} V_{1}^{T}] χ_{j} + ε U_{1} Σ_{1} V_{1}^{T} χ_{j} = U_{1} Σ_{1}^{- 1} V_{1}^{T} u_{j}$ ⁠. Substituting u_j = −∇^THe_j into this expression yields $() ε I + Γ_{1} H Γ_{1} () \nabla χ_{j} = g_{j}^{1}$ ⁠, where $Γ_{1} = U_{1} U_{1}^{T}$ and $g_{j}^{1} = - Γ_{1} H e_{j}$ which is analogous to Eq. (D9). As in Appendix F, the matrix $Γ_{1} = U_{1} U_{1}^{T}$ projects subspaces of $R^{N}$ onto the range of ∇. The antisymmetric matrix Γ₁HΓ₁ has the spectral decomposition $Γ_{1} H Γ_{1} = ı W_{1} \tilde{Λ} W_{1}^{†}$ ⁠, where $\tilde{Λ}$ is a diagonal real-valued matrix and W₁ is a unitary matrix $W_{1}^{†} W_{1} = W_{1} W_{1}^{†} = I$ ⁠, which yields a generalization of (E4). From Γ₁∇ = ∇, we have ⟨∇^TH∇χ_j ⋅ χ_k⟩ = ⟨Γ₁HΓ₁∇χ_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₁ = U hence W₁ = 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₀ P₁], where the columns of the matrices P₀ and P₁ are orthonormal eigenvectors that span the null space and range of Γ, respectively. Since the eigenvalues γ_n associated with the eigenvectors in the matrix P₁ 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₁ span the range of $Γ_{1} = U_{1} U_{1}^{T}$ ⁠, we may take P₁ = U₁ so that P = [P₀ U₁]. Consequently, we can rewrite Eq. (E11) as $U_{1}^{T} H U_{1} = ı R_{11} Λ_{11} R_{11}^{†}$ ⁠. Identifying the notations of Appendix E with this section, we have R₁₁ = R₁ and Λ₁₁ = Λ₁. This equation and Eq. (E12) establish that $\tilde{Λ} = diag (O_{00}, Λ_{1})$ ⁠.

We now establish that the spectral weights associated with each approach are identical. Using the notation of this section, Eq. (E12) yields W₁ = [P₀ U₁R₁]. Consequently, from $\nabla = U_{1} Σ_{1} V_{1}^{T}$ ⁠, $U_{1}^{T} U_{1} = V_{1}^{T} V_{1} = I_{1}$ ⁠, and the formula $Z_{1} = V_{1} Σ_{1}^{- 1} R_{1}$ in (G5), we have that ∇Z₁ = U₁R₁, implying Eq. (G6), which is a generalization of Eq. (F5). Equation (G7) now follows from $Γ_{1}^{T} = Γ_{1}$ ⁠, Γ₁P₀ = O, Γ₁U₁ = U₁, and the formula u_j = −∇^THe_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 $z_{n}^{1} = V_{1} Σ_{1}^{- 1} r_{n}^{1}$ and $λ_{n}^{1}$ ⁠, respectively, where $(λ_{n}^{1}, r_{n}^{1})$ are eigenpairs of matrix $- ı U_{1}^{T} H U_{1}$ ⁠. We emphasize that this matrix is of size K₁ which is more than a factor of d smaller than the matrix −ıΓ₁HΓ₁. 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 $- ı U_{1}^{T} H U_{1}$ ⁠, which is a less costly numerical computation than the computation associated with the generalized eigenvalue problem⁴⁷ 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.

REFERENCES

1.

N. I.

Akhiezer

,

The Classical Moment Problem

(

Oliver & Boyd

,

1965

).

Google Scholar

2.

G. S.

Aslanyan

,

I. L.

Maikov

, and

I. Z.

Filimonova

, “

Simulation of pulverized coal combustion in a turbulent flow

,”

Combust., Explos. Shock Waves

30

(

4

),

448

–

453

(

1994

).

https://doi.org/10.1007/bf00790149

Google Scholar

Crossref

3.

M.

Avellaneda

and

A.

Majda

, “

Stieltjes integral representation and effective diffusivity bounds for turbulent transport

,”

Phys. Rev. Lett.

62

,

753

–

755

(

1989

).

https://doi.org/10.1103/physrevlett.62.753

Google Scholar

Crossref

PubMed

4.

M.

Avellaneda

and

A.

Majda

, “

An integral representation and bounds on the effective diffusivity in passive advection by laminar and turbulent flows

,”

Commun. Math. Phys.

138

,

339

–

391

(

1991

).

https://doi.org/10.1007/bf02099496

Google Scholar

Crossref

5.

M.

Avellaneda

and

M.

Vergassola

, “

Stieltjes integral representation of effective diffusivities in time-dependent flows

,”

Phys. Rev. E

52

(

3

),

3249

–

3251

(

1995

).

https://doi.org/10.1103/physreve.52.3249

Google Scholar

Crossref

6.

G. A.

Baker

and

P. R.

Graves-Morris

,

Padé Approximants

, Encyclopedia of Mathematics and Its Applications (

Cambridge University Press

,

1996

).

Google Scholar

Crossref

7.

A.

Bensoussan

,

J.-L.

Lions

, and

G.

Papanicolaou

,

Asymptotic Analysis for Periodic Structures

(

North-Holland

,

Amsterdam

,

1978

).

Google Scholar

8.

D. J.

Bergman

, “

Exactly solvable microscopic geometries and rigorous bounds for the complex dielectric constant of a two-component composite material

,”

Phys. Rev. Lett.

44

,

1285

–

1287

(

1980

).

https://doi.org/10.1103/physrevlett.44.1285

Google Scholar

Crossref

9.

R.

Bhattacharya

, “

Multiscale diffusion processes with periodic coefficients and an application to solute transport in porous media

,”

Ann. Appl. Probab.

9

(

4

),

951

–

1020

(

1999

).

https://doi.org/10.1214/aoap/1029962863

Google Scholar

Crossref

10.

R. N.

Bhattacharya

,

V. K.

Gupta

, and

H. F.

Walker

, “

Asymptotics of solute dispersion in periodic porous media

,”

SIAM J. Appl. Math.

49

(

1

),

86

–

98

(

1989

).

https://doi.org/10.1137/0149005

Google Scholar

Crossref

11.

L.

Biferale

,

A.

Crisanti

,

M.

Vergassola

, and

A.

Vulpiani

, “

Eddy diffusivities in scalar transport

,”

Phys. Fluids

7

,

2725

–

2734

(

1995

).

https://doi.org/10.1063/1.868651

Google Scholar

Crossref

12.

R. W.

Bilger

,

S. B.

Pope

,

K. N. C.

Bray

, and

J. F.

Driscoll

, “

Paradigms in turbulent combustion research

,”

Proc. Combust. Inst.

30

,

21

–

42

(

2005

).

https://doi.org/10.1016/j.proci.2004.08.273

Google Scholar

Crossref

13.

J.

Bonn

and

R. M.

McLaughlin

, “

Sensitive enhanced diffusivities for flows with fluctuating mean winds: A two-parameter study

,”

J. Fluid Mech.

445

,

345

–

375

(

2001

).

https://doi.org/10.1017/s002211200100564x

Google Scholar

Crossref

14.

J. W.

Demmel

,

Applied Numerical Linear Algebra

(

SIAM

,

1997

).

Google Scholar

Crossref

15.

E.

Di Lorenzo

,

D.

Mountain

,

H. P.

Batchelder

,

N.

Bond

, and

E. E.

Hofmann

, “

Advances in marine ecosystem dynamics from US GLOBEC: The horizontal-advection bottom-up forcing paradigm

,”

Oceanography

26

(

4

),

22

–

33

(

2013

).

https://doi.org/10.5670/oceanog.2013.73

Google Scholar

Crossref

16.

R.

Durrett

,

Probability: Theory and Examples

, 4th ed. (

Cambridge University Press

,

2010

).

Google Scholar

Crossref

17.

A.

Fannjiang

and

G.

Papanicolaou

, “

Convection enhanced diffusion for periodic flows

,”

SIAM J. Appl. Math.

54

(

2

),

333

–

408

(

1994

).

https://doi.org/10.1137/s0036139992236785

Google Scholar

Crossref

18.

A.

Fannjiang

and

G.

Papanicolaou

, “

Convection-enhanced diffusion for random flows

,”

J. Stat. Phys.

88

(

5-6

),

1033

–

1076

(

1997

).

https://doi.org/10.1007/bf02732425

Google Scholar

Crossref

19.

G. B.

Folland

,

Introduction to Partial Differential Equations

(

Princeton University Press

,

Princeton, NJ

,

1995

).

Google Scholar

20.

G. B.

Folland

,

Real Analysis: Modern Techniques and Their Applications

(

Wiley-Interscience

,

New York, NY

,

1999

).

Google Scholar

21.

K. M.

Golden

,

H.

Eicken

,

A. L.

Heaton

,

J.

Miner

,

D.

Pringle

, and

J.

Zhu

, “

Thermal evolution of permeability and microstructure in sea ice

,”

Geophys. Res. Lett.

34

,

L16501

,

https://doi.org/10.1029/2007gl030447

(

2007

).

https://doi.org/10.1029/2007gl030447

Google Scholar

Crossref

22.

K. M.

Golden

and

G.

Papanicolaou

, “

Bounds for effective parameters of heterogeneous media by analytic continuation

,”

Commun. Math. Phys.

90

,

473

–

491

(

1983

).

https://doi.org/10.1007/bf01216179

Google Scholar

Crossref

23.

Y.

Gorb

,

D.

Nam

, and

A.

Novikov

, “

Numerical simulations of diffusion in cellular flows at high Péclet numbers

,”

Discrete Contin. Dyn. Syst., Ser. B

15

(

1

),

75

–

92

(

2011

).

https://doi.org/10.3934/dcdsb.2011.15.75

Google Scholar

Crossref

24.

E. E.

Hofmann

and

E. J.

Murphy

, “

Advection, krill, and Antarctic marine ecosystems

,”

Antarct. Sci.

16

(

04

),

487

–

499

(

2004

).

https://doi.org/10.1017/s0954102004002275

Google Scholar

Crossref

25.

R. A.

Horn

and

C. R.

Johnson

,

Matrix Analysis

(

Cambridge University Press

,

1990

).

Google Scholar

26.

M.

Isichenko

and

J.

Kalda

, “

Statistical topography. II. Two-dimensional transport of passive scalar

,”

J. Nonlinear Sci.

1

(

4

),

375

–

396

(

1991

).

https://doi.org/10.1007/bf02429846

Google Scholar

Crossref

27.

J. D.

Jackson

,

Classical Electrodynamics

(

John Wiley & Sons

,

New York

,

1999

).

Google Scholar

28.

J. P.

Keener

,

Principles of Applied Mathematics: Transformation and Approximation

, Advanced Book Program (

Westview Press

,

Cambridge, MA

,

2000

).

Google Scholar

29.

D.

Khoshnevisan

,

Probability

, Graduate Studies in Mathematics (

American Mathematical Society

,

2007

).

Google Scholar

30.

N.

Kraitzman

,

R.

Hardenbrook

,

N. B.

Murphy

,

E.

Cherkaev

,

J.

Zhu

, and

K. M.

Golden

, “

Bounds on convection enhanced thermal transport in sea ice

” (unpublished).

31.

J. V.

Lukovich

,

D. G.

Babb

, and

D. G.

Barber

, “

On the scaling laws derived from ice beacon trajectories in the southern Beaufort Sea during the International Polar Year—Circumpolar Flaw Lead study, 2007–2008

,”

J. Geophys. Res.: Oceans

116

(

C9

),

C00G07

,

https://doi.org/10.1029/2011jc007049

(

2011

).

https://doi.org/10.1029/2011jc007049

Google Scholar

Crossref

32.

V. I.

Lytle

and

S. F.

Ackley

, “

Heat flux through sea ice in the Western Weddell Sea: Convective and conductive transfer processes

,”

J. Geophys. Res.

101

(

C4

),

8853

–

8868

,

https://doi.org/10.1029/95jc03675

(

1996

).

https://doi.org/10.1029/95jc03675

Google Scholar

Crossref

33.

A.

Majda

and

P. R.

Kramer

,

Simplified Models for Turbulent Diffusion: Theory, Numerical Modelling, and Physical Phenomena

, Physics Reports (

North-Holland

,

1999

).

Google Scholar

34.

A. J.

Majda

and

P. E.

Souganidis

, “

Large scale front dynamics for turbulent reaction-diffusion equations with separated velocity scales

,”

Nonlinearity

7

(

1

),

1

–

30

(

1994

).

https://doi.org/10.1088/0951-7715/7/1/001

Google Scholar

Crossref

35.

D.

McLaughlin

,

G.

Papanicolaou

, and

O.

Pironneau

, “

Convection of microstructure and related problems

,”

SIAM J. Appl. Math.

45

,

780

–

797

(

1985

).

https://doi.org/10.1137/0145046

Google Scholar

Crossref

36.

R. M.

McLaughlin

and

M. G.

Forest

, “

An anelastic, scale-separated model for mixing, with application to atmospheric transport phenomena

,”

Phys. Fluids

11

(

4

),

880

–

892

(

1999

).

https://doi.org/10.1063/1.869967

Google Scholar

Crossref

37.

R. C.

McOwen

,

Partial Differential Equations: Methods and Applications

(

Prentice Hall PTR

,

2003

).

Google Scholar

38.

Y. A.

Melnikov

and

M. Y.

Melnikov

,

Green’s Functions: Construction and Applications

, De Gruyter Studies in Mathematics (

De Gruyter

,

2012

).

Google Scholar

Crossref

39.

G. W.

Milton

, “

Bounds on the complex dielectric constant of a composite material

,”

Appl. Phys. Lett.

37

,

300

–

302

(

1980

).

https://doi.org/10.1063/1.91895

Google Scholar

Crossref

40.

G. W.

Milton

,

Theory of Composites

(

Cambridge University Press

,

Cambridge

,

2002

).

Google Scholar

Crossref

41.

H. K.

Moffatt

, “

Transport effects associated with turbulence with particular attention to the influence of helicity

,”

Rep. Prog. Phys.

46

(

5

),

621

–

664

(

1983

).

https://doi.org/10.1088/0034-4885/46/5/002

Google Scholar

Crossref

42.

N. B.

Murphy

,

E.

Cherkaev

,

C.

Hohenegger

, and

K. M.

Golden

, “

Spectral measure computations for composite materials

,”

Commun. Math. Sci.

13

(

4

),

825

–

862

(

2015

).

https://doi.org/10.4310/cms.2015.v13.n4.a1

Google Scholar

Crossref

43.

N. B.

Murphy

,

E.

Cherkaev

,

J.

Xin

, and

K. M.

Golden

, “

Spectral measures and iterative bounds for effective diffusivity for periodic flows

” (unpublished).

44.

N. B.

Murphy

,

E.

Cherkaev

,

J.

Xin

,

J.

Zhu

, and

K. M.

Golden

, “

Spectral analysis and computation of effective diffusivities in space-time periodic incompressible flows

,”

Ann. Math. Sci. Appl.

2

,

3

–

66

(

2017

).

https://doi.org/10.4310/amsa.2017.v2.n1.a1

Google Scholar

Crossref

45.

N.

Benjamin Murphy

,

E.

Cherkaev

, and

K. M.

Golden

, “

Anderson transition for classical transport in composite materials

,”

Phys. Rev. Lett.

118

,

036401

(

2017

).

https://doi.org/10.1103/physrevlett.118.036401

Google Scholar

Crossref

PubMed

46.

G.

Papanicolaou

and

S.

Varadhan

, “

Boundary value problems with rapidly oscillating coefficients

,” in

Colloquia Mathematica Societatis János Bolyai

, Random Fields (Esztergom, Hungary 1979) (

North-Holland

,

1982

), Vol. 27, pp.

835

–

873

.

Google Scholar

47.

B. N.

Parlett

,

The Symmetric Eigenvalue Problem

(

Prentice-Hall

,

Upper Saddle River, NJ, USA

,

1998

).

Google Scholar

Crossref

48.

G. A.

Pavliotis

, “

Homogenization theory for advection-diffusion equations with mean flow

,” Ph.D. thesis,

Rensselaer Polytechnic Institute Troy

,

New York

,

2002

.

Google Scholar

49.

G. A.

Pavliotis

, “

Asymptotic analysis of the Green–Kubo formula

,”

IMA J. Appl. Math.

75

,

951

–

967

(

2010

).

https://doi.org/10.1093/imamat/hxq039

Google Scholar

Crossref

50.

N.

Peters

,

Turbulent Combustion

, Cambridge Monographs on Mechanics (

Cambridge University Press

,

2000

).

Google Scholar

Crossref

51.

P.

Rampal

,

J.

Weiss

,

D.

Marsan

, and

M.

Bourgoin

, “

Arctic sea ice velocity field: General circulation and turbulent-like fluctuations

,”

J. Geophys. Res.: Oceans

114

(

C10

),

C10014

,

https://doi.org/10.1029/2008jc005227

(

2009

).

https://doi.org/10.1029/2008jc005227

Google Scholar

Crossref

52.

P.

Rampal

,

J.

Weiss

,

D.

Marsan

,

R.

Lindsay

, and

H.

Stern

, “

Scaling properties of sea ice deformation from buoy dispersion analysis

,”

J. Geophys. Res.: Oceans

113

(

C3

),

C03002

,

https://doi.org/10.1029/2007jc004143

(

2008

).

https://doi.org/10.1029/2007jc004143

Google Scholar

Crossref

53.

M. C.

Reed

and

B.

Simon

,

Functional Analysis

(

Academic Press

,

San Diego, CA

,

1980

).

Google Scholar

54.

W.

Rudin

,

Real and Complex Analysis

(

McGraw-Hill

,

New York, NY

,

1987

).

Google Scholar

55.

P. J.

Samson

, “

Atmospheric transport and dispersion of air pollutants associated with vehicular emissions

,” in

Air Pollution, The Automobile, and Public Health

, edited by

A. Y.

Watson

,

R. R.

Bates

, and

D.

Kennedy

(

National Academy Press

,

US

,

1988

), pp.

77

–

97

.

Google Scholar

PubMed

56.

I.

Stakgold

,

Boundary Value Problems of Mathematical Physics

, Classics in Applied Mathematics (

SIAM

,

2000

), Vol. 2, set.

Google Scholar

57.

T.-J.

Stieltjes

, “

Recherches sur les fractions continues

,”

Ann. Fact. Sci. Toulouse

,

4

(

1

),

J1

–

J35

(

1995

).

https://doi.org/10.5802/afst.789

Google Scholar

Crossref

58.

M. H.

Stone

,

Linear Transformations in Hilbert Space

(

American Mathematical Society

,

Providence, RI

,

1964

).

Google Scholar

59.

R. J.

Tabaczynski

, “

Turbulent flows in reciprocating internal combustion engines

,” in

Internal Combustion Engineering: Science and Technology

, edited by

J. H.

Weaving

(

Springer Netherlands

,

1990

), pp.

243

–

285

.

Google Scholar

Crossref

60.

G. I.

Taylor

, “

Eddy motion in the atmosphere

,”

Philos. Trans. Royal Soc. A

215

(

523-537

),

1

–

26

(

1915

).

https://doi.org/10.1098/rsta.1915.0001

Google Scholar

Crossref

61.

G. I.

Taylor

, “

Diffusion by continuous movements

,”

Proc. London Math. Soc.

s2-20

,

196

–

212

(

1922

).

https://doi.org/10.1112/plms/s2-20.1.196

Google Scholar

Crossref

62.

H. J.

Trodahl

,

S. O. F.

Wilkinson

,

M. J.

McGuinness

, and

T. G.

Haskell

, “

Thermal conductivity of sea ice; dependence on temperature and depth

,”

Geophys. Res. Lett.

28

(

7

),

1279

–

1282

,

https://doi.org/10.1029/2000gl012088

(

2001

).

https://doi.org/10.1029/2000gl012088

Google Scholar

Crossref

63.

W. M.

Washington

and

C. L.

Parkinson

,

An Introduction to Three-Dimensional Climate Modeling

(

University Science Books

,

1986

).

Google Scholar

64.

E.

Watanabe

and

H.

Hasumi

, “

Pacific water transport in the western Arctic Ocean simulated by an eddy-resolving coupled sea ice–ocean model

,”

J. Phys. Oceanogr.

39

(

9

),

2194

–

2211

(

2009

).

https://doi.org/10.1175/2009JPO4010.1

Google Scholar

Crossref

65.

F. A.

Williams

,

Turbulent Combustion

, Frontiers in Applied Mathematics (

Society for Industrial and Applied Mathematics (SIAM)

,

Philadelphia, PA

,

1985

), Chap. 3, pp.

97

–

131

.

Google Scholar

66.

J.

Xin

,

An Introduction to Fronts in Random Media

, Surveys and Tutorials in the Applied Mathematical Sciences (

Springer

,

New York

,

2009

).

Google Scholar

Crossref

2020

Author(s)

Spectral analysis and computation for homogenization of advection diffusion processes in steady flows

I. INTRODUCTION

II. HOMOGENIZATION OF THE ADVECTION DIFFUSION EQUATION

III. HILBERT SPACE AND INTEGRAL REPRESENTATIONS

IV. DISCRETE SETTING: SOBOLEV SPACE OF SCALAR FIELDS

V. SPECTRAL COMPUTATIONS OF THE EFFECTIVE DIFFUSIVITY MATRIX

A. Numerical methods

B. Numerical results

1. BC-shear flow

2. Cat’s eye flow

3. Randomly perturbed cat’s eye flow

VI. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: APPENDIX OVERVIEW

APPENDIX B: NOTATION

APPENDIX C: PROPERTIES OF THE LINEAR OPERATOR A

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

1. Functional formulas for the effective diffusivity matrix

2. The analytic continuation method and integral representations of $D^{*}$

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

APPENDIX F: DISCRETE EQUIVALENCE OF THE EFFECTIVE PARAMETER PROBLEMS

APPENDIX G: RANK DEFICIENCY AND A UNIFYING STANDARD EIGENVALUE PROBLEM

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Spectral analysis and computation for homogenization of advection diffusion processes in steady flows

I. INTRODUCTION

II. HOMOGENIZATION OF THE ADVECTION DIFFUSION EQUATION

III. HILBERT SPACE AND INTEGRAL REPRESENTATIONS

IV. DISCRETE SETTING: SOBOLEV SPACE OF SCALAR FIELDS

V. SPECTRAL COMPUTATIONS OF THE EFFECTIVE DIFFUSIVITY MATRIX

A. Numerical methods

B. Numerical results

1. BC-shear flow

2. Cat’s eye flow

3. Randomly perturbed cat’s eye flow

VI. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: APPENDIX OVERVIEW

APPENDIX B: NOTATION

APPENDIX C: PROPERTIES OF THE LINEAR OPERATOR A

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

1. Functional formulas for the effective diffusivity matrix

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

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

APPENDIX F: DISCRETE EQUIVALENCE OF THE EFFECTIVE PARAMETER PROBLEMS

APPENDIX G: RANK DEFICIENCY AND A UNIFYING STANDARD EIGENVALUE PROBLEM

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

2. The analytic continuation method and integral representations of $D^{*}$