EMBED IN ENSEMBLE TO RIGOROUSLY AND ACCURATELY HOMOGENIZE QUASI-PERIODIC MULTI-SCALE HETEROGENEOUS MATERIAL

Abstract

For microscale heterogeneous partial differential equations (PDEs), this article further develops novel theory and methodology for their macroscale mathematical/asymptotic homogenization. This article specifically encompasses the case of quasi-periodic heterogeneity with finite scale separation: no scale separation limit is required. A key innovation herein is to analyse the ensemble of all phase-shifts of the heterogeneity. Dynamical systems theory then frames the homogenization as a slow manifold of the ensemble. Depending upon any perceived scale separation within the quasi-periodic heterogeneity, the homogenization may be done in either one step or two sequential steps: the results are equivalent. The theory not only assures us of the existence and emergence of an exact homogenization at finite scale separation, it also provides a practical systematic method to construct the homogenization to any specified order. For a class of heterogeneities, we show that the macroscale homogenization is potentially valid down to lengths which are just twice that of the microscale heterogeneity! This methodology complements existing well-established results by providing a new rigorous and flexible approach to homogenization that potentially also provides correct macroscale initial and boundary conditions, treatment of forcing and control, and analysis of uncertainty.

1 Introduction

Many engineering and scientific structures have important microscale variations, such as lattice materials [6, 62, 63], the windings in electrical machinery [58], electromagnetism in micro-structured material [20, 44], food drying [70] and slow Stokes flow through porous media [11]. The engineering challenge is to predict and understand the dynamics on a system scale significantly larger than the microscale details. Homogenization, via the recognition of multiple physical length scales, is the common approach [1, 17, 26, 30]. This article further develops an alternative approach to homogenized modelling of the large scale dynamics via a newly established general rigorous theory [52, Section 2.5]. This new complementary approach provides a route to systematic refinements of the homogenization for a finite separation of scales and to novel quantification of the remainder error. The significance of this work is not that it effectively reproduces some classic homogenization results, but that it invokes a powerful dynamical systems framework [12, 16, 19, 32, 42, 69] which provides complementary new insights and powerful practical extensions of previously known aspects of homogenization. Consequently, a background in modern dynamical systems [64, 65] is useful to appreciate the direction developed here for homogenization.

Consider materials with complicated microstructure: we want to model their large-scale dynamics by equations with effective ‘average’ coefficients. Heterogeneous diffusion in one-dimensional (1D) space is the simplest such class of systems: suppose the material, in spatial domain $[0,L]$ , has structure so that ‘material’ or ‘heat’ $u(t,x)$ diffuses according to

(1.1)

where the heterogeneous diffusion coefficient $\kappa (x)$ has some complicated multiscale structure. For definiteness, let us initially suppose the macroscale boundary conditions are simply that u is L-periodic in space x (the end of Section 2 discusses other boundary conditions). Our challenge is to derive, support and characterize the macroscale model PDE

(1.2)

for some effective mean field $U(x,t)$ , some effective macroscale diffusivity K, with some higher order corrections, and possibly an error estimate—both the latter represented by the ellipsis “ $\ldots $ ” in (1.2). The problem we address is not to find approximate solutions to PDE (1.1), but instead to provide new justification, construction and interpretation of the simpler PDE (1.2), and cognate higher-order PDEs, as purely macroscale models of the multiscale PDE (1.1). I also comment throughout on the closely connected problem [1, 3, 18, 25, 61] of wave propagation through heterogeneous material,

(1.3)

and this corresponding effective homogenized PDE model that predicts its macroscale dispersive waves. The homogenization approach for both classes are much the same, so the results here apply to elastic materials and other waves, with some caveats. However, the treatment and theory is clearer and more rigorous for the class of diffusive systems (1.1).

Here we focus on cases where the heterogeneity $\kappa (x)$ is quasi-periodic in space. For an example application, Sections 4 and 5 focus on homogenization for the family of heterogeneous diffusivities $\kappa =1/[1+a_1\cos (2\pi x/\ell _1)+a_2\cos (2\pi x/\ell _2)]$ for arbitrary microscale periods $\ell _1,\ell _2$ with $\ell _1>\ell _2$ . Sections 2 and 3 develop theory and supporting characteristics of homogenization in the general case where the heterogeneity $\kappa $ is expressible as a function of two variables, say $\kappa (x_1,x_2)$ , in which $\kappa $ is $\ell _i$ -periodic in variable $x_i$ . Then the quasi-periodic heterogeneity in (1.1) is $\kappa (x,x)$ . I emphasize that this work does not address the PDE $u_t=\partial _x\{\kappa (x/\epsilon ,x/\epsilon )u_x\}$ : there is no  $\epsilon $ in the analysis here; and consequently, and in contrast to extant well established theory of quasi-periodic homogenization (for example, see [33], [2, Section 7.4]), the framework here does not invoke the infinite length scale separation limit “ $\epsilon \to 0$ ”. Instead, the results herein apply for finite scale separation (quantified by Section 5), and form arbitrarily high-order homogenizations, not just leading order [33] nor just second-order [57].

Clearly, the approach developed here straightforwardly extends to an arbitrary number of incommensurable periods in the structure of the diffusivity $\kappa $ , but for simplicity of introducing the key new ideas, techniques and results, we address here the case of just two periods in $\kappa $ .

Most extant conventional approaches to practically computing a homogenization depend upon the identification of a representative volume element (RVE) [62], or a representative cell. Even in the simple case of periodic media, an uncertain issue is an appropriate size of an RVE [34, 38, 43, 59, 60, 63, p. 2, Section 5.1, Section 2.3, Section 1.2.2, p. 2, p. 13, respectively]. A further issue with such RVE-based arguments is that for quasi-periodic heterogeneities, there is no one RVE [5, Section 2.1]. When the ratio between microscale lengths $\ell _1,\ell _2$ is irrational, then all volume elements, no matter what size is chosen, are different: no volume element is representative. Relatedly, when the microscale lengths $\ell _1,\ell _2$ have a large least common multiple, then any true RVE must be correspondingly large—resulting in onerous analysis or computation. For example, Blanc and Le Bris [9] comment that for quasi-periodic homogenization, “cell problems must be in general set on the whole space”: herein, we show that this “must” is incorrect for our approach. In contrast to “the whole space”, here, Sections 2 and 3 instead develop a novel rigorous systematic approach that uses an ‘RVE’ cell which is the size of the largest microscale length, $\max \{\ell _1,\ell _2\}$ , albeit an RVE of more ‘spatial’ dimensions.

I expect straightforward extensions of the approach developed herein to empower cognate analysis and modelling of the following:

• • cases where the heterogeneity is quasi-periodic in any number of microscale lengths $\ell _1>\ell _2>\ell _3>\cdots $ ;

• • PDEs (1.1) where the field u is itself in a vector space such as the cross-sectional structure of either shear dispersion [40, 41], or of elastic beams [49];

• • multiple macroscale space dimensions as in the modelling of plates and shells [55];

• • PDEs like (1.1) but with macroscale variations in the material parameters as well as microscale variations [5]—sometimes called functionally graded materials;

• • PDEs like (1.1) but with nonlinear effects [38]—in contrast to the experience of Blanc and Le Bris [9], in the framework developed here, nonlinear problems do not involve any difficulty other than extra algebraic details.

Further, Section 6 shows that the approach, analysis and results developed herein also include three-scale multiscale systems where notionally $\ell _2\ll \ell _1\ll L$ , that is, where there is a microscale, a mesoscale and a macroscale [21, 45].

Other analytic techniques such as mathematical homogenization or the method of multiple scales [18, 61] provide a mechanism for deriving effective PDE models such as (1.2). However, such previous underlying theory requires the limit of infinite scale separation between the microscale “cell”-sizes $\ell _i$ and the macroscale domain size L by invoking $\epsilon :=\ell _i/L\to 0$ (for example, [2, 6, 70, p. 5, p. 6 and Section 2.4, respectively] and [4]). In contrast, and arising from developments in dynamical systems, Sections 2 and 3 respectively use phase-shift embedding and slowly varying theory [52, 55] to provide rigorous support to the modelling at finite scale separation: that is, when the macroscale length scale L of interest is larger than the microscales $\ell _1$ and $\ell _2$ , but not infinitely so. In the approach here, not only is there no limit, there is no need for such a defined “ $\epsilon $ ”. Such finite scale separation empowers this approach to capture the dispersion in wave homogenization which Abdulle and Pouchon [1, p. 3] call for in their discussion stating that “new effective models are required”. The phase-shift embedding developed here further illustrates the general principle that macroscale models of microscale structures are often best phrased as ensemble averages [8, 23, 68].

It is well recognized that homogenization is physically effective for reasonably well separated length scales: that the resolved macroscale structures, with length $\ell _M$ say, satisfy $\ell _M\gg \ell _1$ . But what does this mean? Much extant homogenization theory requires the mathematical limit $\ell _1/\ell _M\to 0$ . However, physically, the ratio $\ell _1/\ell _M$ is always finite. By setting homogenization in a new rigorous framework, Section 5 analyses homogenization for finite $\ell _1/\ell _M$ to high-order and presents evidence that, in a class of problems, a quantitative limit to the spatial resolution of homogenization is that $\ell _1/\ell _M\lesssim 1/2$ , equivalently, $\ell _M\gtrsim 2\ell _1$ . Mathematically, a wider scale separation is not necessary for homogenization.

Of course, practically, we usually prefer a wider scale separation because of many complicating factors in construction and use of practical homogenization models. However, in the sense discussed by Section 5, this quantitative, near-sharp, lower bound on the required scale separation indicates how little scale separation is necessary. Indeed, because it applies at finite scale separation, the analysis here encompasses high-contrast media and so also detects features occurring in this particularly interesting class of problems [17, p. 2].

The new approach to homogenization that we develop here creates a powerful dynamical systems framework for addressing other modelling issues in the future. For a time-dependent simulation of homogenized dynamics, one must provide some initial conditions for the homogenized PDE (1.2). These initial conditions are surprisingly nontrivial at finite scale separation, but the new framework developed here comes with well-established methods to derive the correct initial conditions for long-term forecasts [19, 47, 49, 51]. The projection that provides the correct initial conditions, also provides the correct projection of any applied “forcing” [47, Section 7], including the cases of control and system uncertainty. System uncertainty is a challenge for other techniques [4, 38, p. 3 and Section 8, respectively]. Further, the homogenized PDE (1.2) needs macroscale boundary conditions—a “highly challenging problem” according to Shahraki and Guzina [61], and which Cornaggia and Guzina [18] called a “complex problem that is still an active research topic”. Such boundary conditions are nontrivial and it is well recognized that consideration of boundary layers must be accounted for. However, again the new framework developed here comes with a rigorous approach to derive correct boundary conditions for (1.2) to empower accurate predictions at finite scale separation [14–16, 48, 49]. For an example in a different field but cognate to the macroscale modelling developed here, Mercer and Roberts [41] generalized the advection-diffusion ‘homogenization’ model of shear dispersion in a narrow pipe to arbitrary order and to a quantitative bound on the spatial resolution (its Section 2), with guaranteed accurate initial conditions (its Section 3), and correct macroscale boundary conditions at the pipe’s entry and exit (its Section 4). That is, this dynamical systems approach to homogenization potentially provides a rational complete modelling of not just the PDE at finite scale separation, but also material variations, uncertainty, control, forcing, and initial and boundary conditions. In their review of nonlinear homogenization, Matous et al. [38, Section 8] commented that “verification ensures both order of accuracy and consistency … and even mathematical foundations for verification are lacking”—the approach developed here provides the required mathematical foundation.

The new approach developed herein has no need of a variational formulation of the governing equations [4, 6, Section 3 and Section 2.4, respectively], and so applies to a much wider class of problems than many extant homogenization methods.

Throughout this article, in the details of the construction of homogenizations, one may observe that the variables $x_0,x_1,x_2$ appear in much the same sort of algebraic expressions as traditionally obtained by the method of multiple scales. What this parallelism implies is that in suitable circumstances, the algebraic machinations of the method of multiple scales are indeed sound [52, Corollary 12]. What is new here is that we place such algebraic details in a new, more rigorous and more extensive framework. The new framework developed here improves the clarity, precision, power and physical interpretability of homogenization.

2 Phase-shift embedding

Let us embed the specific given physical quasi-periodic PDE (1.1) into a family of PDE problems formed by all phase-shifts of the quasi-periodic microscale. This embedding is a novel and rigorous twist to the concept of an RVE.

As indicated by Figure 1, we create the desired embedding by considering a field ${\mathfrak {u}}(t,x_0,x_1,x_2)$ satisfying the PDE

(2.1)

in the domain ${\mathcal D}:=\{(x_0,x_1,x_2) \mid {x_i}\in [0,\ell _i]\},$ for the macroscale $\ell _0:=L$ , the specified microscale $\ell _1,\ell _2$ and with $\ell _i$ -periodic boundary conditions in ${x_i}$ . We assume $\kappa (x_1,x_2) \geq \kappa _{\min }>0$ . The domain ${\mathcal D}$ in $x_0x_1x_2$ -space (Figure 1) is often termed cylindrical as it is long and thin. I emphasize that this domain has finite aspect ratio: we do not take any limit involving its aspect ratio tending to zero nor to infinity.

Figure 1 indicates that we consider $x_0=x$ . The distinction between $x_0$ and x is that partial derivatives in $x_0$ are done keeping the other $x_i$ constant, whereas partial derivatives in x are done keeping phases $\phi _i$ constant.

Lemma 2.1. For every solution ${\mathfrak {u}}(t,x_0,x_1,x_2)$ of the embedding PDE (2.1) and for every vector of phases , the field (that is, the field ${\mathfrak {u}}$ evaluated on the solid-blue lines in Figure 1) satisfies the heterogeneous diffusion PDE

(2.2)

Hence, the field $u(t,x):={\mathfrak {u}}(t,x,x,x)$ satisfies the given heterogeneous PDE (1.1).

Figure 1 Schematic cylindrical domain ${\mathcal D}$ of the multiscale embedding PDE (2.1) for a field ${\mathfrak {u}}(t,x_0,x_1,x_2)$ (here $L=\ell _0=6$ , $\ell _1=1.62$ , $\ell _2=0.72$ ). We obtain solutions of the heterogeneous diffusion PDE (1.1), or (2.2), on such blue lines as $u(t,x)={\mathfrak {u}}(t,x,x+\phi _1,x+\phi _2)$ for every pair of constant phases $\phi _1$ and $\phi _2$ (here $\phi _1=0.82$ and $\phi _2=0.32$ ), and where the last two arguments of ${\mathfrak {u}}$ are modulo $\ell _1$ and $\ell _2$ , respectively.

Recall that the most common cell boundary conditions assumed for RVEs are periodic, although in the usual homogenization arguments, other boundary conditions are equally as valid in practice albeit giving slightly different results [39]. In contrast, here the $\ell _1,\ell _2$ -periodic cell boundary conditions are not assumed, but arise naturally due to the ensemble of phase-shifts. That is, what previously had to be assumed, here arises naturally.

Proof. Start by considering the left-hand side time derivative in PDE (2.2):

namely the right-hand side of (2.2). Hence, provided PDE (2.1) has boundary conditions of $\ell _i$ -periodicity in $x_i$ , every solution of the embedding PDE (2.1) gives a solution of the original PDE (1.1) for every multi-dimensional phase-shift of the heterogeneity.

In particular, the field satisfies the given heterogeneous PDE (1.1).

Lemma 2.2 (Converse).

Suppose we have a set of solutions of the phase-shifted PDE (2.2)—a set parametrized by the phase vector —and the set depends smoothly upon . Then the field ${\mathfrak {u}}(t,x,x_1,x_2):=u_{(x_1-x,x_2-x)}(t,x)$ satisfies the embedding PDE (2.1).

Proof. First, from the PDE (2.1), consider

Second, since $\phi _i=x_i-x$ , that is, $x_i=x+\phi _i$ , then

Third, hence the right-hand-side of PDE (2.1) becomes

the right-hand side of PDE (2.2). Lastly, since , it follows that ${\mathfrak {u}}:=u_{(x_1-x,x_2-x)}(t,x)$ satisfies the embedding PDE (2.1).

Corresponding results to Lemmas 2.1 and 2.2 hold for the heterogeneous wave PDE $u_{tt}=\partial _x\{\kappa (x)u_x\}$ , with replaced by .

Consequently, PDEs (2.1) and (2.2) are equivalent, and they may provide us with a set of solutions for an ensemble of materials all with the same microscale heterogeneity structure, but with the structural phase of the material shifted through all possibilities. The key difference between PDEs (2.1) and (2.2) is that although PDE (2.2) is heterogeneous in space x (as is (1.1)), the embedding PDE (2.1) is homogeneous in $x_0$ . Because of this homogeneity, Section 3 is empowered to apply an existing rigorous theory for slow variations in space that leads to the desired homogenized model PDE (1.2).

Discussion on translational symmetry. One should look at a macroscale homogenized model such as PDE (1.2) and wonder from where the continuous x-translation symmetry comes. The symmetry is not in the underlying microscale heterogeneous PDE (1.1). In the case where heterogeneity $\kappa (x)$ is periodic, the PDE (1.1) has discrete x-translational symmetry, which via the infinite scale separation limit as $\ell _1\to 0$ , might be viewed as approximately the continuous symmetry in the homogenization. However, in the case of quasi-periodic heterogeneity $\kappa (x)$ , for incommensurate $\ell _1$ and $\ell _2$ , there is no discrete translational symmetry. Instead of such a problematic limit, we here instead create the continuous translational symmetry by considering all together the family of problems formed by all phase shifts of the heterogeneity—the PDEs (2.1) and (2.2). The analysis of the whole family then naturally provides the continuous x-translational symmetry then inherited by the homogenization (1.2).

How does one get predictions correct for a specific heterogeneity? That is, for a given specific phase? The answer is that it is the correct boundary conditions for the homogenization PDE (1.2) that both reflect the specific phase of the heterogeneity at and near a specific boundary, and that also account for boundary layers. Such boundary conditions lie just beyond the scope of this article, but should be derivable from those for the microscale PDE (1.1) by further developing established dynamical systems analysis [14–16, 41, 48, 49].

3 The slow manifold of any-order homogenization

Let us analyse the embedding PDE (2.1) for a useful slow manifold. This slow manifold expresses and supports the emergence of a precise homogenization of the original heterogeneous PDE (1.1). Since the PDEs herein are linear, the slow manifold is more specifically a slow subspace, but I use the term manifold as the same framework and theory immediately also applies to cognate nonlinear systems.

Rigorous theory [52], inspired by earlier more formal arguments [50, 46], establishes how to create a PDE model for the macroscale spatial structure of PDE solutions in cylindrical domains ${\mathcal D}$ like Figure 1. The technique is to base analysis on the case where variations in $x_0$ are approximately negligible, and then treat slow, macroscale, variations in $x_0$ as a regular perturbation [53, Part III] ¹ .

Being “regular”, the perturbation analysis proceeds to any chosen order N in the ‘small’ derivatives in $x_0$ . For two examples, the usual leading order homogenizations are the case $N=2$ , and the so-called second-order homogenizations [5, 18, 31] correspond to the higher-order $N=4$ . Here we consider arbitrary order N. To ensure the required derivatives and operators exist and have requisite properties, we analyse the systems in Sobolov spaces $\mathbb {H}_D:=W^{N+1,2}(D)$ for various spatial domains D.

To establish the basis of a slow manifold, consider PDE (2.1), in the Sobolov space with $\mathbb {H}_{\mathcal D}$ , with neglected:

(3.1)

called the cell problem, to be solved with $\ell _i$ -periodic boundary conditions in ${x_i}$ . The basis applies at each and every $x_0$ . Hence, at every $x_0$ , we consider PDE (3.1) on the Sobolov space $\mathbb {H}_{\mathcal C}$ , where ${\mathcal C}$ is the $x_1x_2$ -cross-section of ${\mathcal D}$ .

Assumption 3.1. Consider the cell eigen-problem for ${\mathfrak {u}}\in \mathbb {H}_{\mathcal {C}}$ ,

with $\ell _i$ -periodic boundary conditions in ${x_i}$ . We assume that $\kappa (x_1,x_2) \geq \kappa _{\min }>0$ is regular enough that the set of eigenvectors are complete in $\mathbb {H}_{\mathcal C}$ .

Equilibria. A family of equilibria of (3.1) is ${\mathfrak {u}}(t,x_1,x_2)=U$ , constant in $x_1,x_2$ , for every U. This family, in $\mathbb {H}_{\mathcal C}$ for every $x_0$ , forms a subspace of equilibria $\mathbb {E}$ for PDE (3.1).

3.1 Spectrum at each equilibrium

We explore the spectrum about every equilibria in $\mathbb {E}$ as the spectrum is crucial in identifying the existence and properties of invariant manifolds/subspaces. Since the cell PDE (3.1) is linear, the perturbation problem is identical at every U, namely (3.1) itself.

Example 3.2. Consider the case where the heterogeneity $\kappa $ is constant. In the vector space $\mathbb {H}_{\mathcal C}$ , and upon defining wavenumbers $k_i =2\pi /\ell _i$ , a complete set of linearly independent eigenfunctions of the right-hand side of (3.1) are $\mathrm{e}^{\mathrm{i}(m k_1 x_1+n k_2 x_2)}$ for every integer $m,n$ . The spectrum of eigenvalues is then $\lambda _{mn}=-\kappa (mk_1+nk_2)^2$ . This spectrum has one zero-eigenvalue and the rest are negative: $\lambda _{mn}\leq -\kappa k_1^2 <0$ (recall $\ell _1>\ell _2$ so $k_1<k_2$ ).

For general heterogeneity $\kappa $ satisfying Assumption 3.1, we here establish that the cell eigenvalue problem associated with (3.1) is self-adjoint, and hence has only real eigenvalues. Denote the right-hand side operator as ${{\mathcal L}_0{\kern-1pt}:={\kern-1pt}(\partial _{x_1}{\kern-1pt}+{\kern-1pt}\partial _{x_2})\{\kappa (x_1,x_2) (\partial _{x_1}{\kern-1pt}+{\kern-1pt}\partial _{x_2})\}}$ to be solved with $\ell _i$ -periodic boundary conditions in ${x_i}$ . Also, we use the inner product $\langle {\mathfrak {u}},{\mathfrak {v}}\rangle :=\iint _{{\mathcal C}} {\mathfrak {u}}\bar {\mathfrak {v}} \,d{x_1}\,d{x_2}$ over the rectangular cross-section $\mathcal {C}:=[0,\ell _1]\times [0,\ell _2]$ , where an overbar denotes complex conjugation. Then, for every pair of functions ${\mathfrak {u}},{\mathfrak {v}}\in \mathbb {H}_{\mathcal C}$ , consider

$$ \begin{align*} \langle{\mathcal L}_0{\mathfrak{u}},{\mathfrak{v}}\rangle& =\iint_{{\mathcal C}} \bar{\mathfrak{v}}(\partial_{x_1}+\partial_{x_2})\{\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})\} \,d{x_1}\,d{x_2} \\[5pt]& =\iint_{{\mathcal C}} \bar{\mathfrak{v}}\partial_{x_1}\{\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})\} + \bar{\mathfrak{v}}\partial_{x_2}\{\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})\} \,d{x_1}\,d{x_2} \\[5pt]& =\int_0^{\ell_2} \bigg\{ [\bar{\mathfrak{v}}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})]_0^{\ell_1} -\int_0^{\ell_1} \bar{\mathfrak{v}}_{x_1}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2}) \,d{x_1} \bigg\}\,d{x_2} \\[5pt]&\quad{} + \int_0^{\ell_1} \bigg\{ [{\mathfrak{v}}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2})]_0^{\ell_2} -\int_0^{\ell_2} \bar{\mathfrak{v}}_{x_2}\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2}) \,d{x_2} \bigg\}\,d{x_1} \\[5pt]& =-\iint_{{\mathcal C}} (\bar{\mathfrak{v}}_{x_1}+\bar{\mathfrak{v}}_{x_2})\kappa ({\mathfrak{u}}_{x_1}+{\mathfrak{u}}_{x_2}) \,d{x_1}\,d{x_2}, \end{align*} $$

as the boundary contributions vanish due to the $\ell _i$ -periodicity. Reverse the above steps with the roles of ${\mathfrak {u}}$ and ${\mathfrak {v}}$ swapped, to then deduce $\langle {\mathcal L}_0{\mathfrak {u}},{\mathfrak {v}}\rangle =\langle {\mathfrak {u}},{\mathcal L}_0{\mathfrak {v}}\rangle $ . Thus, the operator ${\mathcal L}_0$ is self-adjoint.

Consequently, every eigenvalue is real and the eigenvectors are pairwise orthogonal (or may be chosen to be so).

The next step is to establish that all the nonzero eigenvalues of ${\mathcal L}_0$ are negative and bounded away from zero. Consider the eigen-problem ${\mathcal L}_0{\mathfrak {v}}=\lambda {\mathfrak {v}}$ for cross-sectional structures ${\mathfrak {v}}\in \mathbb {H}_{\mathcal C}$ . Apply $\langle \cdot ,{\mathfrak {v}}\rangle $ to the eigen-problem to see that

$$ \begin{align*} \lambda\langle{\mathfrak{v}},{\mathfrak{v}}\rangle &=\iint_{\mathcal C} \bar{\mathfrak{v}}(\partial_{x_1}+\partial_{x_2})\{\kappa ({\mathfrak{v}}_{x_1}+{\mathfrak{v}}_{x_2})\}\,d{x_1}\,d{x_2} \\& =\cdots=-\iint_{\mathcal C} \kappa |{\mathfrak{v}}_{x_1}+{\mathfrak{v}}_{x_2}|^2\,d{x_1}\,d{x_2} \\& \leq-\kappa_{\min}\iint_{\mathcal C} |{\mathfrak{v}}_{x_1}+{\mathfrak{v}}_{x_2}|^2\,d{x_1}\,d{x_2} \,, \end{align*} $$

via integration by parts and the $\ell _i$ -periodic conditions in $x_i$ . To deduce an upper bound on the nonzero eigenvalues $\lambda $ , let us minimize $\iint _{\mathcal C} |{\mathfrak {v}}_{x_1}+{\mathfrak {v}}_{x_2}|^2\,d{x_1}\,d{x_2}$ in $\mathbb {H}_{\mathcal C}$ subject to $\langle {\mathfrak {v}},{\mathfrak {v}}\rangle =\iint _{\mathcal C}|{\mathfrak {v}}|^2\,d{x_1}\,d{x_2}=1$ and that ${\mathfrak {v}}({x_1},x_2)$ has zero mean (being orthogonal to the constant eigenvector corresponding to $\lambda =0$ ). Introduce a Lagrange multiplier $\mu $ , and use calculus of variations to minimize $\iint _{{\mathcal C}} |{\mathfrak {v}}_{x_1}+{\mathfrak {v}}_{x_2}|^2 \,d{x_1}\,d{x_2} +\mu ( \iint _{{\mathcal C}}|{\mathfrak {v}}|^2 \,d{x_1}\,d{x_2}-1 )$ . Via integration by parts, it follows that a minimizer must satisfy $(\partial _{x_1}+\partial _{x_2})^2{\mathfrak {v}}=-\mu {\mathfrak {v}}$ over ${\mathfrak {v}}\in \mathbb {H}_{{\mathcal C}}$ . This is the problem Example 3.2 solves, giving potential minimizers as $\mathrm{e}^{\mathrm{i}(m k_1 x_1+n k_2 x_2)}$ for multiplier $\mu =(mk_1+nk_2)^2$ . The zero mean condition requires at least one of $m,n$ nonzero, and since $\ell _1>\ell _2$ , so minima are thus linear combinations of $\mathrm{e}^{\mathrm{i} m k_1 x_1}$ for $m=\pm 1$ . These give $\iint _{{\mathcal C}} |{\mathfrak {v}}_{x_1}+{\mathfrak {v}}_{x_2}|^2\,d{x_1}\,d{x_2}/\langle {\mathfrak {v}},{\mathfrak {v}}\rangle =k_1^2=4\pi ^2/\ell _1^2$ . Hence, we obtain the upper bound on the nonzero eigenvalues of $\lambda \leq -\beta _1$ for $\beta _1=4\pi ^2\kappa _{\min }/\ell _1^2$ .

In the case of the heterogeneous wave propagation problem, $u_{tt}=\partial _x\{\kappa (x)u_x\}$ , this subsection’s eigenvalue analysis establishes that for the cross-sectional dynamics, there exists a zero eigenvalue separated from pure imaginary eigenvalues of magnitude at least $\sqrt {\beta _1}=2\pi \sqrt {\kappa _{\min }}/\ell _1$ .

3.2 Results of slowly varying theory

Let us comment on the preconditions for the next Proposition 3.3 as listed in [52, Assumption 2]. First, Section 3.1 establishes that the spectrum of ${\mathcal L}_0$ is as required for the Hilbert space $\mathbb {H}_{\mathcal C}$ . Hence, we decompose the space into two closed ${\mathcal L}_0$ -invariant subspaces, the (slow) centre $\mathbb {E}_c^0$ and the stable $\mathbb {E}_s^0$ . Since the operator ${\mathcal L}_0$ is self-adjoint, it generates the required strongly continuous semigroups in the Sobolev space $\mathbb {H}_{\mathcal C}$ [13, Ch. 6]. Hence, the following result holds for homogenization.

Proposition 3.3 (Homogenization exists and emerges).

By [52, Proposition 6], for every chosen truncation order N, in a domain $\mathbb {X}$ of ‘slowly varying solutions’ and after the decay of transients for $\beta '\approx \beta _1$ , the mean field $U(t,x):=(\ell _1\ell _2)^{-1} \iint _{\mathcal C} {\mathfrak {u}}(t,x,x_1,x_2) \,d{x_1}\,d{x_2}$ satisfies a generalized homogeneous PDE

(3.2)

in terms of some deterministic constants $K_n$ , to a quantifiable error [52, (9)–(10) and (23)].

The above mentioned error [52, (23)] involves too many new parameters to meaningfully detail here—but typically, the most important factor at each and every location $x=x_0$ is the $(N+1)$ th spatial derivative of the field ${\mathfrak {u}}$ in a spatial neighbourhood of $x_0$ . The open subset $\mathbb {X}$ of the domain $[0,L]$ is then that part of the domain where the error in (3.2) [52, (23)], dominantly $\partial _x^{N+1}{\mathfrak {u}}$ , is small enough for the modelling purposes at hand.

Also, using the specific mean field defined in Proposition 3.3—the usual “cell average”—is optional: one may instead parametrize the macroscale slow manifold in terms of any reasonable alternative amplitude [53, Section 5.3.3].

A generating function argument [52, Corollaries 12 and 13] establishes that classic, formal, ‘slowly varying’ analyses of the embedding PDE (2.1) are valid. One just needs to formally treat derivatives in the ‘large’ direction, , as ‘small’ in a consistent asymptotic sense [50, 46].

3.3 Wave propagation in heterogeneous media

In the case of the heterogeneous wave propagation problem, $u_{tt}=\partial _x\{\kappa (x)u_x\}$ , I conjecture that a backwards approach [32, 54] to the slowly varying theory here would support a cognate version of Proposition 3.3 about the generalized homogeneous wave PDE $U_{tt}=\sum _{n=0}^N K_n\partial _x^nU$ describing its dispersive macroscale waves.

However, in this case of heterogeneous wave propagation in heterogeneous media, the slow manifold is no longer attractive/emergent (unless there is some significant perturbing dissipative mechanism not included in the ideal wave PDE). Instead, we expect that the slow manifold acts as a ‘centre’ for any fast oscillations arising from the initial conditions. Indeed, conjectured backwards theory [54, Section 2.5] would more precisely assert that there is a system ‘close’ to the given physical PDE $u_{tt}=\partial _x\{\kappa (x)u_x\}$ (close to some specified order in some sense) that precisely possesses the constructed homogeneous slow manifold for the macroscale, and that precisely acts as a centre for all nearby dynamics. In this way, the backwards theory would straightforwardly establish accurate modelling over “time scales ” [1, 24], typically for exponent $\alpha =N$ , the order of truncation.

Although this article addresses homogenization of linear systems, a comment about homogenizing nonlinear wave systems is appropriate here. As geophysical fluid dynamists have understood for decades, for any given wave system, it is common that a nonlinear slow manifold, such as homogenization seeks, cannot actually exist [36]. Nonetheless, backwards theory [32, 54] potentially applies to establish that there exists a wave system ‘close’ to the specified nonlinear wave system, where the ‘close’ systems precisely possess the constructed nonlinear slow manifold and the constructed evolution thereon. Thus, the dynamical systems approach developed here provides a methodology to prevail over subtleties in homogenizing nonlinear wave problems.

4 Example construction of any-order homogenized PDEs

Corollary 13 of [52] established that a procedure [46], previously viewed as formal, is indeed a rigorous method to construct the slow manifold modelling PDEs (3.2). Here we use the computer algebra system (cas) version of the procedure [50], as detailed in generality and examples of the book [53, Part III].

Appendix A lists and documents the computer algebra code. The code constructs the homogenization for the case of heterogeneity,

(4.1) $$ \begin{align} \kappa(x_1,x_2):=1/[1+a_1\cos k_1x_1+a_2\cos k_2x_2], \end{align} $$

in terms of microscale wavenumbers $k_i=2\pi /\ell _i$ for the given microscale periodicities $\ell _1,\ell _2$ . We write approximations to the slow manifold model of the embedding PDE (2.1) in terms of a ‘mean’ field $U(t,x)$ that evolves according to an homogenized PDE of the form (3.2).

Quickly verify an approximation. First, Appendix A.1 verifies that an approximate field is

(4.2a)

such that the mean field U evolves according to the homogenized PDE

(4.2b)

(the coefficient $1$ of the mean diffusivity is the classic harmonic mean of the specific heterogeneous diffusivity (4.1)). To verify, substitute (4.2) into the governing PDE (2.1) and find the PDE’s residual contains only terms in $U_{xxx}$ and higher derivatives. We denote this by saying that the residual is zero to an error . Since the PDE residual is , then theory [52, Proposition 6] assures us that the derived slow manifold and its evolution (4.2) is correct to error .

4.1 Iteration systematically constructs approximations

To construct approximations systematically, we repeatedly compute the residual, which then drives corrections, until the residual is zero to any specified order of error. Theory then assures us that the slow manifold is approximated to the same order of error [52]. For example, the cas iteratively constructs the improved homogenization (the so-called ‘fourth-order’ case $N=6$ ) that

(4.3)

The code of Appendix A.2 constructs this ‘fourth-order’ homogenization (4.3), and executes in approximately 26 seconds. In the case of the heterogeneous wave propagation problem, $u_{tt}=\partial _x\{\kappa (x)u_x\}$ , we expect the macroscale homogenized dispersive wave PDE to have $U_{tt}$ equal to the same right-hand side as the above. Most other rigorous homogenizations invoke a limit ‘ $\epsilon \to 0$ ’ for some defined scale separation parameter such as $\epsilon =\ell _1/L$ . Here such a limit corresponds to the scale separation limit $k_1, k_2\to \infty $ , in which case, the homogenization (4.3) reduces to the classic, simple, diffusion PDE $U_t=U_{xx}$ . However, I emphasize that our homogenized PDEs like (4.3) are established for the finite scale separation of real physical materials: finite $\epsilon ,k_1,k_2$ .

In use, any specific truncation of higher-order homogenization, such as (4.3), may need some asymptotically consistent regularization [7, Section 2]. For an example of regularization, let us suppose we decide on the $N=4$ truncation ${U_t=U_{xx}+\alpha ^2U_{xxxx}}$ with derived parameter . The fourth derivative term, $+\alpha ^2U_{xxxx}$ is undesirably destabilizing. However, to the same order of error, namely , the chosen truncation is asymptotically equivalent to the regularized equation ${U_t= (1-\alpha ^2\partial _x^2)^{-1} U_{xx}}$ . This regularization may be alternatively written as

$$ \begin{align*} (1-\alpha^2\partial_x^2)U_t=U_{xx}\,, \quad\text{or}\quad U_t =\frac1{2\alpha} \mathrm{e}^{-|x|/\alpha}\star U_{xx}\,. \end{align*} $$

The spatially nonlocal convolution $\mathrm{e}^{-|x/\alpha |}\star $ , whether explicit or implicit, stabilizes the regularized chosen truncation while preserving asymptotic consistency [41, cf. Section 2.2]. Such nonlocal PDEs have previously been found desirable in various homogenizations of heterogeneous systems [17, pp. 2–3].

Optional high-order construction. We may set $a_2=0$ to get simpler results for a single periodic component in the heterogeneity, namely ${\kappa (x_1)=1/[1+a_1\cos k_1x_1]}$ . In this case, it becomes feasible to construct the homogenization to high-order, $N=34$ , in spatial derivatives. To construct expressions that are exact in heterogeneity amplitude $a_1$ , we also need to construct to the same order in $a_1$ because the evidence is that here, the expansions in heterogeneity amplitude $a_1$ truncate at just less than the same order as the order of derivatives. For this option (Appendix A.2), the algebra of order $N=34$ takes roughly one minute of compute time.

In general. The code of Appendix A.2 assumes the heterogeneity $\kappa \approx 1$ , specifically, $\kappa =1/[1+\kappa ']$ for some $\kappa '\propto a$ . Then multiplication by $\kappa $ is realized as multiplication by $\sum _{n=0}^{N-1}(-\kappa ')^n$ .

The method [50] is to start the iteration with the leading ‘trivial’ approximation for the field and its evolution: that is, $u\approx u^{(1)}: = U$ such that . Then the nth iteration is given approximations $u^{(n)}$ and $g^{(n)}$ , and starts by evaluating, via the flux $f^{(n)}: =-\kappa (u^{(n)}_x+u^{(n)}_{x_1}+u^{(n)}_{x_2})$ , and then the residual of the embedding PDE (2.1):

$$ \begin{align*} {\operatorname{Res}}^{(n)}:=u^{(n)}_t+f^{(n)}_x+f^{(n)}_{x_1}+f^{(n)}_{x_2}. \end{align*} $$

First update the evolution via the solvability integral that here is the mean over $x_1,x_2$ : $g^{(n+1)}:=g^{(n)}+g'$ , where $g' :=\overline {{\operatorname {Res}}^{(n)}}$ . Second, update the field by a correction $u'$ through using two steps to solve $(\partial _{x_1} +\partial _{x_2}) [ \kappa (u^{\prime }_{x_1} +u^{\prime }_{x_2}) ] ={\operatorname {Res}}^{(n)} +g'$ : step 1 integrates to solve $v_{x_1}+v_{x_2}={\operatorname {Res}}^{(n)} +g'$ with zero mean to avoid unbounded updates; and step 2 integrates again to solve $u^{\prime }_{x_1}+u^{\prime }_{x_2}=v/\kappa $ with zero mean. This second “zero mean” ensures preservation of our free choice that the macroscale field parameter U is to be the local mean of the microscale field u. Then update the field approximation with $u^{(n+1)}: =u^{(n)}+u'$ .

The iterative loop terminates when the PDE residual ${\operatorname {Res}}^{(n)}$ is zero to the specified orders of error. Then the constructed homogenization, such as (4.2) and (4.3), is correct to the same order [52, Proposition 6].

Analogy with machine learning. The recursive iteration used in the above construction was originally developed nearly 30 years ago [50]. Let us draw an analogy between this iteration and machine learning algorithms where an ai learns the generic form of the macroscale evolution from many thousands of simulations [28, 29, 35]. In the algorithm here, each recursive iteration is analogous to a layer in a deep neural network, evaluating residuals is analogous to a nonlinear neuronal function, and the linear update corrections are analogous to using weighted linear combination of outputs of one layer as the inputs of the next layer. However here, the recursion is a ‘smart neural network’ in that both the neurons and the ‘linear weights’ are crafted to the problem at hand using the physics encoded in the PDE operator—they are ‘physics-informed’. Further, being analytic, a single analysis encompasses all ‘data points’ in the state space’s domain (here all ‘slowly varying’ functions in $\mathbb {H}_{{\mathcal D}}$ ), not just a finite sample as in machine learning (possibly a biased sample). Consequently, via such algorithms, I contend that mathematicians have for decades been using such recursive iteration for implementing smart analogues of machine learning. Such analytic learning empowers the physical interpretation, validation and verification required by modern science [10].

4.2 Numerics verify these homogenizations

For a computational verification, let us use the specific instance of the heterogeneity (4.1) with $a_1=a_2=0.4$ , and wavenumbers $k_1=21$ and $k_2=34$ , that is, microscale lengths are $\ell _1=2\pi /k_1\approx 0.30$ and $\ell _2=2\pi /k_2\approx 0.18$ (all nondimensional). Consequently, let the macroscale then refer to lengths ${}\gtrsim 2\ell _1\approx 0.6$ (Section 5), that is, wavenumbers $|k|\lesssim 10$ . We compare the macroscale characteristics of the right-hand side of the homogenized (4.3) with the right-hand side of the original heterogeneous PDE (1.1).

Here the microscale heterogeneity (4.1) (see Figure 2), is near-quasi-periodic with $k_2/k_1=1.6190\approx (\sqrt (5)+1)/2$ , although this heterogeneity does have a large macroscale period $2\pi $ (as plotted). Consequently, numerically, we solve heterogeneous PDE (1.1) on the domain of length $L=2\pi $ with a microgrid of $714$ points, . The “numeric” column of Table 1 gives the resultant eigenvalues of the full heterogeneous PDE (1.1) for macroscale wavenumbers k, $|k|<10$ . These eigenvalues reflect either the decay rate of the diffusion problem (1.1) or the square of frequencies for the corresponding wave problem (1.3). The homogenization (4.3) of this quasi-periodic heterogeneity may be truncated at any chosen order in $\partial _x$ , so we compare its predictions for three truncations: the usual leading order homogenization with truncation to errors ; the ‘second-order’ homogenization with errors ; and the ‘fourth-order’ homogenization with errors . We simply substitute the mode $U=\mathrm{e}^{\lambda t+\mathrm{i}kx}$ into the truncated right-hand side of (4.3) to obtain rates to compare with that of the full problem. Table 1 shows the two expected results: first, for each wavenumber, the accuracy increases with order of truncation; and second, for each order, as the macroscale wavenumber becomes larger (the macroscale length of the mode decreases), the error increases. These expected trends in the errors are shown clearly in Figure 3 that plots the relative errors on log–log axes ² . Table 1 and Figure 3 verify our approach to modelling the homogenized dynamics of quasi-periodic heterogeneous problems.

Figure 2 Example (near) quasi-periodic microscale heterogeneous coefficient (4.1) for the case of ${a_1=a_2=0.4}$ and wavenumbers $k_1=21$ and $k_2=34$ .

Figure 3 Relative errors in the macroscale predictions by the three truncations of the homogenized PDE (4.3). Table 1 gives the data underlying this plot of errors. These are for the quasi-periodic heterogeneity (4.1) with $a_1=a_2=0.4$ and wavenumbers $k_1=21$ and $k_2=34$ .

Figure 3 indicates that the homogenization errors become large for wavenumbers $|k|\gtrsim 10\approx k_1/2$ . This apparent limit of the resolvable wavenumbers matches the estimate obtained in Section 5, the lower bound (5.1), that the macroscale lengths resolved by the homogenization are longer than roughly twice $\ell _1$ .

Table 1 Macroscale eigenvalues as a function of macroscale wavenumber k: the “numeric” column lists the computed eigenvalues of (1.1) for spatial discretization ; the last three columns are the rates predicted for $U=\mathrm{e}^{\mathrm{i} kx}$ by the homogenization (4.3) truncated to errors for orders $n=2,4,6$ , respectively.

5 Spatial resolution of high-order homogenization

An homogenized PDE, such as (3.2) and (4.3), accurately predicts the evolution of the mean field $U(t,x)$  provided that, on the scale of the heterogeneity, the spatial gradients are sufficiently small, that is, provided the solutions are sufficiently slowly varying in space. Matous et al. [38, Section 3.2] phrased the proviso as “the scale of the microstructural fluctuations, [ $\ell _1$ ], must be … much smaller than the macroscopic field fluctuations, $\ell _M$ ” (expressed similarly by Nguyen et al. [43, Section 2.2.2]). The following question arises: what does “sufficiently small”, “sufficiently slowly varying” or “much smaller” actually mean physically? This section provides an indicative quantitative answer.

Here we show that in a prototypical family of homogenization problems, we only need to require that the

(5.1) $$ \begin{align} \text{macroscale lengths } \ell_M\gtrsim2\ell_1. \end{align} $$

This lower bound on the macroscale lengths is approximate and indicative because the actual performance of an homogenization depends upon the details of the heterogeneity $\kappa $ , a chosen truncation N of the homogenized PDE, any chosen specific regularization of the PDE, the chosen desired error of predictions and so on. Recall that this lower bound (5.1) agrees with the growth of errors seen in the example of Figure 3. However here, the numerical coefficient $2$ in (5.1) arises from a well-established quantitative procedure.

Such a result contrasts with the prevailing expectation provided by other homogenization frameworks: for example, Allaire et al. [3, p. 2] comment “no convergence is expected”.

The procedure we invoke is that of high-order asymptotic construction followed by Domb–Sykes plots [22, 67] or Mercer–Roberts plots [40, Appendix] that estimate the distance to the nearest, convergence limiting, singularity in the series approximations. We analyse an homogenized PDE (3.2) and (4.3) by considering its corresponding Fourier transform ( $\partial _x\leftrightarrow \mathrm{i} k$ ): in terms of the Fourier transform $\tilde U(t,k)$ of the mean field $U(t,x)$ . Then the coefficient ${\mathcal K}(k):=\sum _{n=0}^\infty (\mathrm{i}^nK_n)k^n$ is a series in wavenumber k. If we can show convergence for $|k|<k_*$ for some wavenumber bound  $k_*$ , then we deduce that the homogenized PDE is valid for spatial structures of length scale $\ell _M>2\pi /k_*$ . We know that the radius of convergence $k_*$ is the distance to the singularity of ${\mathcal K}(k)$ nearest the origin ( $k=0$ ). However, we only know N terms in the series for ${\mathcal K}(k)$ , so Domb–Sykes and Mercer–Roberts plots estimate the distance  $k_*$ by extrapolating the ‘ratio test’ from finite n to infinity, in the cases of a single singularity or a complex conjugate pair of singularities, respectively. For example, this technique has been used to quantitatively predict the spatial limit on PDEs modelling shear dispersion in channels and pipes [40, 41]. Figure 4 is an example Mercer–Roberts plot for one homogenization case developed here. The extrapolations to $1/n=0$ ( $n\to \infty $ ) in the figure indicate that the convergence is limited by a complex conjugate pair of singularities at a distance in wavenumber space of $k_*\approx 1/\sqrt {3.6} \approx 0.52$ .

Figure 4 Macroscale wavenumbers k of convergence are estimated from Mercer–Roberts plots [40, Appendix] for the series ${\mathcal K}(k)$ in the case of heterogeneity (5.2), here with $a=0.975$ . We use the first 17 terms of the series in $k^2$ . Extrapolation to $1/n=0$ estimates the location of a pair of complex conjugate singularities that limit the radius of convergence.

To obtain such plots and estimates, we need to compute the homogenization series (3.2) and (4.3) to high-order. For the case of quasi-periodic heterogeneity, the algebraic complexity explodes combinatorially with increasing order. Consequently, as a representative example, we explore the accessible case of homogenization with a single microscale period, specifically the heterogeneity

(5.2) $$ \begin{align} \kappa(x):=1/[1+a\cos x]. \end{align} $$

The computer algebra code overviewed in Section 4 constructs the homogenization in this case via the option of setting $a_2=0$ , $a_1=a$ and $k_1=1$ . For algebraic simplicity, in this case, the analysis is nondimensionalized on the microscale length so that nondimensionally, the microscale period $\ell _1=2\pi $ . Executing the code with this option gives that the homogenized PDE (3.2) begins with

(5.3)

As a first step in investigating such high-order homogenization, observe that the $a^2$ terms are evidently all of the form

$$ \begin{align*}\tfrac12a^2 \partial_x^4 (-4\partial_x^2)^{n/2-2}\end{align*} $$

(this form holds up to at least order $n=34$ ). Hence, in Fourier space ( $\partial _x\leftrightarrow {\mathrm{i}} k$ ), the PDE (5.3) becomes the series for ³

The divisor $(1-4k^2)$ in this expression indicates that for small heterogeneity amplitude a, the wavenumbers resolvable by such an homogenization series are limited by the pole singularities at $k=\pm 1/2$ , that is, a bound on the convergence of the series is $k_*=1/2$ . Hence, for small a, the homogenized PDE (5.3) potentially resolves all wavelengths longer than $\ell _*:=2\pi /k_*=4\pi =2\ell _1$ ⁴ .

I conjecture that the singularities at wavenumber $k=\pm 1/2$ are connected to, in wave propagation through heterogeneous media, the well-known phenomenon of spectral gaps, namely frequencies at which no wave can propagate through the underlying medium [17, p. 2].

Figure 5, for small a, confirms this bound on the wavenumbers and wavelengths. The estimates plotted in Figure 5 are obtained from Section 4 constructing the series (5.3) to errors of order $35$ in $\partial _x$ , or equivalently to errors . This cas construction provides us with the first 17 nonzero coefficients in the series in $k^2$ for ${\mathcal K}(k)$ .

• • For heterogeneity amplitude $a\lesssim 0.4$ , these coefficients have the same sign. Consequently, Domb–Sykes plots [22, 67], akin to the top panel of Figure 4, predict convergence limiting singularities at real $k=\pm k_*$ as shown in the left-part of Figure 5.

• • For heterogeneity amplitude $a\gtrsim 0.4$ , the coefficients have a more complicated sign pattern, indicating the convergence limiting singularity is now a complex conjugate pair in the complex $k^2$ -plane. Mercer–Roberts plots [40, Appendix], such as Figure 4, estimate the radius and angle of these convergence limiting singularities as shown in the right-part of Figure 5.

Figure 5 In complex wavenumber space, plot estimates of the radius $k_*$ (circles) and angle $\theta _*$ (crosses) of the convergence limiting singularity of the homogenization (5.3) as a function of the heterogeneity amplitude a in (5.2). The estimates around $a\approx 0.4$ are unreliable (see text).

The two distinct behaviours summarized in Figure 5 are likely due to singularities moving as a function of heterogeneity parameter a. We conjecture that for $a\lesssim 0.3$ , there is a second real singularity for wavenumbers larger than $k_*$ . As a increases, both these singularities move and appear to collide for $a\approx 0.35$ . After collision, they separate and move out into the complex k plane ( $a\gtrsim 0.4$ ). The Domb–Sykes and Mercer–Roberts estimates are both unreliable near such collision, so the estimates for $a\approx 0.3$ to $\approx 0.4$ are unreliable. With this caveat, the evidence of Figure 5 indicates that across all $0<a<1$ , our homogenization can resolve all wavenumbers $|k|<1/2$ . That is, high-order homogenization can resolve wavelengths longer than just twice the wavelength of the underlying periodicity.

In the case of the heterogeneous wave propagation problem, $u_{tt}=\partial _x\{\kappa (x)u_x\}$ , we expect the same restrictions on the spatial resolution because we expect the right-hand side of PDE (4.3) to remain the same for the homogenized waves.

The case of quasi-periodic heterogeneity. In the case of the quasi-periodic heterogeneity (4.1), observe in the homogenization PDE (4.3) that the coefficients in $a_i$ are divided by wavenumbers $k_i$ , and with increasing powers of $k_i$ as the order increases. This suggests that higher-order terms in the series are dominated by the smallest wavenumber $k_i$ , namely $k_1$ . That is, the higher-order terms should be dominated by the longest microscale heterogeneity, here that of $a_1\cos (k_1x)$ . Hence, we expect the limiting wavenumber and wavelengths to be the same as those for the single periodic heterogeneity (5.2). Consequently, expect the bound (5.1) to be valid for this quasi-periodic heterogeneity.

Temporal resolution. The above exploration finds that the spatial resolution of the homogenization appears to vary little with heterogeneity strength a. The temporal resolution is quite a different issue, and it may depend significantly upon a. For example, Figure 4 comes from the case $a=0.975$ , where the microscale heterogeneity varies by a factor of $80$ across a microscale period. In similar scenarios, the time scale of decay to the homogenized model, roughly $1/\beta _1=\ell _1^2/(4\pi ^2\kappa _{\min })$ (end of Section 3.1), may be several orders of magnitude longer than for small heterogeneity. For valid predictions by the homogenization, the macroscale time scale of interest needs to be longer than that of this decay.

Of course, all equilibrium problems satisfy this temporal bound.

6 Three-scale homogenization

There is interest in physical systems possessing a microscale, a mesoscale and a macroscale [21, 45]. Such three-scale systems are encompassed in the framework of Section 2 when $\ell _2\ll \ell _1\ll L$ . That is, the cylindrical domain of Figure 1 is shaped like a physical ruler in that the domain is significantly thinner vertically than its width, and its width is in turn significantly thinner than its length. Although Section 5 finds that these length ratios need to be at least two, they are not necessarily much bigger.

In this scenario, the arguments of Section 2 still apply: we solve the phase-shift embedding PDE (2.1) to obtain solutions of the heterogeneous PDEs (1.1) and (2.2). It is the slow manifold support and analysis that differs in detail.

Because this rigorous homogenization could, in principle, be exact at finite scale separation, we find here the expected result that a slow manifold homogenization of a ‘meso-slow’ manifold homogenization of the microscale system is the same as the slow manifold homogenization of the microscale (Sections 3 and 4). That is, in principle, this approach to homogenization is transitive.

The same transitivity of homogenization would also hold in the case of the heterogeneous wave propagation problem, $u_{tt}=\partial _x\{\kappa (x)u_x\}$ .

6.1 Meso-slow homogenization

To analyse the heterogeneous PDEs (1.1) and (2.2) as a three-scale system, we first establish and create the ‘meso-slow’ manifold homogenization. That is, we use that the domain ${\mathcal D}$ , Figure 1, is thin in $x_2$ , that is, $\ell _2\ll \ell _1,L$ . Consequently, we consider the embedding PDE (2.1) by assuming the spatial variations in both $x_0$ and $x_1$ are on a relatively large length scale. That is, the derivatives of ${\mathfrak {u}}$ and $\kappa $ in both $x_0$ and $x_1$ are treated as small.

This homogenization is an example of that of a functionally graded material [5], because the mesoscale variations in heterogeneity in $x_1$ are on a length scale $\ell _1\gg \ell _2$ . Our slow manifold approach systematically incorporates all such variations within the analysis: for example, giving the fourth-order homogenization (6.3)–(6.4). At orders higher than the leading order, we systematically find higher-order derivatives of the larger scale material structure: for example, the fourth-order coefficient in the upcoming PDE (6.4) involves the third derivative ${\mathfrak {K}}"'(x_1)$ of the mesoscale variations. This approach automatically incorporates effects obtained by the “second-order homogenization” of Anthoine [5].

Akin to Section 3, we here establish the basis of a slow manifold. Consider PDE (2.1), with and neglected, and in the Sobolov space $\mathbb {H}_{[0,\ell _2]}$ with $\ell _2$ -periodic boundary conditions in ${x_2}$ :

(6.1)

The basis established here applies at each and every $x_0,x_1$ .

Equilibria. For the cell PDE (6.1), and for every $x_0,x_1$ , in $\mathbb {H}_{[0,\ell _2]}$ , clearly ${\mathfrak {u}}(t,x_2)={\mathfrak {U}}$ , constant in $x_2$ , forms a subspace of equilibria $\mathbb {E}$ .

Spectrum identifies invariant manifolds. Since the cell PDE (6.1) is linear, the perturbation problem is identical at every equilibria in $\mathbb {E}$ , namely (6.1) itself. We need to characterize the right-hand side operator in (6.1), with $\ell _2$ -periodic boundary conditions, in $\mathbb {H}_{[0,\ell _2]}$ .

For example, in the case where the heterogeneity is constant with respect to $x_2$ , a complete set of linearly independent eigenfunctions of ${\mathfrak {L}}_0$ are $\mathrm{e}^{\mathrm{i} n k_2 x_2}$ for every integer n. The spectrum of eigenvalues at $(x_0,x_1)$ is then $\lambda _{n}(x_1):=-\kappa (x_1)k_2^2n^2$ . This spectrum has one zero-eigenvalue and the rest are negative: $\lambda _{n}(x_1)\leq -\kappa (x_1) k_2^2 <0$ . Hence, for every $x_0,x_1$ , the nonzero eigenvalues satisfy $\lambda _n(x_1)\leq -\kappa _{\min }k_2^2<0$ .

For general heterogeneity, similar deductions to those of Section 3.1 establish that the operator ${\mathfrak {L}}_0$ is self-adjoint. Hence, the eigenvalues are real and eigenfunctions orthogonal. By a similar argument to that of Section 3.1, it follows that the nonzero eigenvalues of the operator are bounded above by $\lambda \leq -\beta _2$ for ${\beta _2:=\kappa _{\min }k_2^2=4\pi ^2\kappa _{\min }/\ell _2^2}$ .

Slowly varying theory. We proceed analogously to Section 3.2, but now with two space dimensions ( $x_0$ and $x_1$ ) in which solutions are slowly varying. Roberts and Bunder [55] developed rigorous slow manifold theory for spatio-temporal systems that have slow variations in multiple space dimensions. The operator ${\mathfrak {L}}_0$ satisfies the preconditions in [55, Assumption 3] with one zero eigenvalue ( $m=1,\alpha =0$ ) and the rest negative. By [55, Proposition 1], and analogous to Proposition 3.3, we are assured that in a regime of slowly varying solutions, the mesoscale field ${\mathfrak {U}}(t,x_0,x_1)$ satisfies a PDE akin to (3.2), but in 2D space, to a quantified error [55, (52)], and upon neglecting exponentially decaying transients (decaying faster than $\mathrm{e}^{-\beta 't}$ for $\beta '\approx \beta _2$ ).

A generating function argument [55, Sections 3.2–3.3] establishes the validity of formal ‘slowly varying’ analysis of the embedding PDE (2.1).

An example mesoscale homogenization. Appendix B lists and documents cas code to construct the mesoscale homogenization for the example case of heterogeneous diffusivity

(6.2) $$ \begin{align} \kappa(x_1,x_2):=1/[1/{\mathfrak{K}}(x_1)+a_2\cos k_2x_2] \quad\text{for }k_2:=2\pi/\ell_2\,, \end{align} $$

and some given mesoscale heterogeneity function ${\mathfrak {K}}(x_1)$ . The code derives that in terms of the mesoscale mean ${\mathfrak {U}}(t,x_0,x_1)$ , and in terms of ${\mathfrak {D}}_x:=\partial _{x_0}+\partial _{x_1}$ , the microscale field (cf. (4.2a)) is

(6.3)

Simultaneously, the code derives that the mesoscale mean field evolves according to the quasi-diffusion PDE

(6.4)

As to be expected, the leading-order effective diffusivity on the mesoscale, ${\mathfrak {K}}(x_1)$ , is the harmonic mean over $x_2$ of the microscale diffusivity (6.2). The theory of this subsection assures us that a truncation of the PDE (6.4) models the dynamics of the embedded PDE (2.1) whenever and wherever the spatial gradients in $x_0$ and $x_1$ are ‘small’ enough on the microscale $\ell _2$ . That is, PDE (6.4) is a mesoscale model.

6.2 Homogenization of the mesoscale model

The mesoscale PDE (6.4) is heterogeneous on the spatial length scale $\ell _1$ of variations in the coefficient ${\mathfrak {K}}(x_1)$ . For the case $\ell _1\ll L$ , we here derive the macroscale homogenization of this mesoscale model. We have to choose some truncation of PDE (6.4) to analyse further, so let us choose the leading order truncation ${\mathfrak {U}}_t={\mathfrak {D}}_x\{ {\mathfrak {K}}(x_1) {\mathfrak {D}}_x{\mathfrak {U}} \}$ as the mesoscale model.

To compare with the one-step homogenization, let us consider PDE (6.4) in the case where the mesoscale diffusivity is

$$ \begin{align*} {\mathfrak{K}}(x_1):=1/[ 1+a_1\cos k_1x_1] \quad\text{for }k_1:=2\pi/\ell_1\,. \end{align*} $$

The theory and construction of the homogenization follow the same procedure as detailed before, so let us just summarize here. First, we treat derivatives as ‘small’. Then the cross-sectional cell PDE is

(6.5)

to be solved with $\ell _1$ -periodic boundary conditions in $x_1$ , and in space ${\mathcal H}_{[0,\ell _1]}$ .

Second, a subspace of equilibria of (6.5) is ${\mathfrak {U}}(t,x_1)=U$ , constant in $x_1$ . The spectrum of the linear operator can be straightforwardly shown to be self-adjoint with one zero eigenvalue, and the rest negative and bounded away from zero, $\lambda \leq -\beta _1$ , for $\beta _1:=k_1^2{{\mathfrak {K}}}_{\min }=4\pi ^2{{\mathfrak {K}}}_{\min }/\ell _1^2$ .

Third, slowly varying theory for cylindrical domains, akin to Section 3.2, applies to assure us there exists a slow manifold homogenization of the mesoscale model in terms of a macroscale mean field $U(t,x_0)$ . Further, the slow manifold emerges exponentially quickly on a time scale of $1/\beta _1$ , and we may approximate the slow manifold via formal ‘slowly varying’ analysis of the mesoscale PDE (6.4).

Fourth, the construction from the mesoscale PDE ${{\mathfrak {U}}}_t={\mathfrak {D}}_x\{ {\mathfrak {K}}(x_1) {{\mathfrak {D}}}_x{\mathfrak {U}} \}$ to the macroscale homogenization is the specific case of Section 4 with $a_2=0$ and no variable $x_2$ . Consequently, corresponding results follow immediately, such as the slow manifold is (4.2a) with $a_2=0$ , and the macroscale homogenization is $U_t=U_{xx}$ to errors .

Here there is no point proceeding to higher order error in $\partial _x$ , because we chose the leading order mesoscale PDE at the start of this subsection. To obtain correct higher-order macroscale homogenizations, we would chose a correspondingly higher-order mesoscale homogenization.

7 Conclusion

This article further develops a novel theory and practice for mathematical/asymptotic homogenization to illuminate and resolve many issues in homogenizing heterogeneous PDEs such as (1.1). Section 2 introduces the technique of analysing an ensemble of all phase shifts, and extends it to multi-periodic cases, including the case of quasi-periodic. Future research may be able to extend the approach to interesting classes of random heterogeneity. The key to this approach is to be able to recast the ensemble as a system which is mathematically homogeneous in the macroscale.

Given the mathematically homogeneous macroscale ensemble, Section 3 applies recent developments of dynamical systems theory to newly establish the homogenization as a rigorous slow manifold of the given PDE system. No “ $\epsilon $ ” is required so the results hold for finite scale separation. No variational principles are required, so the approach applies to a wide range of physical problems, both dissipative and wave-like. The dynamical systems approach provides a framework that future research and applications may exploit to rigorously derive accurate initial conditions, boundary conditions, forcing projections, uncertainty evaluation and error bounds.

This dynamical systems framework provides a very practical way to construct homogenized models, as in the example of Section 4, and as in Appendix A. It illustrates the ease in using computer algebra to perform and check the homogenization (4.2). The method systematically and straightforwardly extends to higher order to derive so-called “second-order homogenization”: indeed, Appendix A.2 simply proceeds further to explicitly construct the “fourth-order homogenization” (4.3).

One advantage in the computer algebra is that, for some classes of heterogeneities, we can proceed to very high order in the analysis. Section 5 computed to 34th-order and uses the resulting generalized homogenization to quantify how ‘much smaller’ the microscale really has to be compared with the macroscale. Remarkably, the evidence is that, ideally, macroscale lengths can be resolved down to just twice the microscale length! Of course, in practice, there are many confounding aspects, but this factor of two is the mathematically ideal lower bound.

The spatio-temporal resolution of homogenizations could be improved, in future research, by multi-modal/multi-zonal homogenization of the heterogeneous microscale system by adapting approaches developed for shear dispersion [56, 69]. Such multi-modal/multi-zonal homogenization would be akin to Timoshenko–Ehrenfest beam theory [66, 71].

Lastly, Section 6 explores how the approach developed here includes the case when there are three separated length scales in the physical problem: a microscale, mesoscale and macroscale. The exploration shows how one can analyse such a system in one step from the combined micro and meso scales to the macroscale, or one can analyse the system in two steps from micro to meso, and then from meso to macro. Because sound modelling is transitive, the resulting homogenizations are the same. Further, if one’s ultimate aim is a macroscale spatial discretization to compute predictions, that is, the aim is not really the homogenized PDE, then instead of the two-step process of first homogenization and second spatial discretization, future research could do it all in one step from microscale heterogeneous PDE to macroscale discretization using a process like that introduced in a shear dispersion example [37].