# High-contrast random systems of PDEs:
homogenisation and spectral theory

Matteo Capoferri Maxwell Institute for Mathematical Sciences, Edinburgh _and_
Department of Mathematics, Heriot-Watt University, Edinburgh EH14 4AS, UK;
<EMAIL_ADDRESS>https://mcapoferri.com.  Mikhail Cherdantsev  School of
Mathematics, Cardiff University, Senghennydd Road, Cardiff CF24 4AG, UK;
<EMAIL_ADDRESS> Igor Velčić  Faculty of Electrical Engineering
and Computing, University of Zagreb, Unska 3, 10000 Zagreb, Croatia;
<EMAIL_ADDRESS>
(3rd September 2024)

###### Abstract

We develop a qualitative homogenisation and spectral theory for elliptic
systems of partial differential equations in divergence form with highly
contrasting (i.e., non uniformly elliptic) random coefficients. The focus of
the paper is on the behaviour of the spectrum as the heterogeneity parameter
tends to zero; in particular, we show that in general one doesn’t have
Hausdorff convergence of spectra. The theoretical analysis is complemented by
several explicit examples, showcasing the wider range of applications and
physical effects of systems with random coefficients, when compared with
systems with periodic coefficients or with scalar operators (both random and
periodic).

Keywords: high-contrast media, random media, stochastic homogenisation,
systems of PDEs.

2020 MSC classes:  primary 74S25, 74A40; secondary 35B27, 74Q15, 35B40, 60H15,
35P05.

###### Contents

  1. 1 Introduction
    1. Structure of the paper
  2. 2 Statement of the problem
    1. 2.1 Probabilistic and geometric setting
    2. 2.2 Our mathematical model
    3. List of notation
  3. 3 Main results
  4. 4 The limiting operator and its spectrum
  5. 5 Spectral convergence
    1. 5.1 Upper bound for the limiting spectrum
    2. 5.2 Lower bound on the limiting spectrum
  6. 6 Examples
    1. 6.1 Example 1: statistical micro-symmetries
    2. 6.2 Example 2: statistical macro-symmetries
    3. 6.3 Example 3: random scaling
  7. 7 Acknowledgements
  8. A Auxiliary materials
    1. A.1 Two-scale convergence
    2. A.2 The Ergodic Theorem
    3. A.3 Higher regularity for the periodic corrector
  9. References

##  1 Introduction

Homogenisation theory addresses the problem of determining effective
properties of heterogeneous materials (_composites_). Over time, this has
become a well-established branch of mathematics, and a large number of
classical monographs are now available. Yet, due to breadth and relevance of
the subject, mathematical and physical understanding of the properties of
composites remains at the forefront of modern research in science and
industry.

A very interesting class of problems — often referred to as _high-contrast
problems_ — is obtained by considering near-resonant inclusions dispersed in a
base medium (often referred to as “matrix”), in such a way that the
macroscopic properties of the material are accounted for by resonances
occurring at the microscopic scale. Mathematically, an important family of
such problems is modelled by second order elliptic operators of the form

$\mathcal{D}^{\varepsilon}=-\nabla\cdot a_{\varepsilon}\nabla\vspace{-2mm}$
(1.1)

defined on bounded or unbounded domains of $\mathbb{R}^{d}$, where the
coefficients $a_{\varepsilon}\in[L^{\infty}(\mathbb{R}^{d})]^{d\times d}$,
describing the physical properties of a material with scale of heterogeneity
$\varepsilon$, are of order $\varepsilon^{2}$ in the inclusions and of order
$1$ in the surrounding base medium. This critical _double-porosity_ type
scaling results in “micro-resonances”, which ultimately gives rise to numerous
novel effects that are not present in the classical (uniformly elliptic)
setting. One such effect, and possibly the most important, is the band-gap
structure of the spectrum, which makes high-contrast problems particularly
relevant from the point of view of applications.

One of the mathematical challenges of working with high-contrast media is that
the family of operators (1.1) is not uniformly elliptic as $\varepsilon\to 0$,
which, in turn, results in a loss of compactness in the limit. In pioneering
work, Allaire [2] showed, in the periodic setting, that the limiting
$\mathcal{D}^{\mathrm{hom}}$, analogue of the homogenised operator, possesses
a _two-scale_ nature, that is, it accounts for both macroscopic and
microscopic properties, coupled in a particular way. Zhikov [41, 42] later
defined and analysed the two-scale limit operator, and established spectral
convergence. By decoupling the macroscopic and microscopic scales in the two-
scale spectral problem $\mathcal{D}^{\mathrm{hom}}u=\lambda u$, one obtains a
spectral problem for the macroscopic component
$\mathcal{D}^{\mathrm{hom}}_{\mathrm{mac}}u_{1}(x)=\beta(\lambda)u_{1}(x)$
with highly nonlinear dependence on the spectral parameter $\lambda$ expressed
via _Zhikov’s $\beta$-function_. The latter, being defined explicitly in terms
of the resonant frequencies of the inclusions, characterises the macroscopic
behaviour of the composite by virtue of capturing length-scale interactions.
In particular, the spectrum of the limit operator is described by

$\sigma(\mathcal{D}^{\mathrm{hom}})=\\{\lambda:\ \beta(\lambda)\geq
0\\}\cup\\{\text{eigenfrequencies of the inclusions}\\},\vspace{-2mm}$

and, as inferred by the explicit representation of $\beta(\lambda)$, has
infinitely many gaps.

The literature in the high-contrast periodic setting is very rich. The seminal
results by Zhikov have been extended in numerous directions, encompassing both
scalar operators and systems of PDEs, rigorously demonstrating new non-trivial
properties, such as, for instance, memory effects and spatial nonlocality [14,
6]. Note that even in the absence of high-contrast in the constituents, by
carefully choosing the frequency of a Bloch wave type solution in relation to
the wavenumber along the fibres one can _artificially_ obtain a high-contrast
operator on the cross-section of the material [17]. Amongst recent
advancements, Cherednichenko, Ershova and Kiselev [12] established a
correspondence between homogenisation limits of micro-resonant periodic
composites and a class of time-dispersive media. Their interesting and new
techniques, based on norm-resolvent estimates, offer a recipe for the
construction of media with prescribed dispersive properties from periodic
composites whose individual components are non-dispersive, opening the way to
the possibility of interpreting high-contrast media as frequency-converting
devices. These scalar results were generalised to the case of systems in [13].
Very recently, significant advancement in the subject was made by Cooper,
Kamotski and Smyshlyaev: in [18] they developed a generic approach to study
abstract family of asymptotically degenerating _periodic_ variational problems
including, while going far beyond, high-contrast highly oscillatory PDEs.
Their asymptotic analysis yields uniform operator-type error estimate, and
their novel approach provides approximations of the spectra of the associated
spectral problems in terms of the spectrum of a certain “bivariate” operator.
The latter arises as an abstract generalisation of the two-scale limiting
operators traditionally obtained in the high-contrast literature.

Since fabricating perfectly periodic composites may be challenging from a
practical point of view, it is natural to ask whether the properties of high-
contrast materials desirable in applications (e.g., the band-gap structure of
the spectrum) persist, in a quantitatively controllable manner, when the
geometry of the inclusions, the coefficients of our operator, or both involve
an element of randomness. Stochastic homogenisation, qualitatively set out in
[30, 34, 31], is nowadays a very lively field of research, with a big
community concerned mainly with uniformly elliptic (non-high contrast)
operators. A most challenging and delicate aspect of the stochastic theory is
to obtain growth estimates in the appropriate norms for the _corrector_ , a
quantity that appears in the first-order term of the asymptotic expansion
ansatz for the homogenisation problem. Comprehensive _quantitative_ results
were first obtained by Gloria, Neukamm and Otto in the early 2010’s in a
series of ground-breaking papers [20, 21, 22] in which, building upon earlier
results by Yurinskii [40] and Naddaf–Spencer [32], they proved that the
corrector is essentially of bounded growth. Similar results were obtained
using a different (variational) argument by Armstrong, Smart, Kuusi and
Mourrat, see e.g. [5], [4] and references in the latter. These developments
prompted an intensive research activity, which resulted in a number of new
elegant ideas whose potential has not yet been fully explored.

For _high-contrast_ operators with random coefficients the picture is
completely different, in that much less is known and even the _qualitative_
theory is still under development. Indeed, moving from high-contrast periodic
to high-contrast stochastic PDEs involves a step change in both conceptual and
technical difficulty. Stochastic two-scale resolvent convergence for high-
contrast scalar operators on bounded domains and Hausdorff convergence of
spectra were established by Cherdantsev, Cherednichenko and Velčić in [9],
building upon earlier results by Zhikov and Piatnitski [44]. While stochastic
two-scale convergence works equally well in unbounded domains, the same
authors showed [10] that in the whole space $\mathbb{R}^{d}$ the lack of
compactness combined with the presence of infinitely many areas with
“atypical” distributions of inclusions allows one to construct quasimodes for
the random operator $\mathcal{D}^{\varepsilon}$ “escaping to infinity”,
namely, weakly converging to zero, as $\varepsilon\to 0$. The latter entails
that generically the spectrum $\sigma(\mathcal{D}^{\mathrm{hom}})$ of the two-
scale limiting operator, characterised by a stochastic version of
$\beta(\lambda)$, is a proper subset of the limiting spectrum
$\lim_{\varepsilon\to 0}\sigma(\mathcal{D}^{\varepsilon})$, in turn
characterised by $\beta_{\infty}(\lambda)$, a “local” modification of the
stochastic $\beta(\lambda)$. In fact, it can be shown that the limiting
spectrum has infinitely many gaps for a range of physically meaningful
examples [10, Section 4.6]. Also in the scalar setting, the authors of the
current paper studied defect modes for operators with random coefficients [8].
Finally, we should mention the paper [3], in which the authors study the large
time behaviour of a Markov process associated with a symmetric diffusion in a
high-contrast random environment and characterise the limit semigroup and the
limit process under the diffusive scaling. Although there are certain
similarities with [10], the setting of [3] is somewhat less general and the
focus quite different.

The properties of high-contrast _systems_ with random coefficients are, to
date, completely unexplored. Already in the periodic setting, high-contrast
systems exhibit additional novel effects, compared with scalar operators. For
example, in the macroscopic component of the limiting spectral problem for
systems the $\beta$-function is replaced by a $\bm{\beta}$-matrix, which can
be thought of as an “effective density”. It has been shown that for high-
contrast periodic systems the $\bm{\beta}$-matrix is frequency dependent, and
may be anisotropic as well as non-positive definite [7, 6, 43], and even
wavenumber dependent [36]. What determines the structure of the spectrum is
the signature of the $\bm{\beta}$-matrix, which no longer has just a binary
nature (positive/negative), thus resulting in a richer (and more complex)
theory. Indeed, on top of “strong” bands (propagation permitted) and gaps
(propagation completely forbidden), one can have the intermediate situation of
“weak” bands. This means that the number of macroscopic propagating modes may
vary with frequency [6] and, in the case of interconnected inclusions, with
direction as well [36].

Due to their multidimensional nature, systems also allow for the intermediate
situation of _partial degeneracy_ (or _partial high contrast_), in which the
coefficients in the inclusions pick up an additional term of order $1$ which
is non-degenerate in certain directions. Examples of this are, for instance,
linear elasticity with inclusions hard in compression and soft in shear or
Maxwell’s equations with high contrast in the electric permittivity. The
introduction of partial high contrast leads to a constrained microscopic
kinematics, which needs to be accounted for in the homogenisation process [16,
11]. Typically, one needs additional assumptions on the coefficients to ensure
the well-posedness of the corrector equation in the space of admissible
microscopic two-scale limit fields. A general theory for high-contrast
partially degenerate PDEs in the periodic setting was developed by Kamotski
and Smyshlyaev in [28].

Given the above rich picture in the periodic setting alongside the recent
developments [9, 10, 8] for operators with random coefficients, it is very
natural and timely to ask what new phenomena are brought about by allowing
high-contrast _systems_ to have _random_ coefficients. Due to the additional
freedom given by the random setting, one could examine scenarios precluded in
the periodic or scalar settings, such as having a mixture of fully degenerate
and partially degenerate inclusions or allowing the direction of partial
degeneracy to vary randomly from inclusion to inclusion, to name just a couple
of examples. Moreover, it may be possible to devise random media with spectral
properties highly desirable in applications, such as eigenvalues embedded into
the weak spectral bands.

In the current paper we carry out the first step of this programme, setting
out a general homogenisation and spectral theory for high-contrast systems
with random coefficients, thus laying firm theoretical foundations for a more
specialised future analysis of the range of problems outlined above. In
particular, random systems with partial degeneracy will be the subject of a
separate paper.

### Structure of the paper

Our paper is structured as follows.

In Section 2 we introduce our geometric and probabilistic settings, and
outline the problem we intend to study. In the end of Section 2 we provide a
list of symbols and notation recurring throughout the paper, with precise
pointers to where each of them is defined, to facilitate the reading.

In Section 3 we summarise the main results of our paper, with forward
references to the key theorems and propositions, and elaborate on the novelty
and differences with respect to the existing literature.

Section 4 is concerned with the examination of the two-scale limiting operator
$\mathcal{A}^{\mathrm{hom}}$. After defining it, we study its spectrum and the
relation thereof with the spectrum of $\mathcal{A}^{\varepsilon}$.

In Section 5 we study the limiting spectrum $\lim_{\varepsilon\to
0}\sigma(\mathcal{A}^{\varepsilon})$. In Subsection 5.1 we prove an upper
bound in terms of a certain set $\mathcal{G}$, which is shown in Subsection
5.2 to be also a lower bound, under additional assumptions on the finiteness
of the range of dependence.

Section 6 is devoted to several explicit examples, illustrating concretely the
role of the various quantities appearing in Sections 4 and 5, and showcasing
the potential of the stochastic setting for high-contrast systems, compared to
the periodic case.

The paper is complemented by one appendix, Appendix A, containing three sub-
appendices with auxiliary technical materials.

##  2 Statement of the problem

###  2.1 Probabilistic and geometric setting

We work in Euclidean space $\mathbb{R}^{d}$, $d\geq 2$, equipped with the
Lebesgue measure. Given a measurable set $A\subset\mathbb{R}^{d}$, we denote
by $|A|$ its Lebesgue measure, by $\overline{A}$ its closure, by
$\mathbbm{1}_{A}$ its characteristic function. We adopt the notation

$\displaystyle\square_{x}^{L}:=[-L/2,L/2]^{d}+x,\qquad x\in\mathbb{R}^{d},\
L>0.$ (2.1)

We define a cut-off function $\eta\in C_{0}^{\infty}(\square_{0}^{1})$ such
that

$0\leq\eta(x)\leq 1\quad\mbox{ and
}\quad\left.\eta\right|_{\square_{0}^{1/2}}=1,$ (2.2)

and put

$\eta_{L}(x):=\eta(x/L).$ (2.3)

Clearly, we have

$|\nabla\eta_{L}(x)|\leq\frac{C}{L},\quad\forall x\in\mathbb{R}^{d}.$ (2.4)

Function spaces with values in $\mathbb{R}^{d}$ will be denoted by bold
letters. For example,

$\displaystyle\mathbf{L}^{2}(\mathbb{R}^{d}):=[L^{2}(\mathbb{R}^{d})]^{d},\qquad\mathbf{W}^{1,p}(\mathbb{R}^{d}):=[W^{1,p}(\mathbb{R}^{d})]^{d},\qquad\mathbf{C}_{0}^{\infty}(\mathbb{R}^{d}):=[C_{0}^{\infty}(\mathbb{R}^{d})]^{d}\,.$
(2.5)

Here $W^{k,p}(\mathbb{R}^{d})$ is the usual Sobolev space of functions
$p$-integrable together with their first $k$ partial derivatives. Similarly,
we shall denote by bold letters vector-valued functions. In the notation for
norms, we will often write
$\|\cdot\|_{\mathbf{L}^{p}}:=\|\cdot\|_{\mathbf{L}^{p}(\mathbb{R}^{d})}$ and
$\|\cdot\|_{\mathbf{W}^{1,p}}:=\|\cdot\|_{\mathbf{W}^{1,p}(\mathbb{R}^{d})}$.
Furthermore, we will write
$\mathbf{H}^{k}(\mathbb{R}^{d}):=\mathbf{W}^{k,2}(\mathbb{R}^{d})$.

Our probability space $(\Omega,\mathcal{F},P)$ is defined as follows.

We define $\Omega$ to be the set of all possible collections of randomly
distributed inclusions in $\mathbb{R}^{d}$, that is, elements
$\omega\in\Omega$ are subsets of $\mathbb{R}^{d}$ satisfying appropriate
geometric conditions. Namely, we essentially require that individual
inclusions are sufficiently regular, of comparable size and that they do not
come too close to one another. These assumptions, rigorously formalised below
as Assumptions 2.1 and 2.2, mainly serve two purposes: on the one hand, they
prevent the formation of clusters of inclusions (_percolation_); on the other
hand, they guarantee the validity of Korn–Sobolev extension theorems with
uniform constants.

######  Assumption 2.1.

There exist constants $\zeta$, $N$ and $M$ such that for all $\omega\in\Omega$
the set $\mathbb{R}^{d}\setminus\omega$ is connected, and $\omega$ can be
written as a disjoint union

$\omega=\bigsqcup_{k\in\mathbb{N}}\omega^{k}$

of sets $\omega^{k}$ (the inclusions) satisfying the following properties.

  1. (a)

For every $k\in\mathbb{N}$ the set $\omega^{k}$ is open and connected.

  2. (b)

There exists a family $\\{\mathcal{B}_{\omega}^{k}\\}_{k\in\mathbb{N}}$ of
bounded open subsets of $\mathbb{R}^{d}$ such that for every $k\in\mathbb{N}$
we have111Here and further on  $\operatorname{diam}A:=\sup_{t,s\in A}|t-s|.$

    1. (i)

$\omega\cap\mathcal{B}_{\omega}^{k}=\omega^{k}$;

    2. (ii)

$\mathcal{B}_{\omega}^{k}\setminus\overline{\omega^{k}}$ is Lipschitz;

    3. (iii)

$\operatorname{diam}\omega^{k}<\frac{1}{2}$ and

$\omega^{k}-D_{\omega}^{k}-d_{1/4}\subset\mathcal{B}_{\omega}^{k}-D_{\omega}^{k}-d_{1/4}\subset[-\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}]^{d},$

where $D_{\omega}^{k}\in\mathbb{R}^{d}$ is defined by

$(D_{\omega}^{k})_{j}:=\inf_{x\in\omega^{k}}x_{j}$ (2.6)

and $d_{1/4}:=(1/4,\ldots,1/4)^{T}\in\mathbb{R}^{d}$.

######  Assumption 2.2.

There exist an integer $N\in\mathbb{N}$, triples
$\\{(P_{1}^{j},P_{2}^{j},P_{3}^{j}),\ j=1,\dots,N\\}$ of (open) Lipschitz
domains in $\mathbb{R}^{d}$, and universal positive constants $c$, $C$ and
$\widetilde{C}$ with the following properties.

  1. (i)

For all $j\in\\{1,\dots,N\\}$ we have

$\overline{P_{1}^{j}}\subset P_{2}^{j},\qquad\overline{P_{2}^{j}}\subset
P_{3}^{j}\,,$

and $P_{2}^{j}\setminus\overline{P_{1}^{j}}$ is Lipschitz.

  2. (ii)

For every $\omega\in\Omega$ and $k\in\mathbb{N}$, there exist

    1. (a)

$j=j(k)\in\\{1,\dots,N\\}$,

    2. (b)

a Lipschitz open neighbourhood $\widetilde{\mathcal{B}_{\omega}^{k}}$ of
$\mathcal{B}_{\omega}^{k}$ and

    3. (c)

a $C^{2}$ diffeomorphism
$\varphi_{k}:P_{1}^{j}\to\widetilde{\mathcal{B}_{\omega}^{k}}$ such that

$\varphi_{k}(P_{3}^{j}\setminus
P_{2}^{j})=\widetilde{\mathcal{B}_{\omega}^{k}}\setminus\mathcal{B}_{\omega}^{k},\qquad\varphi_{k}(P_{2}^{j}\setminus
P_{1}^{j})=\mathcal{B}_{\omega}^{k}\setminus\omega^{k},\qquad\varphi_{k}(P_{1}^{j})=\omega^{k},$

and

$c|x-y|\leq|\varphi_{k}(x)-\varphi_{k}(y)|\leq
C|x-y|,\qquad\|\varphi_{k}\|_{C^{2}(P_{3}^{j})}\leq\widetilde{C}\,.$

######  Remark 2.3.

Assumptions 2.1 and 2.2 warrant a number of remarks.

  1. (i)

Assumption 2.1 requires inclusions to fit within the hypercube
$[-\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}}]^{d}$, after a suitable rigid
translation (see, in particular, item (b)(iii)). The bounds on
$\mathcal{B}_{\omega}^{k}$ and $\omega^{k}$ are merely conventional and, at
the same time, unimportant: one can recover our setting from more general
boundedness assumptions on the inclusions by elementary scaling arguments.

  2. (ii)

The papers [10, 8] feature an analogue of Assumption 2.1 with Lipschitzness
replaced by the stronger assumption of minimal smoothness [37, Chapter VI,
Section 3.3]. Indeed, minimal smoothness was sufficient there to ensure the
uniformity of the constants in the relevant inequalities needed in the proofs.
In the current paper, Lipschitzness is enough for Theorem 2.7, but since we
will be dealing with symmetric gradients — and, hence, Korn-type inequalities
— stronger conditions (such as, e.g., Assumption 2.2) are needed to close the
argument.

  3. (iii)

Let us emphasise that our regularity assumptions on the inclusions are _not_
optimal. We need the regularity properties from Assumption 2.2 in exactly two
places: for the uniformity of Korn’s constant given by Proposition 2.8, and in
the proof of Theorem A.4 on the higher regularity of corrector, which, in
turn, is used in the proof of Theorem 5.7 — see (5.46), (5.57) and (5.67).

As far as the former is concerned, our assumptions are sufficient (although
not necessary) to ensure the uniformity of the Korn constant on the domain.
Note that Assumption 2.2 is a considerably stronger than minimal smoothness
assumption adopted for scalar operators in [10, 8] and, hence, takes us quite
far from optimality. We should like to point out, however, that the study of
the dependence of Korn’s inequality on the domain and its uniformity with
respect to certain classes of domains is an active area of current research
(see, e.g., the relatively recent result in [23] for star-shaped domains); in
this perspective, pursuing optimal conditions for the uniformity of the Korn
constant in our setting goes beyond the scope of the current paper.

As for the latter, in Theorem A.4 one needs a bit more than just a uniform
Korn constant. Indeed, by examining the argument in the proof of Theorem A.4
one realised that a uniform _local_ Korn inequality is required. Here is where
the $C^{2}$ equivalence of the inclusions to a _finite_ number of shapes plays
a role.

We define $\mathcal{F}$ to be the $\sigma$-algebra on $\Omega$ generated by
the mappings $\pi_{q}:\Omega\to\\{0,1\\}$, where $q\in\mathbb{Q}^{d}$ and

$\pi_{q}(\omega):=\mathbbm{1}_{\omega}(q).$

That is, $\mathcal{F}$ is the smallest $\sigma$-algebra that makes the
mappings $\pi_{q}$, $q\in\mathbb{Q}^{d}$, measurable222Under Assumption 2.1,
it is easy to see that this implies that the mappings
$\pi_{x}:\Omega\to\\{0,1\\}$ for $x\in\mathbb{R}^{d}$ are also measurable. We
postulate it only for $q\in\mathbb{Q}^{d}$ to ensure that $\mathcal{F}$ is
countably generated, which in turn implies that the spaces $L^{p}(\Omega)$ are
separable for $1\leq p<\infty$.. There is a natural group action of
$\mathbb{R}^{d}$ on $\Omega$, and hence on $\mathcal{F}$, given by
translations. Namely, for every $y\in\mathbb{R}^{d}$ the translation map
$T_{y}:\Omega\to\Omega$ acts on $\omega\in\Omega$ as

$\omega\mapsto T_{y}\omega=\\{z-y\ |\ z\in\omega\\}\subset\mathbb{R}^{d}.$
(2.7)

We equip $(\Omega,\mathcal{F})$ with a probability measure $P$ assumed to be
invariant under translations, i.e., we assume that $P(T_{y}F)=P(F)$ for every
$F\in\mathcal{F}$ and $y\in\mathbb{R}^{d}$, where $T_{y}F:=\bigcup_{\omega\in
F}T_{y}\omega$. It is not hard to see that $(T_{y})_{y\in\mathbb{R}^{d}}$
satisfies the following properties:

  1. (a)

$T_{y_{1}}\circ T_{y_{2}}=T_{y_{1}+y_{2}}\ $ for all
$y_{1},y_{2}\in\mathbb{R}^{d}$, where $\circ$ stands for composition;

  2. (b)

the map $\mathcal{T}:\mathbb{R}^{d}\times\Omega\to\Omega,\ $ $(y,\omega)\to
T_{y}\omega$ is measurable with respect to the standard $\sigma$-algebra on
the product space induced by $\mathcal{F}$ and the Borel $\sigma$-algebra on
$\mathbb{R}^{d}$.

Finally, we assume the translation group action $(T_{y})_{y\in\mathbb{R}^{d}}$
to be ergodic, i.e. if an element $F\in\mathcal{F}$ satisfies

$P((T_{y}F\cup F)\setminus(T_{y}F\cap F))=0\quad\text{for all}\quad
y\in\mathbb{R}^{d},$

then $P(F)\in\\{0,1\\}$. This will allow us to use the classical Ergodic
Theorem, which we recall as Theorem A.2 in Appendix A.2, for the reader’s
convenience.

We adopt the standard notation

$\displaystyle\mathbb{E}[\overline{\mathbf{f}}]:=\int_{\Omega}\overline{\mathbf{f}}(\omega)\,dP(\omega)$
(2.8)

and we denote by $\mathbf{L}^{p}(\Omega)=[L^{p}(\Omega)]^{d}$ the spaces of
vector-valued $p$-integrable functions in $(\Omega,\mathcal{F},P)$. Since
$\mathcal{F}$ is, clearly, countably generated, $\mathbf{L}^{p}(\Omega)$,
$1\leq p<\infty$, is separable. Given
$\overline{\mathbf{f}}:\Omega\to\mathbb{R}^{d}$, we denote by
$\mathbf{f}(y,\omega):=\overline{\mathbf{f}}(T_{y}\omega)$ its _realisation_
or _stationary extension_. Note that if
$\overline{\mathbf{f}}\in\mathbf{L}^{p}(\Omega)$, then $\mathbf{f}\in
L^{p}_{\mathrm{loc}}(\mathbb{R}^{d};\mathbf{L}^{p}(\Omega))$ [26, Chapter 7].
In the current paper we are mostly concerned with the case $p=2$.

######  Notation 2.4.

Throughout our paper, unless otherwise stated, we will denote with an overline
functions on $\Omega$ and we will remove the overline to denote the
corresponding realisation (stationary extension). We will reserve the letter
$y$ for the (stationary) extension variable. So, for example:

$\displaystyle\overline{\mathbf{f}}=\overline{\mathbf{f}}(\omega),\quad\mathbf{f}=\mathbf{f}(y,\omega):=\overline{\mathbf{f}}(T_{y}\omega),\quad\mathbb{E}[\mathbf{f}]:=\mathbb{E}[\overline{\mathbf{f}}].$
(2.9)

Functions on $\Omega$ may additionally depend on other variables, not
necessarily in a stationary manner. For example, if
$\overline{\mathbf{f}}(x,\omega):\mathbb{R}^{d}\times\Omega\to\mathbb{R}^{d}$,
then $\mathbf{f}(x,y,\omega):=\overline{\mathbf{f}}(x,T_{y}\omega)$ and
$\mathbb{E}[\mathbf{f}]=\mathbb{E}[\overline{\mathbf{f}}(x,\cdot)]$ (note that
the dependence on the non-stationary variable $x$ remains upon taking the
expectation).

Let $\\{e_{j}\\}_{j=1}^{d}$ be the standard basis of $\mathbb{R}^{d}$. We
define the Sobolev spaces $\mathbf{H}^{s}(\Omega)$, $s\in\mathbb{N}$, as

$\mathbf{H}^{s}(\Omega):=\left\\{\overline{\mathbf{f}}\in\mathbf{L}^{2}(\Omega)\
|\ \mathbf{f}\in
H^{s}_{\mathrm{loc}}(\mathbb{R}^{d};\mathbf{L}^{2}(\Omega))\right\\}.$ (2.10)

Moreover, we define

$\mathbf{H}^{\infty}(\Omega):=\bigcap_{s\in\mathbb{N}}\mathbf{H}^{s}(\Omega)$

and

$\mathbf{C}^{\infty}(\Omega):=\left\\{\overline{\mathbf{f}}\in\mathbf{H}^{\infty}(\Omega)\
|\
\overline{\partial^{\bm{\alpha}}\mathbf{f}}\in\mathbf{L}^{\infty}(\Omega)\text{
for every multi-index
}\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{N}_{0}^{d}\right\\}.$
(2.11)

We should like to emphasise that, given
$\overline{\mathbf{f}}\in\mathbf{H}^{s}(\Omega)$, for every multi-index
$\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{N}_{0}^{d}$,
$\sum_{l=1}^{d}\alpha_{l}\leq s$, the quantity
$\partial^{\bm{\alpha}}\mathbf{f}$ is the stationary extension of a vector
field in $\mathbf{L}^{2}(\Omega)$. Accordingly, the quantity
$\overline{\partial^{\bm{\alpha}}\mathbf{f}}$ in the RHS of (2.11) is to be
understood as the random variable whose stationary extension is
$\partial^{\bm{\alpha}}\mathbf{f}$. Here
$\partial^{\bm{\alpha}}:=\partial^{\alpha_{1}}_{y_{1}}\cdots\partial^{\alpha_{d}}_{y_{d}}$.

In view of the above one defines the norm on (2.10) as

$\left\|\overline{\mathbf{f}}\right\|^{2}_{\mathbf{H}^{s}(\Omega)}:=\sum_{|\bm{\alpha}|\leq
s}\left\|\overline{\partial^{\bm{\alpha}}\mathbf{f}}\right\|^{2}_{\mathbf{L}^{2}(\Omega)}.$

######  Remark 2.5.

Observe that the definition of probabilistic Sobolev spaces retraces the
classical one. We refrain from introducing the stationary differential
calculus on $\Omega$ more formally, as it will not be needed in this paper. We
refer the interested reader to [19, Appendices A.2 and A.3] for further
details.

Finally, we define

$\mathbf{L}^{2}_{0}(\Omega):=\left\\{\overline{\mathbf{f}}\in\mathbf{L}^{2}(\Omega)\
|\ \left.\mathbf{f}(\cdot,\omega)\right|_{\mathbb{R}^{d}\setminus\omega}=0\
\text{for all}\ \omega\in\Omega\right\\},$
$\mathbf{H}^{s}_{0}(\Omega):=\left\\{\overline{\mathbf{f}}\in\mathbf{H}^{s}(\Omega)\
|\ \left.\mathbf{f}(\cdot,\omega)\right|_{\mathbb{R}^{d}\setminus\omega}=0\
\text{for all}\ \omega\in\Omega\right\\},$

and

$\mathbf{C}^{\infty}_{0}(\Omega):=\left\\{\overline{\mathbf{f}}\in\mathbf{C}^{\infty}(\Omega)\
|\ \left.\mathbf{f}(\cdot,\omega)\right|_{\mathbb{R}^{d}\setminus\omega}=0\
\text{for all}\ \omega\in\Omega\right\\}.$

The latter are spaces of functions in $\mathbf{L}^{2}(\Omega)$,
$\mathbf{H}^{s}(\Omega)$ and $\mathbf{C}^{\infty}(\Omega)$, respectively,
whose realisations vanish identically outside of the inclusions.

The above function spaces enjoy the following properties: (i)
$\mathbf{H}^{s}(\Omega)$ is a separable Hilbert space; (ii)
$\mathbf{H}^{\infty}(\Omega)$ is dense in $\mathbf{L}^{2}(\Omega)$; (iii)
$\mathbf{C}^{\infty}(\Omega)$ is dense in $\mathbf{L}^{p}(\Omega)$, $1\leq
p<\infty$; (iv) $\mathbf{C}^{\infty}(\Omega)$ is dense in
$\mathbf{H}^{s}(\Omega)$. Furthermore, we have at our disposal the following
result, whose proof may be found in [9, 10].

######  Theorem 2.6.

Under Assumption 2.1, $\mathbf{C}_{0}^{\infty}(\Omega)$ is dense in
$\mathbf{L}_{0}^{2}(\Omega)$ and in $\mathbf{H}^{1}_{0}(\Omega)$ with respect
to $\|\cdot\|_{\mathbf{L}^{2}(\Omega)}$ and
$\|\cdot\|_{\mathbf{H}^{1}(\Omega)}$, respectively.

The underlying motivation to much of our geometric assumptions is the fact
that they are conducive to good “extension properties”.

######  Theorem 2.7 (Extension theorem [33, Lemma 4.1]).

Let $\mathsf{C}_{0}$ be a constant fourth-order tensor satisfying the symmetry
and ellipticity conditions (2.16) and (2.17) below. Under Assumption 2.1, for
every $p\geq 1$ there exists a bounded linear extension operator
$E_{k}:\mathbf{W}^{1,p}(\mathcal{B}_{\omega}^{k}\setminus\omega^{k})\to\mathbf{W}^{1,p}(\mathcal{B}_{\omega}^{k})$
such that for every
$\mathbf{u}\in\mathbf{W}^{1,p}(\mathcal{B}_{\omega}^{k}\setminus\omega^{k})$
the extension
$\widetilde{\mathbf{u}}=E_{k}\mathbf{u}\in\mathbf{W}^{1,p}(\mathcal{B}_{\omega}^{k})$
satisfies

$\widetilde{\mathbf{u}}=\mathbf{u}\qquad\text{in}\qquad\mathcal{B}_{\omega}^{k}\setminus\omega^{k},$
(2.12)
$\operatorname{div}\mathsf{C}_{0}\nabla\widetilde{\mathbf{u}}=0\qquad\text{in}\qquad\omega^{k},$
(2.13)
$\|\widetilde{\mathbf{u}}\|_{\mathbf{W}^{1,p}(\mathcal{B}_{\omega}^{k})}\leq
C\,\left(\|\mathbf{u}\|_{\mathbf{L}^{p}(\mathcal{B}_{\omega}^{k}\setminus\omega^{k})}+\|\operatorname{sym}\nabla\mathbf{u}\|_{\mathbf{L}^{p}(\mathcal{B}_{\omega}^{k}\setminus\omega^{k})}\right),$
(2.14)
$\|\operatorname{sym}\nabla\widetilde{\mathbf{u}}\|_{\mathbf{L}^{p}(\mathcal{B}_{\omega}^{k})}\leq
C\,\|\operatorname{sym}\nabla\mathbf{u}\|_{\mathbf{L}^{p}(\mathcal{B}_{\omega}^{k}\setminus\omega^{k})},$
(2.15)

where the constants $C$ in (2.14), (2.15) depend on $\omega$, $k$, $p$, and
the ellipticity constant of $\mathsf{C}_{0}$.

In plain English, the above extension theorem tells us that we can extend
$\mathbf{H}^{1}$ vector-functions into the inclusions in a
$\mathsf{C}_{0}$-harmonic manner, whilst controlling gradient and symmetric
gradient of the extended function by gradient and symmetric gradient of the
original function.

######  Proposition 2.8.

Under Assumption 2.2, the constants $C$ in (2.14), (2.15) is independent of
$\omega$ and $k$ almost surely.

###### Proof.

The claim is obtained arguing by contradiction and retracing the steps of the
proof [38, Lemma 1] with account of [38, Remark 6]. See also [25] for the
treatment of transformations of Lipschitz domains under sufficiently smooth
mappings. ∎

######  Remark 2.9.

The constants $C$ in (2.14), (2.15) depends on $p$. However, in our paper we
will only use Theorem 2.7 for $p=2$ and for one $p>2$. Therefore, for all
practical purposes we can suppress the dependence of $C$ on $p$, and treat the
constant as if it were uniform with respect to $p$.

###  2.2 Our mathematical model

Let $\mathsf{C}_{j}\in\mathbb{R}^{d^{4}}$, $j=0,1$, be (constant) fourth-order
tensors satisfying the following properties:

  1. (a)

_Symmetry_ :

$(\mathsf{C}_{j})_{\alpha\beta\mu\nu}=(\mathsf{C}_{j})_{\mu\nu\alpha\beta},\qquad(\mathsf{C}_{j})_{\alpha\beta\mu\nu}=(\mathsf{C}_{j})_{\alpha\beta\nu\mu},$
(2.16)

for all $\alpha,\beta,\mu,\nu=1,2,\ldots,d$;

  2. (b)

_Ellipticity_ : there exists $c>0$ such that

$(\mathsf{C}_{j})_{\alpha\beta\mu\nu}\xi_{\alpha\beta}\xi_{\mu\nu}\geq
c\,|\xi|^{2},\qquad\forall\xi\in\mathbb{R}^{d^{2}},\
\xi_{\alpha\beta}=\xi_{\beta\alpha},\qquad j=0,1.$ (2.17)

Here and further on, we adopt Einstein’s summation convention over repeated
indices, unless otherwise stated.

For every $0<\varepsilon<1$, for each random set of inclusions $\omega$ we
partition $\mathbb{R}^{d}$ into two regions, up to a set of measure zero: the
_inclusions_

$\mathcal{I}^{\varepsilon}(\omega):=\varepsilon\omega$ (2.18)

and the _matrix_

$\mathcal{M}^{\varepsilon}(\omega):=\mathbb{R}^{d}\setminus\overline{\mathcal{I}^{\varepsilon}(\omega)}.$
(2.19)

We define the high-contrast tensor field of coefficients as

$\mathsf{C}^{\varepsilon}=\mathsf{C}^{\varepsilon}(x,\omega):=\mathbbm{1}_{\mathbb{R}^{d}\setminus\varepsilon\omega}\,\mathsf{C}_{1}+\,\mathbbm{1}_{\varepsilon\omega}\,\varepsilon^{2}\mathsf{C}_{0}.$
(2.20)

Here the characteristic size of the microstructure is of order $\varepsilon$,
and the coefficients equal to $\mathsf{C}_{1}$ in the stiff matrix and to
$\varepsilon^{2}\,\mathsf{C}_{0}$ in the soft inclusions. This type of scaling
is often referred to as ‘double-porosity’.

######  Definition 2.10.

We define $\mathcal{A}^{\varepsilon}(\omega)$ to be the self-adjoint linear
operator in $\mathbf{L}^{2}(\mathbb{R}^{d})$ associated with the bilinear
form333Here and further on
$\mathsf{C}^{\varepsilon}\nabla\mathbf{u}\,\cdot\nabla\mathbf{v}=(\mathsf{C}^{\varepsilon})_{\alpha\beta\mu\nu}\partial_{\mu}\mathbf{u}_{\nu}\,\partial_{\alpha}\mathbf{v}_{\beta}.$

$\int_{\mathbb{R}^{d}}\mathsf{C}^{\varepsilon}\nabla\mathbf{u}\,\cdot\nabla\mathbf{v}=\int_{\mathbb{R}^{d}}\mathsf{C}^{\varepsilon}(\operatorname{sym}\nabla\mathbf{u})\,\cdot(\operatorname{sym}\nabla\mathbf{v}),\,\qquad\mathbf{u},\mathbf{v}\in\mathbf{H}^{1}(\mathbb{R}^{d}).$
(2.21)

The operator $\mathcal{A}^{\varepsilon}$ is a matrix operator acting on vector
fields.

The overall goal of the paper is to examine the spectral properties of the
operator $\mathcal{A}^{\varepsilon}$ for sufficiently small $\varepsilon$,
through the prism of (two-scale) stochastic homogenisation.

######  Remark 2.11.

Our setting can be generalised in a straightforward manner to cover the
following.

  1. (i)

One can consider operators acting on sections of more general trivial vector
bundles $[L^{2}(\mathbb{R}^{d})]^{m}$, $m\in\mathbb{N}$. When $m=1$, one
recovers the scalar theory.

  2. (ii)

The coefficients $\mathsf{C}_{0}$ and $\mathsf{C}_{1}$ can be allowed to
depend on $\omega$, provided the ellipticity constant in (2.17) is uniform in
$\omega$.

  3. (iii)

Throughout the paper, we set the material density to be equal 1, for
simplicity. Studying the problem at hand in the presence of material density
can be done analogously.

Upgrading our results so as to encompass (i)–(iii) only involves minor
adjustments to the proofs. However, we carry out our arguments in the slightly
more restrictive setting set out above for the sake of clarity, to avoid
inessential technical details.

### List of notation

Symbol | Description
---|---
$\overset{2}{\to}$ ($\stackrel{{\scriptstyle 2}}{{\rightharpoonup}}$) | Strong (resp. weak) stochastic two-scale convergence — Definition 4.3
$\square_{x}^{L}$ | Hypercube of size $L$ centred at $x$, see (2.1)
$\mathbbm{1}_{A}$ | Characteristic function of $A\subset\mathbb{R}^{d}$
$1_{\Omega}$ | Random function (4.3)
$|A|$ | Lebesgue measure of $A\subset\mathbb{R}^{d}$
$\overline{A}$ | Closure of $A\subset\mathbb{R}^{d}$
$\mathcal{A}^{\varepsilon}$ | Definition 2.10
$\mathcal{A}^{\mathrm{hom}}$ | Limiting two-scale operator — see Definition 4.1
$\mathcal{A}_{0}$ | Microscopic part of $\mathcal{A}^{\mathrm{hom}}$ — see Definition 4.1
$\mathcal{A}_{0}^{Q}$ | Operator $\mathcal{A}_{0}$ in $\mathbf{L}^{2}(Q)$ with Dirichlet b.c. — see Definition 4.2
$\mathcal{A}_{1}$ | Macroscopic part of $\mathcal{A}^{\mathrm{hom}}$ — see Definition 4.1
$\overline{\mathbf{b}_{\lambda}}$ | Matrix function (4.14), (4.15)
$\overline{\mathbf{b}_{\lambda}^{(i)}}$ |  $i$-th column of $\overline{\mathbf{b}_{\lambda}}$, see also 4.12 and (4.13)
$\mathcal{B}_{\omega}^{k}$ | Extension domain of $\omega^{k}$ — see Assumption 2.1
$\bm{\beta}(\lambda)$ | Zhikov $\bm{\beta}$-matrix (4.9)
$\beta(\lambda)$ | Largest eigenvalue of $\bm{\beta}(\lambda)$, see (4.21)
$\beta_{\infty}(\lambda)$ | Deterministic function given by Definition 5.1 almost surely
$\mathsf{C}_{0}$, $\mathsf{C}_{1}$ | Symmetric and elliptic fourth-order tensors, see (2.16), (2.17)
$\mathsf{C}^{\varepsilon}$ | Tensor field of coefficients (2.20)
$\mathsf{C}^{\mathrm{hom}}$ | Homogenised tensor field of coefficients (4.4)
$\mathbf{C}_{0}^{\infty}$, $\mathbf{H}^{k}$, $\mathbf{L}^{2}$, $\mathbf{W}^{1,p}$ | Bold letters denote functions spaces with values in $\mathbb{R}^{d}$, see (2.5)
$\mathcal{D}(\mathcal{A})$ | Domain of the opeartor $\mathcal{A}$
$D_{\omega}^{k}$ | Shift vector (2.6)
$\mathbf{e}_{i}$, $i\in\\{1,\dots,d\\}$ | Basis vectors in $\mathbb{R}^{d}$, $[\mathbf{e}_{i}]_{j}=\delta_{ij}$
$\mathbb{E}[\overline{\mathbf{f}}]$ | Expectation (2.8)
$\eta$ and $\eta_{L}$ | Cut-off (2.2) and its scaled version (2.3), respectively
$f(y,\omega):=\overline{f}(T_{y}\omega)$ | Stationary extension 2.9 of the random variable $\overline{f}$
$\mathcal{G}$ | Set (5.3)
$\mathbf{H}$ | Function space (4.1)
$\mathbf{H}^{1}_{\mathrm{per}}(Q)$ | Elements of $\mathbf{H}^{1}$ that are periodic with period $Q$
$I$ | Identity matrix
$\mathcal{I}^{\varepsilon}$ |  $\varepsilon$-scaled inclusions (2.18)
$\mathcal{M}^{\varepsilon}$ |  $\varepsilon$-scaled matrix (2.19)
$\sigma(\mathcal{A})$ | Spectrum of the operator $\mathcal{A}$
$T_{y}$ | Dynamical system (2.7) on $\Omega$ acting by translation
$\mathbf{V}$ | Function space (4.2)
$\bm{\mathcal{V}}^{2}_{\mathrm{pot}}$ | Potential vector fields with zero mean (4.5)
$\omega$ | Collection of inclusions ($\omega\subset\mathbb{R}^{d}$, $\omega\in\Omega$) — see Assumption 2.1
$\omega^{k}$ | Individual inclusion, connected component of $\omega$ — see Assumption 2.1
$(\Omega,\mathcal{F},P)$ | Probability space

##  3 Main results

In this section we summarise, for the convenience of the reader, the main
results of our paper.

Theorem 4.5 establishes that the family of operators
$\mathcal{A}^{\varepsilon}$ converges, as $\varepsilon\to 0$, to a limiting
operator $\mathcal{A}^{\mathrm{hom}}$ (see Definition 4.1). Our first main
result is Theorem 4.11, which characterises the spectrum of the limiting
operator $\mathcal{A}^{\mathrm{hom}}$. The latter is comprised of two parts:
the micro-resonances — i.e., the spectrum of a microscopic limiting operator
$\mathcal{A}_{0}$ — and the nonnegative real numbers $\lambda$ for which the
_matrix_ function $\bm{\beta}(\lambda)$ (4.19) has at least one nonnegative
eigenvalue.

Our second main result, Theorem 5.5, provides an upper bound for the limiting
spectrum $\lim_{\varepsilon\to 0}\sigma(\mathcal{A}^{\varepsilon})$ in terms
of the set $\mathcal{G}$ (5.3), in turn defined by a nonlocal _scalar_
function $\beta_{\infty}(\lambda)$, deterministic almost surely, which is
sensitive to areas of space with atypical arrangements of inclusions.

As it turns out, under the additional assumption of finite-distance
dependence, $\mathcal{G}$ is also an upper bound for the limiting spectrum.
This is the subject of our third main result, Theorem 5.7, which, in
conjunction with Theorem 5.5, yields $\lim_{\varepsilon\to
0}\sigma(\mathcal{A}^{\varepsilon})=\mathcal{G}$. Since, in general,
$\sigma(\mathcal{A}^{\mathrm{hom}})$ is a proper subset of $\mathcal{G}$, this
implies one does not have Hausdorff convergence of spectra, so that the
limiting operator $\mathcal{A}^{\mathrm{hom}}$ captures only partially
information about the spectrum of $\mathcal{A}^{\varepsilon}$ for arbitrary
small but finite $\varepsilon$. This is a distinctive feature of the
stochastic setting, already observed in [10] in the scalar case, which
persists for systems with random coefficients.

Lastly, our fourth main result is Proposition 6.1, which, together with the
other examples from Section 6, showcases the wider range of applications and
physical effects of stochastic systems, when compared with periodic ones, thus
paving the way for the further analysis outlined in Section 1.

##  4 The limiting operator and its spectrum

Consider the ‘two-scale’ spaces

$\mathbf{H}:=\mathbf{L}^{2}(\mathbb{R}^{d})+L^{2}(\mathbb{R}^{d};\mathbf{L}^{2}_{0}(\Omega)),$
(4.1)

which is naturally embedded in
$\mathbf{L}^{2}(\mathbb{R}^{d};\mathbf{L}^{2}(\Omega))$, and its dense
subspace

$\mathbf{V}:=\mathbf{H}^{1}(\mathbb{R}^{d})+L^{2}(\mathbb{R}^{d};\mathbf{H}^{1}_{0}(\Omega))\,.$
(4.2)

The spaces $\mathbf{H}$ and $\mathbf{V}$ are comprised of vector functions of
the form $\mathbf{u}_{1}(x)+\overline{\mathbf{u}}_{0}(x,\omega)$, where the
stationary extension $\mathbf{u}_{0}(x,y,\omega)$ of
$\overline{\mathbf{u}}_{0}(x,\omega)$ vanishes for $y$ outside of the set of
inclusions $\omega$.

Let us recall that a tensor field
$\mathbf{p}\in[\mathbf{L}^{2}_{\mathrm{loc}}(\mathbb{R}^{d})]^{d}$ is said to
be

  * •

_potential_ if there exists
$\bm{\varphi}\in\mathbf{H}^{1}_{\mathrm{loc}}(\mathbb{R}^{d})$ such that
$\mathbf{p}=\nabla\bm{\varphi}$;

  * •

_solenoidal_ if

$\int_{\mathbb{R}^{d}}\mathbf{p}\cdot\nabla\bm{\varphi}=0\qquad\forall\bm{\varphi}\in
C_{0}^{\infty}(\mathbb{R}^{d}).$

See, e.g., [26], for further details. By analogy, one defines a tensor field
$\overline{\bm{p}}\in[\mathbf{L}^{2}(\Omega)]^{d}$ to be potential
(respectively solenoidal) if for a.e. $\omega$ its stationary extension
$y\mapsto\mathbf{p}(y,\omega)\in[\mathbf{L}^{2}_{\mathrm{loc}}(\mathbb{R}^{d})]^{d}$
is such. The classical Weyl’s decomposition theorem holds for tensor fields in
$[\mathbf{L}^{2}(\Omega)]^{d}$ as well, see e.g. [26, Chapter 12.7], see also
[8, Theorem 2.14].

Let $1_{\Omega}:\mathbb{R}^{d}\times\Omega\to\\{0,1\\}$ be the stationary
random function defined in accordance with

$1_{\Omega}(y,\omega):=\mathbbm{1}_{\omega}(y)\,.$ (4.3)

Let $\mathsf{C}^{\mathrm{hom}}$ be the homogenised tensor of coefficients
defined in accordance with

$\mathsf{C}^{\mathrm{hom}}\xi\cdot\xi:=\inf_{\overline{\bm{\psi}}\in\bm{\mathcal{V}}^{2}_{\mathrm{pot}}(\Omega)}\mathbb{E}\left[(1-1_{\Omega})\,\mathsf{C}_{1}(\bm{\psi}+\xi)\cdot(\bm{\psi}+\xi)\right],\quad\forall\xi\in\mathbb{R}^{d^{2}},$
(4.4)

where

$\bm{\mathcal{V}}^{2}_{\mathrm{pot}}(\Omega):=\\{\overline{\bm{\psi}}\in[\mathbf{L}^{2}(\Omega)]^{d}\
|\ \bm{\psi}\text{ is potential and }\
\mathbb{E}[\overline{\bm{\psi}}]=0\\}\,.$ (4.5)

Consider the following bilinear forms:

$\mathbb{E}\left[\mathsf{C}_{0}\nabla_{y}\mathbf{u}\cdot\nabla_{y}\bm{\varphi}\right]\,,\qquad\overline{\mathbf{u}},\
\overline{\bm{\varphi}}\in\mathbf{H}^{1}_{0}(\Omega),$ (4.6)
$\int_{\mathbb{R}^{d}}\mathsf{C}^{\mathrm{hom}}\nabla\mathbf{u}\cdot\nabla\bm{\varphi}\,,\qquad{\mathbf{u}},\
{\bm{\varphi}}\in\mathbf{H}^{1}(\mathbb{R}^{d}),$ (4.7)

and

$\int_{\mathbb{R}^{d}}\mathsf{C}^{\mathrm{hom}}\nabla\mathbf{u}_{1}\cdot\nabla\bm{\varphi}_{1}+\int_{\mathbb{R}^{d}}\mathbb{E}\left[\mathsf{C}_{0}\nabla_{y}\mathbf{u}_{0}\cdot\nabla_{y}\bm{\varphi}_{0}\right]\,,\qquad\mathbf{u}_{1}+\overline{\mathbf{u}_{0}},\
\bm{\varphi}_{1}+\overline{\bm{\varphi}_{0}}\in\mathbf{V}\,.$ (4.8)

By the standard Korn’s inequality and its direct analogue for the space
$\mathbf{H}^{1}_{0}(\Omega)$ we see that the above bilinear forms are
coercive, hence, closed.

######  Definition 4.1 (Limiting two-scale operators).

We define the linear operators $\mathcal{A}_{0}$ in
$\mathbf{L}^{2}_{0}(\Omega)$, $\mathcal{A}_{1}$ in
$\mathbf{L}^{2}(\mathbb{R}^{d})$, and $\mathcal{A}^{\mathrm{hom}}$ in
$\mathbf{H}$ as the self-adjoint operators associated with the bilinear forms
(4.6), (4.7) and (4.8), respectively. We will refer to the operators
$\mathcal{A}_{0}$ and $\mathcal{A}_{1}$ as the micro- and macroscopic parts of
the limiting operator $\mathcal{A}^{\mathrm{hom}}$.

The operator $\mathcal{A}_{0}$ can formally be identified with the operator
$-\operatorname{div}_{y}\mathsf{C}_{0}\nabla_{y}$ acting on the stationary
extensions of the random variables from $\mathrm{dom}(\mathcal{A}_{0})$. For
later use, it is also convenient to introduce a “restriction” of this operator
to a single inclusion $\omega^{k}$ or, more generally, to a bounded domain in
$\mathbb{R}^{d}$.

######  Definition 4.2.

Given a bounded domain $Q\subset\mathbb{R}^{d}$, we denote by
$\mathcal{A}_{0}^{Q}$ the positive definite self-adjoint operator associated
with the bilinear form

$\displaystyle\int_{Q}\mathsf{C}_{0}\,\nabla\mathbf{u}\cdot\nabla\mathbf{v},\qquad\mathbf{u},\mathbf{v}\in\mathbf{H}^{1}_{0}(Q)\,.$

In other word, $\mathcal{A}_{0}^{Q}$ is nothing but the operator
$-\operatorname{div}\mathsf{C}_{0}\nabla$ in $\mathbf{L}^{2}(Q)$ with
Dirichlet boundary conditions.

Let $\mathcal{C}\subset\mathbf{C}^{\infty}(\Omega)$ be a countable dense
family of vector-functions in $\mathbf{L}^{2}(\Omega)$ (recall that the latter
is separable) and let $\Omega_{t}=\Omega_{t}(\mathcal{C})\subset\Omega$ be a
set of probability one such that the claim of the Ergodic Theorem A.2 holds
for all $\omega\in\Omega_{t}$ and $\overline{\mathbf{f}}\in\mathcal{C}$.
Elements of $\Omega_{t}$ are often referred to as _typical_ , hence the
subscript "t".

######  Definition 4.3 (Stochastic two-scale convergence [44]).

Let $\\{\mathbf{u}^{\varepsilon}\\}$ be a bounded sequence in
$\mathbf{L}^{2}(\mathbb{R}^{d})$. We say that $\\{\mathbf{u}^{\varepsilon}\\}$
weakly stochastically two-scale converges to
$\overline{\mathbf{u}}\in\mathbf{L}^{2}(\mathbb{R}^{d}\times\Omega)$ (for a
given $\omega_{0}\in{\Omega_{t}}$) and write
$\mathbf{u}^{\varepsilon}\overset{2}{\rightharpoonup}\overline{\mathbf{u}}$ if

$\lim_{\varepsilon\to
0}\int_{\mathbb{R}^{d}}\mathbf{u}^{\varepsilon}(x)\cdot\varphi(x)\mathbf{f}(x/\varepsilon,\omega_{0})\,dx=\mathbb{E}\left[\int_{\mathbb{R}^{d}}\overline{\mathbf{u}}\cdot\varphi\,\overline{\mathbf{f}}\right]\quad\forall\varphi\in
C^{\infty}_{0}(\mathbb{R}^{d}),\,\overline{\mathbf{f}}\in{\mathcal{C}}.$ (4.9)

We say that $\\{\mathbf{u}^{\varepsilon}\\}$ strongly stochastically two-scale
converges to
$\overline{\mathbf{u}}\in\mathbf{L}^{2}(\mathbb{R}^{d}\times\Omega)$ and write
$\mathbf{u}^{\varepsilon}\overset{2}{\to}\overline{\mathbf{u}}$ if it
satisfies (4.9) and

$\lim_{\varepsilon\to
0}\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}=\|\overline{\mathbf{u}}\|_{\mathbf{L}^{2}(\mathbb{R}^{d}\times\Omega)}.$

For the reader’s convenience, basic properties of two-scale convergence are
summarised in Appendix A.1.

######  Remark 4.4.

Observe that one can always add countably many elements to the family
$\mathcal{C}$ in the above definition, provided one updates the choice of the
set of full measure $\Omega_{t}$ accordingly. See also [44, 24] for further
details. We will tacitly make use of this freedom throughout the paper.
Indeed, one is always able to identify the limiting objects resulting from
two-scale convergence that appear in our paper by testing the relevant weak
identities with countable dense families of (vector-)functions with higher-
than-$L^{2}$ — typically, smooth — regularity.

The operator $\mathcal{A}^{\mathrm{hom}}$ is the stochastic two-scale limit of
$\mathcal{A}^{\varepsilon}$ as $\varepsilon\to 0$ in the sense of the
following theorem.

######  Theorem 4.5.

Let $\lambda<0$ and let
$\\{\mathbf{f}^{\varepsilon}\\}\subset\mathbf{L}^{2}(\mathbb{R}^{d})$,
$\|\mathbf{f}^{\varepsilon}\|_{\mathbf{L}^{2}}\leq C$, be a bounded sequence
of square integrable functions. Suppose
$\mathbf{f}^{\varepsilon}\stackrel{{\scriptstyle
2}}{{\rightharpoonup}}\mathbf{f}$
($\mathbf{f}^{\varepsilon}\overset{2}{\to}\mathbf{f}$). Then the solution
$\mathbf{u}^{\varepsilon}$ of the resolvent problem

$\mathcal{A}^{\varepsilon}\mathbf{u}^{\varepsilon}-\lambda\mathbf{u}^{\varepsilon}=\mathbf{f}^{\varepsilon}$
(4.10)

weakly (strongly) stochastically two-scale converges to some
$\mathbf{u}=\mathbf{u}_{1}+\overline{\mathbf{u}_{0}}\in V$,
$\mathbf{u}^{\varepsilon}\stackrel{{\scriptstyle
2}}{{\rightharpoonup}}\mathbf{u}$
($\mathbf{u}^{\varepsilon}\overset{2}{\to}\mathbf{u}$), almost surely.
Furthermore, $\mathbf{u}$ is the solution of

$\mathcal{A}^{\rm hom}\mathbf{u}-\lambda\mathbf{u}=\mathbf{f}.$ (4.11)

###### Proof.

The argument is standards and we provide only the general sketch for the sake
of completeness.

Testing (4.10) with $\mathbf{u}^{\varepsilon}$ one gets the following energy
bound:

$\|\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathcal{M}^{\varepsilon})}+\varepsilon\|\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathcal{I}^{\varepsilon})}+\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq
C.$

Applying the Extension Theorem 2.7 (modulo a rescaling argument) together with
Korn’s inequality, we obtain the extension
$\widetilde{\mathbf{u}}^{\varepsilon}$ to the whole of $\mathbb{R}^{d}$ of the
restriction $\mathbf{u}^{\varepsilon}|_{\mathcal{M}^{\varepsilon}}$ satisfying
the estimate

$\|\widetilde{\mathbf{u}}^{\varepsilon}\|_{\mathbf{H}^{1}(\mathbb{R}^{d})}\leq
C.$

Moreover, the difference
$\mathbf{u}^{\varepsilon}-\widetilde{\mathbf{u}}^{\varepsilon}\in\mathbf{H}^{1}_{0}(\mathcal{I}^{\varepsilon})$
satisfies (via the two preceding bound, Korn’s and the Poincaré inequalities)

$\varepsilon\|\nabla(\mathbf{u}^{\varepsilon}-\widetilde{\mathbf{u}}^{\varepsilon})\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}+\|\mathbf{u}^{\varepsilon}-\widetilde{\mathbf{u}}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq
C.$

Then, by the properties of stochastic two-scale convergence (see Appendix
A.1), the following convergences (up to extracting a subsequence) take place,

$\displaystyle\mathbf{u}^{\varepsilon}$ $\displaystyle\rightharpoonup$
$\displaystyle\mathbf{u}_{1}\mbox{ weakly in }\mathbf{L}^{2}(\mathbb{R}^{d}),$
$\displaystyle\mathbf{u}^{\varepsilon}-\widetilde{\mathbf{u}}^{\varepsilon}$
$\displaystyle\stackrel{{\scriptstyle 2}}{{\rightharpoonup}}$
$\displaystyle\overline{\mathbf{u}_{0}},$
$\displaystyle\nabla\widetilde{\mathbf{u}}^{\varepsilon}$
$\displaystyle\stackrel{{\scriptstyle 2}}{{\rightharpoonup}}$
$\displaystyle\nabla\mathbf{u}_{1}+\overline{\bm{p}},$
$\displaystyle\varepsilon\nabla(\mathbf{u}^{\varepsilon}-\widetilde{\mathbf{u}}^{\varepsilon})$
$\displaystyle\stackrel{{\scriptstyle 2}}{{\rightharpoonup}}$
$\displaystyle\overline{\nabla_{y}\mathbf{u}_{0}},$

where
$\overline{\bm{p}}\in\mathbf{L}^{2}(\mathbb{R}^{d};\bm{\mathcal{V}}^{2}_{\mathrm{pot}}(\Omega))$.

Testing the equation (4.10) with $\varepsilon
a(x)\bm{\varphi}(x/\varepsilon,\omega)$, $a\in
C^{\infty}_{0}(\mathbb{R}^{d})$, $\overline{\bm{\varphi}}\in\mathcal{C}$, and
passing to the limit as $\varepsilon\to 0$, we see that $\overline{\bm{p}}$ is
the corrector term, i.e. for every $x\in\mathbb{R}^{d}$ the random variable
$\overline{\bm{p}}(x,\cdot)$ is a minimiser of the problem (4.4) with
$\xi=\nabla\mathbf{u}_{1}(x)$. In particular,

$\mathbb{E}\left[\mathbbm{1}_{\mathbb{R}^{d}\setminus\omega}\mathsf{C}_{1}(\nabla\mathbf{u}_{1}+\bm{p})\right]=\mathsf{C}^{\mathrm{hom}}\nabla\mathbf{u}_{1}.$

Then passing to the limit in (4.10) with a test function
$\bm{\varphi}_{1}(x)+a(x)\bm{\varphi}(x/\varepsilon,\omega)$,
$\bm{\varphi}_{1}\in\mathbf{H}^{1}(\mathbb{R}^{d})$, $a\in
C^{\infty}_{0}(\mathbb{R}^{d})$, $\overline{\bm{\varphi}}\in\\{\text{countable
dense subset of }\mathbf{H}^{1}_{0}(\Omega)\\}$, and using the density
argument, one recovers the equation (4.11).

Finally, the strong stochastic two-scale convergence of
$\mathbf{u}^{\varepsilon}$ under the assumption
$\mathbf{f}^{\varepsilon}\overset{2}{\to}\mathbf{f}$ follows from the weak
stochastic two-scale convergence by a standard duality argument, see e.g.
[41]. ∎

The remainder of this section will be devoted to the examination of the
spectrum of $\mathcal{A}^{\mathrm{hom}}$. To this end, for
$\lambda\not\in\sigma(\mathcal{A}_{0})$, let
$\overline{\mathbf{b}^{(i)}_{\lambda}}\in\mathbf{H}^{1}_{0}(\Omega)$,
$i=1,\ldots,d$, be the unique solution of

$\mathcal{A}_{0}\overline{\mathbf{b}^{(i)}_{\lambda}}=\lambda\,\overline{\mathbf{b}^{(i)}_{\lambda}}+\mathbf{e}_{i},$
(4.12)

where $\mathbf{e}_{i}\in\mathbb{R}^{d}$, $(\mathbf{e}_{i})_{j}=\delta_{ij}$.

######  Lemma 4.6.

For every nonnegative $\lambda$, the stationary extension
$\mathbf{b}^{(i)}_{\lambda}$ of $\overline{\mathbf{b}^{(i)}_{\lambda}}$
satisfies

$-\operatorname{div}_{y}\mathsf{C}_{0}\nabla_{y}\,\mathbf{b}^{(i)}_{\lambda}=\lambda\,\mathbf{b}^{(i)}_{\lambda}+\mathbf{e}_{i}\qquad\text{in
}\ \omega\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{almost
surely}}.$ (4.13)

###### Proof.

The proof below relies on a classical mollification argument which can be
found, e.g., in [26, § 7.2, p. 232 ss.]. Let $K\in
C_{0}^{\infty}(\mathbb{R}^{d})$ be a non-negative even function such that
$\int_{\mathbb{R}^{d}}K(x)\,dx=1$, and, for $\delta>0$, set
$K^{\delta}(x):=\delta^{-d}K(\delta^{-1}x)$. For a random variable
$\overline{{f}}\in L^{2}(\Omega)$ (and, analogously, for vector/tensor-valued
random variables) we define its mollification to be

$\overline{{f}}^{\delta}(\omega):=\int_{\mathbb{R}^{d}}K^{\delta}(y)\,{f}(y,\omega)\,dy\,.$

It is easy to check that the mollification $\overline{{f}}^{\delta}$ possesses
the following properties: $\overline{{f}}^{\delta}\in C^{\infty}(\Omega)$,
${f}^{\delta}(y,\omega)=({f}\ast K^{\delta})(y,\omega)$,
$\partial_{y_{i}}{f}^{\delta}={f}\ast(\partial_{y_{i}}K^{\delta})$,
$\mathbb{E}[\overline{{f}}\,\overline{{g}}^{\delta}]=\mathbb{E}[\overline{{f}}^{\delta}\,\overline{{g}}]$,
and
$\mathbb{E}[\overline{{f}}\,\overline{\partial^{\bm{\alpha}}{g}}^{\delta}]=(-1)^{|{\bm{\alpha}}|}\mathbb{E}[\overline{\partial^{\bm{\alpha}}{f}}^{\delta}\,\overline{{g}}]$,
where $*$ stands for convolution in the variable $y$ and $\overline{{g}}\in
L^{2}(\Omega)$.

Let us define the measurable subset $\mathcal{O}\subset\Omega$ as

$\mathcal{O}:=\\{\omega\in\Omega\ |\ 0\in\omega\\}$

and let us introduce the random variable $\mathfrak{d}:\Omega\to\mathbb{R}$
defined as

$\mathfrak{d}(\omega):=\operatorname{dist}(0,\mathbb{R}^{d}\setminus\omega)\,,$

where $0$ is the origin in $\mathbb{R}^{d}$. Observe that the map
$\mathfrak{d}$ is measurable, cf., e.g., [10, Lemma B.3]. Furthermore, for
$\delta>0$ let $\mathcal{O}^{\delta}$ be the (measurable, in view of the
measurability of the translation group action) set

$\mathcal{O}^{\delta}:=\\{\omega\in\mathcal{O}\ |\
\mathfrak{d}(\omega)>\delta\\}\,.$

Clearly, $\mathcal{O}^{c\delta}\subset\mathcal{O}^{\delta}$ for $c>1$.

Let $\overline{\mathbf{v}}$ be an element in $\mathbf{L}^{2}(\Omega)$
vanishing outside $\mathcal{O}^{\delta}$. Then
$\overline{\mathbf{v}}^{\delta}\in\mathbf{C}^{\infty}(\Omega)$ and it vanishes
outside $\mathcal{O}$, therefore it can be used as a test function in the weak
formulation of (4.12) to obtain

$\mathbb{E}\left[\mathsf{C}_{0}\nabla_{y}\mathbf{b}_{\lambda}^{(i)}\cdot\nabla_{y}\mathbf{v}^{\delta}\right]=\mathbb{E}\left[\lambda\,\overline{\mathbf{b}_{\lambda}^{(i)}}\cdot\overline{\mathbf{v}}^{\delta}+(\overline{\mathbf{v}}^{\delta})_{i}\right].$

Using the properties of the mollification listed above we obtain

$\mathbb{E}\left[-\operatorname{div}_{y}\left(\mathsf{C}_{0}\nabla_{y}[\mathbf{b}_{\lambda}^{(i)}]^{\delta}\right)\cdot\mathbf{v}\right]=\mathbb{E}\left[\lambda\,\overline{\mathbf{b}_{\lambda}^{(i)}}^{\delta}\cdot\overline{\mathbf{v}}+\mathbf{e}_{i}^{\delta}\cdot\overline{\mathbf{v}}\right],$

which, since $\overline{\mathbf{v}}$ is arbitrary, in turn implies

$-\overline{\operatorname{div}_{y}\left(\mathsf{C}_{0}\nabla_{y}[\mathbf{b}_{\lambda}^{(i)}]^{\delta}\right)}=\lambda\,\overline{\mathbf{b}_{\lambda}^{(i)}}^{\delta}+\mathbf{e}_{i}^{\delta}\qquad\text{in}\quad\mathcal{O}^{\delta}\,.$

Upon taking the stationary extension of all quantities, the above equation is
valid without the bar in $\omega^{\delta}\subset\mathbb{R}^{d}$ almost surely,
where
$\omega^{\delta}:=\\{x\in\omega\,|\,\operatorname{dist}(x,\mathbb{R}^{d}\setminus\omega)>\delta\\}$.

Let us fix an element $\omega\in\Omega$ and choose
$\mathbf{g}\in\mathbf{C}_{0}^{\infty}(\mathbb{R}^{d})$ compactly supported in
one connected component of $\omega$ (i.e., in a single inclusion). Then for
sufficiently small $\delta$ we have

$\int_{\mathbb{R}^{n}}\mathsf{C}_{0}\nabla\mathbf{b}_{\lambda}^{(i)}\cdot[\nabla\mathbf{g}]^{\delta}=\lambda\int_{\mathbb{R}^{n}}\mathbf{b}_{\lambda}^{(i)}\cdot\mathbf{g}^{\delta}+\int_{\mathbf{R}^{n}}\mathbf{e}_{i}\cdot\mathbf{g}^{\delta}\,.$

Passing to the limit in the above equation as $\delta\to 0$ and using the fact
that $\mathbf{g}$ is smooth, we arrive at

$\int_{\mathbb{R}^{n}}\mathsf{C}_{0}\nabla\mathbf{b}_{\lambda}^{(i)}\cdot\nabla\mathbf{g}=\lambda\int_{\mathbb{R}^{n}}\mathbf{b}_{\lambda}^{(i)}\cdot\mathbf{g}+\int_{\mathbb{R}^{n}}\mathbf{e}_{i}\cdot\mathbf{g}.$

Since $\mathbf{g}$ is arbitrarily and chosen from a dense subspace of
$\mathbf{H}^{1}(\mathbb{R}^{d})$, the latter implies (4.13). ∎

Denote

$\overline{\mathbf{b}_{\lambda}}:=\left(\overline{\mathbf{b}^{(1)}_{\lambda}}\
|\ \overline{\mathbf{b}^{(2)}_{\lambda}}\ |\ \ldots\ |\
\overline{\mathbf{b}^{(d)}_{\lambda}}\right)\in[\mathbf{H}^{1}_{0}(\Omega)]^{d}.$
(4.14)

Clearly, $\overline{\mathbf{b}_{\lambda}}$ satisfies

$\mathcal{A}_{0}\overline{\mathbf{b}_{\lambda}}=\lambda\,\overline{\mathbf{b}_{\lambda}}+I\,.$
(4.15)

Note that, unlike in the periodic setting, the random operator
$\mathcal{A}_{0}$ does not have, in general, discrete spectrum. For example,
in the case we inclusions are allowed to “continuously” change size (e.g.,
being homothetic to a fixed shape), the spectrum of $\mathcal{A}_{0}$ is a
union of intervals. The following statement asserts that
$\sigma(\mathcal{A}_{0})$ can be recovered from a single $\omega$ almost
surely.

######  Theorem 4.7.

We have

$\sigma(\mathcal{A}_{0})=\overline{\bigcup_{k\in\mathbb{N}}\sigma(\mathcal{A}_{0}^{\omega^{k}})}$

for almost every $\omega\in\Omega$.

###### Proof.

The claim is proven analogously to [9, Theorem 5.6], with account of [10,
Remark 3.12]. ∎

Using Theorem 4.7, one can estimate $\mathbf{b}_{\lambda}$ on hypercubes.

######  Lemma 4.8.

There exists a positive constant $C=C_{\lambda,d}$ uniform in $\varepsilon$
such that

$\|\mathbf{b}_{\lambda}(\cdot/\varepsilon,\omega)\|_{[\mathbf{L}^{2}(\square_{z}^{L})]^{d}}+\|\varepsilon\,\nabla\mathbf{b}_{\lambda}(\cdot/\varepsilon,\omega)\|_{[\mathbf{L}^{2}(\square_{z}^{L})]^{d}}\leq
C\,L^{d/2}$ (4.16)

for all $z\in\mathbb{R}^{d}$ and $L>0$ almost surely.

###### Proof.

In each inclusion $\omega^{k}$ the extension
$\mathbf{b}_{\lambda}^{(i)}(\cdot,\omega)$, $i=1,\ldots,d$, solves the problem
(cf. (4.13))

$\mathcal{A}_{0}^{\omega^{k}}\mathbf{b}^{(i)}_{\lambda}-\lambda\,\mathbf{b}^{(i)}_{\lambda}=\mathbf{e}_{i}.$
(4.17)

Since $\lambda\not\in\sigma(\mathcal{A}_{0})$, in view of Theorem 4.7, the
resolvent $(\mathcal{A}_{0}^{\omega^{k}}-\lambda)^{-1}$ is bounded. Hence,
(4.17) implies

$\|\mathbf{b}^{(i)}_{\lambda}\|_{\mathbf{L}^{2}(\omega^{k})}\leq
C\,|\omega^{k}|^{1/2}\,.$ (4.18)

Moreover, taking the inner product of (4.17) with
$\mathbf{b}^{(i)}_{\lambda}$, integrating by parts over $\omega^{k}$, and
invoking Korn’s inequality, one obtains

$\begin{split}\|\nabla\mathbf{b}^{(i)}_{\lambda}\|_{\mathbf{L}^{2}(\omega^{k})}^{2}\leq\lambda\|\mathbf{b}^{(i)}_{\lambda}\|_{\mathbf{L}^{2}(\omega^{k})}^{2}+\|\mathbf{b}^{(i)}_{\lambda}\|_{\mathbf{L}^{2}(\omega^{k})}\,|\omega^{k}|^{1/2}\overset{\eqref{proof
lemma properties b lambda equation 2}}{\leq}C\,|\omega^{k}|\,.\end{split}$

Applying an elementary rescaling argument and observing that the $L^{2}$-norm
of $\mathbf{b}^{(i)}_{\lambda}(\cdot/\varepsilon,\omega)$ and
$\nabla\mathbf{b}^{(i)}_{\lambda}(\cdot/\varepsilon,\omega)$ over any given
box of size $L$ is bounded above by the said $L^{2}$-norms over all inclusions
$\varepsilon\omega^{k}$ fully contained in a box of size $L+2\varepsilon$ (see
Assumption 2.1) we finally arrive at (4.16). ∎

The matrix-function (4.14) allows us to define a stochastic matrix-valued
version of the Zhikov $\beta$-function [41].

######  Definition 4.9.

We define the _Zhikov $\beta$-matrix_ to be the matrix-function

$\bm{\beta}:[0,+\infty)\setminus\sigma(\mathcal{A}_{0})\to\mathrm{Mat}(d;\mathbb{R})\simeq\mathbb{R}^{d^{2}}$

given by

$\bm{\beta}(\lambda):=\mathbb{E}\left[\lambda\,I+\lambda^{2}\overline{\mathbf{b}_{\lambda}}\right],$
(4.19)

where $I$ is the $d\times d$ identity matrix.

######  Lemma 4.10.

The Zhikov $\beta$-matrix (4.19) is symmetric.

###### Proof.

The symmetry of $\mathbb{E}[\overline{\mathbf{b}_{\lambda}}]$ follows from the
self-adjointness of $\mathcal{A}_{0}$, namely,

$\mathbb{E}\left[\overline{\mathbf{b}_{\lambda}^{(i)}}\cdot\mathbf{e}_{j}\right]=\mathbb{E}\left[\overline{\mathbf{b}_{\lambda}^{(i)}}\cdot(\mathcal{A}_{0}-\lambda)\overline{\mathbf{b}_{\lambda}^{(j)}}\right]=\mathbb{E}\left[(\mathcal{A}_{0}-\lambda)\overline{\mathbf{b}_{\lambda}^{(i)}}\cdot\overline{\mathbf{b}_{\lambda}^{(j)}}\right]=\mathbb{E}\left[\mathbf{e}_{i}\cdot\overline{\mathbf{b}_{\lambda}^{(j)}}\right].$

∎

Very much like its scalar counterpart [10, Theorem 4.3], the
$\bm{\beta}$-matrix determines the spectrum of $\mathcal{A}^{\mathrm{hom}}$.

######  Theorem 4.11.

We have

$\sigma(\mathcal{A}^{\mathrm{hom}})=\sigma(\mathcal{A}_{0})\cup\overline{\\{\lambda\
|\ \bm{\beta}(\lambda)<0\\}^{\complement}}\,,$ (4.20)

where $\bm{\beta}(\lambda)<0$ means that the symmetric matrix
$\bm{\beta}(\lambda)$ is negative definite.

######  Remark 4.12.

Observe that if we denote by $\beta$ (non-bold beta) the scalar function

$\beta:[0,+\infty)\setminus\sigma(\mathcal{A}_{0})\to\mathbb{R},\qquad\beta(\lambda):=\max_{\widehat{k}\in\mathbb{S}^{d-1}}\widehat{k}\cdot\bm{\beta}(\lambda)\widehat{k},$
(4.21)

returning, for each $\lambda$ in its domain, the largest eigenvalue of the
matrix $\bm{\beta}(\lambda)$, then

$\overline{\\{\lambda\ |\ \bm{\beta}(\lambda)<0\\}^{\complement}}=\\{\lambda\
|\beta(\lambda)\geq 0\\}.$ (4.22)

The RHS of (4.22) is precisely the expression appearing in the scalar analogue
of Theorem 4.11 (although the two $\beta$’s are, in general, different
functions).

###### Proof of Theorem 4.11.

To begin with, let us show that

$\mathbb{R}\setminus\left(\sigma(\mathcal{A}_{0})\cup\overline{\\{\lambda\ |\
\bm{\beta}(\lambda)<0\\}^{\complement}}\right)\subset\rho(\mathcal{A}^{\mathrm{hom}})\,.$
(4.23)

The resolvent problem (4.11) can be recast as the following coupled system

$\int_{\mathbb{R}^{d}}\mathsf{C}^{\mathrm{hom}}\nabla\mathbf{u}_{1}\cdot\nabla\bm{\varphi}_{1}-\lambda\int_{\mathbb{R}^{d}}(\mathbf{u}_{1}+\mathbb{E}\left[\mathbf{u}_{0}\right])\cdot\bm{\varphi}_{1}\\\
=\int_{\mathbb{R}^{d}}\mathbb{E}\left[\mathbf{f}\right]\cdot\bm{\varphi}_{1},$
(4.24a)
$\mathbb{E}\left[\mathsf{C}_{0}\nabla_{y}\mathbf{u}_{0}(x,\cdot)\cdot\nabla_{y}\bm{\varphi}_{0}(x,\cdot)-\lambda(\mathbf{u}_{1}(x)+\mathbf{u}_{0}(x,\cdot))\cdot\bm{\varphi}_{0}(x,\cdot)\right]\\\
=\mathbb{E}\left[\mathbf{f}(x,\cdot)\cdot\bm{\varphi}_{0}(x,\cdot)\right],\forall
x\in\mathbb{R}^{d},$ (4.24b)

for every $\bm{\varphi}_{1}+\bm{\varphi}_{0}\in\mathbf{V}$. Since
$\lambda\not\in\sigma(\mathcal{A}_{0})$, equation (4.24b) is solvable with the
unique solution given by

$\mathbf{u}_{0}=(\mathcal{A}_{0}-\lambda\operatorname{Id})^{-1}\mathbf{f}+\lambda\,\overline{\mathbf{b}_{\lambda}}\,\mathbf{u}_{1}$
(4.25)

(recall that $\mathbf{b}_{\lambda}$ is the matrix-valued random variable
defined in (4.14)). Substituting (4.25) into (4.24a) and using (4.19) we
obtain

$\mathcal{A}_{1}\mathbf{u}_{1}-\bm{\beta}(\lambda)\mathbf{u}_{1}=\mathbb{E}[\mathbf{f}]+\lambda\mathbb{E}[(\mathcal{A}_{0}-\lambda\operatorname{Id})^{-1}\mathbf{f}].$
(4.26)

Since $\mathcal{A}_{1}$ is non-negative and $\bm{\beta}(\lambda)<0$, the
operator $\mathcal{A}_{1}-\bm{\beta}(\lambda)$ is positive and (4.26) is
uniquely solvable by the Fredholm alternative. This implies (4.23).

Next, let us assume that $\lambda\not\in\sigma(\mathcal{A}_{0})$ and
$\bm{\beta}(\lambda)$ has at least one non-negative eigenvalue. The task is to
show that $\lambda\in\sigma(\mathcal{A}^{\mathrm{hom}})$. To this end, it
suffices to construct a sequence
$\mathbf{u}_{n}\in\mathbf{C}_{0}^{\infty}(\mathbb{R}^{d})$, $n\in\mathbb{N}$,
$\|\mathbf{u}_{n}\|_{\mathbf{L}^{2}}=1$, such that

$\lim_{n\to\infty}\|(\mathcal{A}_{1}-\bm{\beta}(\lambda))\mathbf{u}_{n}\|_{\mathbf{L}^{2}}=0.$
(4.27)

Indeed, in this case the sequence $\mathbf{u}_{n}+\overline{\mathbf{v}_{n}}$,
with
$\overline{\mathbf{v}_{n}}:=\lambda\overline{\mathbf{b}_{\lambda}}\mathbf{u}_{n}\in\mathbf{L}^{2}(\mathbb{R}^{d};\mathbf{H}^{1}_{0}(\Omega))$,
is singular for the limit operator $\mathcal{A}^{\rm hom}$. In order to see
this we first observe that for two functions
$\mathbf{u}\in\mathbf{H}^{1}(\mathbb{R}^{d})$ and
$\overline{\mathbf{v}}\in\mathbf{L}^{2}(\mathbb{R}^{d};\mathbf{H}^{1}_{0}(\Omega))$
such that $\mathcal{A}_{1}\mathbf{u}=\mathbf{f}$ and
$\mathcal{A}_{0}\overline{\mathbf{v}}=\overline{\mathbf{g}}$ one has

$\mathcal{A}^{\rm
hom}\mathbf{u}=(1-\overline{1_{\Omega}})\mathbb{E}[(1-\overline{1_{\Omega}})]^{-1}\,\mathbf{f},$
$\mathcal{A}^{\rm
hom}\overline{\mathbf{v}}=\overline{\mathbf{g}}-(1-\overline{1_{\Omega}})\mathbb{E}[(1-\overline{1_{\Omega}})]^{-1}\mathbb{E}[\overline{\mathbf{g}}].$

Then elementary calculations give

$\lim_{n\to\infty}\|(\mathcal{A}^{\rm
hom}-\lambda)(\mathbf{u}_{n}+\overline{\mathbf{v}_{n}})\|_{\mathbf{L}^{2}(\mathbb{R}^{d}\times\Omega)}=0.$

We seek $\mathbf{u}_{n}$ in the form

$\mathbf{u}_{n}(x)=\frac{\eta_{n}(x)\,e^{ir\,k\cdot
x}\,c}{\|\eta_{n}(x)c\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}},$ (4.28)

where $k$ is a given vector in $\mathbb{S}^{d-1}$, the cut-off function
$\eta_{n}$ is defined in accordance with (2.3), and $c\in\mathbb{R}^{d}$,
$r\in\mathbb{R}$ are constants chosen as follows. First, we choose $r$ in such
a way that444Note that $k\cdot\mathsf{C}\cdot k$ is a $d\times d$ matrix with
matrix elements  $(k\cdot\mathsf{C}^{\mathrm{hom}}\cdot
k)_{\alpha\beta}=\mathsf{C}^{\mathrm{hom}}_{\alpha\mu\nu\beta}k_{\mu}k_{\nu}.$
Due to the properties of the tensor $\mathsf{C}^{\mathrm{hom}}$, the matrix
$k\cdot\mathsf{C}^{\mathrm{hom}}\cdot k$ is clearly symmetric and positive
definite.

$\det(r^{2}\,k\cdot\mathsf{C}^{\mathrm{hom}}\cdot k-\bm{\beta}(\lambda))=0.$

This is clearly possible, because $\beta(\lambda)$ has at least one non-
negative eigenvalue. The vector $c$ is then chosen in the kernel of
$r^{2}\,k\cdot\mathsf{C}^{\mathrm{hom}}\cdot k-\bm{\beta}(\lambda)$, i.e.,

$\left[r^{2}\,k\cdot\mathsf{C}^{\mathrm{hom}}\cdot
k-\bm{\beta}(\lambda)\right]c=0.$ (4.29)

It should be noted that $r$ depends on the choice of $k$, whereas $c$ depends
on the choice of both $k$ and $r$.

Then formulae (2.4), (4.28) and (4.29) imply

$\|(\mathcal{A}_{1}-\bm{\beta}(\lambda))\mathbf{u}_{n}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq\frac{C}{n}.$

which, in turn, gives us (4.27).

It only remains to show that
$\sigma(\mathcal{A}_{0})\subset\sigma(\mathcal{A}^{\mathrm{hom}})$. Let
$\lambda\in\rho(\mathcal{A}^{\mathrm{hom}})$. In order to conclude, by the
inverse mapping theorem it is enough to show that for every
$\overline{\mathbf{g}}\in\mathbf{L}^{2}_{0}(\Omega)$ the equation

$(\mathcal{A}_{0}-\lambda)\overline{\mathbf{v}}=\overline{\mathbf{g}}$ (4.30)

admits a solution.

Since $\lambda\not\in\sigma(\mathcal{A}^{\mathrm{hom}})$, let
$\mathbf{u}=\mathbf{u}_{1}+\overline{\mathbf{u}_{0}}$ the solution to (4.24a),
(4.24b) with $\overline{\mathbf{f}}=h\overline{\mathbf{g}}$, for some $h\in
L^{2}(\mathbb{R}^{d})$. Then, by (4.24b), we have

$(\mathcal{A}_{0}-\lambda)\overline{\mathbf{u}_{0}}=\lambda\mathbf{u}_{1}+h\,\overline{\mathbf{g}},$

and hence555Recall that $\overline{\mathbf{b}_{\lambda}}$ is defined by
(4.14).

$(\mathcal{A}_{0}-\lambda)(\overline{\mathbf{u}_{0}}-\lambda\overline{\mathbf{b}_{\lambda}}\mathbf{u}_{1})=h\,\overline{\mathbf{g}}.$
(4.31)

Equation (4.31) now implies that

$\overline{\mathbf{v}}:=\int_{\mathbb{R}^{d}}\frac{h}{\|h\|_{L^{2}}}(\overline{\mathbf{u}_{0}}-\lambda\overline{\mathbf{b}_{\lambda}}\mathbf{u}_{1})$

is a solution of (4.30), as one can establish by a straightforward
calculation. This completes the proof. ∎

In plain English, the above theorem tells us that the spectrum of
$\mathcal{A}^{\mathrm{hom}}$ is comprised of two parts: (i) those $\lambda$’s
for which $\bm{\beta}(\lambda)$ has at least one non-negative eigenvalue and
(ii) the ‘microscopic resonances’ $\sigma(\mathcal{A}_{0})$.

The following statement is an immediate consequence of Theorem 4.5.

######  Corollary 4.13.

We have

$\sigma(\mathcal{A}^{\mathrm{hom}})\subset\lim_{\varepsilon\to
0}\sigma(\mathcal{A}^{\varepsilon})$ (4.32)

almost surely.

##  5 Spectral convergence

The goal of this section is to characterise the RHS of (4.32).

Recall that the spectrum of the homogenised operator is characterised by
Zhikov’s $\bm{\beta}$-matrix, cf. Definition 4.9 and Theorem 4.11. As it turns
out, the limiting spectrum $\lim_{\varepsilon\to
0}\sigma(\mathcal{A}^{\varepsilon})$ is determined by a similar quantity,
albeit scalar and nonlocal — the $\beta_{\infty}$-function — which “feels” the
atypical local arrangements of inclusions.

######  Definition 5.1.

We define the _$\beta_{\infty}$ -function_ to be the function

$\beta_{\infty}:\mathbb{R}^{+}\backslash\sigma(\mathcal{A}_{0})\times\Omega\to\mathbb{R},$
$\beta_{\infty}(\lambda,\omega):=\liminf_{L\to+\infty}\sup_{x\in\mathbb{R}^{d}}\sup_{k\in\mathbb{S}^{d-1}}\frac{1}{|\square_{x}^{L}|}\int_{\square_{x}^{L}}\left(\lambda+\lambda^{2}\,k\cdot\mathbf{b}_{\lambda}(y,\omega)k\right)dy\,.$
(5.1)

######  Remark 5.2.

Since for every fixed $\lambda$ the random variable
$\overline{\mathbf{b}_{\lambda}}$ is defined up to a set of probability
measure zero (which in general depends on $\lambda$), in order define
$\beta_{\infty}$ as a function of $\lambda$ almost surely we need to work with
an appropriate version (in the probabilistic sense of the word) of the process
$(\overline{\mathbf{b}_{\lambda}})_{\lambda\notin\sigma(\mathcal{A}_{0})}$.
This version can be obtained by using the Kolmogorov–Čentsov Theorem (see e.g.
[29, Theorem 2.8]) and the fact that

$\mathbb{E}\left[|\overline{\mathbf{b}_{\lambda_{1}}}-\overline{\mathbf{b}_{\lambda_{2}}}|^{2}\right]\leq
C(\lambda_{1}-\lambda_{2})^{2},$

whenever one stays away from the spectrum of $\mathcal{A}_{0}$. Note that in
[10] the authors adopted a different approach, defining an appropriate version
of the scalar analogue of $\overline{\mathbf{b}_{\lambda}}$ through
realisations, cf. (4.13).

The following proposition lists some basic properties of $\beta_{\infty}$.
Their proof is a straightforward adaptation of the arguments presented in [10,
Lemma 5.9 and Proposition 5.11] and is not provided here for the sake of
brevity.

######  Proposition 5.3.

  1. (i)

The limit inferior in the right hand side of (5.1) can be replaced by the
limit.

  2. (ii)

The function

$(\lambda,\omega)\mapsto\beta_{\infty}(\lambda,\omega)$

is deterministic almost surely, continuous on
$\mathbb{R}^{+}\backslash\sigma(\mathcal{A}_{0})$, and strictly increasing on
every connected component (interval) thereof.

In view of Proposition 5.3, in what follows we suppress the dependence of
$\beta_{\infty}$ on $\omega$ in our notation and write simply
$\beta_{\infty}(\lambda)$.

######  Remark 5.4.

Note that for $\lambda\geq 0$ the largest eigenvalue $\beta(\lambda)$ of the
matrix $\bm{\beta}(\lambda)$ defined in accordance with (4.21) satisfies the
following bound:

$\begin{split}\beta(\lambda)&=\max_{\widehat{k}\in\mathbb{S}^{d-1}}\widehat{k}\cdot\bm{\beta}(\lambda)\widehat{k}\\\
&=\max_{\widehat{k}\in\mathbb{S}^{d-1}}\lim_{L\to+\infty}\frac{1}{|\square_{x}^{L}|}\int_{\square_{x}^{L}}\left(\lambda+\lambda^{2}\,\widehat{k}\cdot\mathbf{b}_{\lambda}(T_{y}\omega)\widehat{k}\right)\\\
&=\lim_{L\to+\infty}\max_{\widehat{k}\in\mathbb{S}^{d-1}}\frac{1}{|\square_{x}^{L}|}\int_{\square_{x}^{L}}\left(\lambda+\lambda^{2}\,\widehat{k}\cdot\mathbf{b}_{\lambda}(T_{y}\omega)\widehat{k}\right)\\\
&\overset{\eqref{eq definition beta
infinity}}{\leq}\beta_{\infty}(\lambda).\end{split}$ (5.2)

Here $x$ is an arbitrary point in $\mathbb{R}^{d}$, and in writing the second
equality we used the Ergodic Theorem. Formula (5.2) provides a generalisation
to systems of its scalar counterpart [10, Remark 5.7], with the Zhikov
$\beta$-function being replaced by the largest eigenvalue of the
$\bm{\beta}$-matrix (denoted here by $\beta(\lambda)$ by analogy) — cf. Remark
4.12.

A key role in the upcoming analysis will be played by the set

$\mathcal{G}:=\sigma(\mathcal{A}_{0})\cup\overline{\\{\lambda:\beta_{\infty}(\lambda)\geq
0\\}}\,,$ (5.3)

compare with (4.20) and (4.22).

###  5.1 Upper bound for the limiting spectrum

v Our first result provides an upper bound for the limiting spectrum.

######  Theorem 5.5.

We have

$\lim_{\varepsilon\to 0}\sigma(\mathcal{A}^{\varepsilon})\subseteq\mathcal{G}$

almost surely.

###### Proof.

Let $\lambda_{0}\in\lim_{\varepsilon\to 0}\sigma(\mathcal{A}^{\varepsilon})$.
This means that there exist
$\lambda_{\varepsilon}\in\sigma(\mathcal{A}^{\varepsilon})$ and
$\mathbf{u}^{\varepsilon}\in\mathcal{D}(\mathcal{A}^{\varepsilon})$,
$\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}=1$, such that

$\lim_{\varepsilon\to 0}\lambda_{\varepsilon}=\lambda_{0},$
$\|\mathcal{A}^{\varepsilon}\mathbf{u}^{\varepsilon}-\lambda_{\varepsilon}\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}=:\delta_{\varepsilon}=o(1)\quad\text{as}\quad\varepsilon\to
0.$ (5.4)

If $\lambda_{0}\in\sigma(\mathcal{A}_{0})$ there is nothing to prove. Hence,
without loss of generality, let us assume that
$\lambda_{0}\not\in\sigma(\mathcal{A}_{0})$. The task at hand is then to show
that

$\beta_{\infty}(\lambda_{0})\geq 0\,.$ (5.5)

For the future reference we note that (5.4) immediately implies the following
bound on the symmetrised gradient of $\mathbf{u}^{\varepsilon}$:

$\|\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|_{[\mathbf{L}^{2}(\mathcal{M}^{\varepsilon})]^{d}}+\varepsilon\|\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|_{[\mathbf{L}^{2}(\mathcal{I}^{\varepsilon})]^{d}}\leq
C.$ (5.6)

The idea of the proof goes as follows. Intuitively, one can think of the soft
inclusions as micro-resonators, and the local average of the term
$\lambda^{2}\,\mathbf{b}_{\lambda}(y,\omega)$ in (5.1) represents a local
resonant (or anti-resonant, depending on the sign) contribution from said
inclusions. Thus, $\beta_{\infty}(\lambda)$ quantifies the maximal possible
resonant contribution of the soft component on large scales. One expects that
this resonant behaviour is present in the structure of the approximate
eigenfunctions $\mathbf{u}^{\varepsilon}$ introduced above. Therefore, we
devise a procedure which will bring this feature to the surface in the limit
$\varepsilon\to 0$, and relate it to $\beta_{\infty}(\lambda)$. We argue that
one can restrict the attention to large fixed size cubes which carry
significant proportion of the energy of $\mathbf{u}^{\varepsilon}$, at the
same time keeping the error in the quasimode approximation relatively small.
Once one has chosen the cubes appropriately, it only remains to shift such
cubes to the origin in order to utilise a compactness argument. We proceed
with this plan in mind.

Let $L>1$. We claim that there exists $y=y_{L,\varepsilon}\in\mathbb{R}^{d}$
such that

$\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\square_{y}^{L})}>0$ (5.7)

and

$\frac{1}{3^{d+1}\delta_{\varepsilon}^{2}}\|\mathcal{A}^{\varepsilon}\mathbf{u}^{\varepsilon}-\lambda_{\varepsilon}\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\square_{y}^{3L})}+\frac{1}{3^{d+1}}\frac{\|\mathbbm{1}_{\mathcal{M}^{\varepsilon}}\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|^{2}_{[\mathbf{L}^{2}(\square_{y}^{3L})]^{d}}}{\|\mathbbm{1}_{\mathcal{M}^{\varepsilon}}\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|^{2}_{[\mathbf{L}^{2}(\mathbb{R}^{d})]^{d}}}+\frac{1}{3^{d+1}}\|\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\square_{y}^{3L})}\leq\|\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\square_{y}^{L})}\,.$
(5.8)

Indeed,

$\frac{1}{3^{d+1}}\sum_{z\in
L\mathbb{Z}}\left[\frac{\|\mathcal{A}^{\varepsilon}\mathbf{u}^{\varepsilon}-\lambda_{\varepsilon}\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\square_{z}^{3L})}}{\|\mathcal{A}^{\varepsilon}\mathbf{u}^{\varepsilon}-\lambda_{\varepsilon}\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\mathbb{R}^{d})}}+\frac{\|\mathbbm{1}_{\mathcal{M}^{\varepsilon}}\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|^{2}_{[\mathbf{L}^{2}(\square_{z}^{3L})]^{d}}}{\|\mathbbm{1}_{\mathcal{M}^{\varepsilon}}\operatorname{sym}\nabla\mathbf{u}^{\varepsilon}\|^{2}_{[\mathbf{L}^{2}(\mathbb{R}^{d})]^{d}}}+\|\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\square_{z}^{3L})}\right]\\\
=1=\sum_{z\in
L\mathbb{Z}}\|\mathbf{u}^{\varepsilon}\|^{2}_{\mathbf{L}^{2}(\square_{z}^{L})}.$
(5.9)

Then, if we restrict the summation in (5.9) to the lattice points $z$ such
that $\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\square_{z}^{L})}\neq 0$,
(5.8) is an immediate consequence of the equality between real series with
strictly positive terms666Of course, there exists at least one $z\in
L\mathbb{Z}$ such that
$\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\square_{z}^{L})}\neq 0$ because
$\mathbf{u}^{\varepsilon}$ is normalised in $\mathbf{L}^{2}$..

We define a new sequence of vector-functions $\mathbf{w}^{\varepsilon}_{L}$,
normalised in $\mathbf{L}^{2}(\square_{0}^{L})$, as

$\mathbf{w}_{L}^{\varepsilon}(x):=\frac{\mathbf{u}^{\varepsilon}(x+y)}{\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\square_{y}^{L})}}\,.$
(5.10)

Recall that the shift $y$ is $\varepsilon$-dependent. At this point it is
convenient to define the set of shifted inclusions

$\varepsilon\widehat{\omega}=\bigsqcup_{k\in\mathbb{N}}\varepsilon\widehat{\omega}^{k}:=\varepsilon\omega-y,$

and the “shifted” operator $\widehat{\mathcal{A}}^{\varepsilon}$ as the
analogue of $\mathcal{A}^{\varepsilon}$ with shifted coefficients
$\mathsf{C}^{\varepsilon}(x+y,\omega)$. Similarly, we shall denote by
$\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}$ and
$\widehat{\mathbbm{1}}_{\mathcal{I}^{\varepsilon}}$ the shifted characteristic
functions of $\mathcal{M}^{\varepsilon}$ and $\mathcal{I}^{\varepsilon}$.
Then, the definition (5.10) combined with the bounds (5.8) and (5.6) gives us

$\|\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}\operatorname{sym}\nabla\mathbf{w}_{L}^{\varepsilon}\|^{2}_{[\mathbf{L}^{2}(\square_{0}^{3L})]^{d}}+\|\mathbf{w}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\square_{0}^{3L})}\leq
C$ (5.11)

and

$\|\widehat{\mathcal{A}}^{\varepsilon}\mathbf{w}_{L}^{\varepsilon}-\lambda_{\varepsilon}\mathbf{w}_{L}^{\varepsilon}\|_{\mathbf{L}^{2}(\square_{0}^{3L})}\leq
C\,\delta_{\varepsilon}\,.$ (5.12)

Let $\widetilde{\mathbf{u}}^{\varepsilon}$ be the extension of
$\mathbf{u}^{\varepsilon}|_{\mathcal{M}^{\varepsilon}}$ into
$\mathcal{I}^{\varepsilon}$ as per Theorem 2.7. Clearly, the corresponding
extension $\widetilde{\mathbf{w}}_{L}^{\varepsilon}$ of
$\mathbf{w}_{L}^{\varepsilon}$ is given by
$\frac{\widetilde{\mathbf{u}}^{\varepsilon}(x+y)}{\|\mathbf{u}^{\varepsilon}\|_{\mathbf{L}^{2}(\square_{y}^{L})}}$.
We denote
$\mathbf{v}^{\varepsilon}_{L}:=\mathbf{w}^{\varepsilon}_{L}-\widetilde{\mathbf{w}}^{\varepsilon}_{L}$
— this function is supported only on the (shifted) inclusions and carries the
information about the resonant effect of the soft component.

First, we will derive some basic bounds for
$\widetilde{\mathbf{w}}^{\varepsilon}_{L}$ and $\mathbf{v}^{\varepsilon}_{L}$.
Recalling that the extension $\widetilde{\mathbf{w}}^{\varepsilon}$ is chosen
to be $\mathsf{C}_{0}$-harmonic on the set of inclusions, we see that
$\mathbf{v}^{\varepsilon}_{L}$ satisfies the identity

$(-\varepsilon^{2}\operatorname{div}\mathsf{C}_{0}\nabla-\lambda_{\varepsilon})\mathbf{v}^{\varepsilon}_{L}=\lambda_{\varepsilon}\widetilde{\mathbf{w}}^{\varepsilon}_{L}+\left(-\varepsilon^{2}\operatorname{div}\mathsf{C}_{0}\nabla-\lambda_{\varepsilon}\right)\mathbf{w}^{\varepsilon}_{L},\qquad\mbox{in
}\varepsilon\widehat{\omega},$ (5.13)

whose rescaled version reads

$(-\operatorname{div}\mathsf{C}_{0}\nabla-\lambda_{\varepsilon})\mathbf{v}^{\varepsilon}_{L}(\varepsilon
x)=\lambda_{\varepsilon}\widetilde{\mathbf{w}}^{\varepsilon}_{L}(\varepsilon
x)+\left(-\operatorname{div}\mathsf{C}_{0}\nabla-\lambda_{\varepsilon}\right)\mathbf{w}^{\varepsilon}_{L}(\varepsilon
x),\qquad x\in\widehat{\omega}.$ (5.14)

The latter immediately implies the bound

$\|\mathbf{v}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\varepsilon\omega^{k})}\leq\frac{\|\lambda_{\varepsilon}\widetilde{\mathbf{w}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\varepsilon\omega^{k})}+\delta_{\varepsilon}}{\operatorname{dist}(\lambda_{\varepsilon},\sigma(\mathcal{A}_{0}))},\qquad\forall
k=1,2,\ldots.$ (5.15)

(Recall the definition of $\delta_{\varepsilon}$ (5.4).) With account of
(5.11), the extension Theorem 2.7 gives us

$\|\widetilde{\mathbf{w}}^{\varepsilon}_{L}\|_{\mathbf{H}^{1}(\square_{0}^{5L/2})}\leq
C.$ (5.16)

Furthermore, (5.15) and (5.4) imply that
$\widetilde{\mathbf{w}}^{\varepsilon}_{L}$ does not vanish in the limit:

$0<C\leq\|\widetilde{\mathbf{w}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\square_{0}^{2L})}\,.$
(5.17)

Finally, testing (5.13) with $\mathbf{v}^{\varepsilon}_{L}$ on each inclusion
and taking into account (5.15), (5.16) and (5.12) we arrive at

$\varepsilon\|\operatorname{sym}\nabla\mathbf{w}^{\varepsilon}_{L}\|_{[\mathbf{L}^{2}(\square_{0}^{2L})]^{d}}\leq
C.$ (5.18)

Let $\mathbf{w}_{L}^{0}$ and $\mathbf{g}_{\lambda}$ be defined as the weak
limits

$\widetilde{\mathbf{w}}_{L}^{\varepsilon}\rightharpoonup\mathbf{w}_{L}^{0}\neq
0\quad\text{in}\quad\mathbf{H}^{1}(\square_{0}^{5L/2})$ (5.19)

and

$\mathbf{b}_{\lambda_{0}}(\varepsilon^{-1}(\cdot+y),\omega)\rightharpoonup\mathbf{g}_{\lambda_{0}}\quad\text{in}\quad[\mathbf{L}^{2}(\square_{0}^{2L})]^{d},$
(5.20)

respectively.

We claim that

  1. (a)

the limits (5.19) and (5.20) exist (indeed, the first follows from (5.16) and
(5.17), and the second is a consequence of Lemma 4.8);

  2. (b)

we have (recall definition (2.3))

$\lim_{\varepsilon\to
0}\int_{\square_{0}^{2L}}\lambda_{\varepsilon}\mathbf{w}_{L}^{\varepsilon}\cdot\widetilde{\mathbf{w}}_{L}^{\varepsilon}\eta_{2L}\geq-
CL^{-1}$ (5.21)

and

  3. (c)

we have777Note that here the order of the terms matters.

$\mathbf{v}^{\varepsilon}_{L}=\mathbf{w}^{\varepsilon}_{L}-\widetilde{\mathbf{w}}^{\varepsilon}_{L}\rightharpoonup\lambda_{0}\,\mathbf{g}_{\lambda_{0}}\,\mathbf{w}_{L}^{0}\quad\text{in}\quad\mathbf{L}^{2}(\square_{0}^{2L}).$
(5.22)

Suppose this is the case. Then we have

$\begin{split}-\frac{C}{L}&\overset{(b)}{\leq}\lim_{\varepsilon\to
0}\int_{\square_{0}^{2L}}\lambda_{\varepsilon}\mathbf{w}_{L}^{\varepsilon}\cdot\widetilde{\mathbf{w}}_{L}^{\varepsilon}\eta_{2L}\\\
&=\lim_{\varepsilon\to
0}\int_{\square_{0}^{2L}}\lambda_{\varepsilon}\left[\widetilde{\mathbf{w}}_{L}^{\varepsilon}+\mathbf{v}_{L}^{\varepsilon}\right]\cdot\widetilde{\mathbf{w}}_{L}^{\varepsilon}\eta_{2L}\\\
&\overset{(a)+(c)}{=}\int_{\square_{0}^{2L}}\lambda_{0}\left(1+\lambda_{0}\,\mathbf{g}_{\lambda_{0}}\right)\mathbf{w}_{L}^{0}\cdot\mathbf{w}_{L}^{0}\eta_{2L}\\\
&\leq C\,\beta_{\infty}(\lambda_{0}).\end{split}$ (5.23)

Here we have used that (5.19) implies strong convergence in $\mathbf{L}^{2}$.
Furthermore, in the last step we applied the estimate

$\lambda_{0}+\lambda_{0}^{2}\,k\cdot\mathbf{g}_{\lambda_{0}}(x)\,k\leq\beta_{\infty}(\lambda_{0})\quad\text{for
a.e.}\quad x\in\square_{0}^{2L},\ \forall k\in\mathbb{S}^{d-1},$

which easily follows888Indeed, taking into account Proposition 5.3 the
assertion follows by integrating
$\lambda_{0}+\lambda_{0}^{2}\,k\cdot\mathbf{b}_{\lambda_{0}}(\varepsilon^{-1}(\cdot+y),\omega)\,k$
over an arbitrary fixed cube contained in $\square^{2L}_{0}$ and passing to
the limit. from (5.1), to

$\int_{\square_{0}^{2L}}\lambda_{0}\left(1+\lambda_{0}\,\mathbf{g}_{\lambda_{0}}\right)\mathbf{w}_{L}^{0}\cdot\mathbf{w}_{L}^{0}\,\eta_{2L}\\\
=\int_{\square_{0}^{2L}}\mathbbm{1}_{\\{\mathbf{w}^{0}_{L}\neq
0\\}}\left(\lambda_{0}+\lambda_{0}^{2}\frac{\mathbf{w}^{0}_{L}}{|\mathbf{w}^{0}_{L}|}\cdot\mathbf{g}_{\lambda_{0}}\frac{\mathbf{w}^{0}_{L}}{|\mathbf{w}^{0}_{L}|}\right)|\mathbf{w}^{0}_{L}|^{2}\,\eta_{2L}\\\
\leq\beta_{\infty}(\lambda_{0})\int_{\square_{0}^{2L}}|\mathbf{w}^{0}_{L}|^{2}\,\eta_{2L}.$

Now, the integral on the RHS of the latter is bounded uniformly in $L$ due to
(5.16) and (5.19). Since (5.23) holds for an arbitrary $L$ we conclude the
validity of (5.5).

The remainder of the proof will be concerned with demonstrating (b) and (c).

Let us prove (b). Multiplying

$(\widehat{\mathcal{A}}^{\varepsilon}-\lambda_{\varepsilon})\mathbf{w}_{L}^{\varepsilon}=:\mathbf{f}_{L}^{\varepsilon}$
(5.24)

by $\eta_{2L}\widetilde{\mathbf{w}}_{L}^{\varepsilon}$ and integrating over
$\square_{0}^{2L}$, one gets

$\int_{\square_{0}^{2L}}\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}\,\mathsf{C}_{1}\nabla\mathbf{w}_{L}^{\varepsilon}\cdot\nabla(\eta_{2L}\mathbf{w}_{L}^{\varepsilon})+\varepsilon^{2}\int_{\square_{0}^{2L}}\widehat{\mathbbm{1}}_{\mathcal{I}^{\varepsilon}}\,\mathsf{C}_{0}\nabla\mathbf{w}_{L}^{\varepsilon}\cdot\nabla(\eta_{2L}\widetilde{\mathbf{w}}_{L}^{\varepsilon})\\\
=\lambda_{\varepsilon}\int_{\square_{0}^{2L}}\eta_{2L}\,\mathbf{w}_{L}^{\varepsilon}\cdot\widetilde{\mathbf{w}}_{L}^{\varepsilon}+\int_{\square_{0}^{2L}}\eta_{2L}\,\mathbf{f}_{L}^{\varepsilon}\cdot\widetilde{\mathbf{w}}_{L}^{\varepsilon}\,.$
(5.25)

We estimate the second integral on LHS via (5.16) and (5.18) to obtain

$\lim_{\varepsilon\to
0}\left|\varepsilon^{2}\int_{\square_{0}^{2L}}\widehat{\mathbbm{1}}_{\mathcal{I}^{\varepsilon}}\,\mathsf{C}_{0}\nabla\mathbf{w}_{L}^{\varepsilon}\cdot\nabla(\eta_{2L}\widetilde{\mathbf{w}}_{L}^{\varepsilon})\right|=0.$
(5.26)

Furthermore, (5.16) and the properties of $\eta_{2L}$ imply

$\liminf_{\varepsilon\to
0}\int_{\square_{0}^{2L}}\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}\,\mathsf{C}_{1}\nabla\mathbf{w}_{L}^{\varepsilon}\cdot\nabla(\eta_{2L}\mathbf{w}_{L}^{\varepsilon})\geq-\frac{C}{L}\,$

for some $C$ independent of $L$. Finally, (5.12) implies

$\lim_{\varepsilon\to
0}\left|\int_{\square_{0}^{2L}}\eta_{2L}\mathbf{f}_{L}^{\varepsilon}\cdot\widetilde{\mathbf{w}}_{L}^{\varepsilon}\right|=0.$
(5.27)

Combining (5.25) and (5.26)–(5.27) we arrive at (5.21).

Next, let us show (c). Following [10], for the given $L$ we introduce a one-
parameter family of averaging operators
$P^{\varepsilon}:\mathbf{L}^{2}(\mathbb{R}^{d})\to\mathbf{L}^{2}(\mathbb{R}^{d})$
defined by

$P^{\varepsilon}\mathbf{g}(x):=\begin{cases}\frac{1}{|\varepsilon\widehat{\omega}^{k}|}\int_{\varepsilon\widehat{\omega}^{k}}\mathbf{g}(z)\,dz&\text{for
}x\in\varepsilon\widehat{\omega}^{k},\,k\in\mathbb{N},\\\
\mathbf{g}(x)&\text{otherwise}.\end{cases}$

What the operator $P^{\varepsilon}$ does is that on every shifted inclusion
$\varepsilon\widehat{\omega}^{k}$ it replaces the function $\mathbf{g}$ with
its average over said inclusion. Note that the set
$\varepsilon\widehat{\omega}$ changes with $\varepsilon$ in a nontrivial way,
because for different $\varepsilon$’s we get different shifts
$y=y_{L,\varepsilon}$; moreover, $\varepsilon\widehat{\omega}$ comprises
disjoint open sets whose diameters tend to zero as $\varepsilon$ goes to zero.
A slavish adaptation of [10, Lemma E.1] to the case of vector-valued functions
gives us

$\lim_{\varepsilon\to
0}\|P^{\varepsilon}\mathbf{g}-\mathbf{g}\|_{\mathbf{L}^{2}(\square_{0}^{3L})}=0$
(5.28)

for every $\mathbf{g}\in\mathbf{L}^{2}(\mathbb{R}^{d})$.

With account of (5.24), equation (5.13) reads

$-\varepsilon^{2}\operatorname{div}\mathsf{C}_{0}\nabla\mathbf{v}^{\varepsilon}_{L}-\lambda_{\varepsilon}\,\mathbf{v}^{\varepsilon}_{L}=\lambda_{\varepsilon}\widetilde{\mathbf{w}}_{L}^{\varepsilon}+\mathbf{f}_{L}^{\varepsilon}\,.$

Let us decompose $\mathbf{v}^{\varepsilon}_{L}$ as
$\mathbf{v}^{\varepsilon}_{L}=\widehat{\mathbf{v}}^{\varepsilon}_{L}+\mathring{\mathbf{v}}^{\varepsilon}_{L}$,
where $\widehat{\mathbf{v}}^{\varepsilon}_{L}$ and
$\mathring{\mathbf{v}}^{\varepsilon}_{L}$ are defined as solutions of

$-\varepsilon^{2}\operatorname{div}\mathsf{C}_{0}\nabla\widehat{\mathbf{v}}^{\varepsilon}_{L}-\lambda_{0}\,\widehat{\mathbf{v}}^{\varepsilon}_{L}=\lambda_{0}\,P^{\varepsilon}(\widetilde{\mathbf{w}}_{L}^{\varepsilon})\,$
(5.29)

and

$-\varepsilon^{2}\operatorname{div}C_{0}\nabla\mathring{\mathbf{v}}^{\varepsilon}_{L}-\lambda_{\varepsilon}\mathring{\mathbf{v}}^{\varepsilon}_{L}=(\lambda_{\varepsilon}-\lambda_{0})\widehat{\mathbf{v}}^{\varepsilon}_{L}+\lambda_{\varepsilon}\widetilde{\mathbf{w}}^{\varepsilon}_{L}-\lambda_{0}P^{\varepsilon}(\mathbf{w}_{L}^{0})+\mathbf{f}^{\varepsilon}_{L}\,.$
(5.30)

Arguing as in (5.14)–(5.15), one obtains from (5.29) the estimate

$\|\widehat{\mathbf{v}}^{\varepsilon}_{L}\|_{L^{2}(\square_{0}^{2L})}\leq C;$
(5.31)

similarly, with the help of (5.31), (5.28), (5.24), (5.12), the identity
(5.30) implies

$\lim_{\varepsilon\to
0}\|\mathring{\mathbf{v}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\square_{0}^{2L})}=0.$
(5.32)

Finally, we observe that

$\widehat{\mathbf{v}}^{\varepsilon}_{L}=\lambda_{0}\,\mathbf{b}_{\lambda_{0}}(x/\varepsilon+y,\omega)P^{\varepsilon}(\mathbf{w}_{L}^{0})$
(5.33)

solves (5.29) (recall that the shift $y$ depends on $L$ and $\varepsilon$, see
(5.7) and the line above).

Formulae (5.31), (5.32), (5.33) and (5.28) imply that the weak limit of
$\mathbf{v}^{\varepsilon}_{L}=\widehat{\mathbf{v}}^{\varepsilon}_{L}+\mathring{\mathbf{v}}^{\varepsilon}_{L}$
is (5.22), as claimed. This concludes the proof. ∎

###  5.2 Lower bound on the limiting spectrum

In this subsection we will prove a lower bound for the limiting spectrum under
the additional assumption of finite range of dependence. Combined with Theorem
5.5, this will establish equality between $\mathcal{G}$ and the limiting
spectrum.

In order to formulate the assumption on the finite range of dependence we
first need to introduce some notation. Let a set $K\subset\mathbb{R}^{d}$ be
compact. We denote by $\widetilde{\mathcal{F}}_{K}$ the $\sigma$-algebra (in
$\mathbb{R}^{d}$) generated by all compact subsets of $K$ with respect to the
Hausdorff distance and define a set-valued map
$\mathcal{H}_{K}:\Omega\to\mathcal{P}(\mathbb{R}^{d})$ by

$\mathcal{H}_{K}(\omega):=\omega\cap K.$

We also define a $\sigma$-algebra $\mathcal{F}_{K}$ (in $\Omega$) as

$\mathcal{F}_{K}:=\\{\mathcal{H}_{K}^{-1}(F_{K})\,|\,F_{K}\in\widetilde{\mathcal{F}}_{K}\\}.$

######  Assumption 5.6 (Finite range of dependence).

There exists a positive real number $\kappa>0$ such that for every pair of
compact sets $K_{j}\subset\mathbb{R}^{d}$, $j=1,2$, such that
$\operatorname{dist}(K_{1},K_{2})>\kappa$ the $\sigma$-algebras
$\mathcal{F}_{K_{1}}$ and $\mathcal{F}_{K_{2}}$ are independent, where
$\operatorname{dist}$ is the Hausdorff distance.

######  Theorem 5.7.

Under Assumption 5.6 we have

$\mathcal{G}\subseteq\lim_{\varepsilon\to
0}\sigma(\mathcal{A}^{\varepsilon})\,,$

which, in conjunction with Theorem 5.5, implies

$\lim_{\varepsilon\to 0}\sigma(\mathcal{A}^{\varepsilon})=\mathcal{G},$

almost surely.

The proof of Theorem 5.7 relies on the observation that given a
$\lambda\not\in\sigma(\mathcal{A}_{0})$, one can find an arbitrarily large
number of cubes of the same size such that

  1. (i)

these cubes are disjoint and their union is a (larger) cube,

  2. (ii)

the configuration of the set of inclusions in each of these cubes, apart from
the boundary layers of thickness $\kappa$ (the range of dependence), is almost
the same, and

  3. (iii)

the quantity $\beta_{\infty}(\lambda)$ is well approximated by the integral on
the RHS of (5.1) computed for each of these cubes.

The above intuition is formalised by the following proposition, whose proof
was given in [10].

######  Proposition 5.8 ([10, Thms. 5.18 and 5.19]).

Under Assumption 5.6, there exists a set of full measure
$\Omega_{1}\subset\Omega$ such that for every $\omega\in\Omega_{1}$ the
following holds. Suppose we are given $R_{0}>0$, $\delta>0$ and
$\lambda\not\in\sigma(\mathcal{A}_{0})$. Then there exists $R\geq R_{0}$, and
for any given $N\in\mathbb{N}$ there exists a collection of points
$\\{x_{n}\\}_{n=1}^{N^{d}}\subset\mathbb{R}^{d}$ and a point
$\overline{x}\in\mathbb{R}^{d}$ such that

$\operatorname{Int}\square_{x_{j}}^{R+\kappa}\cap\operatorname{Int}\square_{x_{k}}^{R+\kappa}=\emptyset\qquad\forall
j\neq k,\ j,k\in\\{1,\ldots,N^{d}\\},$ (5.34)
$\bigcup_{n=1}^{N^{d}}\square_{x_{n}}^{R+\kappa}=\square_{\overline{x}}^{N(R+\kappa)},$
(5.35) $\operatorname{dist}\left(\mathcal{H}_{\square_{0}^{R}}(\omega-
x_{j}),\mathcal{H}_{\square_{0}^{R}}(\omega-x_{k})\right)<\delta\qquad\forall
j,k\in\\{1,\ldots,N^{d}\\},$ (5.36)
$\left|\beta_{\infty}(\lambda)-\sup_{k\in\mathbb{S}^{d-1}}\frac{1}{|\square_{x_{n}}^{R+\kappa}|}\int_{\square_{x_{n}}^{R+\kappa}}\left(\lambda+\lambda^{2}\,k\cdot\mathbf{b}_{\lambda}(y,\omega)k\right)dy\right|<\delta\qquad\forall
n\in\\{1,\dots,N^{d}\\},$ (5.37)
$\left|\frac{1}{|\square_{x_{n}}^{R+\kappa}|}\int_{\square_{x_{n}}^{R+\kappa}}\left(\lambda\,\mathrm{I}+\lambda^{2}\,\mathbf{b}_{\lambda}(y,\omega)\right)dy-\frac{1}{|\square_{x_{m}}^{R+\kappa}|}\int_{\square_{x_{m}}^{R+\kappa}}\left(\lambda\,\mathrm{I}+\lambda^{2}\,\mathbf{b}_{\lambda}(y,\omega)\right)dy\right|<\delta\\\
\qquad\forall n,m\in\\{1,\dots,N^{d}\\}\,.$ (5.38)

######  Remark 5.9.

Note that the effect of the range of dependence manifests itself explicitly in
the bound (5.36), which is only true when one removes a boundary layer of
width $\kappa$ from our cubes $\square_{x_{n}}^{R+\kappa}$. Indeed,
distributions of inclusions in $\square_{x_{j}}^{R}$ and $\square_{x_{k}}^{R}$
are independent for all $j\neq k$, $j,k\in\\{1,\ldots,N^{d}\\}$ by Assumption
5.6.

###### Proof of Theorem 5.7.

Suppose we are given $\lambda_{0}\in\mathcal{G}$. The task at hand is to show
that

$\lim_{\varepsilon\to
0}\operatorname{dist}(\lambda_{0},\sigma(\mathcal{A}^{\varepsilon}))=0.$
(5.39)

Due to Corollary 4.13 and the continuity of $\beta_{\infty}$, see Proposition
5.3, it is sufficient to prove this for $\lambda_{0}$ such that
$\beta_{\infty}(\lambda_{0})>0$. The proof, rather long and technical, will be
broken into several steps for the sake of clarity.

Step 1: _Approximating $\beta_{\infty}(\lambda_{0})$ on almost-periodic
cubes_.

Let us fix sufficiently small $\delta>0$ and let $R_{0},L>0$ be sufficiently
large positive numbers and suppose we have fixed $\varepsilon>0$. Then by
Proposition 5.8 there exists $R\geq R_{0}$ such that the following holds:
there exist $\overline{x}\in\mathbb{R}^{d}$ and
$\\{x_{k}\\}_{k=1}^{N^{d}}\subset\mathbb{R}^{d}$ for which (5.34)–(5.37) are
satisfied, where $N=N(\varepsilon)$ is chosen to be the smallest positive
integer such that

$\varepsilon N(R+\kappa)\geq L.$ (5.40)

Note that $N(\varepsilon)=O(\varepsilon^{-1})$ as $\varepsilon\to 0$, and the
choice of the centres of the cubes $\overline{x}$, $\\{x_{k}\\}_{k=1}^{N^{d}}$
depends on $\varepsilon$, which we, however, do not reflect in the notation
for the sake of brevity. The above condition (5.40) ensures that

$\square_{\varepsilon\overline{x}}^{L}\subseteq\square_{\varepsilon\overline{x}}^{\varepsilon
N(R+\kappa)}=\bigcup_{k}\square_{\varepsilon{x_{k}}}^{\varepsilon(R+\kappa)}.$
(5.41)

We pick one of the cubes as the approximation basis for the construction that
follows. Since it is not important which cube we use, we pick
$\square_{x_{1}}^{R+\kappa}$. Let

$\bm{\ell}:=\frac{1}{|\square_{x_{1}}^{R+\kappa}|}\int_{\square_{x_{1}}^{R+\kappa}}(\lambda_{0}\,\mathrm{I}+\lambda_{0}^{2}\,\mathbf{b}_{\lambda_{0}}(y,\omega))\,dy\,.$

Then (5.37) can be recast as

$|\beta_{\infty}(\lambda_{0})-\lambda_{\mathrm{max}}(\bm{\ell})|<\delta,$

where $\lambda_{\mathrm{max}}(\bm{\ell})$ denotes the greatest eigenvalue of
the $d\times d$ matrix $\bm{\ell}$. Without loss of generality, we can assume
that $\lambda_{\mathrm{max}}(\bm{\ell})>0$. This can always be achieved by
choosing $\delta$ small enough.

Step 2: _Auxiliary periodic problem_.
Let

$\widetilde{\square}_{x_{1}}^{R+\kappa}:=\square_{x_{1}}^{R+\kappa}\setminus\left(\bigcup_{\omega^{k}\subset\subset\square_{x_{1}}^{R}}\omega^{k}\right)$
(5.42)

be the set obtained from $\square_{x_{1}}^{R+\kappa}$ by removing the
inclusions fully contained within the smaller cube $\square_{x_{1}}^{R}$.

For $j,k\in\\{1,\ldots,d\\}$, let
$\widehat{\mathbf{N}}^{(j,k)}\in\mathbf{H}^{1}_{\mathrm{per}}(\widetilde{\square}_{x_{1}}^{R+\kappa})$
be the unique zero-mean solution to

$\int_{\widetilde{\square}_{x_{1}}^{R+\kappa}}(\mathsf{C}_{1})_{\alpha\beta\mu\nu}\left(\delta_{\mu
j}\delta_{\nu
k}+\nabla_{\mu}\widehat{N}_{\nu}^{(j,k)}\right)\,\nabla_{\alpha}\varphi_{\beta}=0,\qquad\forall\bm{\varphi}\in\mathbf{H}^{1}_{\mathrm{per}}(\widetilde{\square}_{x_{1}}^{R+\kappa})\,.$
(5.43)

The quantity $\widehat{\mathbf{N}}$ is a tensor of order $d^{2}$ valued in
vector functions over $\mathbb{R}^{d}$, and it has the meaning of
_homogenisation corrector_ for the auxiliary periodic problem obtained by
covering $\mathbb{R}^{d}$ with copies of the perforated cubical domain
$\widetilde{\square}_{x_{1}}^{R+\kappa}$. Let us first extend
$\widehat{\mathbf{N}}^{(j,k)}$ to the whole of $\square_{x_{1}}^{R+\kappa}$ by
Theorem 2.7, and then to the whole of $\mathbb{R}^{d}$ by periodicity. With
slight abuse of notation, we denote the resulting object
$\widehat{\mathbf{N}}^{(j,k)}\in\mathbf{H}^{1}_{\mathrm{per}}(\square_{x_{1}}^{R+\kappa})$
by the same symbol.

We define the periodic homogenised tensor
$\widehat{\mathsf{C}}^{\mathrm{hom}}$ via

$\widehat{\mathsf{C}}^{\mathrm{hom}}\xi:=\frac{1}{|\square_{x_{1}}^{R+\kappa}|}\int_{\widetilde{\square}_{x_{1}}^{R+\kappa}}\mathsf{C}_{1}(\xi+\xi_{jk}\nabla\widehat{\mathbf{N}}^{(j,k)})\qquad\forall\xi\in\mathbb{R}^{d^{2}},\
\xi=\operatorname{sym}\xi\,.$ (5.44)

Step 3: _Estimates for the periodic corrector_.
Formula (5.43) immediately gives us

$\int_{\square_{x_{1}}^{R+\kappa}}\mathsf{C}_{1}\,\nabla\widehat{\mathbf{N}}^{(j,k)}\cdot\nabla\widehat{\mathbf{N}}^{(j,k)}=-\int_{\square_{x_{1}}^{R+\kappa}}(\mathsf{C}_{1})_{\alpha\beta
jk}\,\nabla_{\alpha}\widehat{N}_{\beta}^{(j,k)},\qquad
j,k\in\\{1,\ldots,d\\},$

(in the above formula there is no summation over $j$ and $k$). The latter and
(2.17) imply

$\|\operatorname{sym}\nabla\widehat{\mathbf{N}}^{(j,k)}\|_{[\mathbf{L}^{2}(\square_{x_{1}}^{R+\kappa})]^{d}}\leq
C(R+\kappa)^{\frac{d}{2}}.$ (5.45)

The Korn and Poincaré inequalities then give us

$\|\widehat{\mathbf{N}}^{(j,k)}\|_{\mathbf{L}^{2}(\square_{x_{1}}^{R+\kappa})}\leq
C(R+\kappa)^{\frac{d}{2}+1}.$ (5.46)

The above estimates are rather straightforward. However, further on in the
proof, in order to close the estimates, we will need higher (than
$\mathbf{H}^{1}$) regularity of the periodic corrector. The key technical
result in this regard is stated and proved in Appendix A.3 in the form of
Theorem A.4, which states that there exists a constant $C>0$ and $p>2$, both
independent of $R$, such that

$\|\nabla\widehat{\mathbf{N}}^{(j,k)}\|_{[\mathbf{L}^{p}(\square_{x_{1}}^{R+\kappa})]^{d}}\leq
C(R+\kappa)^{\frac{d}{p}}.$ (5.47)

Step 4: _Construction of “quasimodes”_.
Let $k_{0}\in\mathbb{S}^{d-1}$ be given and let us define a $d\times d$ matrix
$\widetilde{\mathsf{C}}=\widetilde{\mathsf{C}}(k_{0})$ component-wise as

$\widetilde{\mathsf{C}}_{\alpha\beta}:=(k_{0})_{\mu}\,\widehat{\mathsf{C}}^{\mathrm{hom}}_{\alpha\mu\nu\beta}\,(k_{0})_{\nu}\,,\qquad\alpha,\beta\in\\{1,\dots,d\\}\,.$
(5.48)

Observe that $\widetilde{\mathsf{C}}$ is symmetric and positive definite. In
the wave equation context the direct analogue of $k_{0}$ has the meaning of
the direction of propagation of a flat wave, cf. (5.49) below. It’s particular
choice is unimportant, and it won’t affect the final conclusion.

Let us choose $\varpi=\varpi(k_{0})\in\mathbb{R}$ and
$A=A(k_{0})\in\mathbb{S}^{d-1}$ in accordance with

$\varpi^{2}=\lambda_{\mathrm{max}}(\widetilde{\mathsf{C}}^{-1/2}\,\bm{\ell}\,\widetilde{\mathsf{C}}^{-1/2}),$
$A\in\ker(\varpi^{2}\widetilde{\mathsf{C}}-\bm{\ell}),$

and let us put

$\mathbf{u}(x):=A\cos(\varpi\,k_{0}\cdot x).$ (5.49)

A direct calculation shows that formulae (5.48)–(5.49) imply

$-\operatorname{div}\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla\mathbf{u}=\bm{\ell}\,\mathbf{u}\,.$
(5.50)

Let $\eta_{L,\varepsilon}(x):=\eta(x/L+\varepsilon\overline{x})$, where the
cut-off $\eta$ is defined as in (2.2). Then the normalised vector-functions

$\mathbf{u}_{L,\varepsilon}:=\frac{\eta_{L,\varepsilon}\mathbf{u}}{\|\eta_{L,\varepsilon}\mathbf{u}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}}$

are approximate solutions of (5.50), that is, they satisfy

$\|\operatorname{div}\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla\mathbf{u}_{L,\varepsilon}+\bm{\ell}\,\mathbf{u}_{L,\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq\frac{C}{L}.$
(5.51)

Note that the constant $C$ in the above estimate is independent of
$\varepsilon$. The vector-functions $\mathbf{u}_{L,\varepsilon}$ can be
thought of as “quasimodes” as $L\to\infty$, where quotation marks are needed
because the “spectral parameter” here is a matrix. Also note that (cf. (5.41))

$\operatorname{supp}\mathbf{u}_{L,\varepsilon}\subseteq\square_{\varepsilon\overline{x}}^{L}\subseteq\square_{\varepsilon\overline{x}}^{\varepsilon
N(R+\kappa)}\,.$ (5.52)

Out of $\mathbf{u}_{L,\varepsilon}$ one can construct candidates for
“quasimodes” for the operator $\mathcal{A}^{\varepsilon}$ as follows. Let
$\mathbf{b}_{\lambda_{0}}^{\varepsilon}(x):=\mathbf{b}_{\lambda_{0}}(x/\varepsilon,\omega)$
(recall that $\mathbf{b}_{\lambda_{0}}$ is defined by (4.14)) and define

$\mathbf{u}_{L}^{\varepsilon}:=(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\mathbf{u}_{L,\varepsilon}.$
(5.53)

While the functions $\mathbf{u}_{L}^{\varepsilon}$ are in the domain of the
form (2.21), they are not in the domain of $\mathcal{A}^{\varepsilon}$.
Therefore, we define

$\widehat{\mathbf{u}}^{\varepsilon}_{L}:=(\lambda_{0}+1)(\mathcal{A}^{\varepsilon}+1)^{-1}\mathbf{u}_{L}^{\varepsilon}$

so that

$(\mathcal{A}^{\varepsilon}-\lambda_{0})\widehat{\mathbf{u}}^{\varepsilon}_{L}=(\lambda_{0}+1)(\mathbf{u}^{\varepsilon}_{L}-\widehat{\mathbf{u}}^{\varepsilon}_{L}).$
(5.54)

The task at hand reduces to show that for every $\sigma>0$ one can choose
$\varepsilon_{0}>0$, $\delta>0$ and $L>0$ such that

$\|\mathbf{u}^{\varepsilon}_{L}-\widehat{\mathbf{u}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}<\sigma\qquad\forall\varepsilon<\varepsilon_{0}.$
(5.55)

The latter, combined with (5.54) and the fact that, clearly,

$\|\mathbf{u}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\geq C>0$

uniformly in $\varepsilon$ and $L$, gives us (5.39).

Step 6: _Adding the corrector_.
Rather than estimating
$\|\mathbf{u}^{\varepsilon}_{L}-\widehat{\mathbf{u}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}$
directly, it is convenient to estimate
$\|\mathbf{u}^{\varepsilon}_{LC}-\widehat{\mathbf{u}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}$
instead, where

$\mathbf{u}^{\varepsilon}_{LC}(x):=\mathbf{u}_{L}^{\varepsilon}(x)+\varepsilon\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}(x)\widehat{\mathbf{N}}^{(j,k)}(x/\varepsilon)$
(5.56)

is a modification of $\mathbf{u}_{L}^{\varepsilon}$ by adding the
homogenisation corrector. (Recall that $\widehat{\mathbf{N}}^{(j,k)}$ is not
the actual stochastic corrector for the operator $\mathcal{A}^{\varepsilon}$,
but the auxiliary periodic corrector introduced in Step 2.) Doing so will
yield the result we are after because, on account of the estimate (5.46) for
the corrector, we have

$\left|\|\mathbf{u}^{\varepsilon}_{L}-\widehat{\mathbf{u}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}-\|\mathbf{u}^{\varepsilon}_{LC}-\widehat{\mathbf{u}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\right|\leq\|\mathbf{u}^{\varepsilon}_{L}-\mathbf{u}^{\varepsilon}_{LC}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq
C\varepsilon(R+\kappa).$ (5.57)

Step 7: _Key bilinear form estimate_.
Let us denote by $a^{\varepsilon}$ the bilinear form associated with the
operator $\mathcal{A}^{\varepsilon}+\mathrm{Id}$:

$a^{\varepsilon}(\mathbf{u},\mathbf{v}):=\int_{\mathbb{R}^{d}}\left(\mathsf{C}^{\varepsilon}\nabla\mathbf{u}\,\cdot\nabla\mathbf{v}+\mathbf{u}\cdot\mathbf{v}\right)\,.$

Following a strategy proposed in [27] (which exploits, in turn, a general idea
found, e.g., in [39]), we claim that the task at hand reduces to proving the
key estimate

$|a^{\varepsilon}(\mathbf{u}^{\varepsilon}_{LC}-\widehat{\mathbf{u}}^{\varepsilon}_{L},\mathbf{v})|\leq
C\,\mathcal{E}(\varepsilon,\delta,L)\sqrt{a^{\varepsilon}(\mathbf{v},\mathbf{v})}\qquad\forall\mathbf{v}\in\mathbf{H}^{1}(\mathbb{R}^{d})$
(5.58)

for some ‘error’ $\mathcal{E}(\varepsilon,\delta,L)$ satisfying

$\lim_{\delta\to 0}\lim_{L\to+\infty}\lim_{\varepsilon\to
0}\mathcal{E}(\varepsilon,\delta,L)=0.$ (5.59)

Suppose (5.58) holds. Then for
$\mathbf{v}=\mathbf{u}^{\varepsilon}_{LC}-\widehat{\mathbf{u}}^{\varepsilon}_{L}$
one obtains

$\|\mathbf{u}^{\varepsilon}_{LC}-\widehat{\mathbf{u}}^{\varepsilon}_{L}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq\mathcal{E}(\varepsilon,\delta,L).$
(5.60)

Formulae (5.60), (5.59) and (5.57) imply (5.55). The remainder of the proof is
devoted to proving (5.58).

Given an arbitrary $\mathbf{v}\in\mathbf{H}^{1}(\mathbb{R}^{d})$, let
$\widetilde{\mathbf{v}}^{\varepsilon}$ be the extension of
$\mathbf{v}|_{\mathcal{M}^{\varepsilon}}$ into $\mathcal{I}^{\varepsilon}$ via
Theorem 2.7 and put
$\mathbf{v}_{0}^{\varepsilon}=\mathbf{v}-\widetilde{\mathbf{v}}^{\varepsilon}$.

By the properties of the extension, we have:

$\|\mathbf{v}_{0}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq
C\varepsilon\|\operatorname{sym}\nabla\mathbf{v}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\,,$
(5.61)
$\|\operatorname{sym}\nabla\mathbf{v}_{0}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq
C\,\|\operatorname{sym}\nabla\widetilde{\mathbf{v}}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}$
(5.62)

and

$\|\widetilde{\mathbf{v}}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}+\|\operatorname{sym}\nabla\widetilde{\mathbf{v}}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}+\varepsilon\|\operatorname{sym}\nabla\mathbf{v}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\leq
C\sqrt{a^{\varepsilon}(\mathbf{v},\mathbf{v})}\,.$ (5.63)

Let us consider the quantity

$a^{\varepsilon}(\mathbf{u}^{\varepsilon}_{LC},\mathbf{v})=\int_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}\nabla\mathbf{u}^{\varepsilon}_{LC}\cdot\nabla\mathbf{v}+\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\mathsf{C}_{0}\nabla\mathbf{u}^{\varepsilon}_{LC}\cdot\nabla\mathbf{v}+\int_{\mathbb{R}^{d}}\mathbf{u}^{\varepsilon}_{LC}\cdot\mathbf{v}.$
(5.64)

Step 8: _Estimates: part 1_.
By adding and subtracting
$\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla\mathbf{u}_{L,\varepsilon}\cdot\nabla\widetilde{\mathbf{v}}$,
the first term in the RHS of (5.64) can be equivalently recast as

$\int_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}\nabla\mathbf{u}^{\varepsilon}_{LC}\cdot\nabla\mathbf{v}=\int_{\mathcal{M}^{\varepsilon}}\left[\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla\mathbf{u}_{L,\varepsilon}+(\mathbbm{1}_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}(\delta_{j\cdot}\delta_{k\cdot}+\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon))-\widehat{\mathsf{C}}^{\mathrm{hom}}\delta_{j\cdot}\delta_{k\cdot})\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\right.\\\
\left.+\varepsilon\,\mathsf{C}_{1}(\nabla\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j})\otimes\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon)\right]\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}.$
(5.65)

In view of (5.46), the elementary estimate

$\|\nabla^{2}\mathbf{u}_{L,\varepsilon}\|_{[L^{\infty}(\mathbb{R}^{d})]^{d^{3}}}\leq
CL^{-d/2}$ (5.66)

and (5.40), the last term in the RHS of (5.65) can be estimated as

$\left|\int_{\mathcal{M}^{\varepsilon}}\left[\varepsilon\,\mathsf{C}_{1}(\nabla\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j})\otimes\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon)\right]\cdot\nabla\widetilde{\mathbf{v}}\right|\leq
C\varepsilon
R\,\|\operatorname{sym}\nabla\widetilde{\mathbf{v}}^{\varepsilon}\|_{\mathbf{L}(\mathbb{R}^{d})}\,.$
(5.67)

Next, let us denote by $\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}$ the
characteristic function of the set $\widetilde{\square}_{x_{1}}^{R+\kappa}$,
cf. (5.42), extended by periodicity to the whole of $\mathbb{R}^{d}$, and let
us rewrite

$\mathbbm{1}_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}(\delta_{j\cdot}\delta_{k\cdot}+\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon))-\widehat{\mathsf{C}}^{\mathrm{hom}}\delta_{j\cdot}\delta_{k\cdot}\\\
=(\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}})\mathsf{C}_{1}(\delta_{j\cdot}\delta_{k\cdot}+\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon))+\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}(\delta_{j\cdot}\delta_{k\cdot}+\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon))-\widehat{\mathsf{C}}^{\mathrm{hom}}\delta_{j\cdot}\delta_{k\cdot}\,.$
(5.68)

In view of (5.52), we will need to estimate
$\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}$
in $\square_{\varepsilon\overline{x}}^{L}$. We claim that, for $R$ large
enough, we have

$\int_{\square_{\varepsilon\overline{x}}^{L}}|\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}|\leq
C\delta L^{d}.$ (5.69)

The inequality (5.69) follows from:

  1. (i)

Formula (5.41);

  2. (ii)

The estimate

$\int_{\square_{\varepsilon x_{j}}^{\varepsilon
R}}|\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}|\leq
C\delta(\varepsilon\,R)^{d},\qquad j\in\\{1,\ldots,N^{d}\\},$

which is a consequence of (5.36) and the fact that the surface area of
inclusions in each $\square_{\varepsilon x_{j}}^{\varepsilon R}$,
$j\in\\{1,\ldots,N^{d}\\}$, is bounded above by the volume of the cube, up to
a constant (this follows immediately from Assumption 2.2);

  3. (iii)

The fact that

$\frac{\left|\square_{\varepsilon
x_{j}}^{\varepsilon(R+\kappa)}\setminus\square_{\varepsilon
x_{j}}^{\varepsilon R}\right|}{|\square_{\varepsilon\overline{x}}^{L}|}$

can be made arbitrarily small by choosing $R$ (at the start of the proof)
sufficiently large.

When substituting (5.68) into (5.65), the first term on the RHS of (5.68)
gives rise to two contributions:
$(\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}})\mathsf{C}_{1}\nabla\mathbf{u}_{L,\varepsilon}\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}$
and
$(\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}})\mathsf{C}_{1}(\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon)\,\partial_{k}(\mathbf{u}_{L,\varepsilon})_{j})\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}$.
The former can be estimated via (5.66) and (5.69) as

$\left|\int_{\mathbb{R}^{d}}(\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}})\mathsf{C}_{1}\nabla\mathbf{u}_{L,\varepsilon}\cdot\nabla\widetilde{\mathbf{v}}\right|\leq
C\,\delta^{\nicefrac{{1}}{{2}}}\|\operatorname{sym}\nabla\widetilde{\mathbf{v}}^{\varepsilon}\|_{[\mathbf{L}^{2}(\mathbb{R}^{d})]^{2}}\,.$
(5.70)

Combining again (5.66), (5.69) and the higher regularity of the corrector
(5.47) by means of Hölder’s inequality, and applying Korn’s inequality to
$\nabla\widetilde{\mathbf{v}}^{\varepsilon}$, we can estimate the latter as

$\left|\int_{\mathbb{R}^{d}}(\mathbbm{1}_{\mathcal{M}^{\varepsilon}}-\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}})\mathsf{C}_{1}(\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon)\,\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j})\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}\right|\leq
C\,\delta^{\frac{p-2}{2p}}\,\|\operatorname{sym}\nabla\widetilde{\mathbf{v}}^{\varepsilon}\|_{[\mathbf{L}^{2}(\mathbb{R}^{d})]^{2}}\,$
(5.71)

for some $p>2$.

We are left to deal with the last two terms on the RHS of (5.68). To this end,
we observe that, for fixed $\beta$, the vector field

$\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}(\mathsf{C}_{1})_{\alpha\beta\mu\nu}(\delta_{j\mu}\delta_{k\nu}+\nabla_{\mu}N^{(j,k)}_{\nu})-(\widehat{\mathsf{C}}^{\mathrm{hom}})_{\alpha\beta\mu\nu}\,\delta_{j\mu}\delta_{k\nu}$

(whose components are labelled by the index $\alpha$) is divergence free —
because of (5.43) — and has zero mean — because of (5.44). Hence, by a
classical construction, see e.g. [26], there exist $d^{3}$ skew-symmetric,
zero-mean matrix functions
$\mathbf{G}_{\cdot\,\cdot;\beta}^{(j,k)}\in[\mathbf{H}^{1}_{\mathrm{per}}(\square_{x_{1}}^{R+\kappa})]^{2}$,
$\beta,j,k\in\\{1,\dots,d\\}$, such that

$\partial_{\gamma}\mathbf{G}_{\gamma\alpha;\beta}^{(j,k)}=\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}(\mathsf{C}_{1})_{\alpha\beta\mu\nu}(\delta_{j\mu}\delta_{k\nu}+\nabla_{\mu}N^{(j,k)}_{\nu}(\cdot/\varepsilon))-(\widehat{\mathsf{C}}^{\mathrm{hom}})_{\alpha\beta\mu\nu}\,\delta_{j\mu}\delta_{k\nu}\,.$
(5.72)

Observe that the RHS of (5.72) is symmetric in the pair of indices $\alpha$
and $\beta$. Moreover, the explicit construction of
$\mathbf{G}_{\gamma\alpha;\beta}^{(j,k)}$ implies

$\|\mathbf{G}_{\cdot\,\cdot;\beta}^{(j,k)}\|_{[\mathbf{L}^{2}(\square_{x_{1}}^{R+\kappa})]^{2}}\\\
\leq
CR\|\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}(\mathsf{C}_{1})_{\cdot\,\cdot\mu\nu}(\delta_{j\mu}\delta_{k\nu}+\nabla_{\mu}N^{(j,k)}_{\nu}(\cdot/\varepsilon)-(\widehat{\mathsf{C}}^{\mathrm{hom}})_{\cdot\,\cdot\mu\nu}\,\delta_{j\mu}\delta_{k\nu}\|_{[\mathbf{L}^{2}(\square_{x_{1}}^{R+\kappa})]^{2}}\\\
\overset{\eqref{proof part2 equation 9}}{\leq}CR^{\nicefrac{{d}}{{2}}+1}\,.$
(5.73)

In view of (5.72), integrating by parts and using the antisymmetry of
$\mathbf{G}_{\cdot\,\cdot;\beta}^{(j,k)}$, we get

$\int_{\mathbb{R}^{d}}\left[(\widehat{\mathbbm{1}}_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}(\delta_{j\cdot}\delta_{k\cdot}+\nabla\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon))-\widehat{\mathsf{C}}^{\mathrm{hom}}\delta_{j\cdot}\delta_{k\cdot})\,\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\right]\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}\\\
=\int_{\mathbb{R}^{d}}\left[\varepsilon\,\partial_{\gamma}\mathbf{G}_{\gamma\alpha;\beta}^{(j,k)}\,\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\right]\,\partial_{\alpha}\widetilde{\mathbf{v}}^{\varepsilon}_{\beta}\\\
=-\int_{\mathbb{R}^{d}}\varepsilon\,\mathbf{G}_{\gamma\alpha;\beta}^{(j,k)}\,\partial_{\gamma}\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\,\partial_{\alpha}\widetilde{\mathbf{v}}^{\varepsilon}_{\beta}-\int_{\mathbb{R}^{d}}\underset{=0}{\underbrace{\varepsilon\,\mathbf{G}_{\gamma\alpha;\beta}^{(j,k)}\,\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\,\partial_{\gamma}\partial_{\alpha}\widetilde{\mathbf{v}}^{\varepsilon}_{\beta}}}\\\
=-\frac{1}{2}\int_{\mathbb{R}^{d}}\varepsilon\,\mathbf{G}_{\gamma\alpha;\beta}^{(j,k)}\,\partial_{\gamma}\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\,\left(\partial_{\alpha}\widetilde{\mathbf{v}}^{\varepsilon}_{\beta}+\partial_{\beta}\widetilde{\mathbf{v}}^{\varepsilon}_{\alpha}\right)\,.$
(5.74)

Hence, by combining (5.67), (5.68), (5.70), (5.71), (5.74), (5.66), (5.73),
and (5.63), we obtain that (5.65) can be rewritten as

$\int_{\mathcal{M}^{\varepsilon}}\mathsf{C}_{1}\nabla\mathbf{u}^{\varepsilon}_{LC}\cdot\nabla\mathbf{v}=\int_{\mathbb{R}^{d}}\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla\mathbf{u}_{L,\varepsilon}\cdot\nabla\mathbf{v}+\mathcal{R}_{1},$
(5.75)

where

$|\mathcal{R}_{1}|\leq
C\left[R\,\varepsilon+\delta^{\frac{p-2}{2p}}\right]\|\operatorname{sym}\nabla\widetilde{\mathbf{v}}^{\varepsilon}\|_{[\mathbf{L}^{2}(\mathbb{R}^{d})]^{2}}\,.$

Here we are assuming, without loss of generality, that $\delta<1$, so that
$\delta^{\frac{1}{2}}<\delta^{\frac{p-2}{2p}}$ for $p>2$.

Step 8: _Estimates: part 2_.
Let us now examine the second term on the RHS of (5.64).

We have

$\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\mathsf{C}_{0}\nabla\mathbf{u}^{\varepsilon}_{LC}\cdot\nabla\mathbf{v}=\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\mathsf{C}_{0}\left((1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\nabla\mathbf{u}_{L,\varepsilon}\right)\cdot\nabla\mathbf{v}\\\
+\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\lambda_{0}[\mathsf{C}_{0}]_{\alpha\beta\mu\nu}\,\partial_{\mu}[\mathbf{b}^{\varepsilon}_{\lambda_{0}}]_{\nu\gamma}\,[\mathbf{u}_{L,\varepsilon}]_{\gamma}\,\partial_{\alpha}[\mathbf{v}^{\varepsilon}_{0}]_{\beta}+\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\mathsf{C}_{0}\nabla(\varepsilon\partial_{k}[\mathbf{u}_{L,\varepsilon}]_{j}\widehat{\mathbf{N}}^{(j,k)}(\cdot/\varepsilon))\cdot\nabla\mathbf{v}\,.$
(5.76)

Now, integrating by parts and taking into account (4.17), the second
contribution to the RHS can be recast as

$\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\lambda_{0}\mathsf{C}_{0}\,\nabla\mathbf{b}^{\varepsilon}_{\lambda_{0}}\,\mathbf{u}_{L,\varepsilon}\cdot\nabla(\mathbf{v}^{\varepsilon}_{0}+\widetilde{\mathbf{v}}^{\varepsilon})\\\
=\int_{\mathcal{I}^{\varepsilon}}\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\,\mathbf{u}_{L,\varepsilon}\cdot\mathbf{v}^{\varepsilon}_{0}-\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\lambda_{0}[\mathsf{C}_{0}]_{\alpha\beta\mu\nu}\,\partial_{\mu}[\mathbf{b}^{\varepsilon}_{\lambda_{0}}]_{\nu\gamma}\,\partial_{\alpha}[\mathbf{u}_{L,\varepsilon}]_{\gamma}\,[\mathbf{v}^{\varepsilon}_{0}]_{\beta}\\\
+\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\lambda_{0}\mathsf{C}_{0}\left(\nabla\mathbf{b}^{\varepsilon}_{\lambda_{0}}\,\mathbf{u}_{L,\varepsilon}\right)\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}\,.$
(5.77)

Substituting (5.77) into (5.76) and using (5.45), (5.46), (4.16), and
(5.61)–(5.63) we arrive at

$\varepsilon^{2}\int_{\mathcal{I}^{\varepsilon}}\mathsf{C}_{0}\nabla\mathbf{u}^{\varepsilon}_{LC}\cdot\nabla\mathbf{v}=\int_{\mathcal{I}^{\varepsilon}}\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\,\mathbf{u}_{L,\varepsilon}\cdot\mathbf{v}^{\varepsilon}_{0}+\mathcal{R}_{2},$
(5.78)

where

$|\mathcal{R}_{2}|\leq
C(\varepsilon+\varepsilon^{2}\,R)\,\sqrt{a^{\varepsilon}(\mathbf{v},\mathbf{v})}\,.$

Step 9: _Back to the bilinear form_.
Substituting (5.75) and (5.78) into (5.64), we obtain

$a^{\varepsilon}(\mathbf{u}^{\varepsilon}_{LC},\mathbf{v})=\int_{\mathbb{R}^{d}}\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla\mathbf{u}_{L,\varepsilon}\cdot\nabla\widetilde{\mathbf{v}}^{\varepsilon}+\int_{\mathbb{R}^{d}}\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\,\mathbf{u}_{L,\varepsilon}\cdot\mathbf{v}^{\varepsilon}_{0}+\int_{\mathbb{R}^{d}}\mathbf{u}^{\varepsilon}_{LC}\cdot\mathbf{v}+\mathcal{R}_{3}\,,$
(5.79)

where

$|\mathcal{R}_{3}|\leq
C\left[(1+R)\varepsilon+\delta^{\frac{p-2}{2p}}\right]\,\sqrt{a^{\varepsilon}(\mathbf{v},\mathbf{v})}\,.$
(5.80)

With (5.51) in mind, and recalling (5.56) and (5.53), we rewrite (5.79) as

$a^{\varepsilon}(\mathbf{u}^{\varepsilon}_{LC},\mathbf{v})=\int_{\mathbb{R}^{d}}\left(-\operatorname{div}\widehat{\mathsf{C}}^{\mathrm{hom}}\nabla-\bm{\ell}\right)\mathbf{u}_{L,\varepsilon}\cdot\widetilde{\mathbf{v}}^{\varepsilon}\\\
+\int_{\mathbb{R}^{d}}\left[\bm{\ell}-\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\right]\,\mathbf{u}_{L,\varepsilon}\cdot\widetilde{\mathbf{v}}^{\varepsilon}\\\
+\int_{\mathbb{R}^{d}}\left[\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\mathbf{u}_{L,\varepsilon}+\mathbf{u}^{\varepsilon}_{LC}\right]\cdot\mathbf{v}+\mathcal{R}_{3}\,.$
(5.81)

By combining formulae (5.81), (5.80), (5.51), and using (5.56) and (5.54), we
obtain

$|a^{\varepsilon}(\mathbf{u}^{\varepsilon}_{LC}-\widehat{\mathbf{u}}^{\varepsilon}_{L},\mathbf{v})|\leq
C\left[(1+R)\varepsilon+\delta^{\frac{p-2}{2p}}+L^{-1}\right]\,\sqrt{a^{\varepsilon}(\mathbf{v},\mathbf{v})}\\\
+\left|\int_{\mathbb{R}^{d}}\left[\bm{\ell}-\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\right]\,\mathbf{u}_{L,\varepsilon}\cdot\widetilde{\mathbf{v}}^{\varepsilon}\right|\,.$
(5.82)

Therefore, the remaining task is to estimate the second term on the RHS of the
latter inequality.

Step 10: _Estimates: part 3_.
The argument consists in utilising (5.51) while replacing
$\mathbf{u}_{L,\varepsilon}\otimes\widetilde{\mathbf{v}}^{\varepsilon}$ with
its average values on the cubes of size $\varepsilon(R+\kappa)$, cf. (5.41),
and using Poincaré’s inequality to control the error.

$\left|\int_{\mathbb{R}^{d}}\left(\bm{\ell}-\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\right)\,\mathbf{u}_{L,\varepsilon}\cdot\widetilde{\mathbf{v}}^{\varepsilon}\right|\\\
\leq\sum_{i=1}^{N^{d}}\left|\int_{\varepsilon\square_{\varepsilon
x_{i}}^{R+\kappa}}\left(\bm{\ell}-\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\right)\cdot\fint_{\varepsilon\square_{\varepsilon
x_{i}}^{R+\kappa}}\mathbf{u}_{L,\varepsilon}\otimes\widetilde{\mathbf{v}}^{\varepsilon}\right|\\\
+\left|\sum_{i=1}^{N^{d}}\int_{\varepsilon\square_{\varepsilon
x_{i}}^{R+\kappa}}\left(\bm{\ell}-\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\right)\cdot\left(\mathbf{u}_{L,\varepsilon}\otimes\widetilde{\mathbf{v}}^{\varepsilon}-\fint_{\varepsilon\square_{\varepsilon
x_{i}}^{R+\kappa}}\mathbf{u}_{L,\varepsilon}\otimes\widetilde{\mathbf{v}}^{\varepsilon}\right)\,\right|\\\
\overset{\eqref{proposition cubes equaiton
5}}{\leq}C\,\delta\,\|\widetilde{\mathbf{v}}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}\\\
+C\varepsilon(R+\kappa)\|\bm{\ell}-\lambda_{0}(1+\lambda_{0}\mathbf{b}^{\varepsilon}_{\lambda_{0}})\|_{\mathbf{L}^{2}(\square_{\varepsilon\overline{x}}^{\varepsilon
N(R+\kappa)})}\left\|\nabla\left(\mathbf{u}_{L,\varepsilon}\otimes\widetilde{\mathbf{v}}^{\varepsilon}\right)\right\|_{\mathbf{L}^{2}(\square_{\varepsilon\overline{x}}^{\varepsilon
N(R+\kappa)})}\\\ \overset{\eqref{lemma properties b lambda equation
1}}{\leq}C\left[\varepsilon(R+\kappa)+\delta\right]\left(\|\widetilde{\mathbf{v}}^{\varepsilon}\|_{\mathbf{L}^{2}(\mathbb{R}^{d})}+\|\operatorname{sym}\widetilde{\nabla}\mathbf{v}^{\varepsilon}\|_{[\mathbf{L}^{2}(\mathbb{R}^{d})]^{d}}\right)\,.$
(5.83)

In the last step we also used the Korn inequality in
$\square_{\varepsilon\overline{x}}^{\varepsilon N(R+\kappa)}$, and then the
domain monotonicity of the $L^{2}$-norm.

Step 11: _Conclusion_.
By combining (5.82), (5.83) and (5.63) we arrive at (5.58) with

$\mathcal{E}(\varepsilon,\delta,L)=(1+R)\varepsilon+\delta^{\frac{p-2}{2p}}+L^{-1}\,,$

which clearly satisfies (5.59).

This concludes the proof. ∎

##  6 Examples

In this section we will discuss some explicit examples, to showcase the
different properties of $\bm{\beta}(\lambda)$ and $\beta_{\infty}(\lambda)$,
as well as to compare and contrast the stochastic and the periodic settings.

###  6.1 Example 1: statistical micro-symmetries

Let us work in Euclidean space $\mathbb{R}^{3}$ equipped with Cartesian
coordinates $x_{j}$, $j\in\\{1,2,3\\}$.

Let us assume, for simplicity, that for almost every $\omega\in\Omega$ each
inclusion $\omega^{k}$ belongs to a finite set of shapes, up to rigid
translations and rotations, and that the tensor $\mathsf{C}_{0}$ is isotropic,
i.e. it acts on 2-tensors as (see, e.g., [15])

$\mathsf{C}_{0}\xi=c_{1}(\textrm{tr}\xi)I+c_{2}\operatorname{sym}\xi\,$ (6.1)

for some $c_{1},c_{2}\in\mathbb{R}$.

Let us denote by $\omega^{0}$ the connected component of $\omega$ containing
the origin of $\mathbb{R}^{3}$. In the case when $0\notin\omega$ we set
$\omega^{0}=\emptyset$. It was shown in [10, Lemma B.4]999The probabilistic
setup of [10] is slightly different from that of the current paper, but it is
equivalent to it. We refer the reader to [8, Remark 2.18] for a detailed
comparison of the two. that the set-valued mapping
$P_{0}:\omega\mapsto\omega^{0}$ is measurable, when the target space is
endowed with the $\sigma$-algebra $\mathcal{H}$ generated by the topology on
$[-2,2]^{d}$ induced by the Hausdorff distance. Observe that, by assumption,
the mapping $P_{0}$ takes values in a finite set of shapes, up to rigid
translations and rotations.

Let us consider the following two types of discrete isometries of
$\mathbb{R}^{3}$.

  1. (i)

Reflections $\\{S_{j}\\}_{j=1}^{3}$: $S_{j}$ describes the reflection across
the plane $\\{x_{j}=0\\}$.

  2. (ii)

Rotations $\\{R_{j}^{\pi/2}\\}_{j=1}^{3}$: $R_{j}^{\pi/2}$ describes the
counterclockwise rotation by $\pi/2$ about the $x_{j}$ coordinate axis.

Suppose that the inclusions possess the following additional statistical
micro-symmetries:

$\displaystyle P(\omega^{0}\in\mathfrak{H})=P(\omega^{0}\in
S_{j}\mathfrak{H}),$
$\displaystyle\qquad\forall\mathfrak{H}\in\mathcal{H},\forall j\in\\{1,2,3\\}$
(6.2)

and

$\displaystyle P(\omega^{0}\in\mathfrak{H})=P(\omega^{0}\in
R^{\pi/2}_{j}\mathfrak{H}),$
$\displaystyle\qquad\forall\mathfrak{H}\in\mathcal{H},\forall
j\in\\{1,2,3\\}.$ (6.3)

Here $S_{j}\mathfrak{H}$ and $R^{\pi/2}_{j}\mathfrak{H}$ denote the sets
obtained by applying the appropriate symmetry on each element of
$\mathfrak{H}$. The prefix “micro” refers to the fact that (6.2) and (6.3) are
concerned with statistical properties of individual inclusions, rather than
with the collective properties of the set $\omega$. In plain English, the
above conditions encode the fact that, on average, we have as many inclusions
containing the origin with a given spatial orientation as we have with
reflected — in the case of (6.2) — or with $\pi/2$-rotated — in the case of
(6.3) — about any of the three axes.

For $\lambda\not\in\sigma(\mathcal{A}_{0})$ and $c\in\mathbb{R}^{3}$, let
$\mathbf{b}^{c}_{\lambda,\omega_{0}}\in H_{0}^{1}(\omega^{0};\mathbb{R}^{3})$
be the unique solution to

$-\operatorname{div}\mathsf{C}_{0}\nabla\mathbf{b}^{c}_{\lambda,\omega^{0}}-\lambda\,\mathbf{b}^{c}_{\lambda,\omega^{0}}=c.$

Note that by Lemma 4.6 we have

$\overline{\mathbf{b}^{(i)}_{\lambda}}(\omega)=\mathbf{b}^{e_{i}}_{\lambda,\omega^{0}}(0)\,.$
(6.4)