The structure of fluctuations in homogenization

# The structure of fluctuations in stochastic homogenization

Mitia Duerinckx Antoine Gloria  and  Felix Otto Université Libre de Bruxelles (ULB)
Brussels, Belgium
Sorbonne Université, UMR 7598, Laboratoire Jacques-Louis Lions, Paris, France & Université Libre de Bruxelles (ULB), Brussels, Belgium Max Planck Institute for Mathematics in the Sciences
Leipzig, Germany
###### Abstract.

Four quantities are fundamental in homogenization of elliptic systems in divergence form and in its applications: the field and the flux of the solution operator (applied to a general deterministic right-hand side), and the field and the flux of the corrector. Homogenization is the study of the large-scale properties of these objects. For random coefficients, these quantities fluctuate and their fluctuations are a priori unrelated. Depending on the law of the coefficient field, and in particular on the decay of its correlations on large scales, these fluctuations may display different scalings and different limiting laws (if any). In this contribution, we identify a fifth and crucial intrinsic quantity, a random 2-tensor field, which we refer to as the homogenization commutator. In the model framework of discrete linear elliptic equations in divergence form with independent and identically distributed coefficients, we show what we believe to be a general principle, namely that the homogenization commutator drives at leading order the fluctuations of each of the four other quantities in probability, which reveals the pathwise structure of fluctuations in stochastic homogenization. In addition, we show in this framework that the (rescaled) homogenization commutator converges in law to a (2-tensor) Gaussian white noise, the distribution of which is thus characterized by some 4-tensor, and we analyze to which precision this tensor can be extracted from the representative volume element method. All these results are optimally quantified and hold in any dimension. This constitutes the first complete theory of fluctuations in stochastic homogenization. Extensions to the (non-symmetric) continuum setting are discussed, while details are postponed to forthcoming work.

## 1. Introduction

This article constitutes the first part of a series of works that develops a theory of fluctuations in stochastic homogenization of elliptic (non-necessarily symmetric) systems. In this first part, using an elementary approach, we provide a complete picture of our theory (with optimal error estimates in terms of convergence rates) in the model framework of discrete elliptic equations with independent and identically distributed (iid) conductances. Links to the literature are discussed in Section 1.3 below, while the extension to more general situations is postponed to forthcoming work and is shortly described in Section 1.4.

### 1.1. General overview

Although in the sequel we shall focus on the case of discrete equations, we use non-symmetric continuum notation in this introduction. Let \boldsymbol{a} be a stationary and ergodic random coefficient field on \mathbb{R}^{d} that satisfies the boundedness and ellipticity properties

for some \lambda>0. For all \varepsilon>0 we set \boldsymbol{a}_{\varepsilon}:=\boldsymbol{a}(\frac{\cdot}{\varepsilon}), and for all deterministic vector fields f\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, we consider the random family (u_{\varepsilon})_{\varepsilon>0} of unique Lax-Milgram solutions in \mathbb{R}^{d} (which, in the rest of this article, means the unique weak solutions in \dot{H}^{1}(\mathbb{R}^{d})) of the rescaled problems

 \displaystyle-D\cdot\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}\,=\,D\cdot f, (1.1)

where D denotes the continuum gradient (while the notation \nabla is reserved in the sequel for the discrete gradient). It is known since the pioneering work of Papanicolaou and Varadhan [PapaVara] and of Kozlov [Kozlov-79] that, almost surely, u_{\varepsilon} converges weakly (in \dot{H}^{1}(\mathbb{R}^{d})) as \varepsilon\downarrow 0 to the unique Lax-Milgram solution \bar{u} in \mathbb{R}^{d} of

 \displaystyle-D\cdot\bar{\boldsymbol{a}}D\bar{u}\,=\,D\cdot f, (1.2)

where \bar{\boldsymbol{a}} is a deterministic and constant matrix that only depends on \boldsymbol{a}. More precisely, for any direction e\in\mathbb{R}^{d}, the projection \bar{\boldsymbol{a}}e is the expectation of the flux of the corrector in the direction e,

 \bar{\boldsymbol{a}}e\,=\,\mathbb{E}\left[\boldsymbol{a}(D\phi_{e}+e)\right],

where the corrector \phi_{e} is the unique (up to a random additive constant) almost-sure solution of the corrector equation in \mathbb{R}^{d},

 -D\cdot\boldsymbol{a}(D\phi_{e}+e)\,=\,0,

in the class of functions the gradient of which is stationary and has finite second moment. We denote by \phi=(\phi_{i})_{i=1}^{d} the vector field the entries of which are the correctors \phi_{i} in the canonical directions e_{i} of \mathbb{R}^{d}. Note that the convergence of Du_{\varepsilon} to D\bar{u} in \operatorname{L}^{2}(\mathbb{R}^{d})^{d} is only weak since Du_{\varepsilon} typically displays spatial oscillations at scale \varepsilon, which are not captured by the limit D\bar{u}. These oscillations are however well-described by those of the corrector field D\phi(\frac{\cdot}{\varepsilon}) through the two-scale expansion (we systematically use Einstein’s summation rule on repeated indices)

 \displaystyle Du_{\varepsilon}\,\approx\,(D\phi_{i}(\tfrac{\cdot}{\varepsilon}% )+e_{i})D_{i}\bar{u}, (1.3)

in the sense that Du_{\varepsilon}-(D\phi_{i}(\frac{\cdot}{\varepsilon})+e_{i})D_{i}\bar{u} converges strongly to zero in \operatorname{L}^{2}(\mathbb{R}^{d})^{d}. In the random setting, this theory of oscillations was recently optimally quantified in [GNO2], [GNO-quant, Theorem 3], [GO4, Corollary 3], and [AKM-book, Chapter 6].

As opposed to periodic homogenization, which boils down to the sole understanding of the (spatial) oscillations of Du_{\varepsilon}, the stochastic setting involves the (random) fluctuations of Du_{\varepsilon} on top of its oscillations. More precisely, whereas oscillations are concerned with the (almost sure) lack of strong compactness for Du_{\varepsilon} in \operatorname{L}^{2}(\mathbb{R}^{d})^{d}, fluctuations are concerned with the leading-order probabilistic behavior of weak-type expressions of the form \int g\cdot Du_{\varepsilon} for g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}. Let us emphasize that in the case of a weakly correlated coefficient field \boldsymbol{a} the error in the two-scale expansion (1.3) is of order \varepsilon in \operatorname{L}^{2}(\mathbb{R}^{d})^{d} (or \varepsilon|\!\log\varepsilon|^{\frac{1}{2}} for d=2) while fluctuations of Du_{\varepsilon} display the central limit theorem (CLT) scaling \varepsilon^{\frac{d}{2}}, so that (1.3) is not expected to be accurate in that scaling. This was indeed first checked in dimension d\geq 3 by Gu and Mourrat [GuM, Section 3.2] (see also the last item in Remarks 2.12 below for d=2), who further argue that accuracy in (1.3) in the fluctuation scaling cannot even be reached by the use of higher-order correctors. The corrector field D\phi is therefore the driving quantity for oscillations but a priori not for fluctuations.

In the present article, we develop a complete theory of fluctuations in stochastic homogenization in line with the known theory of oscillations, and our main achievement is the identification of the driving quantity for fluctuations. The key in our theory consists in focusing on the homogenization commutator of the solution \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}}Du_{\varepsilon} and in studying its relation to the (standard) homogenization commutator \Xi:=(\Xi_{i})_{i=1}^{d} defined by

This stationary random (non-symmetric) 2-tensor field \Xi enjoys the following three crucial properties, which lead to our complete theory of fluctuations:

1. First and most importantly, the two-scale expansion of the homogenization commutator of the solution

 \displaystyle\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}}% Du_{\varepsilon}-\mathbb{E}\left[\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-% \bar{\boldsymbol{a}}Du_{\varepsilon}\right]\,\approx\,\Xi_{i}(\tfrac{\cdot}{% \varepsilon})D_{i}\bar{u} (1.5)

is (generically) accurate in the fluctuation scaling in the sense of

 \displaystyle\hskip 31.298031pt\mathbb{E}\left[\Big{|}\int_{\mathbb{R}^{d}}g% \cdot\big{(}\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}}% Du_{\varepsilon}-\mathbb{E}\left[\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-% \bar{\boldsymbol{a}}Du_{\varepsilon}\right]\big{)}-\int_{\mathbb{R}^{d}}g\cdot% \Xi_{i}(\tfrac{\cdot}{\varepsilon})D_{i}\bar{u}\Big{|}^{2}\right]^{\frac{1}{2}% }\\ \displaystyle\leq\leavevmode\nobreak\ o(1)\,\mathbb{E}\left[\Big{|}\int_{% \mathbb{R}^{d}}g\cdot\Xi_{i}(\tfrac{\cdot}{\varepsilon})D_{i}\bar{u}\Big{|}^{2% }\right]^{\frac{1}{2}}, (1.6)

where o(1)\downarrow 0 as \varepsilon\downarrow 0, for all g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}. Let us emphasize again that this property is nontrivial and is due to the form of the commutator.

2. Second, both the fluctuations of the field Du_{\varepsilon} and of the flux \boldsymbol{a}_{\varepsilon}Du_{\varepsilon} can be recovered through deterministic projections of those of the homogenization commutator \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}}Du_{\varepsilon}, which shows that no information is lost by passing to the homogenization commutator. More precisely,

 \displaystyle        \int_{\mathbb{R}^{d}}g\cdot D(u_{\varepsilon}-\mathbb{E}% \left[u_{\varepsilon}\right])=-\int_{\mathbb{R}^{d}}(\bar{\mathcal{P}}_{H}^{*}% g)\cdot\big{(}\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}% }Du_{\varepsilon}-\mathbb{E}\big{[}\boldsymbol{a}_{\varepsilon}Du_{\varepsilon% }-\bar{\boldsymbol{a}}Du_{\varepsilon}\big{]}\big{)}, (1.7) \displaystyle        \int_{\mathbb{R}^{d}}g\cdot(\boldsymbol{a}_{\varepsilon}% Du_{\varepsilon}-\mathbb{E}\left[\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}% \right])=\int_{\mathbb{R}^{d}}(\bar{\mathcal{P}}_{L}^{*}g)\cdot\big{(}% \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}}Du_{% \varepsilon}-\mathbb{E}\big{[}\boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-% \bar{\boldsymbol{a}}Du_{\varepsilon}\big{]}\big{)},

in terms of the Helmholtz and Leray projections in \operatorname{L}^{2}(\mathbb{R}^{d})^{d},

where \bar{\boldsymbol{a}}^{*} denotes the transpose of \bar{\boldsymbol{a}}. In addition, the fluctuations of the field D\phi(\tfrac{\cdot}{\varepsilon}) and of the flux \boldsymbol{a}_{\varepsilon}D\phi(\tfrac{\cdot}{\varepsilon}) of the corrector are clearly determined by those of the standard commutator \Xi(\tfrac{\cdot}{\varepsilon}). Indeed, the definition of \Xi leads to -D\cdot\bar{\boldsymbol{a}}D\phi_{i}=D\cdot\Xi_{i} and \boldsymbol{a}(D\phi_{i}+e_{i})-\bar{\boldsymbol{a}}e_{i}=\Xi_{i}+\bar{% \boldsymbol{a}}D\phi_{i}, to the effect of D\phi_{i}=-\bar{\mathcal{P}}_{H}\Xi_{i} and \boldsymbol{a}(D\phi_{i}+e_{i})-\bar{\boldsymbol{a}}e_{i}=(\operatorname{Id}-% \bar{\boldsymbol{a}}\bar{\mathcal{P}}_{H})\Xi_{i} in the stationary sense, and hence, formally,

 \displaystyle\int_{\mathbb{R}^{d}}F:D\phi(\tfrac{\cdot}{\varepsilon}) \displaystyle= \displaystyle-\int_{\mathbb{R}^{d}}\bar{\mathcal{P}}_{H}^{*}F:\Xi(\tfrac{\cdot% }{\varepsilon}), \displaystyle\int_{\mathbb{R}^{d}}F:\big{(}\boldsymbol{a}_{\varepsilon}(D\phi(% \tfrac{\cdot}{\varepsilon})+\operatorname{Id})-\bar{\boldsymbol{a}}\big{)} \displaystyle= \displaystyle\int_{\mathbb{R}^{d}}\mathcal{P}_{L}^{*}F:\Xi(\tfrac{\cdot}{% \varepsilon}), (1.9)

where \bar{\mathcal{P}}_{H}^{*} and \bar{\mathcal{P}}_{L}^{*} act on the second index of the tensor field F. A suitable sense to these identities is given as part of Corollary 2.4 below.
Let us highlight the pathwise structure of fluctuations revealed here. Combined with (1.6), identities (1.7) and (1.9) imply that the fluctuations of Du_{\varepsilon}, \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}, D\phi(\tfrac{\cdot}{\varepsilon}), and \boldsymbol{a}_{\varepsilon}D\phi(\tfrac{\cdot}{\varepsilon}) are determined at leading order by those of \Xi(\tfrac{\cdot}{\varepsilon}) in a strong norm in probability. This almost sure (“pathwise” in the language of SPDE) relation thus reduces the leading-order fluctuations of all quantities of interest to those of the sole homogenization commutator \Xi in a pathwise sense. Besides its theoretical importance, this pathwise structure is bound to affect multi-scale computing and uncertainty quantification in an essential way.

3. Third, the standard homogenization commutator \Xi is an approximately local function of the coefficients \boldsymbol{a}, which allows to infer the large-scale behavior of \Xi from the large-scale behavior of \boldsymbol{a} itself. This locality is best seen when formally computing partial derivatives of \Xi with respect to \boldsymbol{a}: letting \phi^{*} denote the corrector associated with the pointwise transpose coefficient field \boldsymbol{a}^{*}, and letting \sigma^{*} denote the corresponding flux corrector (cf. (3.1) below), we obtain (cf. (3.18))

 \displaystyle\hskip 28.452756pt\frac{\partial}{\partial\boldsymbol{a}(x)}\Xi_{% ij}=(D\phi_{j}^{*}+e_{j})\cdot\frac{\partial\boldsymbol{a}}{\partial% \boldsymbol{a}(x)}(D\phi_{i}+e_{i})\\ \displaystyle-D\cdot\bigg{(}\phi_{j}^{*}\frac{\partial\boldsymbol{a}}{\partial% \boldsymbol{a}(x)}(D\phi_{i}+e_{i})\bigg{)}-D\cdot\bigg{(}(\phi_{j}^{*}% \boldsymbol{a}+\sigma^{*}_{j})\frac{\partial D\phi_{i}}{\partial\boldsymbol{a}% (x)}\bigg{)}. (1.10)

In view of \frac{\partial\boldsymbol{a}}{\partial\boldsymbol{a}(x)}=\delta(\cdot-x), the first right-hand side term reveals an exactly local dependence upon \boldsymbol{a}. The second term is exactly local as well, but since it is written in divergence form its contribution is negligible when integrating on large scales. The only non-local effect comes from the last term due to \frac{\partial D\phi}{\partial\boldsymbol{a}}, which is given by the mixed derivative of the Green’s function for -D\cdot\boldsymbol{a}D and thus is expected to have only borderline integrable decay. However, it also appears inside a divergence, hence is negligible on large scales. In fact, in this paper, the accuracy in (1.5) is established relying on a similar representation of \frac{\partial}{\partial\boldsymbol{a}}(\boldsymbol{a}_{\varepsilon}Du_{% \varepsilon}-\bar{\boldsymbol{a}}Du_{\varepsilon}-\Xi_{i}(\tfrac{\cdot}{% \varepsilon})D_{i}\bar{u}), cf. Lemma 3.2.

We quickly comment on the form of the homogenization commutator. As well-known in applications, homogenization is the rigorous version of averaging fields and fluxes in a consistent way. This is made precise in the very definition of H-convergence by Murat and Tartar [MuratTartar], which requires both weak convergence of the fields Du_{\varepsilon}\rightharpoonup D\bar{u} and of the fluxes \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}\rightharpoonup\bar{\boldsymbol{a}% }D\bar{u} in \operatorname{L}^{2}(\mathbb{R}^{d})^{d}, to the effect of

 \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\bar{\boldsymbol{a}}Du_{% \varepsilon}\rightharpoonup 0.

This weak convergence of the homogenization commutator is the mathematical formulation of the so-called Hill-Mandel relation in mechanics [H63, H72]. Applied to the corrector, this justifies the definition of the (standard) homogenization commutator \Xi, cf. (1.4), which is seen as a natural and intrinsic measure of the accuracy of homogenization for large-scale averages.

###### Remark 1.1.

If a suitable rescaling \varepsilon^{-\frac{\beta}{2}}\Xi(\tfrac{\cdot}{\varepsilon}) of the homogenization commutator converges in law to some random linear functional \Gamma, then combining (1.6) with identities (1.7) and (1.9) leads to the joint convergence in law

 \displaystyle\bigg{(}\varepsilon^{-\frac{\beta}{2}}\int_{\mathbb{R}^{d}}F:\Xi(% \tfrac{\cdot}{\varepsilon})\leavevmode\nobreak\ ,\leavevmode\nobreak\ % \varepsilon^{-\frac{\beta}{2}}\int_{\mathbb{R}^{d}}g\cdot(Du_{\varepsilon}-% \mathbb{E}\left[Du_{\varepsilon}\right])\leavevmode\nobreak\ ,\leavevmode% \nobreak\ \varepsilon^{-\frac{\beta}{2}}\int_{\mathbb{R}^{d}}g\cdot(% \boldsymbol{a}_{\varepsilon}Du_{\varepsilon}-\mathbb{E}\left[\boldsymbol{a}_{% \varepsilon}Du_{\varepsilon}\right])\,, \displaystyle                                           \varepsilon^{-\frac{% \beta}{2}}\int_{\mathbb{R}^{d}}F:D\phi(\tfrac{\cdot}{\varepsilon})\leavevmode% \nobreak\ ,\leavevmode\nobreak\ \varepsilon^{-\frac{\beta}{2}}\int_{\mathbb{R}% ^{d}}F:\big{(}\boldsymbol{a}_{\varepsilon}(D\phi(\tfrac{\cdot}{\varepsilon})+% \operatorname{Id})-\bar{\boldsymbol{a}}\big{)}\bigg{)} \displaystyle     \,\to\,\Big{(}\Gamma(F)\leavevmode\nobreak\ ,\leavevmode% \nobreak\ \Gamma\big{(}\bar{\mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g% \big{)}\leavevmode\nobreak\ ,\leavevmode\nobreak\ -\Gamma(\bar{\mathcal{P}}_{H% }f\otimes\bar{\mathcal{P}}_{L}^{*}g)\leavevmode\nobreak\ ,\leavevmode\nobreak% \ -\Gamma\big{(}\bar{\mathcal{P}}_{H}^{*}F\big{)}\leavevmode\nobreak\ ,% \leavevmode\nobreak\ \Gamma(\bar{\mathcal{P}}_{L}^{*}F)\Big{)}, (1.11)

thus establishing the non-trivial fact that the limiting joint law is degenerate. This almost sure relation between the marginals is precisely the manifestation of the pathwise structure (I)–(II) above. For a weakly correlated coefficient field \boldsymbol{a} we expect from property (III) that \Xi(\frac{\cdot}{\varepsilon}) displays the CLT scaling \beta=d and that its rescaling converges in law to a Gaussian white noise \Gamma, so that the convergence in (1.11) then leads to the known (or expected) scaling limit results for the different quantities of interest in stochastic homogenization. We however emphasize that the main novelty of the present contribution does not rely in such convergence results in themselves, but rather in the mechanism that leads to them, which is summarized in items (I)–(III) above. ∎

### 1.2. Main results

In order to present our complete theory of fluctuations and address items (I)–(III), we place ourselves in the simplest setting possible and consider discrete elliptic equations with iid conductances, which we think of as the prototype for weakly correlated coefficient fields. Although conceptually simpler on the stochastic side, the discrete setting has some technical inconveniences on the deterministic side, including a discretization error.

Our main results take on the following guise. Precise notation and assumptions are postponed to Section 2, as well as many remarks and corollaries. While items (i) and (ii) below (together with the non-degeneracy in (iv)) imply property (I) in the form (1.6) with the optimal rate o(1)\simeq_{f,g}\varepsilon\mu_{d}(\tfrac{1}{\varepsilon})^{\frac{1}{2}}, items (iii) and (iv) are manifestations of property (III). Regarding the proofs, items (i) and (ii) are established using a Poincaré inequality in the probability space, which holds in the iid setting (cf. Lemma 3.1), item (iii) using a second-order Poincaré inequality due to Chatterjee [C1, C2] (cf. Lemma 4.1), and item (iv) using (an iid version of) the so-called Helffer-Sjöstrand representation formula for covariances [HS-94, Sjostrand-96, NS2] (cf. Lemma 5.1). Apart from these simplifying tools of the discrete iid setting (which can be either extended or avoided in more general contexts, cf. Section 1.4), the proofs only rely on arguments that extend to the continuum setting and to the case of systems. At a technical level, we make strong use of the (pathwise) large-scale weighted Calderón-Zygmund theory for the operator -\nabla\cdot\boldsymbol{a}\nabla (cf. [Armstrong-Daniel-16, GNO-reg]).

###### Theorem 1.

Consider the discrete iid setting, and assume that the law is non-degenerate. Then the following hold for all \varepsilon>0,

1. CLT scaling: for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d},

 \mathbb{E}\left[\Big{|}\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}F:\Xi(% \tfrac{\cdot}{\varepsilon})\Big{|}^{2}\right]^{\frac{1}{2}}\,\lesssim_{F}\,1.
2. Pathwise structure (with optimal error estimates): for all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, letting u_{\varepsilon} and \bar{u} denote the solutions of (the discrete version of) (1.1) and of (1.2),

 \displaystyle\hskip 17.071654pt\mathbb{E}\left[\Big{|}\varepsilon^{-\frac{d}{2% }}\int_{\mathbb{R}^{d}}g\cdot\big{(}\boldsymbol{a}_{\varepsilon}\nabla_{% \varepsilon}u_{\varepsilon}-\bar{\boldsymbol{a}}\nabla_{\varepsilon}u_{% \varepsilon}-\mathbb{E}\left[\boldsymbol{a}_{\varepsilon}\nabla_{\varepsilon}u% _{\varepsilon}-\bar{\boldsymbol{a}}\nabla_{\varepsilon}u_{\varepsilon}\right]% \big{)}-\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\Xi_{i}(\tfrac{% \cdot}{\varepsilon})D_{i}\bar{u}\Big{|}^{2}\right]^{\frac{1}{2}}\\ \displaystyle\lesssim_{f,g}\,\varepsilon\mu_{d}(\tfrac{1}{\varepsilon})^{\frac% {1}{2}}, (1.12)

where we set for all r>0,

 \displaystyle\mu_{d}(r):=\begin{cases}r&:\leavevmode\nobreak\ \leavevmode% \nobreak\ d=1,\\ \log(2+r)&:\leavevmode\nobreak\ \leavevmode\nobreak\ d=2,\\ 1&:\leavevmode\nobreak\ \leavevmode\nobreak\ d>2.\end{cases} (1.13)
3. Asymptotic normality (with nearly optimal rate): for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d},

 \displaystyle\delta_{\mathcal{N}}\bigg{(}\varepsilon^{-\frac{d}{2}}\int_{% \mathbb{R}^{d}}F:\Xi(\tfrac{\cdot}{\varepsilon})\bigg{)}\,\lesssim_{F}\,% \varepsilon^{\frac{d}{2}}\log(2+\tfrac{1}{\varepsilon}),

where for a random variable X\in\operatorname{L}^{2}(\Omega) its distance to normality is defined by

 \displaystyle\delta_{\mathcal{N}}(X):=\operatorname{d}_{\operatorname{W}}\left% ({\frac{X}{\mathrm{Var}\left[X\right]^{\frac{1}{2}}}},{\mathcal{N}}\right)+% \operatorname{d}_{\operatorname{K}}\left({\frac{X}{\mathrm{Var}\left[X\right]^% {\frac{1}{2}}}},{\mathcal{N}}\right), (1.14)

with \mathcal{N} a standard Gaussian random variable and with \operatorname{d}_{\operatorname{W}}\left({\cdot},{\cdot}\right) and \operatorname{d}_{\operatorname{K}}\left({\cdot},{\cdot}\right) the Wasserstein and Kolmogorov metrics.

4. Convergence of the covariance structure (with optimal rate): there exists a non-degenerate symmetric 4-tensor \mathcal{Q} such that for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d},

 \Big{|}\mathrm{Var}\left[\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}F:\Xi(% \tfrac{\cdot}{\varepsilon})\right]-\int_{\mathbb{R}^{d}}F:\mathcal{Q}\,F\Big{|% }\,\lesssim_{F}\,\varepsilon\mu_{d}(\tfrac{1}{\varepsilon})^{\frac{1}{2}}.

In particular, combined with item (iii), this yields the convergence in law of \varepsilon^{-\frac{d}{2}}\Xi(\tfrac{\cdot}{\varepsilon}) to a 2-tensor Gaussian white noise \Gamma with covariance structure \mathcal{Q}, and the (discrete version of the) joint convergence result (1.11) follows.

In the estimates above, and in the rest of this paper, \lesssim_{\gamma} (with possibly a subscript \gamma) stands for \leq up to a multiplicative constant C=C(d,\lambda,\gamma) that only depends on d, \lambda, and \gamma (through a suitable norm of \gamma, should \gamma be a function). ∎

This fluctuation theory is complemented by the following characterization of the fluctuation tensor \mathcal{Q} by periodization in law. This characterization comes in form of a representative volume element (RVE) method, of which we give the optimal error estimate. In particular, comparing with the results for the RVE approximation \bar{\boldsymbol{a}}_{L,N} of the homogenized coefficients \bar{\boldsymbol{a}} (cf. [GO1, GO2, GNO1]), and choosing N\simeq L^{d} below, we may conclude that an RVE approximation for \mathcal{Q} with accuracy O(L^{-\frac{d}{2}}) (up to logarithmic corrections) is extracted at the same cost as an RVE approximation for \bar{\boldsymbol{a}} with accuracy O(L^{-d}). Precise assumptions and notation are again postponed to Section 2.

###### Theorem 2.

Consider the discrete iid setting. Define

 \displaystyle\bar{\boldsymbol{a}}_{L}e_{i}:=\fint_{Q_{L}}\boldsymbol{a}_{L}(% \nabla\phi_{L,i}+e_{i}), (1.15)

in terms of the L-periodized coefficient field \boldsymbol{a}_{L} and corrector \phi_{L}. Then the fluctuation tensor \mathcal{Q} defined in Theorem 1(iv) satisfies

 \displaystyle\mathcal{Q}=\lim_{L\uparrow\infty}\mathrm{Var}\left[L^{\frac{d}{2% }}\bar{\boldsymbol{a}}_{L}^{*}\right]. (1.16)

In addition, considering iid realizations (\boldsymbol{a}_{L}^{(n)})_{n=1}^{N} of \boldsymbol{a}_{L} and setting \bar{\boldsymbol{a}}_{L}^{(n)}:=\bar{\boldsymbol{a}}_{L}(\boldsymbol{a}_{L}^{(% n)}), we define the RVE approximation as the square of the sample standard deviation

 \displaystyle\mathcal{Q}_{L,N}:=\frac{L^{d}}{N-1}\sum_{n=1}^{N}\big{(}\bar{% \boldsymbol{a}}_{L}^{(n)}-\bar{\boldsymbol{a}}_{L,N}\big{)}^{*}\otimes\big{(}% \bar{\boldsymbol{a}}_{L}^{(n)}-\bar{\boldsymbol{a}}_{L,N}\big{)}^{*},\qquad% \bar{\boldsymbol{a}}_{L,N}:=\frac{1}{N}\sum_{n=1}^{N}\bar{\boldsymbol{a}}_{L}^% {(n)}, (1.17)

and for all L,N\geq 2 there holds

### 1.3. Relation to previous works

Uniform moment bounds on \varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla_{\varepsilon}u_{\varepsilon} were first obtained by Conlon and Naddaf [CN] and by the second author [G] in weaker forms, and were established in their optimal form in any dimension d\geq 2 by Marahrens and the third author [MaO] for discrete elliptic equations with iid conductances; see also [GMa] for continuum scalar equations. The convergence in law of \varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla_{\varepsilon}u_{\varepsilon} and \varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}F:\nabla\phi(\frac{\cdot}{% \varepsilon}) was established in dimensions d>2 for a discrete elliptic equation with iid Gaussian conductances by Mourrat and the last author [MO], Mourrat and Nolen [MN], and Gu and Mourrat [GuM], where in addition the pathwise structure in the form (1.11) was partially discovered. These works were based on the Helffer-Sjöstrand representation formula [HS-94, Sjostrand-96, NS2] and on tools introduced by the last two authors [GO1, GO2] (inspired by the unpublished work by Naddaf and Spencer [NS]), the last two authors and Neukamm [GNO1, GNO2], and Marahrens and the last author [MaO]. We also mention related work on the fluctuations of \bar{\boldsymbol{a}}_{L} (cf. (1.15)) by Nolen [N, Nolen-16] (who was first to use Chatterjee’s second-order Poincaré inequalities [C1, C2] in stochastic homogenization), the second author and Nolen [GN], Rossignol [R], and Biskup, Salvi, and Wolff [BSW]. Note that the present contribution originates in the attempt to upgrade the results in [GN] into a functional CLT for an energy density.

A variational quantity related to the homogenization commutator can also be traced back to the work [AS] by Armstrong and Smart. The stationary version of this quantity, which is key to study fluctuations and essentially coincides with the standard homogenization commutator \Xi that we define here, was independently introduced by Armstrong, Kuusi, and Mourrat in [AKM2, Definition 5.3 and Paragraph 8.1] around the same time as the first version of the present work. As here it contains the relation (1.9) to the corrector, and as in [GO4] it contains the convergence to white noise under the assumption of finite range of dependence. A partial version of the pathwise structure was conjectured by Armstrong, Gu, and Mourrat in [GuM, MourratFun] based on a string of heuristic arguments within the variational and renormalization perspective of [AS, AKM2], but the crucial validity of the two-scale expansion of the homogenization commutator (cf. property (I)) does not appear in [GuM, MourratFun, AKM2] nor in any subsequent work.

Summing up, the pathwise structure of fluctuations revealed in this contribution and the underlying mechanism described in (I)–(III) are made precise and rigorous here for the first time — in any setting.

### 1.4. Extensions to the continuum setting and robustness of results

In the case of (non-symmetric) continuum systems, we may extend our fluctuation theory in two different directions, which are postponed to the companion articles [DGO2, DFGO1, DFGO2].

A first extension concerns the case of a coefficient field with finite range of dependence, for which no functional inequality is satisfied in general. In [GO4], the last two authors showed the convergence in law of \varepsilon^{-\frac{d}{2}}\Xi(\tfrac{\cdot}{\varepsilon}) to a Gaussian white noise (albeit without convergence rate) by combining a semi-group approach with the approximate locality of \Xi with respect to the coefficient field and with the assumption of finite range of dependence of the coefficients. A similar result is achieved in [AKM2]. The proof of the validity of the two-scale expansion (1.12) of the homogenization commutator in the CLT scaling is more involved and will be presented in [DFGO2] based on the semi-group approach of [GO4].

A second extension concerns the case of a coefficient field with strong correlations. In [DGO2, DFGO1], we consider the model framework of a coefficient field given by (the image by a local Lipschitz function of) a Gaussian field that has an algebraically decaying (but not necessarily integrable) covariance function c, say at some fixed (yet arbitrary) rate c(x)\simeq(1+|x|)^{-\beta} parametrized by \beta>0. For such coefficient fields, we establish in [DGO2] the validity of the two-scale expansion (1.12) of the homogenization commutator in the suitable fluctuation scaling. The proof relies on a weighted version of a Poincaré inequality in the probability space (cf. [DG1, DG2]), (pathwise) large-scale Calderón-Zygmund theory for -\nabla\cdot\boldsymbol{a}\nabla (cf. [GNO-reg]), and moment bounds on the corrector (cf. [GNO-quant]). This illustrates the surprising robustness of the pathwise structure with respect to the large-scale behavior of the homogenization commutator. Indeed, in dimension d=1 (in which case the quantities under investigation are simpler and explicit), two typical behaviors have been identified in terms of scaling limit of the homogenization commutator \Xi, depending on the parameter \beta (cf. [BGMP-08]):

• For \beta>d=1: The commutator \Xi displays the CLT scaling and \varepsilon^{-\frac{d}{2}}\Xi(\frac{\cdot}{\varepsilon}) converges to a Gaussian white noise (Gaussian fluctuations, local limiting covariance structure).

• For 0<\beta<d=1: The suitable rescaling \varepsilon^{-\frac{\beta}{2}}\Xi(\frac{\cdot}{\varepsilon}) converges up to a subsequence to a fractional Gaussian field (Gaussian fluctuations, nonlocal limiting covariance structure, potentially no uniqueness of the limit).

(Note that a different, non-Gaussian behavior may also occur, cf. [Gu-Bal-12, LNZH-17].) In particular, the pathwise result holds in these two examples whereas the rescaled homogenization commutator \Xi does not necessarily converge to white noise or may even not converge at all. The identification of the scaling limit of the homogenization commutator in higher dimensions is addressed in [DFGO1] for the whole range of values of \beta>0, where we investigate the consequences of the locality of \Xi with respect to the coefficient field, combining techniques developed in [GNO-reg, GNO-quant] with Malliavin calculus versions of the Helffer-Sjöstrand representation formula and of a second-order Poincaré inequality. This work extends [BGMP-08] to dimensions d\geq 2.

## 2. Main results

In this section, we introduce notation and assumptions, we state precise versions of the main results (and in particular make explicit the norms of the test functions in the estimates), and we discuss various corollaries.

### 2.1. Notation and assumptions

We start by introducing the discrete iid framework in which our main results are established. We consider a random conductance problem on the integer lattice \mathbb{Z}^{d}, and denote by \{e_{i}\}_{i=1}^{d} the canonical basis of \mathbb{R}^{d}. We regard \mathbb{Z}^{d} as a graph with (unoriented) edge set \mathcal{B}=\{(x,z)\in\mathbb{Z}^{d}\times\mathbb{Z}^{d}:|x-z|=1\}. For edges (x,z)\in\mathcal{B}, we also write x\sim z. We define the set of conductances \{a(b)\}_{b\in\mathcal{B}} by \Omega=[\lambda,1]^{\mathcal{B}} for some fixed 0<\lambda\leq 1. We endow \Omega with the \sigma-algebra generated by cylinder sets and with a probability measure \mathbb{P}. We denote by \mathbb{E}\left[\cdot\right], \mathrm{Var}\left[\cdot\right], and \operatorname{Cov}\left[{\cdot};{\cdot}\right] the associated expectation, variance, and covariance. A realization a\in\Omega is by definition a countable collection \{a(b)\}_{b\in\mathcal{B}} of conductances. A random field u:\mathbb{R}^{d}\times\Omega\to\mathbb{R} is said to be stationary if it is shift-covariant, in the sense of u(x,a(\cdot-z))=u(x-z,a) for all x,z\in\mathbb{R}^{d} and a\in\Omega. In this contribution, we focus on the case when the probability measure \mathbb{P} is a product measure, that is, when the conductances \{a(b)\}_{b\in\mathcal{B}} are iid random variables, and we shall make use of available functional inequalities in this product probability space.

Let \nabla denote the forward discrete gradient (u:\mathbb{Z}^{d}\to\mathbb{R})\mapsto(\nabla u:\mathbb{Z}^{d}\to\mathbb{R}^{d}) defined componentwise by \nabla_{i}u(x)=u(x+e_{i})-u(x) for 1\leq i\leq d, and let \nabla^{*} denote the backward discrete gradient (u:\mathbb{Z}^{d}\to\mathbb{R})\mapsto(\nabla^{*}u:\mathbb{Z}^{d}\to\mathbb{R}% ^{d}) defined componentwise by \nabla_{i}^{*}u(x)=u(x)-u(x-e_{i}) for all 1\leq i\leq d. The operator -\nabla^{*}\cdot is thus the adjoint of \nabla on \ell^{2}(\mathbb{Z}^{d}), and we consider the elliptic operator -\nabla^{*}\cdot\boldsymbol{a}\nabla with coefficients

 \boldsymbol{a}:x\mapsto\boldsymbol{a}(x):=\mathrm{diag}\left[a(x,x+e_{1}),% \dots,a(x,x+e_{d})\right],

acting on functions u:\mathbb{Z}^{d}\to\mathbb{R} as

 -\nabla^{*}\cdot\boldsymbol{a}\nabla u(x)\,:=\,\sum_{z:z\sim x}a(x,z)(u(x)-u(z% )).

In order to state the standard qualitative homogenization result [Ku, Kozlov-85] for the corresponding discrete elliptic equation, we consider for all \varepsilon>0 the rescaled operator -\nabla_{\varepsilon}^{*}\cdot\boldsymbol{a}_{\varepsilon}\nabla_{\varepsilon}, where \boldsymbol{a}_{\varepsilon}(\cdot):=\boldsymbol{a}(\frac{\cdot}{\varepsilon}), and where \nabla_{\varepsilon} and \nabla^{*}_{\varepsilon} act on functions u_{\varepsilon}:\mathbb{Z}_{\varepsilon}^{d}:=\varepsilon\mathbb{Z}^{d}\to% \mathbb{R}, and are defined componentwise by \nabla_{\varepsilon,i}u_{\varepsilon}(x)=\varepsilon^{-1}(u_{\varepsilon}(x+% \varepsilon e_{i})-u_{\varepsilon}(x)) and \nabla_{\varepsilon,i}^{*}u_{\varepsilon}(x)=\varepsilon^{-1}(u_{\varepsilon}(% x)-u_{\varepsilon}(x-\varepsilon e_{i})) for all i. We shall also let \nabla_{\varepsilon} and \nabla^{*}_{\varepsilon} act on continuous functions u:\mathbb{R}^{d}\to\mathbb{R}, so that \nabla_{\varepsilon}u and \nabla_{\varepsilon}^{*}u are continuous functions as well. If u\in C^{1}(\mathbb{R}^{d}), then \nabla_{\varepsilon}u(x) and \nabla_{\varepsilon}^{*}u(x) converge to the continuum gradient Du(x) for all x\in\mathbb{R}^{d} as \varepsilon\downarrow 0. In what follows, for all m\geq 1, we systematically extend maps v:\mathbb{Z}^{d}\to\mathbb{R}^{m} to piecewise constant maps \mathbb{R}^{d}\to\mathbb{R}^{m} (still denoted by v) by setting v|_{Q(x)}:=v(x) for all x\in\mathbb{Z}^{d} (where Q(x):=x+[-\frac{1}{2},\frac{1}{2})^{d} is the unit cube centered at x), and we use this notation e.g. for \boldsymbol{a}:\mathbb{Z}^{d}\to\mathbb{R}^{d\times d} (but also for \phi:\mathbb{Z}^{d}\to\mathbb{R}^{d} and \Xi:\mathbb{Z}^{d}\to\mathbb{R}^{d\times d} defined below). This systematic extension of functions from the lattice \mathbb{Z}^{d} to \mathbb{R}^{d} allows to state all discrete results in a form that would hold mutatis mutandis in the continuum setting. In addition, although in the discrete setting it is more natural to consider a symmetric coefficient field \boldsymbol{a}, we use non-symmetric notation in the statement of the results in view of the non-symmetric continuum setting, and we denote by \boldsymbol{a}^{*} the pointwise transpose field associated with \boldsymbol{a}.

Qualitative stochastic homogenization [Ku, Kozlov-85] ensures that, for all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, almost surely, the unique Lax-Milgram solutions u_{\varepsilon} and v_{\varepsilon} in \mathbb{R}^{d} of111 These equations are understood as follows: for all x\in Q the function u_{\varepsilon}(\varepsilon x+\cdot) on \mathbb{Z}_{\varepsilon}^{d} is the solution of the discrete elliptic equation with coefficient \boldsymbol{a}_{\varepsilon} and with right-hand side \nabla_{\varepsilon}^{*}\cdot f(\varepsilon x+\cdot). This definition allows to state results in a form that holds in the continuum setting, and in terms of norms of the right-hand side that needn’t embed into the space of continuous functions.

 -\nabla_{\varepsilon}^{*}\cdot\boldsymbol{a}_{\varepsilon}\nabla_{\varepsilon}% u_{\varepsilon}\,=\,\nabla_{\varepsilon}^{*}\cdot f,\qquad-\nabla_{\varepsilon% }^{*}\cdot\boldsymbol{a}_{\varepsilon}^{*}\nabla_{\varepsilon}v_{\varepsilon}% \,=\,\nabla_{\varepsilon}^{*}\cdot f, (2.1)

converge weakly as \varepsilon\downarrow 0 to the unique Lax-Milgram solutions \bar{u} and \bar{v} in \mathbb{R}^{d} of the (continuum) elliptic equations

respectively, where \bar{\boldsymbol{a}} is the homogenized matrix characterized by

 \bar{\boldsymbol{a}}e_{i}\,=\,\mathbb{E}\left[\boldsymbol{a}(\nabla\phi_{i}+e_% {i})\right], (2.3)

for all 1\leq i\leq d, where \phi_{i} is the so-called corrector in direction e_{i}. It is defined, for almost every realization \boldsymbol{a}, as the unique solution in \mathbb{Z}^{d} of

 -\nabla^{*}\cdot\boldsymbol{a}(\nabla\phi_{i}+e_{i})\,=\,0, (2.4)

with \nabla\phi_{i} stationary and having finite second moment, and with \phi_{i}(0)=0. We then set \phi:=(\phi_{i})_{i=1}^{d}. Note that \overline{(\boldsymbol{a}^{*})}=(\bar{\boldsymbol{a}})^{*}. For symmetric coefficient fields, \boldsymbol{a}^{*}=\boldsymbol{a} and \bar{\boldsymbol{a}}^{*}=\bar{\boldsymbol{a}}.

We consider the fluctuations of the field \nabla u_{\varepsilon} and of the flux \boldsymbol{a}_{\varepsilon}\nabla u_{\varepsilon}, as encoded in the random linear functionals I_{1}^{\varepsilon}:(f,g)\mapsto I_{1}^{\varepsilon}(f,g) and I_{2}^{\varepsilon}:(f,g)\mapsto I_{2}^{\varepsilon}(f,g) defined for all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d} by

 \displaystyle I_{1}^{\varepsilon}(f,g) \displaystyle:= \displaystyle\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla_{% \varepsilon}(u_{\varepsilon}-\mathbb{E}\left[u_{\varepsilon}\right]), \displaystyle I_{2}^{\varepsilon}(f,g) \displaystyle:= \displaystyle\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\big{(}% \boldsymbol{a}_{\varepsilon}\nabla_{\varepsilon}u_{\varepsilon}-\mathbb{E}% \left[\boldsymbol{a}_{\varepsilon}\nabla_{\varepsilon}u_{\varepsilon}\right]% \big{)}.

We further encode the fluctuations of the corrector field \nabla\phi and flux \boldsymbol{a}(\nabla\phi+\operatorname{Id}) in the random linear functionals J_{1}^{\varepsilon}:F\mapsto J_{1}^{\varepsilon}(F) and J_{2}^{\varepsilon}:F\mapsto J_{2}^{\varepsilon}(F) defined for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} by

 \displaystyle J_{1}^{\varepsilon}(F) \displaystyle:= \displaystyle\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}F(x):\nabla\phi(% \tfrac{x}{\varepsilon})\,dx, \displaystyle J_{2}^{\varepsilon}(F) \displaystyle:= \displaystyle\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}F(x):\big{(}% \boldsymbol{a}_{\varepsilon}(x)(\nabla\phi(\tfrac{x}{\varepsilon})+% \operatorname{Id})-\bar{\boldsymbol{a}}\big{)}\,dx.

As explained above, a crucial role is played by the (standard) homogenization commutator, which in the present discrete setting takes the form \Xi:=(\Xi_{i})_{i=1}^{d} with

and by the error in the two-scale expansion of the homogenization commutator of the solution. These quantities are encoded in the random linear functionals J_{0}^{\varepsilon}:F\mapsto J_{0}^{\varepsilon}(F) and E^{\varepsilon}:(f,g)\mapsto E^{\varepsilon}(f,g) defined for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} and all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d} by

 \displaystyle J_{0}^{\varepsilon}(F) \displaystyle:= \displaystyle\!\!\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}F(x):\Xi(% \tfrac{x}{\varepsilon})\,dx, \displaystyle E^{\varepsilon}(f,g) \displaystyle:= \displaystyle\!\!\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\big{(}% \boldsymbol{a}_{\varepsilon}\nabla_{\varepsilon}u_{\varepsilon}-\bar{% \boldsymbol{a}}\nabla_{\varepsilon}u_{\varepsilon}-\mathbb{E}\left[\boldsymbol% {a}_{\varepsilon}\nabla_{\varepsilon}u_{\varepsilon}-\bar{\boldsymbol{a}}% \nabla_{\varepsilon}u_{\varepsilon}\right]\big{)}-\varepsilon^{-\frac{d}{2}}% \int_{\mathbb{R}^{d}}g\cdot\Xi_{i}(\tfrac{\cdot}{\varepsilon})D_{i}\bar{u}.

Since the case d=1 is much simpler and well-understood [Gu-16], we shall only focus in the sequel on dimensions d\geq 2.

We first recall the following uniform boundedness result for J_{0}^{\varepsilon}, establishing the CLT scaling for the fluctuations of the homogenization commutator (cf. Theorem 1(i)). Although essentially contained in the main result of the first contribution [GO1] of the second and third authors to the field, a short proof with up-to-date tools is included for completeness in Section 3. Note that the norm of the test function is substantially weaker than \operatorname{L}^{1}(\mathbb{R}^{d}) in terms of integrability and is thus compatible with the behavior of Helmholtz projections of smooth and compactly supported functions, which is necessary for the pathwise result of Corollary 2.4 below.

###### Proposition 2.1.

Let d\geq 2, let \mathbb{P} be a product measure, and set w_{1}(z):=1+|z|. For all \varepsilon>0 and all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} we have for all 0<p-1\ll 1 and all \alpha>d\frac{p-1}{4p},

 \mathbb{E}\left[|J_{0}^{\varepsilon}(F)|^{2}\right]^{\frac{1}{2}}+\mathbb{E}% \left[|J_{1}^{\varepsilon}(F)|^{2}\right]^{\frac{1}{2}}+\mathbb{E}\left[|J_{2}% ^{\varepsilon}(F)|^{2}\right]^{\frac{1}{2}}\,\lesssim_{\alpha,p}\,\|w_{1}^{2% \alpha}F\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})}.

Above, and in the rest of this article, \ll stands for \leq up to a small enough multiplicative constant C=C(d)>0 that only depends on d. ∎

### 2.2. Pathwise structure

Our first main result establishes the smallness of the rescaled error E^{\varepsilon} in the two-scale expansion of the homogenization commutator (cf. Theorem 1(ii)), which is the key to the pathwise structure (1.11). As for Proposition 2.1, the proof relies on the Poincaré inequality in the probability space that is satisfied for iid coefficients. From a technical point of view, we exploit the large-scale Calderón-Zygmund theory for the operator -\nabla^{*}\cdot\boldsymbol{a}\nabla in the form developed in [GNO-reg].

###### Proposition 2.2.

Let d\geq 2, let \mathbb{P} be a product measure, let \mu_{d} be defined in (1.13), and set w_{1}(z):=1+|z|. For all \varepsilon>0 and all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d} we have for all 0<p-1\ll 1 and all \alpha>d\frac{p-1}{4p},

 \mathbb{E}\left[|E^{\varepsilon}(f,g)|^{2}\right]^{\frac{1}{2}}\,\lesssim_{% \alpha,p}\,\varepsilon\mu_{d}(\tfrac{1}{\varepsilon})^{\frac{1}{2}}\Big{(}\|f% \|_{\operatorname{L}^{4}(\mathbb{R}^{d})}\|w_{1}^{\alpha}Dg\|_{\operatorname{L% }^{4p}(\mathbb{R}^{d})}+\|g\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}\|w_{1}^{% \alpha}Df\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}\Big{)}.\qed
###### Remark 2.3.

For simplicity the estimates in Propositions 2.1 and 2.2 above are stated and proved for second moments only, but the same arguments yield similar estimates for all algebraic (and even stretched exponential) moments (cf. [DGO2, DFGO2]). ∎

In view of identity (1.7) (which indeed holds in the discrete setting up to a higher-order discretization error), the above result implies that the large-scale fluctuations of I_{1}^{\varepsilon} and I_{2}^{\varepsilon} are driven by those of J_{0}^{\varepsilon} in a pathwise sense. Identity (1.9) (which again holds up to a discretization error) yields a similar pathwise result for J_{1}^{\varepsilon} and J^{\varepsilon}_{2}.

###### Corollary 2.4.

Let d\geq 2, let \mathbb{P} be a product measure, let \bar{\mathcal{P}}_{H}, \bar{\mathcal{P}}_{H}^{*}, and \bar{\mathcal{P}}_{L}^{*} be defined in (1.8), let \mu_{d} be defined in (1.13), and set w_{1}(z):=1+|z|. For all \varepsilon>0, all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, and all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d}, we have for all 0<p-1\ll 1 and all \alpha>d\frac{p-1}{4p},

 \displaystyle\mathbb{E}\left[|I_{1}^{\varepsilon}(f,g)-J_{0}^{\varepsilon}(% \bar{\mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g)|^{2}\right]^{\frac{1}% {2}}+\mathbb{E}\left[|I_{2}^{\varepsilon}(f,g)+J_{0}^{\varepsilon}(\bar{% \mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{L}^{*}g)|^{2}\right]^{\frac{1}{2}}% \\ \displaystyle\lesssim_{\alpha,p}\varepsilon\mu_{d}(\tfrac{1}{\varepsilon})^{% \frac{1}{2}}\,\Big{(}\|f\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}\|w_{1}^{% \alpha}Dg\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}+\|g\|_{\operatorname{L}^{4% }(\mathbb{R}^{d})}\|w_{1}^{\alpha}Df\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}% \Big{)}, (2.6)

and

 \displaystyle\mathbb{E}\left[|J_{1}^{\varepsilon}(F)+J_{0}^{\varepsilon}(\bar{% \mathcal{P}}_{H}^{*}F)|^{2}\right]^{\frac{1}{2}}+\mathbb{E}\left[|J_{2}^{% \varepsilon}(F)-J_{0}^{\varepsilon}(\bar{\mathcal{P}}_{L}^{*}F)|^{2}\right]^{% \frac{1}{2}}\lesssim_{\alpha,p}\,\varepsilon\|w_{1}^{2\alpha}DF\|_{% \operatorname{L}^{2p}(\mathbb{R}^{d})}, (2.7)

where by definition we have \bar{\mathcal{P}}_{H}f=-D\bar{u} and \bar{\mathcal{P}}_{H}^{*}g=-D\bar{v}. In particular, we give meaning to J_{0}^{\varepsilon}(\bar{\mathcal{P}}_{H}^{*}F) and J_{0}^{\varepsilon}(\bar{\mathcal{P}}_{L}^{*}F) in \operatorname{L}^{2}(\Omega) for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d}, even when \bar{\mathcal{P}}_{H}^{*}F and \bar{\mathcal{P}}_{L}^{*}F do not have integrable decay. ∎

###### Remark 2.5.

For all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, we may also consider the unique Lax-Milgram solutions u_{\varepsilon}^{\circ} and v_{\varepsilon}^{\circ} in \mathbb{R}^{d} of

which, almost surely, converge weakly as \varepsilon\downarrow 0 to the unique Lax-Milgram solutions \bar{u}^{\circ} and \bar{v}^{\circ} in \mathbb{R}^{d} of

respectively. Similar considerations as in the proof of Proposition 2.2 and Corollary 2.4 then lead to a pathwise result for the fluctuations of the random linear functionals I_{3}^{\varepsilon}:(f,g)\mapsto I_{3}^{\varepsilon}(f,g) and I_{4}^{\varepsilon}:(f,g)\mapsto I_{4}^{\varepsilon}(f,g) defined for all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d} by

 \displaystyle I_{3}^{\varepsilon}(f,g) \displaystyle:= \displaystyle\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla_{% \varepsilon}(u_{\varepsilon}^{\circ}-\mathbb{E}\left[u_{\varepsilon}^{\circ}% \right]), \displaystyle I_{4}^{\varepsilon}(f,g) \displaystyle:= \displaystyle\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\big{(}% \boldsymbol{a}_{\varepsilon}(\nabla_{\varepsilon}u_{\varepsilon}^{\circ}+f)-% \mathbb{E}\left[\boldsymbol{a}_{\varepsilon}(\nabla_{\varepsilon}u_{% \varepsilon}^{\circ}+f)\right]\big{)},

that takes the form

 \displaystyle\mathbb{E}\left[\big{|}I_{3}^{\varepsilon}(f,g)+J_{0}^{% \varepsilon}\big{(}\bar{\mathcal{P}}_{L}f\otimes\bar{\mathcal{P}}_{H}^{*}g\big% {)}\big{|}^{2}\right]^{\frac{1}{2}}+\mathbb{E}\left[\big{|}I_{4}^{\varepsilon}% (f,g)-J_{0}^{\varepsilon}\big{(}\bar{\mathcal{P}}_{L}f\otimes\bar{\mathcal{P}}% _{L}^{*}g\big{)}\big{|}^{2}\right]^{\frac{1}{2}}\lesssim_{f,g}\varepsilon\mu_{% d}(\tfrac{1}{\varepsilon})^{\frac{1}{2}},

where by definition \bar{\mathcal{P}}_{H}f=-D\bar{u}, \bar{\mathcal{P}}_{L}f=D\bar{u}^{\circ}+f, \bar{\mathcal{P}}_{H}^{*}g=-D\bar{v}, and \bar{\mathcal{P}}_{L}^{*}g=D\bar{v}^{\circ}+g. ∎

Incidentally, as a consequence of our analysis, combining the two-scale expansion of the homogenization commutator (1.5) with identity (1.7), we obtain a new (nonlocal) two-scale expansion for the solution \nabla_{\varepsilon}u_{\varepsilon} that is not only accurate at order 1 for the strong \operatorname{L}^{2}(\mathbb{R}^{d}) topology but also at the order of the fluctuation scaling for the weak \operatorname{L}^{2}(\mathbb{R}^{d}) topology — in contrast to the usual two-scale expansion (1.3) (cf. [GuM]). (The second estimate below is a reformulation of Proposition 2.2, whereas the first estimate is a corollary of [GNO-quant, Theorem 3].)

###### Corollary 2.6.

Let d\geq 2, let \mathbb{P} be a product measure, and let \mu_{d} be defined in (1.13). For all \varepsilon>0 and all f\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, we set

 r_{\varepsilon}(f):=\nabla_{\varepsilon}u_{\varepsilon}-\Big{(}\underbrace{% \mathbb{E}\left[\nabla_{\varepsilon}u_{\varepsilon}\right]+\nabla_{\varepsilon% }(-\nabla_{\varepsilon}^{*}\cdot\bar{\boldsymbol{a}}\nabla_{\varepsilon})^{-1}% \nabla_{\varepsilon}^{*}\cdot\big{(}\Xi_{i}(\tfrac{\cdot}{\varepsilon})D_{i}% \bar{u}\big{)}}_{\displaystyle\text{nonlocal two-scale expansion of $\nabla_{% \varepsilon}u_{\varepsilon}$}}\Big{)}.

Then this (nonlocal) two-scale expansion correctly captures

• the spatial oscillations of \nabla_{\varepsilon}u_{\varepsilon} in a strong norm: for all f\in C^{\infty}_{c}(\mathbb{R}^{d})^{d},

 \mathbb{E}\left[\|r_{\varepsilon}(f)\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}^% {2}\right]^{\frac{1}{2}}\,\lesssim_{f}\,\varepsilon\mu_{d}(\tfrac{1}{% \varepsilon})^{\frac{1}{2}};
• the random fluctuations of \nabla_{\varepsilon}u_{\varepsilon} in the CLT scaling: for all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d},

 \mathbb{E}\left[\Big{|}\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot r% _{\varepsilon}(f)\Big{|}^{2}\right]^{\frac{1}{2}}\,\lesssim_{f,g}\,\varepsilon% \mu_{d}(\tfrac{1}{\varepsilon})^{\frac{1}{2}}.\qed

### 2.3. Approximate normality

We turn to the normal approximation result for the homogenization commutator (cf. Theorem 1(iii)), which states that the fluctuations of \varepsilon^{-\frac{d}{2}}\Xi(\frac{\cdot}{\varepsilon}) are asymptotically Gaussian (up to a non-degeneracy condition that is elucidated in Proposition 2.9 below). The approach is inspired by previous works by Nolen [N, Nolen-16], based on a second-order Poincaré inequality à la Chatterjee [C1, LRP-15], which is key to optimal convergence rates. Since such functional inequalities are not easily amenable to the use of large-scale Calderón-Zygmund theory for the operator -\nabla^{*}\cdot\boldsymbol{a}\nabla, we rather have to exploit optimal annealed estimates on mixed gradients of the Green’s function [MaO], which leads to an additional \log(2+\frac{1}{\varepsilon}) factor in the rate below (we do not know whether this is optimal). The proof exploits the approximate locality of the homogenization commutator \Xi.

###### Proposition 2.7.

Let d\geq 2, let \mathbb{P} be a product measure, let \mu_{d} and \delta_{\mathcal{N}} be defined in (1.13) and (1.14), and set w_{1}(z):=1+|z|. For all \varepsilon>0 and all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} with \mathrm{Var}\left[J_{0}^{\varepsilon}(F)\right]>0, we have for all \alpha>0,

 \displaystyle\delta_{\mathcal{N}}(J_{0}^{\varepsilon}(F))\,\lesssim_{\alpha}\,% \varepsilon^{\frac{d}{2}}\,\frac{\|F\|_{\operatorname{L}^{3}(\mathbb{R}^{d})}^% {3}+\|w_{1}^{\alpha}DF\|_{\operatorname{L}^{3}(\mathbb{R}^{d})}^{3}}{\mathrm{% Var}\left[J_{0}^{\varepsilon}(F)\right]^{\frac{3}{2}}} \displaystyle                                           +\varepsilon^{\frac{d}% {2}}\log(2+\tfrac{1}{\varepsilon})\,\frac{\|w_{1}^{\alpha}F\|_{\operatorname{L% }^{4}(\mathbb{R}^{d})}^{2}+\|w_{1}^{\alpha}DF\|_{\operatorname{L}^{4}(\mathbb{% R}^{d})}^{2}}{\mathrm{Var}\left[J_{0}^{\varepsilon}(F)\right]}.\qed
###### Remark 2.8.

In the case of iid conductances that are (smooth local transformations of) Gaussian random variables, a nicer version of a second-order Poincaré inequality is available (cf. [C2, Theorem 2.2]), which in addition controls the total variation distance. It allows to use large-scale Calderón-Zygmund theory and so to avoid Green’s functions, which leads to the optimal rate \varepsilon^{\frac{d}{2}} in the above estimate (i.e. without the spurious logarithmic factor). This argument is given in the continuum setting in the forthcoming work [DFGO1]. ∎

### 2.4. Covariance structure

Since J_{0}^{\varepsilon} is asymptotically Gaussian, it remains to identify the limit of its covariance structure (cf. Theorem 1(iv)). The following shows that the limiting covariance is that of a (tensorial) white noise with some non-degenerate covariance tensor \mathcal{Q}. The convergence rate in (2.8) below is new in any dimension and is expected to be optimal. The proof crucially relies on the approximate locality of the homogenization commutator and on (an iid version of) the Helffer-Sjöstrand representation formula for the variance [HS-94, Sjostrand-96, NS2], which is a stronger tool than the Poincaré inequality. As for the pathwise result, the proof exploits the large-scale Calderón-Zygmund theory for the operator -\nabla^{*}\cdot\boldsymbol{a}\nabla.

###### Proposition 2.9.

Let d\geq 2, let \mathbb{P} be a product measure, let \mu_{d} be defined in (1.13), and set w_{1}(z):=1+|z|.

1. There exists a symmetric222Since \mathcal{Q} is a (limiting) covariance, it is of course symmetric in the sense of \mathcal{Q}_{ijkl}=\mathcal{Q}_{klij}. If the coefficients \boldsymbol{a} are symmetric, then it has the additional symmetry \mathcal{Q}_{ijkl}=\mathcal{Q}_{jikl} (hence also \mathcal{Q}_{ijkl}=\mathcal{Q}_{ijlk}). 4-tensor \mathcal{Q} such that for all \varepsilon>0 and all F,G\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} we have for all 0<p-1\ll 1 and all \alpha>d\frac{p-1}{4p},

 \displaystyle\hskip 28.452756pt\bigg{|}\operatorname{Cov}\left[{J_{0}^{% \varepsilon}(F)};{J_{0}^{\varepsilon}(G)}\right]-\int_{\mathbb{R}^{d}}F(x):% \mathcal{Q}\,G(x)dx\bigg{|}\,\lesssim_{\alpha,p}\,\varepsilon\mu_{d}(\tfrac{1}% {\varepsilon})^{\frac{1}{2}}\\ \displaystyle\times\big{(}\|F\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}+\|w_{1}% ^{2\alpha}DF\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})}\big{)}\big{(}\|G\|_{% \operatorname{L}^{2}(\mathbb{R}^{d})}+\|w_{1}^{2\alpha}DG\|_{\operatorname{L}^% {2p}(\mathbb{R}^{d})}\big{)}. (2.8)

Moreover, for all 1\leq i,j,k,l\leq d and all \delta>0, we have for all L\geq 1,

 \displaystyle\bigg{|}\mathcal{Q}_{ijkl}-\int_{Q_{2L}}\frac{|Q_{L}\cap(x+Q_{L})% |}{|Q_{L}|}\operatorname{Cov}\left[{\Xi_{ij}(x)};{\Xi_{kl}(0)}\right]dx\bigg{|% }\,\lesssim_{\delta}\,L^{\delta-\frac{1}{2}}, (2.9)

where Q_{L}:=[\frac{L}{2},\frac{L}{2})^{d} denotes the cube of sidelength L centered at the origin.

2. If in addition \mathbb{P} is nontrivial, then this effective fluctuation tensor \mathcal{Q} is non-degenerate in the sense that (e\otimes e):\mathcal{Q}\,(e\otimes e)>0 for all e\in\mathbb{R}^{d}\setminus\{0\}.∎

###### Remarks 2.10.

• When applying a covariance inequality (cf. Lemma 5.1 below) to the argument of the limit in the Green-Kubo formula (2.9), we end up with the bound

 \int_{Q_{2L}}\frac{|Q_{L}\cap(x+Q_{L})|}{|Q_{L}|}|\operatorname{Cov}\left[{\Xi% _{ij}(x)};{\Xi_{kl}(0)}\right]\!|\,dx\,\lesssim\,\log L,

which is sharp. The main difficulty to characterize the limiting covariance structure is that, as usual for Green-Kubo formulas, the covariance of the homogenization commutator \Xi is not integrable, and cancellations have to be unravelled.

• The optimal rate (2.8) for the convergence of the covariance structure of J_{0}^{\varepsilon} owes to the very local structure of the homogenization commutator \Xi. Combined with the pathwise result of Corollary 2.4, it carries over to I_{1}^{\varepsilon}, I_{2}^{\varepsilon}, I_{3}^{\varepsilon}, I_{4}^{\varepsilon}, J_{1}^{\varepsilon}, and J_{2}^{\varepsilon}. In [MO, GuM], the Gaussian Helffer-Sjöstrand representation formula for the variance [HS-94, Sjostrand-96, NS2] was already used in order to prove the convergence of the covariance structure of I_{1}^{\varepsilon} and J_{1}^{\varepsilon} for d>2, but the obtained convergence rate was suboptimal in every dimension.

• The non-degeneracy property (ii) already follows from [GN, Proposition 2.1] (modulo the identification (1.16)); see also [MO, Remark 2.3] in the Gaussian case. ∎

The combination of Propositions 2.7 and 2.9 leads to a complete scaling limit result for J_{0}^{\varepsilon}, and proves the convergence in law to a Gaussian white noise. As in Proposition 2.7, we do not know whether the full logarithmic factor is optimal for d=2.

###### Corollary 2.11.

Let d\geq 2, let \mathbb{P} be a product measure, and let \mu_{d} be defined as in (1.13). Let \mathcal{Q} be the 4-tensor defined in Proposition 2.9(i), and let \Gamma denote the 2-tensor Gaussian white noise with covariance tensor \mathcal{Q}, that is, the Gaussian random linear functional with zero expectation \mathbb{E}\left[\Gamma(F)\right]=0 and with covariance structure \operatorname{Cov}\left[{\Gamma(F)};{\Gamma(G)}\right]=\int_{\mathbb{R}^{d}}F:% \mathcal{Q}\,G for all F,G\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d}. Then for all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} the random variable J_{0}^{\varepsilon}(F) converges in law to \Gamma(F), and for \int_{\mathbb{R}^{d}}F:\mathcal{Q}\,F\neq 0 there holds

 (\operatorname{d}_{\operatorname{W}}+\operatorname{d}_{\operatorname{K}})\!% \left({J_{0}^{\varepsilon}(F)\!}\,,\,{\!\Gamma(F)}\right)\leavevmode\nobreak\ % \lesssim_{F}\leavevmode\nobreak\ \varepsilon\mu_{d}(\tfrac{1}{\varepsilon}).\qed
###### Remarks 2.12.

• Combined with the pathwise result of Corollary 2.4, this result leads to a proof of the joint convergence (1.11) and implies in particular quantitative versions of the known scaling limit results for I_{1}^{\varepsilon} and J_{1}^{\varepsilon}: For all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d} and all F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d} the random variables I_{1}^{\varepsilon}(f,g) and J_{1}^{\varepsilon}(F) converge in law to \Gamma(\bar{\mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g) and -\Gamma(\bar{\mathcal{P}}_{H}^{*}F), respectively. Moreover, whenever \int_{\mathbb{R}^{d}}(\bar{\mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g)% :\mathcal{Q}\,(\bar{\mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g)\neq 0 and \int_{\mathbb{R}^{d}}\bar{\mathcal{P}}_{H}^{*}F:\mathcal{Q}\,\bar{\mathcal{P}}% _{H}^{*}F\neq 0, we have

 \displaystyle(\operatorname{d}_{\operatorname{W}}+\operatorname{d}_{% \operatorname{K}})\!\left({I_{1}^{\varepsilon}(f,g)}\,,\,{\Gamma(\bar{\mathcal% {P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g)}\right) \displaystyle\lesssim_{f,g} \displaystyle\varepsilon\mu_{d}(\tfrac{1}{\varepsilon}), \displaystyle(\operatorname{d}_{\operatorname{W}}+\operatorname{d}_{% \operatorname{K}})\!\left({J_{1}^{\varepsilon}(F)}\,,\,{-\Gamma(\bar{\mathcal{% P}}_{H}^{*}F)}\right) \displaystyle\lesssim_{F} \displaystyle\varepsilon\mu_{d}(\tfrac{1}{\varepsilon}).

This extends and unifies [GN, MO, MN, GuM], and yields the first scaling limit results in the critical dimension d=2. Convergence rates are new in any dimension, and are optimal (at least for d>2).

• SPDE representation for the scaling limit of the solution operator. The scaling limit result for I_{1}^{\varepsilon} above indicates that \varepsilon^{-\frac{d}{2}}\nabla_{\varepsilon}(u_{\varepsilon}-\mathbb{E}\left% [u_{\varepsilon}\right]) (seen as a random linear functional) converges in law to the solution DU in \mathbb{R}^{d} of

 -D\cdot\bar{\boldsymbol{a}}DU=D\cdot(\Gamma_{i}D_{i}\bar{u}).

This justifies a posteriori the conclusion (although not the strategy of proof) of the heuristics by Armstrong, Gu, and Mourrat [GuM] in dimensions d\geq 2. (See also [Gu-16] for a rigorous treatment of the easier case of dimension d=1.)

• Scaling limit of the corrector. The scaling limit result for J_{1}^{\varepsilon} above shows that the rescaled corrector field \varepsilon^{-\frac{d}{2}}D\phi(\frac{\cdot}{\varepsilon}) (seen as a random linear functional) converges in law to D(-D\cdot\bar{\boldsymbol{a}}D)^{-1}D\cdot\Gamma, that is, to the gradient of a variant of the so-called Gaussian free field. This variant involves both \bar{\boldsymbol{a}} and \mathcal{Q}. As pointed out in [MO], it is easily checked in Fourier space that this variant does not coincide in general with the standard Gaussian free field (unless the compatibility condition \mathcal{Q}_{ijkl}=\eta_{ik}\bar{\boldsymbol{a}}_{lj} is satisfied for some matrix \eta, which however does not even hold in elementary examples, see e.g. [Gu-Mourrat-17, Section 3] and [D-thesis, equation (5.35)]). This variant of the Gaussian free field is studied in [Gu-Mourrat-17], where it is shown to be Markovian only in the standard case. In the critical dimension d=2, since the whole-space Gaussian free field is not well-defined (only its gradient is), this implies the non-existence of stationary correctors.

• Gu and Mourrat’s observation. With the above results at hand, we recover the observation by Gu and Mourrat [GuM] that the usual two-scale expansion (1.3) of u_{\varepsilon} is not accurate in the fluctuation scaling. The above indeed shows that the fluctuations of \varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla_{\varepsilon}(u_{% \varepsilon}-\mathbb{E}\left[u_{\varepsilon}\right]) and of \varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla\phi_{i}(\frac{% \cdot}{\varepsilon})\nabla_{\varepsilon,i}\bar{u} are asymptotically given by \Gamma(\bar{\mathcal{P}}_{H}f\otimes\bar{\mathcal{P}}_{H}^{*}g) and by \Gamma(\bar{\mathcal{P}}_{H}^{*}((\bar{\mathcal{P}}_{H}f)\otimes g)), respectively, and therefore do not coincide. ∎

### 2.5. Approximation of the fluctuation tensor

We finally turn to the representative volume element (RVE) approximation of \mathcal{Q} (cf. Theorem 2). Indeed, the Green-Kubo formula (2.9) for the fluctuation tensor \mathcal{Q} is of no practical use in applications since it requires to solve the corrector equation on the whole space and for every realization of the random coefficient field. It is therefore natural to seek a suitable RVE approximation. It consists in introducing an artificial period L>0 and in considering an L-periodized coefficient field \boldsymbol{a}_{L}, typically given by a suitable periodization in law (cf. [EGMN-15]). In the present iid setting, we simply define \boldsymbol{a}_{L}(x+Ly):=\boldsymbol{a}(x) for all y\in\mathbb{Z}^{d} and x\in Q_{L}:=[-\frac{L}{2},\frac{L}{2})^{d}. Note that the map \boldsymbol{a}\mapsto\boldsymbol{a}_{L} on \Omega pushes forward the measure \mathbb{P} to a measure \mathbb{P}_{L} concentrated on L-periodic coefficients, so that we may view \boldsymbol{a}_{L} as an element of \Omega_{L}=[\lambda,1]^{\mathcal{B}_{L}} equipped with the product measure \mathbb{P}_{L}=\pi^{\otimes\mathcal{B}_{L}}, where \mathcal{B}_{L}:=\{(x,x+e_{i}):x\in Q_{L}\cap\mathbb{Z}^{d},1\leq i\leq d\}. We then define the L-periodized corrector \phi_{L,i} in the direction e_{i} as the unique L-periodic solution in Q_{L}\cap\mathbb{Z}^{d} of

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}_{L}(\nabla\phi_{L,i}+e_{i})=0, (2.10)

satisfying \sum_{z\in Q_{L}\cap\mathbb{Z}^{d}}\phi_{L,i}(z)=0, and we set \phi_{L}:=(\phi_{L,i})_{i=1}^{d} (which we implicitly extend as usual into a periodic piecewise constant map on \mathbb{R}^{d}). The spatial average of the flux,

 \bar{\boldsymbol{a}}_{L}e_{i}:=\fint_{Q_{L}}\boldsymbol{a}_{L}(\nabla\phi_{L,i% }+e_{i}),

is then an RVE approximation for the homogenized coefficient \bar{\boldsymbol{a}}e_{i}=\mathbb{E}\left[\boldsymbol{a}(\nabla\phi_{i}+e_{i})\right]. The optimal numerical analysis of this approximation was originally performed in [GO1, GO2, GNO1], where it was established that for all L\geq 2 there holds

 \displaystyle|\mathrm{Var}\left[\bar{\boldsymbol{a}}_{L}\right]\!|^{\frac{1}{2% }}\lesssim L^{-\frac{d}{2}},\qquad|\mathbb{E}\left[\bar{\boldsymbol{a}}_{L}% \right]-\bar{\boldsymbol{a}}|\lesssim L^{-d}\log^{d}L. (2.11)

In Theorem 2, we claim that the fluctuation tensor \mathcal{Q} coincides with the limit of the rescaled variance of \bar{\boldsymbol{a}}_{L}^{*}. In addition, this characterization naturally leads to an RVE approximation \mathcal{Q}_{L,N} for \mathcal{Q}, of which we obtain the optimal error estimate.

###### Remarks 2.13.

• Definition (1.17) for \mathcal{Q}_{L,N} is equivalent to

 \displaystyle\mathcal{Q}_{L,N}=\frac{L^{d}}{N-1}\sum_{n=1}^{N}\Big{(}\fint_{Q_% {L}}\Xi_{L,N}^{(n)}\Big{)}\otimes\Big{(}\fint_{Q_{L}}\Xi_{L,N}^{(n)}\Big{)}, (2.12)

where

 \displaystyle\Xi_{L,N,i}^{(n)}:=\boldsymbol{a}_{L}^{(n)}(\nabla\phi_{L,i}^{(n)% }+e_{i})-\bar{\boldsymbol{a}}_{L,N}(\nabla\phi_{L,i}^{(n)}+e_{i}),

with the obvious notation \nabla\phi_{L}^{(n)}:=\nabla\phi_{L}(\boldsymbol{a}_{L}^{(n)}). Since by stationarity

 \int_{Q_{L}}\operatorname{Cov}\left[{\Xi_{L,N}(x)};{\Xi_{L,N}(0)}\right]dx=L^{% d}\,\mathrm{Var}\left[\fint_{Q_{L}}\Xi_{L,N}\right],

formula (2.12) is in the spirit of the Green-Kubo formula (2.9).

• In (2.11) the standard deviation |\mathrm{Var}\left[\bar{\boldsymbol{a}}_{L}\right]\!|^{\frac{1}{2}} of the RVE approximation for \bar{\boldsymbol{a}} is seen to be O(L^{\frac{d}{2}}) times larger than the systematic error |\mathbb{E}\left[\bar{\boldsymbol{a}}_{L}\right]-\bar{\boldsymbol{a}}| (up to a logarithmic correction). In practice, we rather use \bar{\boldsymbol{a}}_{L,N} as an approximation for \bar{\boldsymbol{a}},

since in the regime N\simeq L^{d} the standard deviation becomes of the same order as the systematic error O(L^{-d}). Combining this with the estimates in Theorem 2, since \mathcal{Q}_{L,N} is extracted at no further cost than \bar{\boldsymbol{a}}_{L,N} itself, we may infer that an RVE approximation for \mathcal{Q} with accuracy O(L^{-\frac{d}{2}}) is extracted at the same cost as an RVE approximation for \bar{\boldsymbol{a}} with accuracy O(L^{-d}).

• In [N, Nolen-16, GN] (see also [R, BSW]), the fluctuations of the RVE approximation \bar{\boldsymbol{a}}_{L} for the homogenized coefficient \bar{\boldsymbol{a}} was investigated. Combined with the characterization (1.16) of the limit of the rescaled variance, the main result in [GN] takes on the following guise, for all L\geq 2 and all N\geq 1,

 \sup_{e\in\mathbb{R}^{d}\setminus\{0\}}\,(\operatorname{d}_{\operatorname{W}}+% \operatorname{d}_{\operatorname{K}})\!\left({N^{\frac{1}{2}}L^{\frac{d}{2}}% \frac{e\cdot(\bar{\boldsymbol{a}}_{L,N}-\bar{\boldsymbol{a}})e}{(e\otimes e:% \mathcal{Q}:e\otimes e)^{\frac{1}{2}}}}\,,\,{\mathcal{N}}\right)\lesssim N^{-% \frac{1}{2}}L^{-\frac{d}{2}}\log^{d}L.\qed

## 3. Pathwise structure

Henceforth we place ourselves in the discrete setting of Section 2. In the present section, we establish the pathwise result stated in Proposition 2.2, that is, the main novelty of this contribution. Similar estimates also lead to the CLT scaling result stated in Proposition 2.1, and we further deduce Corollary 2.4.

### 3.1. Structure of the proof and auxiliary results

The main tool that we use to prove Propositions 2.1 and 2.2 is the following Poincaré inequality (or spectral gap estimate) in the probability space, which holds for any product measure \mathbb{P} on \Omega (see e.g. [GO1, Lemma 2.3] for a proof). Let us first fix some notation. Let X=X(\boldsymbol{a}) be a random variable on \Omega, that is, a measurable function of (a(b))_{b\in\mathcal{B}}. We choose an iid copy \boldsymbol{a}^{\prime} of \boldsymbol{a},333Although we are then working on a product probability space \Omega\times\Omega, we use for simplicity the same notation \mathbb{P} (and \mathbb{E}) for the product probability measure (and expectation), that is, with respect to both \boldsymbol{a} and \boldsymbol{a}^{\prime}. and for all b\in\mathcal{B} we denote by \boldsymbol{a}^{b} the random field that coincides with \boldsymbol{a} on all edges b^{\prime}\neq b and with \boldsymbol{a}^{\prime} on edge b. In particular, \boldsymbol{a} and \boldsymbol{a}^{b} have the same law. We use the abbreviation X^{b}=X(\boldsymbol{a}^{b}) and define the difference operator \Delta_{b}X:=X-X^{b}, which we call the (Glauber) vertical derivative at edge b.

###### Lemma 3.1 (e.g. [GO1]).

Let \mathbb{P} be a product measure. For all X=X(\boldsymbol{a})\in\operatorname{L}^{2}(\Omega) we have

 \mathrm{Var}\left[X\right]\,\leq\,\frac{1}{2}\,\mathbb{E}\left[\sum_{b\in% \mathcal{B}}|\Delta_{b}X|^{2}\right].\qed

Next to the corrector \phi, we need to recall the notion of flux corrector \sigma, which was recently introduced in [GNO-reg] in the continuum stochastic setting (see also [GuM, Lemma 4.4] and [Nguyen-thesis, Proposition III.2.2] for its subsequent introduction in the discrete setting) and was crucially used in [GNO-quant, GO4]. It allows to put the equation for the two-scale homogenization error in divergence form (cf. (3.29)). Let \sigma=(\sigma_{ijk})_{i,j,k=1}^{d} be the 3-tensor defined as the unique solution in \mathbb{Z}^{d} of

 -\triangle\sigma_{ijk}\,:=\,-\nabla^{*}\cdot\nabla\sigma_{ijk}\,=\,\nabla_{j}q% _{ik}-\nabla_{k}q_{ij}, (3.1)

with \nabla\sigma stationary and having finite second moment, and with \sigma(0)=0, where q_{i} denotes the flux of the corrector

Note that for all i the 2-tensor field \sigma_{i}:=(\sigma_{ijk})_{j,k=1}^{d} is skew-symmetric, that is,

 \sigma_{ijk}=-\sigma_{ikj}, (3.3)

and is shown to satisfy

 \nabla^{*}\cdot\sigma_{i}\,:=\,e_{j}\nabla_{k}^{*}\sigma_{ijk}\,=\,q_{i}. (3.4)

Although considering a symmetric coefficient field, we use non-symmetric notation in view of the extension to the continuum setting, and we denote by \phi^{*} and \sigma^{*} the corrector and flux corrector associated with the pointwise transpose coefficient field \boldsymbol{a}^{*}. For symmetric coefficient fields, \phi^{*}=\phi and \sigma^{*}=\sigma.

We now describe the string of arguments that leads to Proposition 2.2. We start with a suitable decomposition of the vertical derivative of E^{\varepsilon}(f,g), which is key to the proof. Note that we rather consider a suitable version E^{\varepsilon}_{0}(f,g) of E^{\varepsilon}(f,g), which only coincides with E^{\varepsilon}(f,g) up to some minor discretization error (in the continuum setting \bar{u}_{\varepsilon} and \bar{v}_{\varepsilon} would simply coincide with \bar{u}(\varepsilon\cdot) and \bar{v}(\varepsilon\cdot)). In the proofs, it is convenient to rescale all quantities down to scale 1.

###### Lemma 3.2.

For all \varepsilon>0 and all f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}, setting f_{\varepsilon}:=f(\varepsilon\cdot) and g_{\varepsilon}:=g(\varepsilon\cdot), we denote by \bar{u}_{\varepsilon} and \bar{v}_{\varepsilon} the unique Lax-Milgram solutions in \mathbb{R}^{d} of

 \displaystyle-\nabla^{*}\cdot\bar{\boldsymbol{a}}\nabla\bar{u}_{\varepsilon}=% \nabla^{*}\cdot(\varepsilon f_{\varepsilon}),\qquad-\nabla^{*}\cdot\bar{% \boldsymbol{a}}^{*}\nabla\bar{v}_{\varepsilon}=\nabla^{*}\cdot(\varepsilon g_{% \varepsilon}), (3.5)

and we define

 \displaystyle E^{\varepsilon}_{0}(f,g)\,:=\,\varepsilon^{\frac{d}{2}-1}\int_{% \mathbb{R}^{d}}g_{\varepsilon}\cdot\big{(}\boldsymbol{a}\nabla(u_{\varepsilon}% (\varepsilon\cdot))-\bar{\boldsymbol{a}}\nabla(u_{\varepsilon}(\varepsilon% \cdot))-\mathbb{E}\left[\boldsymbol{a}\nabla(u_{\varepsilon}(\varepsilon\cdot)% )-\bar{\boldsymbol{a}}\nabla(u_{\varepsilon}(\varepsilon\cdot))\right]\big{)}% \\ \displaystyle-\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}g_{\varepsilon}% \cdot\Xi_{i}\nabla_{i}\bar{u}_{\varepsilon}, (3.6)

as well as the two-scale expansion error w_{f,\varepsilon}:=u_{\varepsilon}(\varepsilon\cdot)-(1+\phi_{i}\nabla_{i})% \bar{u}_{\varepsilon}. Then we have for all b\in\mathcal{B},

 \displaystyle\Delta_{b}E_{0}^{\varepsilon}(f,g)=\varepsilon^{\frac{d}{2}-1}% \int_{\mathbb{R}^{d}}g_{\varepsilon,j}(\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b% }\boldsymbol{a}(\nabla w^{b}_{f,\varepsilon}+\phi_{i}^{b}\nabla\nabla_{i}\bar{% u}_{\varepsilon})\\ \displaystyle+\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(% \cdot+e_{k})\nabla_{k}g_{\varepsilon,j}e_{k}\cdot\Delta_{b}\boldsymbol{a}% \nabla(u_{\varepsilon}^{b}(\varepsilon\cdot))\\ \displaystyle-\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(% \cdot+e_{k})\nabla_{k}(g_{\varepsilon,j}\nabla_{i}\bar{u}_{\varepsilon})e_{k}% \cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i})\\ \displaystyle+\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\nabla r_{% \varepsilon}\cdot\Delta_{b}\boldsymbol{a}\nabla(u_{\varepsilon}^{b}(% \varepsilon\cdot))-\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\nabla R_{% \varepsilon,i}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}), (3.7)

where the auxiliary fields r_{\varepsilon} and R_{\varepsilon}=(R_{\varepsilon,i})_{i=1}^{d} are the unique Lax-Milgram solutions in \mathbb{R}^{d} of

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}^{*}\nabla r_{\varepsilon} \displaystyle= \displaystyle\nabla_{l}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})\boldsymbol{a}_{kl}% \nabla_{k}g_{\varepsilon,j}+\sigma_{jkl}^{*}(\cdot-e_{k})\nabla_{k}^{*}g_{% \varepsilon,j}\big{)}, (3.8) \displaystyle-\nabla^{*}\cdot\boldsymbol{a}^{*}\nabla R_{\varepsilon,i} \displaystyle= \displaystyle\nabla_{l}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})\boldsymbol{a}_{kl}% \nabla_{k}(g_{\varepsilon,j}\nabla_{i}\bar{u}_{\varepsilon})+\sigma_{jkl}^{*}(% \cdot-e_{k})\nabla_{k}^{*}(g_{\varepsilon,j}\nabla_{i}\bar{u}_{\varepsilon})% \big{)}. (3.9)

By the spectral gap estimate of Lemma 3.1, the desired pathwise result (2.6) follows from a suitable estimate of the sum over \mathcal{B} of the squares of the right-hand side terms in (3.7). For that purpose, we make crucial use of the following moment bounds for the extended corrector (\phi,\sigma) and its gradient. (These bounds are a variation of [GNO1] and are the discrete versions of a result in [GNO-quant], the proof of which extends to the discrete setting considered here.)

###### Lemma 3.3 ([GNO1, GNO-quant]).

Let d\geq 2, let \mathbb{P} be a product measure, and let \mu_{d} be defined in (1.13). For all q<\infty and all z\in\mathbb{Z}^{d} we have

 \mathbb{E}\left[|\phi(z)|^{q}\right]^{\frac{1}{q}}+\mathbb{E}\left[|\sigma(z)|% ^{q}\right]^{\frac{1}{q}}\,\lesssim_{q}\,\mu_{d}(|z|)^{\frac{1}{2}},

and

 \mathbb{E}\left[|\nabla\phi(z)|^{q}\right]^{\frac{1}{q}}+\mathbb{E}\left[|% \nabla\sigma(z)|^{q}\right]^{\frac{1}{q}}\,\lesssim_{q}\,1.\qed

An additional crucial ingredient is the following large-scale weighted Calderón-Zygmund estimate for the operator -\nabla^{*}\cdot\boldsymbol{a}\nabla. (A proof in the continuum setting was originally given in the first version of this article, see now [GNO-reg, Corollary 5] — the adaptation to the discrete setting is straightforward since it is solely based on the energy and Caccioppoli estimates.)

###### Lemma 3.4 ([GNO-reg]).

Let d\geq 1 and let \mathbb{P} be a product measure. There exists a \frac{1}{8}-Lipschitz stationary random field r_{*}\geq 1 on \mathbb{R}^{d} with \mathbb{E}\left[r_{*}^{q}\right]\lesssim_{q}1 for all q<\infty, such that the following holds almost surely: For \varepsilon>0, 2\leq p<\infty, and 0\leq\gamma<d(p-1), for any (sufficiently fast) decaying scalar field w and vector field h related in \mathbb{R}^{d} by

 -\nabla^{*}\cdot\boldsymbol{a}\nabla w=\nabla^{*}\cdot h,

we have

 \int_{\mathbb{R}^{d}}\big{(}1+\varepsilon(|x|+r_{*}(0))\big{)}^{\gamma}\Big{(}% \fint_{B_{*}(x)}|\nabla w|^{2}\Big{)}^{\frac{p}{2}}dx\,\lesssim_{\gamma,p}\,% \int_{\mathbb{R}^{d}}\big{(}1+\varepsilon(|x|+r_{*}(0))\big{)}^{\gamma}|h(x)|^% {p}dx

with the short-hand notation B_{*}(x):=B_{r_{*}(x)}(x). ∎

### 3.2. Proof of Proposition 2.1

We focus on J_{0}^{\varepsilon}, while the proof is similar for J_{1}^{\varepsilon} and J_{2}^{\varepsilon}. Let F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d}, and set F_{\varepsilon}:=F(\varepsilon\cdot). We split the proof into two steps: we start by giving a suitable representation formula for the vertical derivative \Delta_{b}J_{0}^{\varepsilon}(F), and then apply the spectral gap estimate.

Step 1. Representation formula for \Delta_{b}J_{0}^{\varepsilon}(F):

 \displaystyle\Delta_{b}J_{0}^{\varepsilon}(F)=\varepsilon^{\frac{d}{2}}\int_{% \mathbb{R}^{d}}F_{\varepsilon,ij}e_{j}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi% _{i}^{b}+e_{i})+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\nabla s_{% \varepsilon,i}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}), (3.10)

where the auxiliary field s_{\varepsilon} is the unique Lax-Milgram solution in \mathbb{R}^{d} of

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}^{*}\nabla s_{\varepsilon,i}=\nabla% ^{*}\cdot\big{(}F_{\varepsilon,ij}(\boldsymbol{a}-\bar{\boldsymbol{a}})e_{j}% \big{)}. (3.11)

By definition of the homogenization commutator,

 \Delta_{b}J_{0}^{\varepsilon}(F)=\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}% }F_{\varepsilon,ij}e_{j}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}% )+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}e_{j}\cdot(% \boldsymbol{a}-\bar{\boldsymbol{a}})\nabla\Delta_{b}\phi_{i}.

By definition (3.11) of s_{\varepsilon,i}, we find

 \Delta_{b}J_{0}^{\varepsilon}(F)=\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}% }F_{\varepsilon,ij}e_{j}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}% )-\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\nabla s_{\varepsilon,i}\cdot% \boldsymbol{a}\nabla\Delta_{b}\phi_{i}.

Using then the vertical derivative of the corrector equation (2.4) in the form

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}\phi_{i}=\nabla^{*}% \cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}), (3.12)

the claim (3.10) follows.

Step 2. Conclusion.

For b\in\mathcal{B} we use the notation b=(z_{b},z_{b}+e_{b}). Inserting the representation formula (3.10) in the spectral gap estimate of Lemma 3.1, and noting that |\Delta_{b}\boldsymbol{a}(x)|\lesssim\mathds{1}_{Q(z_{b})}(x), we obtain

 \displaystyle\mathrm{Var}\left[J_{0}^{\varepsilon}(F)\right]\,\lesssim\,% \varepsilon^{d}\sum_{b\in\mathcal{B}}\mathbb{E}\left[|\nabla\phi^{b}(z_{b})+% \operatorname{Id}|^{2}\right]\int_{Q(z_{b})}|F_{\varepsilon}|^{2}\\ \displaystyle+\varepsilon^{d}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}|\nabla% \phi^{b}(z_{b})+\operatorname{Id}|^{2}\int_{Q(z_{b})}|\nabla s_{\varepsilon}|^% {2}\right],

and hence, by Lemma 3.3 in the form \mathbb{E}\left[|\nabla\phi^{b}|^{2}\right]=\mathbb{E}\left[|\nabla\phi|^{2}% \right]\lesssim 1,

 \displaystyle\mathrm{Var}\left[J_{0}^{\varepsilon}(F)\right]\,\lesssim\,% \varepsilon^{d}\|F_{\varepsilon}\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}^{2}+% \varepsilon^{d}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}|\nabla\phi^{b}(z_{b})+% \operatorname{Id}|^{2}\int_{Q(z_{b})}|\nabla s_{\varepsilon}|^{2}\right]. (3.13)

It remains to estimate the last right-hand side term. Using equation (3.12) in the form -\nabla^{*}\cdot\boldsymbol{a}^{b}\nabla(\phi^{b}-\phi)=\nabla^{*}\cdot(% \boldsymbol{a}^{b}-\boldsymbol{a})(\nabla\phi+\operatorname{Id}), an energy estimate yields

 \displaystyle|\nabla(\phi^{b}-\phi)(z_{b})|^{2}\leq\int_{\mathbb{R}^{d}}|% \nabla(\phi^{b}-\phi)|^{2}\lesssim\int_{\mathbb{R}^{d}}|\boldsymbol{a}^{b}-% \boldsymbol{a}|^{2}|\nabla\phi+\operatorname{Id}|^{2}\lesssim|\nabla\phi(z_{b}% )+\operatorname{Id}|^{2}, (3.14)

so that |\nabla\phi^{b}(z_{b})+\operatorname{Id}|\lesssim|\nabla\phi(z_{b})+% \operatorname{Id}|. Further estimating in (3.13) integrals over unit cubes by integrals over balls at scale r_{*} (cf. Lemma 3.4), smuggling in a power \alpha\frac{p-1}{p} of the weight w_{\varepsilon}(z):=1+\varepsilon|z|, and applying Hölder’s inequality in space with exponent p, we deduce for all p>1,

 \displaystyle\varepsilon^{d}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}|\nabla% \phi^{b}(z_{b})+\operatorname{Id}|^{2}\int_{Q(z_{b})}|\nabla s_{\varepsilon}|^% {2}\right]\,\lesssim\,\varepsilon^{d}\,\mathbb{E}\left[\int_{\mathbb{R}^{d}}|% \nabla\phi(z)+\operatorname{Id}|^{2}\Big{(}\int_{Q_{2}(z)}|\nabla s_{% \varepsilon}|^{2}\Big{)}dz\right]\\ \displaystyle\lesssim\,\varepsilon^{d}\,\mathbb{E}\bigg{[}\bigg{(}\int_{% \mathbb{R}^{d}}|\nabla\phi(z)+\operatorname{Id}|^{\frac{2p}{p-1}}r_{*}(z)^{% \frac{dp}{p-1}}w_{\varepsilon}(z)^{-\alpha}dz\bigg{)}^{\frac{p-1}{p}}\\ \displaystyle\times\bigg{(}\int_{\mathbb{R}^{d}}w_{\varepsilon}(z)^{\alpha(p-1% )}\Big{(}\fint_{B_{*}(z)}|\nabla s_{\varepsilon}|^{2}\Big{)}^{p}dz\bigg{)}^{% \frac{1}{p}}\bigg{]}.

Applying Hölder’s inequality in the probability space, using Lemmas 3.3 and 3.4 in the form \mathbb{E}\big{[}|\nabla\phi+\operatorname{Id}|^{q}+r_{*}^{q}\big{]}\lesssim_{% q}1 for all q<\infty, and noting that \int_{\mathbb{R}^{d}}w_{\varepsilon}(z)^{-\alpha}dz\lesssim_{\alpha}% \varepsilon^{-d} provided \alpha>d, we obtain for all p>1 and all \alpha>d,

 \displaystyle\varepsilon^{d}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}|\nabla% \phi^{b}(z_{b})+\operatorname{Id}|^{2}\int_{Q(z_{b})}|\nabla s_{\varepsilon}|^% {2}\right]\\ \displaystyle\lesssim_{\alpha,p}\,\varepsilon^{\frac{d}{p}}\,\mathbb{E}\left[% \int_{\mathbb{R}^{d}}w_{\varepsilon}(z)^{\alpha(p-1)}\Big{(}\fint_{B_{*}(z)}|% \nabla s_{\varepsilon}|^{2}\Big{)}^{p}dz\right]^{\frac{1}{p}}. (3.15)

By large-scale weighted Calderón-Zygmund theory (cf. Lemma 3.4) applied to equation (3.11) for s_{\varepsilon} with \alpha(p-1)<d(2p-1), using again the moment bounds on r_{*}, and rescaling spatial integrals, we deduce for all 0<p-1\ll 1 and all 0<\alpha-d\ll 1,

 \displaystyle\varepsilon^{d}\mathbb{E}\left[\sum_{b\in\mathcal{B}}|\nabla\phi^% {b}(z_{b})+\operatorname{Id}|^{2}\int_{Q(z_{b})}|\nabla s_{\varepsilon}|^{2}\right] \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{\frac{d}{p}}\,\mathbb{E}\left[r_{*}(0)^{\alpha(p-1)}% \int_{\mathbb{R}^{d}}w_{\varepsilon}^{\alpha(p-1)}|F_{\varepsilon}|^{2p}\right% ]^{\frac{1}{p}} (3.16) \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{\frac{d}{p}}\|w_{\varepsilon}^{\alpha\frac{p-1}{2p}}% F_{\varepsilon}\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})}^{2}.

Inserting this into (3.13) and rescaling spatial integrals, we deduce for all 0<p-1\ll 1 and all 0<\alpha-d\ll 1,

 \mathrm{Var}\left[J_{0}^{\varepsilon}(F)\right]\,\lesssim_{\alpha,p}\,\|F\|_{% \operatorname{L}^{2}(\mathbb{R}^{d})}^{2}+\|w_{1}^{\alpha\frac{p-1}{2p}}F\|_{% \operatorname{L}^{2p}(\mathbb{R}^{d})}^{2}.

Further using Hölder’s inequality in the form

 \|F\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}\leq\Big{(}\int_{\mathbb{R}^{d}}w_% {1}^{-\alpha}\Big{)}^{\frac{p-1}{2p}}\Big{(}\int_{\mathbb{R}^{d}}w_{1}^{\alpha% (p-1)}|F|^{2p}\Big{)}^{\frac{1}{2p}}\lesssim_{\alpha,p}\|w_{1}^{\alpha\frac{p-% 1}{2p}}F\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})},

the conclusion follows (after replacing the exponent \alpha\frac{p-1}{2p} by 2\alpha).

### 3.3. Proof of Lemma 3.2

We split the proof into two steps. To simplify notation, in this proof (and only in this proof), we write u:=u_{\varepsilon}(\varepsilon\cdot).

Step 1. Representation formula for \Delta_{b}((\boldsymbol{a}-\bar{\boldsymbol{a}})\nabla u):

 \displaystyle\Delta_{b}\big{(}e_{j}\cdot(\boldsymbol{a}-\bar{\boldsymbol{a}})% \nabla u\big{)}=(\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b}\boldsymbol{a}\nabla u% ^{b}-\nabla_{k}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})e_{k}\cdot\Delta_{b}% \boldsymbol{a}\nabla u^{b}\big{)}\\ \displaystyle-\nabla_{k}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})e_{k}\cdot% \boldsymbol{a}\nabla\Delta_{b}u\big{)}-\nabla_{k}\big{(}\sigma_{jkl}^{*}(\cdot% -e_{k})\nabla_{l}\Delta_{b}u\big{)}. (3.17)

In particular, replacing x\mapsto u(x) by x\mapsto\phi_{i}(x)+x_{i}, we deduce the following discrete version of (1.10),

 \displaystyle\Delta_{b}\Xi_{ij}=(\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b}% \boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i})-\nabla_{k}^{*}\big{(}\phi_{j}^{*}(% \cdot+e_{k})e_{k}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i})\big{)% }\\ \displaystyle-\nabla_{k}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})\boldsymbol{a}_{kl% }\nabla_{l}\Delta_{b}\phi_{i}\big{)}-\nabla_{k}\big{(}\sigma_{jkl}^{*}(\cdot-e% _{k})\nabla_{l}\Delta_{b}\phi_{i}\big{)}. (3.18)

Using the definition (3.4) of \sigma_{j}^{*} in the form (\boldsymbol{a}^{*}-\bar{\boldsymbol{a}}^{*})e_{j}=-\boldsymbol{a}^{*}\nabla% \phi_{j}^{*}+\nabla^{*}\cdot\sigma_{j}^{*}, we find

 \displaystyle\Delta_{b}\big{(}e_{j}\cdot(\boldsymbol{a}-\bar{\boldsymbol{a}})% \nabla u\big{)} \displaystyle= \displaystyle e_{j}\cdot\Delta_{b}\boldsymbol{a}\nabla u^{b}+e_{j}\cdot(% \boldsymbol{a}-\bar{\boldsymbol{a}})\nabla\Delta_{b}u (3.19) \displaystyle= \displaystyle e_{j}\cdot\Delta_{b}\boldsymbol{a}\nabla u^{b}+(\nabla^{*}\cdot% \sigma_{j}^{*})\cdot\nabla\Delta_{b}u-\nabla\phi_{j}^{*}\cdot\boldsymbol{a}% \nabla\Delta_{b}u.

On the one hand, using the following discrete version of the Leibniz rule, for all \chi_{1},\chi_{2}:\mathbb{Z}^{d}\to\mathbb{R},

 \displaystyle\nabla_{l}^{*}(e_{l}\chi_{1}(\cdot+e_{l})\chi_{2})=\chi_{2}\nabla% \chi_{1}+\chi_{1}\nabla^{*}\chi_{2}, (3.20)

we obtain

 (\nabla^{*}\cdot\sigma_{j}^{*})\cdot\nabla\Delta_{b}u\,=\,\nabla_{l}^{*}\big{(% }\sigma_{jkl}^{*}\nabla_{k}\Delta_{b}u(\cdot+e_{l})\big{)}-\sigma_{jkl}^{*}% \nabla_{k}\nabla_{l}\Delta_{b}u,

so that the skew-symmetry (3.3) of \sigma_{j} leads to

 \displaystyle(\nabla^{*}\cdot\sigma_{j}^{*})\cdot\nabla\Delta_{b}u\,=\,-\nabla% _{k}^{*}\big{(}\sigma_{jkl}^{*}\nabla_{l}\Delta_{b}u(\cdot+e_{k})\big{)}\,=\,-% \nabla_{k}\big{(}\sigma_{jkl}^{*}(\cdot-e_{k})\nabla_{l}\Delta_{b}u\big{)}. (3.21)

On the other hand, using the vertical derivative of equation (2.1) in the form -\nabla^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}u=\nabla^{*}\cdot\Delta_{b}% \boldsymbol{a}\nabla u^{b}, the discrete Leibniz rule (3.20) yields

 \displaystyle\nabla\phi_{j}^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}u \displaystyle= \displaystyle-\phi_{j}^{*}\nabla^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}u+% \nabla_{k}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})e_{k}\cdot\boldsymbol{a}\nabla% \Delta_{b}u\big{)} (3.22) \displaystyle= \displaystyle-\nabla\phi_{j}^{*}\cdot\Delta_{b}\boldsymbol{a}\nabla u^{b}+% \nabla_{k}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})e_{k}\cdot\Delta_{b}\boldsymbol{% a}\nabla u^{b}\big{)} \displaystyle                                       +\nabla_{k}^{*}\big{(}\phi% _{j}^{*}(\cdot+e_{k})e_{k}\cdot\boldsymbol{a}\nabla\Delta_{b}u\big{)}.

Inserting (3.21) and (3.22) into (3.19), the claim (3.17) follows.

Step 2. Conclusion.

Integrating identities (3.17) and (3.18) with the test functions g_{\varepsilon} and \nabla\bar{u}_{\varepsilon}\otimes g_{\varepsilon}, respectively, and integrating by parts, we obtain by definition of E^{\varepsilon}_{0},

 \displaystyle\Delta_{b}E_{0}^{\varepsilon}(f,g) \displaystyle= \displaystyle\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}g_{\varepsilon,j}% (\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b}\boldsymbol{a}\big{(}\nabla u^{b}-(% \nabla\phi_{i}^{b}+e_{i})\nabla_{i}\bar{u}_{\varepsilon}\big{)} \displaystyle+\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(% \cdot+e_{k})\nabla_{k}g_{\varepsilon,j}e_{k}\cdot\Delta_{b}\boldsymbol{a}% \nabla u^{b} \displaystyle-\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(% \cdot+e_{k})\nabla_{k}(g_{\varepsilon,j}\nabla_{i}\bar{u}_{\varepsilon})e_{k}% \cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}) \displaystyle+\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\big{(}\phi_{j}^% {*}(\cdot+e_{k})\boldsymbol{a}_{kl}\nabla_{k}g_{\varepsilon,j}+\sigma_{jkl}^{*% }(\cdot-e_{k})\nabla_{k}^{*}g_{\varepsilon,j}\big{)}\nabla_{l}\Delta_{b}u \displaystyle-\varepsilon^{\frac{d}{2}-1}\int_{\mathbb{R}^{d}}\big{(}\phi_{j}^% {*}(\cdot+e_{k})\boldsymbol{a}_{kl}\nabla_{k}(g_{\varepsilon,j}\nabla_{i}\bar{% u}_{\varepsilon})+\sigma_{jkl}^{*}(\cdot-e_{k})\nabla_{k}^{*}(g_{\varepsilon,j% }\nabla_{j}\bar{u}_{\varepsilon})\big{)}\nabla_{l}\Delta_{b}\phi_{i}.

The first right-hand side term is reformulated using the definition of w_{f,\varepsilon} in the form \nabla u^{b}-(\nabla\phi_{i}^{b}+e_{i})\nabla_{i}\bar{u}_{\varepsilon}=\nabla w% ^{b}_{f,\varepsilon}+\phi_{i}^{b}\nabla\nabla_{i}\bar{u}_{\varepsilon}. It remains to post-process the last two right-hand side terms. Using equation (3.8) for r_{\varepsilon} and using the vertical derivative of equation (2.1) for u_{\varepsilon} in the form -\nabla^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}u=\nabla^{*}\cdot\Delta_{b}% \boldsymbol{a}\nabla u^{b}, we find

 \displaystyle\int_{\mathbb{R}^{d}}\big{(}\phi_{j}^{*}(\cdot+e_{k})\boldsymbol{% a}_{kl}\nabla_{k}g_{\varepsilon,j}+\sigma_{jkl}^{*}(\cdot-e_{k})\nabla_{k}^{*}% g_{\varepsilon,j}\big{)}\nabla_{l}\Delta_{b}u\\ \displaystyle=-\int_{\mathbb{R}^{d}}\nabla r_{\varepsilon}\cdot\boldsymbol{a}% \nabla\Delta_{b}u=\int_{\mathbb{R}^{d}}\nabla r_{\varepsilon}\cdot\Delta_{b}% \boldsymbol{a}\nabla u^{b}.

Similarly, equations (3.9) and (3.12) lead to

 \displaystyle\int_{\mathbb{R}^{d}}\big{(}\phi_{j}^{*}(\cdot+e_{k})\boldsymbol{% a}_{kl}\nabla_{k}(g_{\varepsilon,j}\nabla_{i}\bar{u}_{\varepsilon})+\sigma_{% jkl}^{*}(\cdot-e_{k})\nabla_{k}^{*}(g_{\varepsilon,j}\nabla_{i}\bar{u}_{% \varepsilon})\big{)}\nabla_{l}\Delta_{b}\phi_{i}\\ \displaystyle=-\int_{\mathbb{R}^{d}}\nabla R_{\varepsilon,i}\cdot\boldsymbol{a% }\nabla\Delta_{b}\phi_{i}=\int_{\mathbb{R}^{d}}\nabla R_{\varepsilon,i}\cdot% \Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i}),

and the conclusion follows.

### 3.4. Proof of Proposition 2.2

Using the representation formula (3.7), and recalling that for symmetric coefficients we have (\phi^{*},\sigma^{*})=(\phi,\sigma), the spectral gap estimate of Lemma 3.1 leads to

 \displaystyle\mathrm{Var}\left[E_{0}^{\varepsilon}(f,g)\right]\,\lesssim\,T_{1% }^{\varepsilon}+T_{2}^{\varepsilon}+T_{3}^{\varepsilon}+T_{4}^{\varepsilon}+T_% {5}^{\varepsilon}, (3.23)

where we have set

 \displaystyle T_{1}^{\varepsilon} \displaystyle:= \displaystyle\sum_{b\in\mathcal{B}}\mathbb{E}\left[\Big{(}\varepsilon^{\frac{d% }{2}-1}\int_{\mathbb{R}^{d}}g_{\varepsilon,j}(\nabla\phi_{j}+e_{j})\cdot\Delta% _{b}\boldsymbol{a}(\nabla w^{b}_{f,\varepsilon}+\phi_{i}^{b}\nabla\nabla_{i}% \bar{u}_{\varepsilon})\Big{)}^{2}\right], \displaystyle T_{2}^{\varepsilon} \displaystyle:= \displaystyle\sum_{b\in\mathcal{B}}\mathbb{E}\left[\Big{(}\varepsilon^{\frac{d% }{2}-1}\int_{\mathbb{R}^{d}}\phi_{j}(\cdot+e_{k})\nabla_{k}g_{\varepsilon,j}e_% {k}\cdot\Delta_{b}\boldsymbol{a}\nabla(u_{\varepsilon}^{b}(\varepsilon\cdot))% \Big{)}^{2}\right], \displaystyle T_{3}^{\varepsilon} \displaystyle:= \displaystyle\sum_{b\in\mathcal{B}}\mathbb{E}\left[\Big{(}\varepsilon^{\frac{d% }{2}-1}\int_{\mathbb{R}^{d}}\phi_{j}(\cdot+e_{k})\nabla_{k}(g_{\varepsilon,j}% \nabla_{i}\bar{u}_{\varepsilon})e_{k}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_% {i}^{b}+e_{i})\Big{)}^{2}\right], \displaystyle T_{4}^{\varepsilon} \displaystyle:= \displaystyle\sum_{b\in\mathcal{B}}\mathbb{E}\left[\Big{(}\varepsilon^{\frac{d% }{2}-1}\int_{\mathbb{R}^{d}}\nabla r_{\varepsilon}\cdot\Delta_{b}\boldsymbol{a% }\nabla(u_{\varepsilon}^{b}(\varepsilon\cdot))\Big{)}^{2}\right], \displaystyle T_{5}^{\varepsilon} \displaystyle:= \displaystyle\sum_{b\in\mathcal{B}}\mathbb{E}\left[\Big{(}\varepsilon^{\frac{d% }{2}-1}\int_{\mathbb{R}^{d}}\nabla R_{\varepsilon,i}\cdot\Delta_{b}\boldsymbol% {a}(\nabla\phi_{i}^{b}+e_{i})\Big{)}^{2}\right],

with the auxiliary fields r_{\varepsilon} and R_{\varepsilon} defined in (3.8) and in (3.9). The conclusion of Proposition 2.2 is a consequence of the following five estimates: for all 0<p-1\ll 1 and all 0<\alpha-d\ll 1,

 \displaystyle T_{1}^{\varepsilon} \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{2}\mu_{d}(\tfrac{1}{\varepsilon})\,\|g\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|w_{1}^{\alpha\frac{p-1}{4p}}\mu_{d}% (|\cdot|)^{\frac{1}{2}}Df\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}^{2}, (3.24) \displaystyle T_{2}^{\varepsilon} \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{2}\mu_{d}(\tfrac{1}{\varepsilon})\,\|f\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|\mu_{d}(|\cdot|)^{\frac{1}{2}}Dg\|_% {\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}, (3.25) \displaystyle T_{3}^{\varepsilon} \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{2}\mu_{d}(\tfrac{1}{\varepsilon})\Big{(}\|f\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|\mu_{d}(|\cdot|)^{\frac{1}{2}}Dg\|_% {\operatorname{L}^{4}(\mathbb{R}^{d})}^{2} (3.26) \displaystyle                                           +\|g\|_{\operatorname{% L}^{4}(\mathbb{R}^{d})}^{2}\|\mu_{d}(|\cdot|)^{\frac{1}{2}}Df\|_{\operatorname% {L}^{4}(\mathbb{R}^{d})}^{2}\Big{)}, \displaystyle T_{4}^{\varepsilon} \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{2}\mu_{d}(\tfrac{1}{\varepsilon})\,\|f\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|w_{1}^{\alpha\frac{p-1}{4p}}\mu_{d}% (|\cdot|)^{\frac{1}{2}}Dg\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}^{2}, (3.27) \displaystyle T_{5}^{\varepsilon} \displaystyle\lesssim_{\alpha,p} \displaystyle\varepsilon^{2}\mu_{d}(\tfrac{1}{\varepsilon})\Big{(}\|f\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|w_{1}^{\alpha\frac{p-1}{4p}}\mu_{d}% (|\cdot|)^{\frac{1}{2}}Dg\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}^{2} (3.28) \displaystyle                                   +\|g\|_{\operatorname{L}^{4}(% \mathbb{R}^{d})}^{2}\|w_{1}^{\alpha\frac{p-1}{4p}}\mu_{d}(|\cdot|)^{\frac{1}{2% }}Df\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}^{2}\Big{)}.

We split the proof into three steps: we prove the above five estimates in the first two steps, and conclude in the last step by controlling the discretization error.

Step 1. Equation for the two-scale expansion error w_{f,\varepsilon} on \mathbb{R}^{d}:

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}\nabla w_{f,\varepsilon}=\nabla_{l}% ^{*}\Big{(}\sigma_{jkl}(\cdot-e_{k})\nabla_{k}^{*}\nabla_{j}\bar{u}_{% \varepsilon}+\phi_{j}(\cdot+e_{k})\boldsymbol{a}_{lk}\nabla_{k}\nabla_{j}\bar{% u}_{\varepsilon}\Big{)}. (3.29)

This is the discrete counterpart of similar identities in [GNO-reg, GNO-quant].

Using equations (2.1) and (3.5) in the form -\nabla^{*}\cdot\boldsymbol{a}\nabla(u_{\varepsilon}(\varepsilon\cdot))=-% \nabla^{*}\cdot\bar{\boldsymbol{a}}\nabla\bar{u}_{\varepsilon}, and using the following discrete version of the Leibniz rule, for all \chi_{1},\chi_{2}:\mathbb{Z}^{d}\to\mathbb{R},

 \displaystyle\nabla(\chi_{1}\chi_{2})=\chi_{1}\nabla\chi_{2}+e_{l}\chi_{2}(% \cdot+e_{l})\nabla_{l}\chi_{1}, (3.30)

we obtain

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}\nabla w_{f,\varepsilon}=-\nabla^{*% }\cdot\boldsymbol{a}\nabla(u_{\varepsilon}(\varepsilon\cdot)-\bar{u}_{% \varepsilon}-\phi_{j}\nabla_{j}\bar{u}_{\varepsilon})\\ \displaystyle=-\nabla^{*}\cdot\bar{\boldsymbol{a}}\nabla\bar{u}_{\varepsilon}+% \nabla^{*}\cdot\boldsymbol{a}\nabla\bar{u}_{\varepsilon}+\nabla^{*}\cdot(% \boldsymbol{a}\nabla\phi_{j}\nabla_{j}\bar{u}_{\varepsilon})+\nabla^{*}\cdot(% \boldsymbol{a}e_{k}\phi_{j}(\cdot+e_{k})\nabla_{k}\nabla_{j}\bar{u}_{% \varepsilon}).

Rearranging the terms and using the definition (3.4) of \sigma_{j}, this turns into

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}\nabla w_{f,\varepsilon} \displaystyle= \displaystyle\nabla^{*}\cdot\big{(}(\boldsymbol{a}(\nabla\phi_{j}+e_{j})-\bar{% \boldsymbol{a}}e_{j})\nabla_{j}\bar{u}_{\varepsilon}\big{)}+\nabla^{*}\cdot(% \boldsymbol{a}e_{k}\phi_{j}(\cdot+e_{k})\nabla_{k}\nabla_{j}\bar{u}_{% \varepsilon}) \displaystyle= \displaystyle\nabla^{*}\cdot\big{(}(\nabla^{*}\cdot\sigma_{j})\nabla_{j}\bar{u% }_{\varepsilon}\big{)}+\nabla^{*}\cdot(\boldsymbol{a}e_{k}\phi_{j}(\cdot+e_{k}% )\nabla_{k}\nabla_{j}\bar{u}_{\varepsilon}).

Using again the discrete Leibniz rule (3.30) and the skew-symmetry (3.3) of \sigma_{j}, we find

 \displaystyle\nabla^{*}\cdot\big{(}(\nabla^{*}\cdot\sigma_{j})\nabla_{j}\bar{u% }_{\varepsilon}\big{)} \displaystyle= \displaystyle\nabla_{k}^{*}(\nabla_{l}^{*}\sigma_{jkl}\nabla_{j}\bar{u}_{% \varepsilon})=\underbrace{\nabla_{k}^{*}\nabla_{l}^{*}\sigma_{jkl}\nabla_{j}% \bar{u}_{\varepsilon}}_{=0}+\nabla_{l}^{*}\sigma_{jkl}(\cdot-e_{k})\nabla_{k}^% {*}\nabla_{j}\bar{u}_{\varepsilon} \displaystyle= \displaystyle\nabla_{l}^{*}(\sigma_{jkl}(\cdot-e_{k})\nabla_{k}^{*}\nabla_{j}% \bar{u}_{\varepsilon})-\underbrace{\sigma_{jkl}(\cdot-e_{k}-e_{l})\nabla_{k}^{% *}\nabla_{l}^{*}\nabla_{j}\bar{u}_{\varepsilon}}_{=0},

and the conclusion (3.29) follows.

Step 2. Proof of estimates (3.24)–(3.28).

We start with the first term T_{1}^{\varepsilon}. For b\in\mathcal{B} we use the notation b=(z_{b},z_{b}+e_{b}). Since |\Delta_{b}\boldsymbol{a}(x)|\lesssim\mathds{1}_{Q(z_{b})}(x), the Cauchy-Schwarz inequality yields

 \displaystyle T_{1}^{\varepsilon} \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\sum_{b\in\mathcal{B}}\mathbb{E}\left[|\nabla% \phi(z_{b})+\operatorname{Id}|^{2}\Big{(}\int_{Q(z_{b})}|g_{\varepsilon}||% \nabla w^{b}_{f,\varepsilon}+\phi_{i}^{b}\nabla\nabla_{i}\bar{u}_{\varepsilon}% |\Big{)}^{2}\right] \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}|\nabla% \phi(z_{b})+\operatorname{Id}|^{4}\Big{(}\int_{Q(z_{b})}|g_{\varepsilon}|^{2}% \Big{)}^{2}\right]^{\frac{1}{2}} \displaystyle                                           \times\,\mathbb{E}% \left[\sum_{b\in\mathcal{B}}\Big{(}\int_{Q(z_{b})}|\nabla w^{b}_{f,\varepsilon% }+\phi_{i}^{b}\nabla\nabla_{i}\bar{u}_{\varepsilon}|^{2}\Big{)}^{2}\right]^{% \frac{1}{2}},

and hence, using the moment bounds of Lemma 3.3 and the exchangeability of (\boldsymbol{a},\boldsymbol{a}^{b}),

 \displaystyle T_{1}^{\varepsilon}\lesssim\varepsilon^{d-2}\,\|g_{\varepsilon}% \|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\,\mathbb{E}\left[\sum_{b\in% \mathcal{B}}\Big{(}\int_{Q(z_{b})}|\nabla w_{f,\varepsilon}+\phi_{i}\nabla% \nabla_{i}\bar{u}_{\varepsilon}|^{2}\Big{)}^{2}\right]^{\frac{1}{2}}.

We argue as in (3.15): We rewrite the second right-hand side factor as a norm of averages at the scale r_{*}, smuggle in a suitable power of the weight w_{\varepsilon}, and apply Hölder’s inequality, so that for all p>1 and all \alpha>d,

 \displaystyle T_{1}^{\varepsilon}\lesssim_{\alpha,p}\varepsilon^{\frac{d}{2}(1% +\frac{1}{p})-2}\,\|g_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{% 2}\\ \displaystyle\times\,\mathbb{E}\left[\int_{\mathbb{R}^{d}}w_{\varepsilon}(z)^{% \alpha(p-1)}\Big{(}\fint_{B_{*}(z)}|\nabla w_{f,\varepsilon}|^{2}\Big{)}^{2p}% dz+\int_{\mathbb{R}^{d}}w_{\varepsilon}^{\alpha(p-1)}|\phi|^{4p}|\nabla^{2}% \bar{u}_{\varepsilon}|^{4p}\right]^{\frac{1}{2p}}. (3.31)

By large-scale weighted Calderón-Zygmund theory (cf. Lemma 3.4) applied to equation (3.29) for w_{f,\varepsilon}, we deduce for all 0<p-1\ll 1 and all 0<\alpha-d\ll 1,

 \displaystyle T_{1}^{\varepsilon}\lesssim_{\alpha,p}\varepsilon^{\frac{d}{2}(1% +\frac{1}{p})-2}\,\|g_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{% 2}\\ \displaystyle\times\,\mathbb{E}\left[r_{*}(0)^{\alpha(p-1)}\int_{\mathbb{R}^{d% }}w_{\varepsilon}^{\alpha(p-1)}\Big{(}|\sigma|^{4p}+|\phi|^{4p}+\sum_{k=1}^{d}% |\phi(\cdot+e_{k})|^{4p}\Big{)}|\nabla^{2}\bar{u}_{\varepsilon}|^{4p}\right]^{% \frac{1}{2p}}.

By the moment bounds of Lemmas 3.3 and 3.4, this yields

 \displaystyle T_{1}^{\varepsilon}\lesssim_{\alpha,p}\varepsilon^{\frac{d}{2}(1% +\frac{1}{p})-2}\,\|g_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{% 2}\|w_{\varepsilon}^{\alpha\frac{p-1}{4p}}\mu_{d}(|\cdot|)^{\frac{1}{2}}\nabla% ^{2}\bar{u}_{\varepsilon}\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}^{2}.

We then apply the standard weighted Calderón-Zygmund theory to the discrete constant-coefficient equation (3.5) for \bar{u}_{\varepsilon} (cf. Lemma 3.4 with r_{*}=1), note that for all \chi,\zeta\in C^{\infty}_{c}(\mathbb{R}^{d}) and all q<\infty the inequality |\nabla(\zeta(\varepsilon\cdot))|\leq\varepsilon\int_{0}^{1}|D_{k}\zeta(% \varepsilon(\cdot+te_{k}))|dt leads to

 \displaystyle\int_{\mathbb{R}^{d}}\chi|\nabla(\zeta(\varepsilon\cdot))|^{q}% \leq\varepsilon^{q}\int_{\mathbb{R}^{d}}\Big{(}\sup_{B(x)}|\chi|\Big{)}|D\zeta% (\varepsilon x)|^{q}dx\leq\varepsilon^{q-d}\int_{\mathbb{R}^{d}}\Big{(}\sup_{B% (\frac{x}{\varepsilon})}|\chi|\Big{)}|D\zeta(x)|^{q}dx, (3.32)

rescale the integrals, estimate \mu_{d}(|\frac{\cdot}{\varepsilon}|)\leq\mu_{d}(\frac{1}{\varepsilon})\mu_{d}(% |\cdot|), and the conclusion (3.24) follows.

We turn to the second term T_{2}^{\varepsilon}. Since |\Delta_{b}\boldsymbol{a}(x)|\lesssim\mathds{1}_{Q(z_{b})}(x), the Cauchy-Schwarz inequality yields

 \displaystyle T_{2}^{\varepsilon}\lesssim\varepsilon^{d-2}\,\mathbb{E}\left[% \sum_{b\in\mathcal{B}}|\phi(z_{b}+e_{k})|^{2}\Big{(}\int_{Q(z_{b})}|\nabla_{k}% g_{\varepsilon}|^{2}\Big{)}\Big{(}\int_{Q(z_{b})}|\nabla(u_{\varepsilon}^{b}(% \varepsilon\cdot))|^{2}\Big{)}\right].

We bound the second local integral by an integral at the scale r_{*}^{b}, set B_{*}^{b}(z):=B_{r_{*}^{b}(z)}(z), apply the Cauchy-Schwarz inequality, appeal to the moment bounds of Lemmas 3.3 and 3.4, and obtain

 \displaystyle T_{2}^{\varepsilon}\lesssim\varepsilon^{d-2}\,\|\mu_{d}(|\cdot|)% ^{\frac{1}{2}}\nabla g_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^% {2}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}\Big{(}\fint_{B_{*}^{b}(z_{b})\cup Q% (z_{b})}|\nabla(u_{\varepsilon}^{b}(\varepsilon\cdot))|^{2}\Big{)}^{2}\right]^% {\frac{1}{2}},

which, by exchangeability of (\boldsymbol{a},\boldsymbol{a}^{b}), takes the form

 \displaystyle T_{2}^{\varepsilon}\lesssim\varepsilon^{d-2}\,\|\mu_{d}(|\cdot|)% ^{\frac{1}{2}}\nabla g_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^% {2}\,\mathbb{E}\left[\int_{\mathbb{R}^{d}}\Big{(}\fint_{B_{*}(z)}|\nabla(u_{% \varepsilon}(\varepsilon\cdot))|^{2}\Big{)}^{2}dz\right]^{\frac{1}{2}}.

By large-scale (unweighted) Calderón-Zygmund theory (cf. Lemma 3.4) applied to equation (2.1) for u_{\varepsilon} in the form -\nabla^{*}\cdot\boldsymbol{a}\nabla(u_{\varepsilon}(\varepsilon\cdot))=\nabla% ^{*}\cdot(\varepsilon f_{\varepsilon}), we deduce

 \displaystyle T_{2}^{\varepsilon}\lesssim\varepsilon^{d}\,\|\mu_{d}(|\cdot|)^{% \frac{1}{2}}\nabla g_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2% }\,\|f_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2},

and the conclusion (3.25) follows similarly as above.

The proof of (3.26) for T_{3}^{\varepsilon} is more direct. Indeed, using the moment bounds of Lemma 3.3, decomposing \nabla(g_{\varepsilon,i}\nabla u_{\varepsilon})=\nabla g_{\varepsilon,i}% \otimes\nabla u_{\varepsilon}+g_{\varepsilon,i}(\cdot+e_{k})e_{k}\otimes\nabla% _{k}\nabla u_{\varepsilon}, and applying the Cauchy-Schwarz inequality, we find

 \displaystyle T_{3}^{\varepsilon} \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\sum_{k=1}^{d}\mathbb{E}\left[\sum_{b\in\mathcal% {B}}|\phi(z_{b}+e_{k})|^{2}|\nabla\phi^{b}(z_{b})+\operatorname{Id}|^{2}\Big{(% }\int_{Q(z_{b})}|\nabla(g_{\varepsilon}\nabla\bar{u}_{\varepsilon})|\Big{)}^{2% }\right] \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\,\int_{\mathbb{R}^{d}}\mu_{d}(|\cdot|)|\nabla(g% _{\varepsilon}\nabla\bar{u}_{\varepsilon})|^{2} \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\Big{(}\|\nabla\bar{u}_{\varepsilon}\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|\mu_{d}(|\cdot|)\nabla g_{% \varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}+\|g_{\varepsilon}\|_% {\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|\mu_{d}(|\cdot|)\nabla^{2}\bar{u}_% {\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\Big{)},

and the conclusion (3.26) follows as above.

We turn to the fourth term T_{4}^{\varepsilon}: We smuggle in a power \alpha\frac{p-1}{2p} of the weight w_{\varepsilon}, apply Hölder’s inequality with exponents (\frac{2p}{p-1},2p,2), appeal to the moment bounds of Lemma 3.4, use the exchangeability of (\boldsymbol{a},\boldsymbol{a}^{b}), and obtain for all p>1 and all \alpha>d,

 \displaystyle T_{4}^{\varepsilon} \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\,\mathbb{E}\left[\sum_{b\in\mathcal{B}}r_{*}(z_% {b})^{d}r_{*}^{b}(z_{b})^{d}\Big{(}\fint_{B_{*}(z_{b})\cup Q(z_{b})}|\nabla r_% {\varepsilon}|^{2}\Big{)}\Big{(}\fint_{B_{*}^{b}(z_{b})\cup Q(z_{b})}|\nabla(u% _{\varepsilon}^{b}(\varepsilon\cdot))|^{2}\Big{)}dz\right] \displaystyle\lesssim \displaystyle\varepsilon^{\frac{d}{2}(1+\frac{1}{p})-2}\,\mathbb{E}\left[\int_% {\mathbb{R}^{d}}w_{\varepsilon}(z)^{\alpha(p-1)}\Big{(}\fint_{B_{*}(z)}|\nabla r% _{\varepsilon}|^{2}\Big{)}^{2p}dz\right]^{\frac{1}{2p}} \displaystyle                                                            % \times\,\mathbb{E}\left[\int_{\mathbb{R}^{d}}\Big{(}\fint_{B_{*}(z)}|\nabla(u_% {\varepsilon}(\varepsilon\cdot))|^{2}\Big{)}^{2}dz\right]^{\frac{1}{2}}.

By the large-scale weighted Calderón-Zygmund theory (cf. Lemma 3.4) applied to equation (3.8) for r_{\varepsilon} and to equation (2.1) for u_{\varepsilon}, and the moment bounds of Lemma 3.3, we deduce for all 0<p-1\ll 1 and all 0<\alpha-d\ll 1,

 \displaystyle T_{4}^{\varepsilon}\lesssim\varepsilon^{\frac{d}{2}(1+\frac{1}{p% })}\,\|w_{\varepsilon}^{\alpha\frac{p-1}{4p}}\mu_{d}(|\cdot|)^{\frac{1}{2}}% \nabla g_{\varepsilon}\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}^{2}\|f_{% \varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2},

and the conclusion (3.27) follows as before.

Finally, we turn to the last term T_{5}^{\varepsilon}: We use (3.14) in the form |\nabla\phi^{b}(z_{b})+\operatorname{Id}|\lesssim|\nabla\phi(z_{b})+% \operatorname{Id}|, smuggle in a power \alpha\frac{p-1}{2p} of the weight w_{\varepsilon}, apply Hölder’s inequality with exponents (\frac{2p}{p-1},\frac{2p}{p+1}), appeal to the moment bounds of Lemmas 3.3 and 3.4, and therefore obtain for all p>1 and all \alpha>d,

 \displaystyle T_{5}^{\varepsilon} \displaystyle\lesssim \displaystyle\varepsilon^{d-2}\,\mathbb{E}\left[\int_{\mathbb{R}^{d}}|\nabla% \phi(z)+\operatorname{Id}|^{2}\Big{(}\int_{Q_{2}(z)}|\nabla R_{\varepsilon}|^{% 2}\Big{)}dz\right] \displaystyle\lesssim \displaystyle\varepsilon^{\frac{d}{2}(1+\frac{1}{p})-2}\,\mathbb{E}\left[\int_% {\mathbb{R}^{d}}w_{\varepsilon}(z)^{\alpha\frac{p-1}{p+1}}\Big{(}\fint_{B_{*}(% z)}|\nabla R_{\varepsilon}|^{2}\Big{)}^{\frac{2p}{p+1}}dz\right]^{\frac{p+1}{2% p}}.

By the large-scale weighted Calderón-Zygmund theory (cf. Lemma 3.4) applied to equation (3.9) for R_{\varepsilon}, and the moment bounds of Lemma 3.3, we deduce for all 0<p-1\ll 1 and all 0<\alpha-d\ll 1,

 \displaystyle T_{5}^{\varepsilon}\lesssim\varepsilon^{\frac{d}{2}(1+\frac{1}{p% })-2}\,\|w_{\varepsilon}^{\alpha\frac{p-1}{4p}}\mu_{d}(|\cdot|)^{\frac{1}{2}}% \nabla(g_{\varepsilon}\nabla\bar{u}_{\varepsilon})\|_{\operatorname{L}^{\frac{% 4p}{p+1}}(\mathbb{R}^{d})}^{2}.

Decomposing \nabla(g_{\varepsilon,i}\nabla u_{\varepsilon})=\nabla g_{\varepsilon,i}% \otimes\nabla u_{\varepsilon}+g_{\varepsilon,i}(\cdot+e_{k})e_{k}\otimes\nabla% _{k}\nabla u_{\varepsilon} and suitably applying Hölder’s inequality with exponents (\frac{p+1}{p},p+1), the conclusion (3.28) follows as before.

Step 3. Conclusion.

Inserting estimates (3.24)–(3.28) into (3.23) yields for all 0<p-1\ll 1 and all \alpha>d\frac{p-1}{4p},

 \displaystyle\|E^{\varepsilon}_{0}(f,g)\|_{\operatorname{L}^{2}(\Omega)}\,% \lesssim_{\alpha,p}\,\varepsilon\mu_{d}(\tfrac{1}{\varepsilon})^{\frac{1}{2}}% \\ \displaystyle\times\Big{(}\|f\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}\|w_{1}^% {\alpha}Dg\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})}+\|g\|_{\operatorname{L}^{% 4}(\mathbb{R}^{d})}\|w_{1}^{\alpha}Df\|_{\operatorname{L}^{4p}(\mathbb{R}^{d})% }\Big{)}. (3.33)

It remains to deduce the corresponding result for E^{\varepsilon}(f,g), and deal with the discretization error. In terms of \tilde{u}_{\varepsilon}:=\bar{u}_{\varepsilon}(\frac{\cdot}{\varepsilon}) and \tilde{v}_{\varepsilon}:=\bar{v}_{\varepsilon}(\frac{\cdot}{\varepsilon}), equations (3.5) take the form

 \displaystyle-\nabla_{\varepsilon}^{*}\cdot\bar{\boldsymbol{a}}\nabla_{% \varepsilon}\tilde{u}_{\varepsilon}=\nabla_{\varepsilon}^{*}\cdot f,\qquad-% \nabla_{\varepsilon}^{*}\cdot\bar{\boldsymbol{a}}^{*}\nabla_{\varepsilon}% \tilde{v}_{\varepsilon}=\nabla_{\varepsilon}^{*}\cdot g. (3.34)

The definitions of E^{\varepsilon} and E_{0}^{\varepsilon} then lead to the relation

 \displaystyle E^{\varepsilon}(f,g)\,=\,E^{\varepsilon}_{0}(f,g)+J_{0}^{% \varepsilon}\big{(}(\nabla_{\varepsilon}\tilde{u}_{\varepsilon}-D\bar{u})% \otimes g\big{)}, (3.35)

where J_{0}^{\varepsilon}((\nabla_{\varepsilon}\tilde{u}_{\varepsilon}-D\bar{u})% \otimes g) is a discretization error. By Proposition 2.1 (and Cauchy-Schwarz’ inequality), it is enough to establish for all 1<p<\infty and all 0\leq\alpha<d\frac{p-1}{p},

 \displaystyle\|w_{1}^{\alpha}(\nabla_{\varepsilon}\tilde{u}_{\varepsilon}-D% \bar{u})\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}\lesssim_{\alpha,p}% \varepsilon\|w_{1}^{\alpha}Df\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}. (3.36)

For that purpose, we note that \bar{u} is an approximate solution of the discrete equation (3.34). Indeed, integrating equation (2.2) for \bar{u} on a unit cube yields for all x\in\mathbb{R}^{d}

 \displaystyle 0=\int_{[-1,0)^{d}}D\cdot(\bar{\boldsymbol{a}}D\bar{u}+f)(x+% \varepsilon y)dy=\nabla_{\varepsilon}^{*}\cdot(\bar{\boldsymbol{a}}D\bar{u}+f)% (x)+\nabla_{\varepsilon}^{*}\cdot T_{\varepsilon}(x), (3.37)

where the error term T_{\varepsilon} is given by T_{\varepsilon}(x):=e_{i}\int_{S_{i}}((\bar{\boldsymbol{a}}D\bar{u}+f)_{i}(x+% \varepsilon y)-(\bar{\boldsymbol{a}}D\bar{u}+f)_{i}(x))dy in terms of S_{i}:=\{y\in[-1,0]^{d}:y_{i}=0\}, and satisfies for all 1\leq p<\infty and all 0\leq\alpha<\infty,

 \displaystyle\|w_{1}^{\alpha}T_{\varepsilon}\|_{\operatorname{L}^{p}(\mathbb{R% }^{d})}\lesssim_{\alpha}\varepsilon\|w_{1}^{\alpha}D(\bar{\boldsymbol{a}}D\bar% {u}+f)\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}\lesssim\varepsilon\|w_{1}^{% \alpha}Df\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}+\varepsilon\|w_{1}^{\alpha}% D^{2}\bar{u}\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}. (3.38)

Comparing equations (3.34) and (3.37), the difference \bar{u}-\tilde{u}_{\varepsilon} satisfies

 -\nabla_{\varepsilon}^{*}\cdot\bar{\boldsymbol{a}}\nabla_{\varepsilon}(\bar{u}% -\tilde{u}_{\varepsilon})=\nabla_{\varepsilon}^{*}\cdot T_{\varepsilon}-\nabla% _{\varepsilon}^{*}\cdot\bar{\boldsymbol{a}}(\nabla_{\varepsilon}\bar{u}-D\bar{% u}).

Hence, using the standard weighted Calderón-Zygmund theory applied to this discrete constant-coefficient equation, we obtain for all 1<p<\infty and all 0\leq\alpha<d\frac{p-1}{p},

 \displaystyle\|w_{1}^{\alpha}\nabla_{\varepsilon}(\bar{u}-\tilde{u}_{% \varepsilon})\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}\,\lesssim_{\alpha,p}\,% \|w_{1}^{\alpha}T_{\varepsilon}\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}+\|w_{% 1}^{\alpha}(\nabla_{\varepsilon}\bar{u}-D\bar{u})\|_{\operatorname{L}^{p}(% \mathbb{R}^{d})}.

Since the second right-hand side term is bounded by \varepsilon\|w_{1}^{\alpha}D^{2}\bar{u}\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}, estimate (3.38) yields

 \displaystyle\|w_{1}^{\alpha}(\nabla_{\varepsilon}\tilde{u}_{\varepsilon}-D% \bar{u})\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}\,\lesssim_{\alpha,p}\,% \varepsilon\|w_{1}^{\alpha}Df\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}+% \varepsilon\|w_{1}^{\alpha}D^{2}\bar{u}\|_{\operatorname{L}^{p}(\mathbb{R}^{d}% )}.

The claim (3.36) then follows from the standard weighted Calderón-Zygmund theory applied to the constant-coefficient equation (2.2) for \bar{u}.

### 3.5. Proof of Corollary 2.4

We start with the proof of (2.6) for I_{1}^{\varepsilon}. By integration by parts, equations (3.34) and (2.1) for \tilde{v}_{\varepsilon}, \tilde{u}_{\varepsilon}, and u_{\varepsilon} lead to

 \displaystyle\int g\cdot\nabla_{\varepsilon}(u_{\varepsilon}-\tilde{u}_{% \varepsilon})\leavevmode\nobreak\ \lx@stackrel{{\scriptstyle\eqref{eq:def-ueps% -tilde}}}{{=}}\leavevmode\nobreak\ -\int\nabla_{\varepsilon}\tilde{v}_{% \varepsilon}\cdot\bar{\boldsymbol{a}}\nabla_{\varepsilon}(u_{\varepsilon}-% \tilde{u}_{\varepsilon}) \displaystyle\lx@stackrel{{\scriptstyle\eqref{eq:def-ueps-tilde}}}{{=}} \displaystyle-\int\nabla_{\varepsilon}\tilde{v}_{\varepsilon}\cdot f-\int% \nabla_{\varepsilon}\tilde{v}_{\varepsilon}\cdot\bar{\boldsymbol{a}}\nabla_{% \varepsilon}u_{\varepsilon} \displaystyle\lx@stackrel{{\scriptstyle\eqref{e.def-ueps}}}{{=}} \displaystyle\int\nabla_{\varepsilon}\tilde{v}_{\varepsilon}\cdot(\boldsymbol{% a}_{\varepsilon}\nabla_{\varepsilon}u_{\varepsilon}-\bar{\boldsymbol{a}}\nabla% _{\varepsilon}u_{\varepsilon}).

Subtracting the expectation of both sides yields a discrete version of identity (1.7). In terms of J_{0}^{\varepsilon}, I_{1}^{\varepsilon}, and E_{0}^{\varepsilon} (cf. Section 2.1 and (3.6)), this takes on the following guise,

 \displaystyle I_{1}^{\varepsilon}(f,g)-J_{0}^{\varepsilon}(D\bar{u}\otimes D% \bar{v})=J_{0}^{\varepsilon}(\nabla_{\varepsilon}\tilde{u}_{\varepsilon}% \otimes\nabla_{\varepsilon}\tilde{v}_{\varepsilon}-D\bar{u}\otimes D\bar{v})+E% _{0}^{\varepsilon}(f,\nabla_{\varepsilon}\tilde{v}_{\varepsilon}). (3.39)

Using (3.36) and the standard weighted Calderón-Zygmund theory applied to the constant-coefficient equations (2.2) and (3.34), the conclusion (2.6) for I_{1}^{\varepsilon} follows from (3.33) together with Proposition 2.1.

We turn to the proof of (2.6) for I_{2}^{\varepsilon}. By definition of J_{0}^{\varepsilon}, I_{1}^{\varepsilon}, I_{2}^{\varepsilon}, and E^{\varepsilon} (cf. Section 2.1), we find

 \displaystyle I_{2}^{\varepsilon}(f,g)=E^{\varepsilon}(f,g)+I_{1}^{\varepsilon% }(f,\bar{\boldsymbol{a}}^{*}g)+J_{0}^{\varepsilon}(D\bar{u}\otimes g).

Inserting identities (3.35) and (3.39) (with g replaced by \bar{\boldsymbol{a}}^{*}g and thus \bar{v} replaced by the solution \bar{v}^{\circ} of -D\cdot\bar{\boldsymbol{a}}^{*}D\bar{v}^{\circ}=D\cdot\bar{\boldsymbol{a}}^{*}g, so that \bar{\mathcal{P}}_{L}^{*}g=D\bar{v}^{\circ}+g), the conclusion (2.6) follows similarly as for I_{1}^{\varepsilon}.

We now turn to the proof of (2.7). Let \mathcal{S}(\mathbb{R}^{d}) denote the Schwartz space of rapidly decaying functions, and consider the subspace \mathcal{K}_{\varepsilon}:=\{g\in\mathcal{S}(\mathbb{R}^{d})^{d}:\bar{v}_{% \varepsilon}\in\mathcal{S}(\mathbb{R}^{d})\}, cf. (3.5). Given some fixed \chi\in C^{\infty}_{c}(\mathbb{R}^{d}), set \chi_{L}:=\chi(L\cdot) for L\geq 1. For g\in\mathcal{K}_{\varepsilon}, we compute by integration by parts, using equation (2.4) for \phi_{j} and equation (3.5) for \bar{v}_{\varepsilon}, together with the discrete Leibniz rule (3.30),

 \displaystyle\int_{\mathbb{R}^{d}}\chi_{L}\nabla\bar{v}_{\varepsilon}\cdot\Xi_% {i}\,=\,\int_{\mathbb{R}^{d}}\chi_{L}\nabla\bar{v}_{\varepsilon}\cdot\big{(}% \boldsymbol{a}(\nabla\phi_{i}+e_{i})-\bar{\boldsymbol{a}}(\nabla\phi_{i}+e_{i}% )\big{)} \displaystyle     \lx@stackrel{{\scriptstyle\eqref{e.corr}}}{{=}} \displaystyle-\int_{\mathbb{R}^{d}}\nabla(\bar{v}_{\varepsilon}\chi_{L})\cdot% \bar{\boldsymbol{a}}\nabla\phi_{i}-\int_{\mathbb{R}^{d}}\bar{v}_{\varepsilon}(% \cdot+e_{j})\nabla_{j}\chi_{L}\,\Xi_{ij} \displaystyle     \lx@stackrel{{\scriptstyle\eqref{e.def-utildeeps}}}{{=}} \displaystyle\varepsilon\int_{\mathbb{R}^{d}}\chi_{L}g_{\varepsilon}\cdot% \nabla\phi_{i}+\varepsilon\int_{\mathbb{R}^{d}}\phi_{i}(\cdot+e_{j})g_{% \varepsilon,j}\nabla_{j}\chi_{L}+\int_{\mathbb{R}^{d}}\phi_{i}(\cdot+e_{j})% \nabla_{j}\chi_{L}e_{j}\cdot\bar{\boldsymbol{a}}\nabla\bar{v}_{\varepsilon} \displaystyle             -\int_{\mathbb{R}^{d}}\bar{v}_{\varepsilon}(\cdot+e_% {j})\nabla_{j}\chi_{L}e_{j}\cdot\bar{\boldsymbol{a}}\nabla\phi_{i}-\int_{% \mathbb{R}^{d}}\bar{v}_{\varepsilon}(\cdot+e_{j})\nabla_{j}\chi_{L}\,\Xi_{ij}.

For fixed \varepsilon and g\in\mathcal{K}_{\varepsilon}, using the moment bounds of Lemma 3.3 and the rapid decay at infinity of g and \bar{v}_{\varepsilon}, we may pass to the limit L\uparrow\infty in both sides in \operatorname{L}^{2}(\Omega), and we deduce almost surely

 \displaystyle\int_{\mathbb{R}^{d}}\nabla\bar{v}_{\varepsilon}:\Xi_{j}\,=\,% \varepsilon\int_{\mathbb{R}^{d}}g_{\varepsilon}\cdot\nabla\phi_{j},

that is, after rescaling,

 \displaystyle J_{1}^{\varepsilon}(e_{j}\otimes g)=J_{0}^{\varepsilon}(e_{j}% \otimes\nabla_{\varepsilon}\tilde{v}_{\varepsilon}). (3.40)

We now argue that for all \varepsilon>0 this almost-sure identity can be extended to hold in \operatorname{L}^{2}(\Omega) for all g\in C_{c}^{\infty}(\mathbb{R}^{d})^{d}. First note that Proposition 2.1 combined with the standard weighted Calderón-Zygmund theory for the constant-coefficient equation (3.34) yields for all 0<p-1\ll 1 and all d\frac{p-1}{4p}<\alpha<d\frac{2p-1}{4p},

 \mathbb{E}\left[|J_{0}^{\varepsilon}(e_{j}\otimes\nabla_{\varepsilon}\tilde{v}% _{\varepsilon})|^{2}\right]^{\frac{1}{2}}\,\lesssim_{\alpha,p}\,\|w_{1}^{2% \alpha}\nabla_{\varepsilon}\tilde{v}_{\varepsilon}\|_{\operatorname{L}^{2p}(% \mathbb{R}^{d})}\,\lesssim_{\alpha,p}\,\|w_{1}^{2\alpha}g\|_{\operatorname{L}^% {2p}(\mathbb{R}^{d})},

 \mathbb{E}\left[|J_{1}^{\varepsilon}(e_{j}\otimes g)|^{2}\right]^{\frac{1}{2}}% \,\lesssim\,\|w_{1}^{2\alpha}g\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})}.

Hence, it suffices to check the following density result: for all test functions g\in C_{c}^{\infty}(\mathbb{R}^{d})^{d} there exist a sequence (g_{n})_{n} of elements of \mathcal{K}_{\varepsilon} such that \|w_{1}^{2\alpha}(g_{n}-g)\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})}\to 0 holds for some 0<p-1\ll 1 and some \alpha>d\frac{p-1}{4p}. Let g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d} be fixed. Up to a convolution argument on large scales, we may already assume that the Fourier transform \hat{g} has compact support, say contained in B_{R}. Since the (continuum) Fourier symbol of the discrete Helmholtz projection \nabla_{\varepsilon}(\nabla_{\varepsilon}^{*}\cdot\bar{a}\nabla_{\varepsilon})% ^{-1}\nabla_{\varepsilon}^{*}\cdot is bounded and smooth outside of the dual lattice (\frac{2\pi}{\varepsilon}\mathbb{Z})^{d}, a function g_{n}\in\mathcal{S}(\mathbb{R}^{d})^{d} actually belongs to {\mathcal{K}}_{\varepsilon} whenever its Fourier transform \hat{g}_{n} vanishes in a neighborhood of (\frac{2\pi}{\varepsilon}\mathbb{Z})^{d}. Choosing \chi\in C^{\infty}_{c}(\mathbb{R}^{d}) with \chi=1 in a neighborhood of 0, and defining

 \chi_{n}\,:=\,1-\sum_{z\in(\frac{2\pi}{\varepsilon}\mathbb{Z})^{d}}\chi(n(% \cdot-z)),

the function g_{n}\in\mathcal{S}(\mathbb{R}^{d})^{d} defined by \hat{g}_{n}:=\chi_{n}\hat{g} thus belongs to \mathcal{K}_{\varepsilon}. For p\geq 1, setting q:=\frac{2p-1}{2p}, since \hat{g} is compactly supported in B_{R}, the Hausdorff-Young inequality leads to

 \displaystyle\|w_{1}^{2\alpha}(g_{n}-g)\|_{\operatorname{L}^{2p}(\mathbb{R}^{d% })}\,\leq\,\|(\chi_{n}-1)\hat{g}\|_{W^{2\alpha,q}(\mathbb{R}^{d})}\,\lesssim_{% \alpha}\,\|\chi_{n}-1\|_{W^{2\alpha,q}(B_{R})}\|\hat{g}\|_{W^{2\alpha,\infty}(% \mathbb{R}^{d})}\\ \displaystyle\lesssim\,\|\chi_{n}-1\|_{W^{2\alpha,q}(B_{R})}\|w_{1}^{2\alpha}g% \|_{\operatorname{L}^{1}(\mathbb{R}^{d})}.

For 2\alpha<\frac{d}{q}=d\frac{2p-1}{2p}, reflecting the fact that the Sobolev space W^{2\alpha,q}(\mathbb{R}^{d}) fails to embed into the space of continuous functions, there holds \chi_{n}\to 1 in W^{2\alpha,q}_{\operatorname{loc}}(\mathbb{R}^{d}) as n\uparrow\infty, and hence \|w_{1}^{2\alpha}(g_{n}-g)\|_{\operatorname{L}^{2p}(\mathbb{R}^{d})}\to 0. This establishes the claimed density result, and we conclude that identity (3.40) can be extended in \operatorname{L}^{2}(\Omega) to all g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}. The estimate (2.7) for J_{1}^{\varepsilon} then follows from the discretization error estimate (3.36) together with Proposition 2.1. The estimate (2.7) for J_{2}^{\varepsilon} is obtained in a similar way.

## 4. Asymptotic normality

We turn to the normal approximation result for the homogenization commutator \Xi as stated in Proposition 2.7.

### 4.1. Structure of the proof and auxiliary results

The main tool to prove Proposition 2.7 is the following suitable form of a second-order Poincaré inequality à la Chatterjee [C1]. Based on Stein’s method, it can be shown to hold for any product measure \mathbb{P} on \Omega. (The proof follows from [C1, Theorem 2.2] and from [LRP-15, Theorem 4.2] in the case of the Wasserstein and of the Kolmogorov metric, respectively, combined with the spectral gap estimate of Lemma 3.1.) Let us first fix some more notation. Let X=X(\boldsymbol{a}) be a random variable on \Omega, that is, a measurable function of (a(b))_{b\in\mathcal{B}}. For all E\subset\mathcal{B} we denote by \boldsymbol{a}^{E} the random field that coincides with \boldsymbol{a} on all edges b\notin E and with the iid copy \boldsymbol{a}^{\prime} on all edges b\in E. In particular, \boldsymbol{a} and \boldsymbol{a}^{E} always have the same law. We use the abbreviation X^{E}:=X(\boldsymbol{a}^{E}) and define \Delta_{b}X^{E}:=X^{E}-X^{E\cup\{b\}}. As before, we write for simplicity X^{b}:=X^{\{b\}}, and similarly X^{b,b^{\prime}}:=X^{\{b,b^{\prime}\}}. In particular, \Delta_{b}\Delta_{b^{\prime}}X=X-X^{b}-X^{b^{\prime}}+X^{b,b^{\prime}}.

###### Lemma 4.1 ([C1, LRP-15]).

Let \mathbb{P} be a product measure and let \delta_{\mathcal{N}} be defined in (1.14). For all X=X(\boldsymbol{a})\in\operatorname{L}^{2}(\Omega), we have

 \displaystyle\delta_{\mathcal{N}}(X)\,\lesssim\,\frac{1}{\mathrm{Var}\left[X% \right]^{\frac{3}{2}}}\sum_{b\in\mathcal{B}_{L}}\mathbb{E}\left[|\Delta_{b}X|^% {6}\right]^{\frac{1}{2}} \displaystyle                                   +\frac{1}{\mathrm{Var}\left[X% \right]}\bigg{(}\sum_{b\in\mathcal{B}}\Big{(}\sum_{e^{\prime}\in\mathcal{B}}% \mathbb{E}\left[|\Delta_{b^{\prime}}X|^{4}\right]^{\frac{1}{4}}\mathbb{E}\left% [|\Delta_{b}\Delta_{b^{\prime}}X|^{4}\right]^{\frac{1}{4}}\Big{)}^{2}\bigg{)}^% {\frac{1}{2}}.\qed

In addition, we make crucial use of the following optimal annealed estimate on the mixed gradient of the Green’s function, first proved by Marahrens and the third author [MaO].

###### Lemma 4.2 ([MaO]).

Let d\geq 2 and let \mathbb{P} be a product measure. For all y\in\mathbb{Z}^{d} there exists a function \nabla G(\cdot,y) that is the unique decaying solution in \mathbb{Z}^{d} of

 -\nabla^{*}\cdot\boldsymbol{a}\nabla G(\cdot,y)=\delta(\cdot-y).

It satisfies the following moment bound: for all q<\infty and all x,y\in\mathbb{Z}^{d},

 \displaystyle\mathbb{E}\left[|\nabla\nabla G(x,y)|^{q}\right]^{\frac{1}{q}}\,% \lesssim_{q}\,(1+|x-y|)^{-d},

where \nabla\nabla denotes the mixed second gradient. ∎

### 4.2. Proof of Proposition 2.7

Let F\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d}, and set F_{\varepsilon}:=F(\varepsilon\cdot). By Lemma 4.1, it is enough to estimate the following two contributions,

 \displaystyle K_{1}^{\varepsilon}:=\,\sum_{b\in\mathcal{B}}\mathbb{E}\left[|% \Delta_{b}I_{0}^{\varepsilon}(F)|^{6}\right]^{\frac{1}{2}},\quad K_{2}^{% \varepsilon}:=\sum_{b\in\mathcal{B}}\left(\sum_{b^{\prime}\in\mathcal{B}}% \mathbb{E}\left[|\Delta_{b^{\prime}}I_{0}^{\varepsilon}(F)|^{4}\right]^{\frac{% 1}{4}}\mathbb{E}\left[|\Delta_{b}\Delta_{b^{\prime}}I_{0}^{\varepsilon}(F)|^{4% }\right]^{\frac{1}{4}}\right)^{2}.

We split the proof into three steps: we start with an auxiliary estimate, and then estimate K_{1}^{\varepsilon} and K_{2}^{\varepsilon} separately.

Step 1. Auxiliary estimate: for all \zeta\in C^{\infty}_{c}(\mathbb{R}^{d}), all 1\leq p<\infty, and all r\geq 0,

 \displaystyle\int_{\mathbb{R}^{d}}\log^{r}(2+|z|)\bigg{(}\int_{\mathbb{R}^{d}}% \frac{|\zeta(x)|}{(1+|x-z|)^{d}}dx\bigg{)}^{p}dz\leavevmode\nobreak\ \lesssim_% {p,r}\leavevmode\nobreak\ \int_{\mathbb{R}^{d}}\log^{p+r}(2+|x|)\,|\zeta(x)|^{% p}\,dx. (4.1)

Let \alpha>0 be fixed. Smuggling in a power \alpha\frac{p-1}{p} of the weight 1+|x|, and applying Hölder’s inequality with exponent p, we find

 \displaystyle\int_{\mathbb{R}^{d}}\log^{r}(2+|z|)\bigg{(}\int_{\mathbb{R}^{d}}% \frac{|\zeta(x)|}{(1+|x-z|)^{d}}dx\bigg{)}^{p}dz\\ \displaystyle\leq\int_{\mathbb{R}^{d}}\log^{r}(2+|z|)\bigg{(}\int_{\mathbb{R}^% {d}}\frac{(1+|x|)^{\alpha(p-1)}|\zeta(x)|^{p}}{(1+|x-z|)^{d}}dx\bigg{)}\bigg{(% }\int_{\mathbb{R}^{d}}\frac{dx}{(1+|x-z|)^{d}(1+|x|)^{\alpha}}\bigg{)}^{p-1}dz.

The last integral is controlled by C(d,\alpha)\frac{\log(2+|z|)}{(1+|z|)^{\alpha}}, hence by Fubini’s theorem

 \displaystyle\int_{\mathbb{R}^{d}}\log^{r}(2+|z|)\bigg{(}\int_{\mathbb{R}^{d}}% \frac{|\zeta(x)|}{(1+|x-z|)^{d}}dx\bigg{)}^{p}dz \displaystyle\lesssim_{\alpha,p} \displaystyle\int_{\mathbb{R}^{d}}\frac{\log^{p+r-1}(2+|z|)}{(1+|z|)^{\alpha(p% -1)}}\bigg{(}\int_{\mathbb{R}^{d}}\frac{(1+|x|)^{\alpha(p-1)}|\zeta(x)|^{p}}{(% 1+|x-z|)^{d}}dx\bigg{)}dz \displaystyle= \displaystyle\int_{\mathbb{R}^{d}}(1+|x|)^{\alpha(p-1)}|\zeta(x)|^{p}\bigg{(}% \int_{\mathbb{R}^{d}}\frac{\log^{p+r-1}(2+|z|)}{(1+|x-z|)^{d}(1+|z|)^{\alpha(p% -1)}}dz\bigg{)}dx.

Since the last integral is controlled by C(d,p,r,\alpha)\frac{\log^{p+r}(2+|x|)}{(1+|x|)^{\alpha(p-1)}}, the conclusion (4.1) follows.

Step 2. Proof of

 K_{1}^{\varepsilon}\,\lesssim\,\varepsilon^{\frac{d}{2}}\big{(}\|F\|_{% \operatorname{L}^{3}(\mathbb{R}^{d})}^{3}+\|\log(2+|\cdot|)\,\mu_{d}(|\cdot|)^% {\frac{1}{2}}DF\|_{\operatorname{L}^{3}(\mathbb{R}^{d})}^{3}\big{)}.

After integration by parts, the representation formula for the vertical derivative \Delta_{b}\Xi in (3.18) leads to

 \displaystyle\Delta_{b}J_{0}^{\varepsilon}(F)=\varepsilon^{\frac{d}{2}}\int_{% \mathbb{R}^{d}}F_{\varepsilon,ij}(\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b}% \boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i})\\ \displaystyle+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(\cdot% +e_{k})\nabla_{k}F_{\varepsilon,ij}e_{k}\cdot\Delta_{b}\boldsymbol{a}(\nabla% \phi_{i}^{b}+e_{i})\\ \displaystyle+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\big{(}\phi_{j}^{*% }(\cdot+e_{k})\boldsymbol{a}_{kl}\nabla_{k}F_{\varepsilon,ij}+\sigma_{jkl}^{*}% (\cdot-e_{k})\nabla_{k}^{*}F_{\varepsilon,ij}\big{)}\nabla_{l}\Delta_{b}\phi_{% i}.

For b=(z_{b},z_{b}+e_{b}), the Green representation formula applied to equation (3.12) takes the following form, for all x\in\mathbb{Z}^{d},

 \displaystyle\nabla\Delta_{b}\phi_{i}(x)\,=\,-\nabla\nabla G(x,z_{b})\Delta_{b% }\boldsymbol{a}(z_{b})(\nabla\phi_{i}^{b}(z_{b})+e_{i}). (4.2)

Inserted into the above representation formula for \Delta_{b}J_{0}^{\varepsilon}(F), and combined with |\Delta_{b}\boldsymbol{a}(x)|\lesssim\mathds{1}_{Q(z_{b})}(x) and the moment bounds of Lemmas 3.3 and 4.2, it yields for all q<\infty,

 \displaystyle\mathbb{E}\left[|\Delta_{b}I_{0}^{\varepsilon}(F)|^{q}\right]^{% \frac{1}{q}}\,\lesssim_{q}\,\varepsilon^{\frac{d}{2}}\int_{Q(z_{b})}|F_{% \varepsilon}|+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\frac{\mu_{d}(|x|)% ^{\frac{1}{2}}|\nabla F_{\varepsilon}(x)|}{(1+|x-z_{b}|)^{d}}dx. (4.3)

Summing the cube of this estimate over b\in\mathcal{B} for q=6, and using (4.1) with p=3, r=0, and \zeta=\mu_{d}(|\cdot|)^{\frac{1}{2}}\nabla F_{\varepsilon}, we obtain

 \displaystyle K_{1}^{\varepsilon} \displaystyle\lesssim \displaystyle\varepsilon^{\frac{3d}{2}}\int_{\mathbb{R}^{d}}|F_{\varepsilon}|^% {3}+\varepsilon^{\frac{3d}{2}}\int_{\mathbb{R}^{d}}\bigg{(}\int_{\mathbb{R}^{d% }}\frac{\mu_{d}(|x|)^{\frac{1}{2}}|\nabla F_{\varepsilon}(x)|}{(1+|x-z|^{d})}% dx\bigg{)}^{3}dz \displaystyle\lesssim \displaystyle\varepsilon^{\frac{3d}{2}}\int_{\mathbb{R}^{d}}|F_{\varepsilon}|^% {3}+\varepsilon^{\frac{3d}{2}}\int_{\mathbb{R}^{d}}\log^{3}(2+|\cdot|)\,\mu_{d% }(|\cdot|)^{\frac{3}{2}}|\nabla F_{\varepsilon}|^{3}.

Rescaling the integrals and using (3.32), the conclusion follows.

Step 3. Proof of

 \displaystyle K_{2}^{\varepsilon}\lesssim\varepsilon^{d}\log^{2}(2+\tfrac{1}{% \varepsilon})\Big{(}\|F\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{4}+\|\log(2+% |\cdot|)F\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{4}\Big{)}\\ \displaystyle+\varepsilon^{d+2}\log^{4}(2+\tfrac{1}{\varepsilon})\mu_{d}(% \tfrac{1}{\varepsilon})\Big{(}\|\log(2+|\cdot|)F\|_{\operatorname{L}^{4}(% \mathbb{R}^{d})}^{4}+\|\log^{2}(2+|\cdot|)\mu_{d}(|\cdot|)^{\frac{1}{2}}DF\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{4}\Big{)}\\ \displaystyle+\varepsilon^{d+4}\log^{6}(2+\tfrac{1}{\varepsilon})\mu_{d}(% \tfrac{1}{\varepsilon})^{2}\|\log^{2}(2+|\cdot|)\mu_{d}(|\cdot|)^{\frac{1}{2}}% DF\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{4}.

We need to iterate the vertical derivative and estimate \Delta_{b}\Delta_{b^{\prime}}I_{0}^{\varepsilon}(F). By definition of the homogenization commutator, we find

 \displaystyle\Delta_{b}\Delta_{b^{\prime}}J_{0}^{\varepsilon}(F)=\varepsilon^{% \frac{d}{2}}\Delta_{b}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}e_{j}\cdot\Delta_% {b^{\prime}}\boldsymbol{a}(\nabla\phi_{i}^{b^{\prime}}+e_{i})+\varepsilon^{% \frac{d}{2}}\Delta_{b}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}e_{j}\cdot(% \boldsymbol{a}-\bar{\boldsymbol{a}})\nabla\Delta_{b^{\prime}}\phi_{i} (4.4) \displaystyle= \displaystyle\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}e% _{j}\cdot\Delta_{b}\Delta_{b^{\prime}}\boldsymbol{a}(\nabla\phi^{b,b^{\prime}}% _{i}+e_{i})+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}e_% {j}\cdot\big{(}\Delta_{b}\boldsymbol{a}\nabla\Delta_{b^{\prime}}\phi_{i}^{b}+% \Delta_{b^{\prime}}\boldsymbol{a}\nabla\Delta_{b}\phi_{i}^{b^{\prime}}\big{)} \displaystyle                                                    +\varepsilon^% {\frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}e_{j}\cdot(\boldsymbol{a}-% \bar{\boldsymbol{a}})\nabla\Delta_{b}\Delta_{b^{\prime}}\phi_{i}.

In order to avoid additional logarithmic factors, we need to suitably rewrite the last right-hand side term, and we argue similarly as in the proof of (3.18). Using the definition (3.4) of \sigma_{j}^{*} in the form (\boldsymbol{a}^{*}-\bar{\boldsymbol{a}}^{*})e_{j}=-\boldsymbol{a}^{*}\nabla% \phi_{j}^{*}+\nabla^{*}\cdot\sigma_{j}^{*}, applying the discrete Leibniz rule (3.20), and using the skew-symmetry (3.3) of \sigma_{i}, we obtain

 \displaystyle e_{j}\cdot(\boldsymbol{a}-\bar{\boldsymbol{a}})\nabla\Delta_{b}% \Delta_{b^{\prime}}\phi_{i}=(\nabla^{*}\cdot\sigma_{j}^{*})\cdot\nabla\Delta_{% b}\Delta_{b^{\prime}}\phi_{i}-\nabla\phi_{j}^{*}\cdot\boldsymbol{a}\nabla% \Delta_{b}\Delta_{b^{\prime}}\phi_{i}\\ \displaystyle=-\nabla_{k}\big{(}\sigma_{jkl}^{*}(\cdot-e_{k})\nabla_{l}\Delta_% {b}\Delta_{b^{\prime}}\phi_{i}\big{)}-\nabla_{k}^{*}\big{(}\phi_{j}^{*}(\cdot+% e_{k})e_{k}\cdot\boldsymbol{a}\nabla\Delta_{b}\Delta_{b^{\prime}}\phi_{i}\big{% )}+\phi_{j}^{*}\nabla^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}\Delta_{b^{\prime}% }\phi_{i}.

The vertical derivative of equation (3.12) takes the form

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}\nabla\Delta_{b}\Delta_{b^{\prime}}% \phi_{i}=\nabla^{*}\cdot\Delta_{b^{\prime}}\boldsymbol{a}\nabla\Delta_{b}\phi_% {i}^{b^{\prime}}+\nabla^{*}\cdot\Delta_{b}\Delta_{b^{\prime}}\boldsymbol{a}(% \nabla\phi_{i}^{b,b^{\prime}}+e_{i})+\nabla^{*}\cdot\Delta_{b}\boldsymbol{a}% \nabla\Delta_{b^{\prime}}\phi_{i}^{b}, (4.5)

which, combined with the above, yields

 \displaystyle e_{j}\cdot(\boldsymbol{a}-\bar{\boldsymbol{a}})\nabla\Delta_{b}% \Delta_{b^{\prime}}\phi_{i}=-\nabla_{k}\big{(}\sigma_{jkl}^{*}(\cdot-e_{k})% \nabla_{l}\Delta_{b}\Delta_{b^{\prime}}\phi_{i}\big{)}-\nabla_{k}^{*}\big{(}% \phi_{j}^{*}(\cdot+e_{k})e_{k}\cdot\boldsymbol{a}\nabla\Delta_{b}\Delta_{b^{% \prime}}\phi_{i}\big{)}\\ \displaystyle-\phi_{j}^{*}\nabla^{*}\cdot\Delta_{b^{\prime}}\boldsymbol{a}% \nabla\Delta_{b}\phi_{i}^{b^{\prime}}-\phi_{j}^{*}\nabla^{*}\cdot\Delta_{b}% \Delta_{b^{\prime}}\boldsymbol{a}(\nabla\phi_{i}^{b,b^{\prime}}+e_{i})-\phi_{j% }^{*}\nabla^{*}\cdot\Delta_{b}\boldsymbol{a}\nabla\Delta_{b^{\prime}}\phi_{i}^% {b}.

Inserting this into (4.4), integrating by parts, and applying the discrete Leibniz rule (3.30), we obtain the following representation formula

 \displaystyle\Delta_{b}\Delta_{b^{\prime}}J_{0}^{\varepsilon}(F)=\varepsilon^{% \frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}(\nabla\phi_{j}^{*}+e_{j})% \cdot\Delta_{b}\Delta_{b^{\prime}}\boldsymbol{a}(\nabla\phi^{b,b^{\prime}}_{i}% +e_{i})\\ \displaystyle+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}% (\nabla\phi_{j}^{*}+e_{j})\cdot\big{(}\Delta_{b}\boldsymbol{a}\nabla\Delta_{b^% {\prime}}\phi_{i}^{b}+\Delta_{b^{\prime}}\boldsymbol{a}\nabla\Delta_{b}\phi_{i% }^{b^{\prime}}\big{)}\\ \displaystyle+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\big{(}\sigma_{jkl% }^{*}(\cdot-e_{k})\nabla_{k}^{*}F_{\varepsilon,ij}+\phi_{j}^{*}(\cdot+e_{k})% \boldsymbol{a}_{kl}\nabla_{k}F_{\varepsilon,ij}\big{)}\nabla_{l}\Delta_{b}% \Delta_{b^{\prime}}\phi_{i}\\ \displaystyle+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(\cdot% +e_{k})\nabla_{k}F_{\varepsilon,ij}e_{k}\cdot\big{(}\Delta_{b}\boldsymbol{a}% \nabla\Delta_{b^{\prime}}\phi_{i}^{b}+\Delta_{b^{\prime}}\boldsymbol{a}\nabla% \Delta_{b}\phi_{i}^{b^{\prime}}\big{)}\\ \displaystyle+\varepsilon^{\frac{d}{2}}\int_{\mathbb{R}^{d}}\phi_{j}^{*}(\cdot% +e_{k})\nabla_{k}F_{\varepsilon,ij}e_{k}\cdot\Delta_{b}\Delta_{b^{\prime}}% \boldsymbol{a}(\nabla\phi_{i}^{b,b^{\prime}}+e_{i}). (4.6)

We need to estimate the moment of each right-hand side term. Fix momentarily b=(z_{b},z_{b}+e_{b}) and b^{\prime}=(z_{b^{\prime}},z_{b^{\prime}}+e_{b^{\prime}}). Applying Lemmas 3.3 and 4.2 to the Green representation formula (4.2) for \nabla\Delta_{b}\phi, we find for all q<\infty,

 \mathbb{E}\left[|\nabla\Delta_{b}\phi(x)|^{q}\right]^{\frac{1}{q}}\lesssim_{q}% (1+|x-z_{b}|)^{-d}.

We then turn to the second vertical derivatives. We obviously have \Delta_{b}\Delta_{b^{\prime}}\boldsymbol{a}=\mathds{1}_{b=b^{\prime}}\Delta_{b% }\boldsymbol{a}. Next, the Green representation formula applied to equation (4.5) yields

 \displaystyle\nabla\Delta_{b}\Delta_{b^{\prime}}\phi_{j}(x)=-\nabla\nabla G(x,% z_{b^{\prime}})\cdot\Delta_{b^{\prime}}\boldsymbol{a}(z_{b^{\prime}})\nabla% \Delta_{b}\phi_{j}^{b^{\prime}}(z_{b^{\prime}})-\nabla\nabla G(x,z_{b})\cdot% \Delta_{b}\boldsymbol{a}(z_{b})\nabla\Delta_{b^{\prime}}\phi_{j}^{b}(z_{b})\\ \displaystyle-\mathds{1}_{b=b^{\prime}}\nabla\nabla G(x,z_{b})\cdot\Delta_{b}% \boldsymbol{a}(z_{b})(\nabla\phi_{j}^{b}+e_{j}),

so that, for all q<\infty, Lemmas 3.3 and 4.2 lead to

 \displaystyle\mathbb{E}\left[|\nabla\Delta_{b}\Delta_{b^{\prime}}\phi(x)|^{q}% \right]^{\frac{1}{q}}\,\lesssim_{q}\,(1+|z_{b}-z_{b^{\prime}}|)^{-d}\big{(}(1+% |x-z_{b^{\prime}}|)^{-d}+(1+|x-z_{b}|)^{-d}\big{)}.

Inserting these estimates into (4.6), we obtain

 \displaystyle\mathbb{E}\left[|\Delta_{b}\Delta_{b^{\prime}}J_{0}^{\varepsilon}% (F)|^{4}\right]^{\frac{1}{4}}\lesssim\frac{\varepsilon^{\frac{d}{2}}}{(1+|z_{b% }-z_{b^{\prime}}|)^{d}}\bigg{(}\int_{Q(z_{b})}|F_{\varepsilon}|+\int_{Q(z_{b^{% \prime}})}|F_{\varepsilon}|\\ \displaystyle+\int_{\mathbb{R}^{d}}\frac{\mu_{d}(|x|)^{\frac{1}{2}}|\nabla F_{% \varepsilon}(x)|}{(1+|x-z_{b^{\prime}}|)^{d}}dx+\int_{\mathbb{R}^{d}}\frac{\mu% _{d}(|x|)^{\frac{1}{2}}|\nabla F_{\varepsilon}(x)|}{(1+|x-z_{b}|)^{d}}dx\bigg{% )}.

Combining this with (4.3) and with the definition of K_{2}^{\varepsilon}, with the short-hand notation

 I(\zeta)(z):=\int_{\mathbb{R}^{d}}\frac{|\zeta(x)|}{(1+|x-z|)^{d}}dx,

and G_{\varepsilon}:=\mu_{d}(|\cdot|)^{\frac{1}{2}}\nabla F_{\varepsilon}, we deduce

 \displaystyle K_{2}^{\varepsilon}\lesssim\varepsilon^{2d}\int_{\mathbb{R}^{d}}% \Big{(}|F_{\varepsilon}|^{2}|I(F_{\varepsilon})|^{2}+|I(|F_{\varepsilon}|^{2})% |^{2}+|I(F_{\varepsilon})|^{2}|I(G_{\varepsilon})|^{2}+|I(F_{\varepsilon}I(G_{% \varepsilon}))|^{2}\\ \displaystyle+|F_{\varepsilon}|^{2}|I(I(G_{\varepsilon}))|^{2}+|I(|I(G_{% \varepsilon})|^{2})|^{2}+|I(G_{\varepsilon})|^{2}|I(I(G_{\varepsilon}))|^{2}% \Big{)}.

By the Cauchy-Schwarz inequality and a multiple use of (4.1) in the form

 \|\log^{r}(2+|\cdot|)\,I(\zeta)\|_{\operatorname{L}^{p}(\mathbb{R}^{d})}% \lesssim_{p,r}\|\log^{r+1}(2+|\cdot|)\,\zeta\|_{\operatorname{L}^{p}(\mathbb{R% }^{d})},

we are led to

 \displaystyle K_{2}^{\varepsilon}\lesssim\varepsilon^{2d}\Big{(}\|F_{% \varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|\log(2+|\cdot|)F_{% \varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}+\|\log^{\frac{1}{2}}% (2+|\cdot|)F_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{4}\\ \displaystyle+\|\log(2+|\cdot|)F_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb% {R}^{d})}^{2}\|\log(2+|\cdot|)G_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{% R}^{d})}^{2}+\|F_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\|% \log^{2}(2+|\cdot|)G_{\varepsilon}\|_{\operatorname{L}^{4}(\mathbb{R}^{d})}^{2% }\\ \displaystyle+\|\log^{\frac{3}{2}}(2+|\cdot|)G_{\varepsilon}\|_{\operatorname{% L}^{4}(\mathbb{R}^{d})}^{4}+\|\log(2+|\cdot|)G_{\varepsilon}\|_{\operatorname{% L}^{4}(\mathbb{R}^{d})}^{2}\|\log^{2}(2+|\cdot|)G_{\varepsilon}\|_{% \operatorname{L}^{4}(\mathbb{R}^{d})}^{2}\Big{)}.

Inserting the definition of G_{\varepsilon}, rescaling the integrals, and using (3.32), the conclusion follows.

## 5. Covariance structure

In this section, we turn to the limiting covariance structure of the homogenization commutator, as stated in Proposition 2.9.

### 5.1. Structure of the proof and auxiliary results

The main tool to prove Proposition 2.9 is the following stronger version of the spectral gap estimate of Lemma 3.1, which gives an identity (rather than a bound) for the variance of a random variable in terms of its variations. This is an iid version of the so-called Helffer-Sjöstrand representation formula [HS-94, Sjostrand-96] (see also [NS2, MO]), which holds for any product measure \mathbb{P} on \Omega. A proof is included for completeness in Subsection 5.2 below. It is more conveniently formulated in terms of \tilde{\Delta}_{b}X:=X-\mathbb{E}_{a(b)}[X], where \mathbb{E}_{a(b)}[\cdot]:=\mathbb{E}\left[\left.\cdot\,\right|\,\!(a(b^{\prime% }))_{b^{\prime}\neq b}\right] denotes the expectation with respect to the random variable a(b) only. This is a natural variant of the vertical derivative \Delta_{b} and satisfies \mathbb{E}\big{[}|\tilde{\Delta}_{b}X|^{2}\big{]}=\frac{1}{2}\mathbb{E}\big{[}% |\Delta_{b}X|^{2}\big{]}. Note that by definition \tilde{\Delta}_{b}X=\mathbb{E}_{a^{b}(b)}[\Delta_{b}X], where \mathbb{E}_{a^{b}(b)}[\cdot] denotes the expectation with respect to the random variable a^{b}(b) only.

###### Lemma 5.1.

Let \mathbb{P} be a product measure. For all X=X(\boldsymbol{a})\in\operatorname{L}^{2}(\Omega) we have

 \mathrm{Var}\left[X\right]=\sum_{b\in\mathcal{B}}\mathbb{E}\left[(\tilde{% \Delta}_{b}X)\,\mathcal{T}\,(\tilde{\Delta}_{b}X)\right],

where \mathcal{T}:=(\sum_{b\in\mathcal{B}}\tilde{\Delta}_{b}\tilde{\Delta}_{b})^{-1} is a self-adjoint positive operator on \operatorname{L}^{2}(\Omega)/\mathbb{R}:=\{X\in\operatorname{L}^{2}(\Omega):% \mathbb{E}\left[X\right]=0\} with operator norm bounded by 1. In particular, it implies the following covariance inequality: for all X,Y\in\operatorname{L}^{2}(\Omega) we have

 \operatorname{Cov}\left[{X};{Y}\right]\,\leq\,\frac{1}{2}\sum_{b\in\mathcal{B}% }\mathbb{E}\left[|\Delta_{b}X|^{2}\right]^{\frac{1}{2}}\mathbb{E}\left[|\Delta% _{b}Y|^{2}\right]^{\frac{1}{2}}.\qed

The proof of Proposition 2.9(i) below further implies that the effective fluctuation tensor \mathcal{Q} is given by the following formula, with the notation b_{n}:=(0,e_{n}),

 \displaystyle\mathcal{Q}_{ijkl} \displaystyle:= \displaystyle\sum_{n=1}^{d}\mathbb{E}\left[\big{(}M_{ij}^{n}-\mathbb{E}\big{[}% M_{ij}^{n}\big{]}\big{)}\,\mathcal{T}\,\big{(}M_{kl}^{n}-\mathbb{E}\big{[}M_{% kl}^{n}\big{]}\big{)}\right], (5.1) \displaystyle M_{ij}^{n} \displaystyle:= \displaystyle\mathbb{E}_{a^{b_{n}}(b_{n})}\Big{[}(a(b_{n})-a^{b_{n}}(b_{n}))% \big{(}e_{n}\cdot(\nabla\phi_{j}^{*}(0)+e_{j})\big{)}\big{(}e_{n}\cdot(\nabla% \phi_{i}^{b_{n}}(0)+e_{i})\big{)}\Big{]},

in terms of the abstract operator \mathcal{T} defined above. Although not convenient for numerical approximation of \mathcal{Q}, this formula allows to easily deduce the non-degeneracy result contained in Proposition 2.9(ii). In addition, this is key to the proof of Theorem 2 on the RVE method.

### 5.2. Proof of Lemma 5.1

We start with some observations on the difference operator \tilde{\Delta}_{b} on \operatorname{L}^{2}(\Omega). For all X,Y\in\operatorname{L}^{2}(\Omega), by exchangeability of (\boldsymbol{a},\boldsymbol{a}^{b}), we find

 \displaystyle\mathbb{E}\big{[}X\tilde{\Delta}_{b}Y\big{]}=\mathbb{E}\big{[}XY% \big{]}-\mathbb{E}\big{[}X\mathbb{E}_{a(b)}[Y]\big{]}=\mathbb{E}\big{[}XY\big{% ]}-\mathbb{E}\big{[}\mathbb{E}_{a(b)}[X]\mathbb{E}_{a(b)}[Y\big{]}\\ \displaystyle=\mathbb{E}\big{[}XY\big{]}-\mathbb{E}\big{[}Y\mathbb{E}_{a(b)}[X% ]\big{]}=\mathbb{E}\big{[}Y\tilde{\Delta}_{b}X\big{]},

so that \tilde{\Delta}_{b} is symmetric on \operatorname{L}^{2}(\Omega). In addition, we easily compute, for all b,b^{\prime}\in\mathcal{B},

With these observations at hand, we now turn to the study of the (densely defined) operator \mathcal{S}:=\sum_{b\in\mathcal{B}}\tilde{\Delta}_{b}\tilde{\Delta}_{b} on \operatorname{L}^{2}(\Omega). More precisely, we consider the space \operatorname{L}^{2}(\Omega)/\mathbb{R}:=\{X\in\operatorname{L}^{2}(\Omega):% \mathbb{E}\left[X\right]=0\} of mean-zero square-integrable random variables, and we show that \mathcal{S} is an essentially self-adjoint, non-negative operator on \operatorname{L}^{2}(\Omega)/\mathbb{R} with dense image. First, since \mathbb{E}\big{[}\tilde{\Delta}_{b}X\big{]}=0 for all b\in\mathcal{B} and X\in\operatorname{L}^{2}(\Omega), the image \text{Im}\,\mathcal{S} is clearly contained in \operatorname{L}^{2}(\Omega)/\mathbb{R}. Second, for all X\in\operatorname{L}^{2}(\Omega) in the domain of \mathcal{S}, we compute

 \mathbb{E}\left[X\mathcal{S}X\right]=\sum_{b\in\mathcal{B}}\mathbb{E}\big{[}|% \tilde{\Delta}_{b}X|^{2}\big{]}\geq 0,

which shows that \mathcal{S} is non-negative. Third, if X\in\operatorname{L}^{2}(\Omega)/\mathbb{R} in the domain of \mathcal{S} is orthogonal to the image \text{Im}\,\mathcal{S}, then we deduce

 0=\mathbb{E}\left[X\mathcal{S}X\right]=\sum_{b\in\mathcal{B}}\mathbb{E}\big{[}% |\tilde{\Delta}_{b}X|^{2}\big{]},

so that \tilde{\Delta}_{b}X=0 almost surely for all b\in\mathcal{B}, which implies that X is constant.

These properties of \mathcal{S} allow us to define (densely) the inverse \mathcal{T}:=\mathcal{S}^{-1} as an essentially self-adjoint, non-negative operator on \operatorname{L}^{2}(\Omega)/\mathbb{R}. Finally, the spectral gap of Lemma 3.1 implies, for all X\in\operatorname{L}^{2}(\Omega)/\mathbb{R} in the domain of \mathcal{S},

 \|X\|_{\operatorname{L}^{2}(\Omega)}^{2}=\mathrm{Var}\left[X\right]\leq\sum_{b% \in\mathcal{B}}\mathbb{E}\big{[}|\tilde{\Delta}_{b}X|^{2}\big{]}=\mathbb{E}% \left[X\mathcal{S}X\right]\leq\|X\|_{\operatorname{L}^{2}(\Omega)}\|\mathcal{S% }X\|_{\operatorname{L}^{2}(\Omega)},

and hence \|X\|_{\operatorname{L}^{2}(\Omega)}\leq\|\mathcal{S}X\|_{\operatorname{L}^{2}% (\Omega)}, which implies that \mathcal{T}=\mathcal{S}^{-1} on \operatorname{L}^{2}(\Omega)/\mathbb{R} has operator norm bounded by 1.

It remains to establish the representation formula for the variance. By density, it suffices to prove it for all X\in\text{Im}\,\mathcal{S}. Writing X=\mathcal{S}Y for some Y\in\operatorname{L}^{2}(\Omega)/\mathbb{R}, we decompose

 \displaystyle\mathrm{Var}\left[X\right]=\mathbb{E}\left[X\mathcal{S}Y\right]=% \sum_{b\in\mathcal{B}}\mathbb{E}\big{[}\tilde{\Delta}_{b}X\tilde{\Delta}_{b}Y% \big{]}=\sum_{b\in\mathcal{B}}\mathbb{E}\big{[}(\tilde{\Delta}_{b}X)(\tilde{% \Delta}_{b}\mathcal{T}X)\big{]}.

Since the commutation relations (5.2) ensure that \tilde{\Delta}_{b}\mathcal{S}=\mathcal{S}\tilde{\Delta}_{b} holds on the domain of \mathcal{S} in \operatorname{L}^{2}(\Omega), we deduce \tilde{\Delta}_{b}\mathcal{T}=\mathcal{T}\tilde{\Delta}_{b} on \operatorname{L}^{2}(\Omega)/\mathbb{R}, and the above then leads to the desired representation

 \displaystyle\mathrm{Var}\left[X\right]=\sum_{b\in\mathcal{B}}\mathbb{E}\big{[% }(\tilde{\Delta}_{b}X)\mathcal{T}(\tilde{\Delta}_{b}X)\big{]}.

### 5.3. Proof of Proposition 2.9(i)

By polarization and linearity, it is enough to prove (2.8) with F=G\in C^{\infty}_{c}(\mathbb{R}^{d})^{d\times d}. We thus need to establish the convergence of the variance

 \nu_{\varepsilon}:=\mathrm{Var}\left[\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R% }^{d}}F:\Xi(\tfrac{\cdot}{\varepsilon})\right]=\mathrm{Var}\left[\varepsilon^{% \frac{d}{2}}\int_{\mathbb{R}^{d}}F_{\varepsilon}:\Xi\right],

where we have set F_{\varepsilon}:=F(\varepsilon\cdot). We split the proof into two steps.

Step 1. Proof of (2.8).

The Helffer-Sjöstrand representation of Lemma 5.1 applied to the variance \nu_{\varepsilon} yields

 \displaystyle\nu_{\varepsilon} \displaystyle=\varepsilon^{d}\sum_{b\in\mathcal{B}}\mathbb{E}\left[\Big{(}% \tilde{\Delta}_{b}\int_{\mathbb{R}^{d}}F_{\varepsilon}:\Xi\Big{)}\,\mathcal{T}% \,\Big{(}\tilde{\Delta}_{b}\int_{\mathbb{R}^{d}}F_{\varepsilon}:\Xi\Big{)}% \right]. (5.3)

We now appeal to (3.18) in the form

 \displaystyle\Delta_{b}\int_{\mathbb{R}^{d}}F_{\varepsilon}:\Xi=\int_{\mathbb{% R}^{d}}F_{\varepsilon,ij}(\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b}\boldsymbol{% a}(\nabla\phi_{i}^{b}+e_{i})\\ \displaystyle+\int_{\mathbb{R}^{d}}\phi_{j}^{*}(\cdot+e_{k})\nabla_{k}F_{% \varepsilon,ij}e_{k}\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+e_{i})+% \int_{\mathbb{R}^{d}}\nabla h_{\varepsilon,i}\cdot\Delta_{b}\boldsymbol{a}(% \nabla\phi_{i}^{b}+e_{i}),

where the auxiliary field h_{\varepsilon,i} is the unique Lax-Milgram solution in \mathbb{R}^{d} of

 \displaystyle-\nabla^{*}\cdot\boldsymbol{a}^{*}\nabla h_{\varepsilon,i}=\nabla% _{l}^{*}\big{(}\phi_{j}^{*}(\cdot+e_{k})\boldsymbol{a}_{kl}\nabla_{k}F_{% \varepsilon,ij}+\sigma_{jkl}^{*}(\cdot-e_{k})\nabla_{k}^{*}F_{\varepsilon,ij}% \big{)}. (5.4)

Recalling that \tilde{\Delta}_{b}X=\mathbb{E}_{a^{b}(b)}[\Delta_{b}X], inserting this representation formula into (5.3), extracting the first term U_{\varepsilon} defined below, and using that \mathcal{T} on \operatorname{L}^{2}(\Omega)/\mathbb{R} has operator norm bounded by 1, we find

 \displaystyle|\nu_{\varepsilon}-\varepsilon^{d}U_{\varepsilon}|\leq{% \varepsilon^{d}}\sum_{b\in\mathcal{B}}\Big{(}S_{\varepsilon}^{b}T_{\varepsilon% }^{b}+\frac{1}{2}(T_{\varepsilon}^{b})^{2}\Big{)}, (5.5)

where for convenience we define

 \displaystyle U_{\varepsilon} \displaystyle:= \displaystyle\sum_{b\in\mathcal{B}}\mathbb{E}\left[\,(V_{\varepsilon}^{b}-% \mathbb{E}\big{[}V_{\varepsilon}^{b}\big{]})\,\mathcal{T}\,(V_{\varepsilon}^{b% }-\mathbb{E}\big{[}V_{\varepsilon}^{b}\big{]})\,\right], \displaystyle V_{\varepsilon}^{b} \displaystyle:= \displaystyle\mathbb{E}_{a^{b}(b)}\Big{[}\int_{\mathbb{R}^{d}}F_{\varepsilon,% ij}(\nabla\phi_{j}^{*}+e_{j})\cdot\Delta_{b}\boldsymbol{a}(\nabla\phi_{i}^{b}+% e_{i})\Big{]},

while for all b\in\mathcal{B} the error terms are given by

 \displaystyle S_{\varepsilon}^{b} \displaystyle:= \displaystyle\mathbb{E}\left[\Big{(}\int_{\mathbb{R}^{d}}|\Delta_{b}% \boldsymbol{a}||\nabla\phi^{*}+\operatorname{Id}||\nabla\phi^{b}+\operatorname% {Id}||F_{\varepsilon}|\Big{)}^{2}\right]^{\frac{1}{2}},

and by T_{\varepsilon}^{b}:=T_{\varepsilon,1}^{b}+T_{\varepsilon,2}^{b} with

 \displaystyle T_{\varepsilon,1}^{b} \displaystyle:= \displaystyle\sum_{k=1}^{d}\mathbb{E}\left[\Big{(}\int_{\mathbb{R}^{d}}|\Delta% _{b}\boldsymbol{a}||\phi^{*}(\cdot+e_{k})||\nabla\phi^{b}+\operatorname{Id}||% \nabla F_{\varepsilon}|\Big{)}^{2}\right]^{\frac{1}{2}}, \displaystyle T_{\varepsilon,2}^{b} \displaystyle:= \displaystyle\mathbb{E}\left[\Big{(}\int_{\mathbb{R}^{d}}|\Delta_{b}% \boldsymbol{a}||\nabla\phi^{b}+\operatorname{Id}||\nabla h_{\varepsilon}|\Big{% )}^{2}\right]^{\frac{1}{2}}.

We start with the analysis of U_{\varepsilon}. Writing \Delta_{b}\boldsymbol{a}(x)=(a(b)-a^{b}(b))\mathds{1}_{Q(z_{b})}(x)e_{b}% \otimes e_{b} for b=(z_{b},z_{b}+e_{b}), we may compute

 \displaystyle V_{\varepsilon}^{b}=\Big{(}\int_{Q(z_{b})}F_{\varepsilon,ij}\Big% {)}\,\mathbb{E}_{a^{b}(b)}\Big{[}(a(b)-a^{b}(b))\big{(}e_{b}\cdot(\nabla\phi_{% j}^{*}(z_{b})+e_{j})\big{)}\big{(}e_{b}\cdot(\nabla\phi_{i}^{b}(z_{b})+e_{i})% \big{)}\Big{]},

so that, by stationarity,

 \displaystyle\varepsilon^{d}U_{\varepsilon}=\mathcal{Q}_{ijkl}\,\varepsilon^{d% }\sum_{z\in\mathbb{Z}^{d}}\Big{(}\int_{Q(z)}F_{\varepsilon,ij}\Big{)}\Big{(}% \int_{Q(z)}F_{\varepsilon,kl}\Big{)},

where the coefficient \mathcal{Q}_{ijkl} is defined in (5.1) above. Since \mathcal{T} on \operatorname{L}^{2}(\Omega)/\mathbb{R} has operator norm bounded by 1, the moment bounds of Lemma 3.3 yield

 \displaystyle|\mathcal{Q}_{ijkl}|\,\lesssim\,\sum_{n=1}^{d}\mathbb{E}\left[|% \nabla\phi^{*}+\operatorname{Id}|^{2}|\nabla\phi^{b_{n}}+\operatorname{Id}|^{2% }\right]\,\lesssim\,1.

We may then estimate the discretization error

 \displaystyle\Big{|}\varepsilon^{d}U_{\varepsilon}-\mathcal{Q}_{ijkl}\,\int_{% \mathbb{R}^{d}}F_{ij}F_{kl}\Big{|} \displaystyle= \displaystyle\Big{|}\varepsilon^{d}U_{\varepsilon}-\mathcal{Q}_{ijkl}\,% \varepsilon^{d}\int_{\mathbb{R}^{d}}F_{\varepsilon,ij}F_{\varepsilon,kl}\Big{|} (5.6) \displaystyle\lesssim \displaystyle\varepsilon^{d}\sum_{z\in\mathbb{Z}^{d}}\int_{Q(z)}\Big{|}F_{% \varepsilon}(x)-\int_{Q(z)}F_{\varepsilon}\Big{|}^{2}dx \displaystyle\lesssim \displaystyle\varepsilon^{d}\int_{\mathbb{R}^{d}}|DF_{\varepsilon}|^{2}=% \varepsilon^{2}\int_{\mathbb{R}^{d}}|DF|^{2}.

We now turn to the estimate of the right-hand side of (5.5). Using |\Delta_{b}\boldsymbol{a}(x)|\lesssim\mathds{1}_{Q(z_{b})}(x) and the moment bounds of Lemma 3.3, we obtain

 \displaystyle S_{\varepsilon}^{b}\,\lesssim\,\mathbb{E}\left[|\nabla\phi^{*}+% \operatorname{Id}|^{2}|\nabla\phi^{b}+\operatorname{Id}|^{2}\right]^{\frac{1}{% 2}}\int_{Q(z_{b})}|F_{\varepsilon}|\,\lesssim\,\int_{Q(z_{b})}|F_{\varepsilon}|.

Hence, by the Cauchy-Schwarz inequality,

 \displaystyle\sum_{b\in\mathcal{B}}S_{\varepsilon}^{b}T_{\varepsilon}^{b}\,% \lesssim\,\sum_{b\in\mathcal{B}}T_{\varepsilon}^{b}\int_{Q(z_{b})}|F_{% \varepsilon}| \displaystyle\lesssim\,\|F_{\varepsilon}\|_{\operatorname{L}^{2}(\mathbb{R}^{d% })}\Big{(}\sum_{b\in\mathcal{B}}(T_{\varepsilon}^{b})^{2}\Big{)}^{\frac{1}{2}} \displaystyle\lesssim\,\varepsilon^{-\frac{d}{2}}\|F\|_{\operatorname{L}^{2}(% \mathbb{R}^{d})}\Big{(}\sum_{b\in\mathcal{B}}(T_{\varepsilon}^{b})^{2}\Big{)}^% {\frac{1}{2}}, (5.7)

and it remains to estimate

 \sum_{b\in\mathcal{B}}(T_{\varepsilon}^{b})^{2}\leq 2\sum_{b\in\mathcal{B}}(T_% {\varepsilon,1}^{b})^{2}+2\sum_{b\in\mathcal{B}}(T_{\varepsilon,2}^{b})^{2}.

First, using |\Delta_{b}\boldsymbol{a}(x)|\lesssim\mathds{1}_{Q(z_{b})}(x) and the moment bounds of Lemma 3.3, we find

 \displaystyle\varepsilon^{d}\sum_{b\in\mathcal{B}}(T_{\varepsilon,1}^{b})^{2}% \,\lesssim_{\alpha,p}\,\varepsilon^{d}\|\mu_{d}(|\cdot|)^{\frac{1}{2}}\nabla F% _{\varepsilon}\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}^{2}. (5.8)

Second, arguing as in the proof of Proposition 2.1 (cf. (3.16)), using the large-scale weighted Calderón-Zygmund theory (cf. Lemma 3.4) applied to equation (5.4) for h_{\varepsilon}, we obtain for all 0<p-1\ll 1 and all \alpha>d,

 \displaystyle\varepsilon^{d}\sum_{b\in\mathcal{B}}(T_{\varepsilon,2}^{b})^{2}% \,\lesssim_{\alpha,p}\,\varepsilon^{\frac{d}{p}}\|w_{\varepsilon}^{\alpha\frac% {p-1}{2p}}\mu_{d}(|\cdot|)^{\frac{1}{2}}\nabla F_{\varepsilon}\|_{% \operatorname{L}^{2p}(\mathbb{R}^{d})}^{2}. (5.9)

Rescaling the integrals and using (3.32) and Hölder’s inequality, we find

 \displaystyle\varepsilon^{d}\sum_{b\in\mathcal{B}}(T_{\varepsilon}^{b})^{2}\,% \lesssim_{\alpha,p}\,\varepsilon^{2}\mu_{d}(\tfrac{1}{\varepsilon})\|w_{1}^{% \alpha\frac{p-1}{2p}}\mu_{d}(|\cdot|)^{\frac{1}{2}}DF\|_{\operatorname{L}^{2p}% (\mathbb{R}^{d})}^{2}.

and the conclusion (2.8) follows.

Step 2. Proof of the Green-Kubo formula (2.9).

In order to establish (2.9), it suffices to repeat the argument of Step 1 with the test function F=\mathds{1}_{Q}\,e_{i}\otimes e_{j} (hence F_{\varepsilon}=\mathds{1}_{\frac{1}{\varepsilon}Q}\,e_{i}\otimes e_{j}), for some fixed 1\leq i,j\leq d. Lemma 5.1 again leads to (5.5), and we briefly indicate how to analyze the different terms in the present setting. First, the estimate (5.6) is replaced by the following (no summation over repeated indices),

 \displaystyle|\varepsilon^{d}U_{\varepsilon}-\mathcal{Q}_{ijij}| \displaystyle\lesssim \displaystyle\varepsilon^{d}\sum_{z\in\mathbb{Z}^{d}}\int_{Q}\Big{(}\mathds{1}% _{\frac{1}{\varepsilon}Q}(z+x)-\int_{Q}\mathds{1}_{\frac{1}{\varepsilon}Q}(z+y% )dy\Big{)}^{2}dx \displaystyle\leq \displaystyle\varepsilon^{d}\sum_{z\in\mathbb{Z}^{d}}\mathds{1}_{(z+Q)\cap% \partial(\frac{1}{\varepsilon}Q)\neq\varnothing}\leavevmode\nobreak\ \lesssim% \leavevmode\nobreak\ \varepsilon.

Second, the estimate (5.7) remains unchanged. Third, using estimates (5.8) and (5.9), and noting that |\nabla F_{\varepsilon}|\lesssim\mathds{1}_{A_{\varepsilon}} with A_{\varepsilon}:=B+\partial Q_{\frac{1}{\varepsilon}} and that w_{\varepsilon}\lesssim 1 and \mu_{d}(|\cdot|)\lesssim\mu_{d}(\frac{1}{\varepsilon}) on A_{\varepsilon}, we deduce

 \displaystyle\sum_{b\in\mathcal{B}}(T_{\varepsilon,1}^{b})^{2} \displaystyle\lesssim \displaystyle\mu_{d}(\tfrac{1}{\varepsilon})|A_{\varepsilon}|\,\lesssim\,% \varepsilon^{1-d}\mu_{d}(\tfrac{1}{\varepsilon}), \displaystyle\sum_{b\in\mathcal{B}}(T_{\varepsilon,2}^{b})^{2} \displaystyle\lesssim \displaystyle\varepsilon^{-d\frac{p-1}{p}}\mu_{d}(\tfrac{1}{\varepsilon})|A_{% \varepsilon}|^{\frac{1}{p}}\,\lesssim\,\varepsilon^{\frac{1}{p}-d}\mu_{d}(% \tfrac{1}{\varepsilon}),

and the conclusion (2.9) follows.

### 5.4. Proof of Proposition 2.9(ii)

The following proof of the non-degeneracy of \mathcal{Q} is based on the Helffer-Sjöstrand representation formula (see also [MO, Remark 2.3]), and constitutes a shorter alternative to the corresponding proof in [GN]. Given a fixed direction e\in\mathbb{R}^{d}\setminus\{0\}, and letting \phi_{e} denote the corrector in this direction, we may write, in view of formula (5.1) (with \phi_{e}^{*}=\phi_{e} by symmetry of the coefficients),

 \displaystyle(e\otimes e):\mathcal{Q}\,(e\otimes e)=\sum_{n=1}^{d}\mathbb{E}% \big{[}(e\cdot M^{n}e)\mathcal{T}(e\cdot M^{n}e)\big{]}, (5.10) \displaystyle e\cdot M^{n}e:=\mathbb{E}_{a^{b_{n}}(b_{n})}\Big{[}(a(b_{n})-a^{% b_{n}}(b_{n}))\big{(}e_{n}\cdot(\nabla\phi_{e}(0)+e)\big{)}\big{(}e_{n}\cdot(% \nabla\phi_{e}^{b_{n}}(0)+e)\big{)}\Big{]},

since the exchangeability of (\boldsymbol{a},\boldsymbol{a}^{b_{n}}) indeed yields \mathbb{E}\big{[}e\cdot M^{n}e\big{]}=0 for all n. We start with a suitable reformulation of e\cdot M^{n}e. Considering the difference of the corrector equation (2.4) for \phi_{e} and \phi_{e}^{b_{n}} in the form -\nabla^{*}\cdot\boldsymbol{a}^{b_{n}}\nabla(\phi_{e}^{b_{n}}-\phi_{e})=\nabla% ^{*}\cdot(\boldsymbol{a}^{b_{n}}-\boldsymbol{a})(\nabla\phi_{e}+e), an integration by parts yields

 \displaystyle\int_{\mathbb{R}^{d}}\nabla(\phi_{e}^{b_{n}}-\phi_{e})\cdot% \boldsymbol{a}^{b_{n}}\nabla(\phi_{e}^{b_{n}}-\phi_{e})=-\int_{\mathbb{R}^{d}}% \nabla(\phi_{e}^{b_{n}}-\phi_{e})\cdot(\boldsymbol{a}^{b_{n}}-\boldsymbol{a})(% \nabla\phi_{e}+e)\\ \displaystyle=(a(b_{n})-a^{b_{n}}(b_{n}))(e_{n}\cdot\nabla(\phi_{e}^{b_{n}}-% \phi_{e})(0))(e_{n}\cdot(\nabla\phi_{e}(0)+e)).

Hence, by definition of e\cdot M^{n}e,

 \displaystyle e\cdot M^{n}e=\mathbb{E}_{a^{b_{n}}(b_{n})}\bigg{[}\int_{\mathbb% {R}^{d}}\nabla(\phi^{b_{n}}_{e}-\phi_{e})\cdot\boldsymbol{a}^{b_{n}}\nabla(% \phi_{e}^{b_{n}}-\phi_{e})\bigg{]}\\ \displaystyle+(a(b_{n})-\mathbb{E}\left[a(b_{n})\right])(e_{n}\cdot(\nabla\phi% _{e}(0)+e))^{2}. (5.11)

We now argue by contradiction. If (e\otimes e):\mathcal{Q}\,(e\otimes e)=0, then by formula (5.10) and by the non-negativity of \mathcal{T} we would have \mathbb{E}\big{[}(e\cdot M^{n}e)\mathcal{T}(e\cdot M^{n}e)\big{]}=0 for all n. Let 1\leq n\leq d be momentarily fixed. Recalling that \mathcal{T}=\mathcal{S}^{-1} with \mathcal{S}=\sum_{b\in\mathcal{B}}\tilde{\Delta}_{b}\tilde{\Delta}_{b}, this would imply

 0=\mathbb{E}\left[(\mathcal{T}(e\cdot M^{n}e))\mathcal{S}(\mathcal{T}(e\cdot M% ^{n}e))\right]=\sum_{b\in\mathcal{B}}\mathbb{E}\left[\big{|}\tilde{\Delta}_{b}% \mathcal{T}(e\cdot M^{n}e)\big{|}^{2}\right],

hence \mathcal{T}(e\cdot M^{n}e)=0, and thus e\cdot M^{n}e=0 almost surely. Formula (5.11) would then imply

 \displaystyle(a(b_{n})-\mathbb{E}\left[a(b_{n})\right])(e_{n}\cdot(\nabla\phi_% {e}(0)+e))^{2}\\ \displaystyle=-\mathbb{E}_{a^{b_{n}}(b_{n})}\bigg{[}\int_{\mathbb{R}^{d}}% \nabla(\phi^{b_{n}}_{e}-\phi_{e})\cdot\boldsymbol{a}^{b_{n}}\nabla(\phi_{e}^{b% _{n}}-\phi_{e})\bigg{]}, (5.12)

almost surely. Since the law of a(b_{n}) is non-degenerate, the event a(b_{n})>\mathbb{E}\left[a(b_{n})\right] occurs with a positive probability. Conditioning on this event, the left-hand side in (5.12) is non-negative, and the non-positivity of the right-hand side would then imply that both sides vanish, that is,

almost surely. Since the integrand in this last expectation is non-negative, we would deduce that the event a(b_{n})>\mathbb{E}\left[a(b_{n})\right] entails e_{n}\cdot(\nabla\phi_{e}(0)+e)=0 and \nabla\phi_{e}(0)=\nabla\phi_{e}^{b_{n}}(0), and thus also e_{n}\cdot(\nabla\phi_{e}^{b_{n}}(0)+e)=0 almost surely. Since this last event is independent of a(b_{n}), hence of the conditioning event, we would deduce unconditionally e_{n}\cdot(\nabla\phi_{e}^{b_{n}}(0)+e)=0 almost surely. By exchangeability of (\boldsymbol{a},\boldsymbol{a}^{b_{n}}), this means e_{n}\cdot(\nabla\phi_{e}(0)+e)=0 almost surely. As this holds for any n, we would conclude \nabla\phi_{e}(0)+e=0 almost surely, and taking the expectation would lead to a contradiction.

## 6. Approximation of the fluctuation tensor

In this section, we analyze the RVE method for the approximation of the fluctuation tensor \mathcal{Q} as stated in Theorem 2.

### 6.1. Structure of the proof and auxiliary results

The estimate on the standard deviation is obtained similarly as the CLT scaling in Proposition 2.1, noting that the large-scale Calderón-Zygmund result of Lemma 3.4 also holds for the periodized operator -\nabla^{*}\cdot\boldsymbol{a}_{L}\nabla on Q_{L}.444The only issue concerns the corresponding moment bound \mathbb{E}\left[r_{*,L}^{q}\right]\lesssim_{q}1 for all q<\infty, which by definition of r_{*,L} (cf. [GNO-reg]) is a consequence of a sup-bound based on the version of Lemma 3.3 for the periodized correctors (\phi_{L},\sigma_{L}) (cf. [GNO1]). The characterization (1.16) of \mathcal{Q} and the estimate of the systematic error are deduced as corollaries of formula (5.1) for the fluctuation tensor \mathcal{Q}, together with the following crucial estimates on the periodized corrector \phi_{L}. (The first estimate on \nabla\phi_{L} is stated as such in [GNO1, Proposition 1], and the second estimate follows from a decomposition of the difference \nabla\phi_{L}-\nabla\phi via massive approximation of the corrector and Richardson extrapolation, applying [GN, Lemma 2.8 and estimate (2.68)], and optimizing the mass.)

###### Lemma 6.1 ([GNO1, GN]).

Let d\geq 2 and let \mathbb{P} be a product measure. For all L\geq 2 and all q<\infty we have

### 6.2. Proof of Theorem 2

We split the proof into two steps: we first estimate the variance of the RVE approximation, and then we turn to the characterization (1.16) of \mathcal{Q} and to the systematic error of the RVE approximation.

Step 1. Proof of |\mathrm{Var}\left[\mathcal{Q}_{L,N}\right]|^{\frac{1}{2}}\lesssim N^{-\frac{1% }{2}}.

Since the realizations \bar{\boldsymbol{a}}_{L}^{(n)} are iid copies of \bar{\boldsymbol{a}}_{L}, the definition (1.17) of \mathcal{Q}_{L,N} leads after straightforward computations to

 \displaystyle\mathrm{Var}\left[\mathcal{Q}_{L,N}\right]=N^{-1}\mathrm{Var}% \left[\big{(}L^{\frac{d}{2}}\bar{\boldsymbol{a}}_{L}^{*}-\mathbb{E}\big{[}L^{% \frac{d}{2}}\bar{\boldsymbol{a}}_{L}^{*}\big{]}\big{)}^{\otimes 2}\right],

and hence,

 \displaystyle|\mathrm{Var}\left[\mathcal{Q}_{L,N}\right]\!|\lesssim N^{-1}\,% \mathbb{E}\left[\big{|}L^{\frac{d}{2}}(\bar{\boldsymbol{a}}_{L}-\mathbb{E}% \left[\bar{\boldsymbol{a}}_{L}\right])\big{|}^{4}\right].

Arguing as in [GNO1, Lemma 2], the spectral gap estimate of Lemma 3.1 is seen to imply the following inequality: for all X=X(\boldsymbol{a})\in\operatorname{L}^{4}(\Omega),

 \displaystyle\mathbb{E}\left[(X-\mathbb{E}\left[X\right])^{4}\right]\,\leq\,4% \,\mathbb{E}\left[\bigg{(}\sum_{b\in\mathcal{B}}|\Delta_{b}X|^{2}\bigg{)}^{2}% \right].

Applying this inequality to (each component of) X=\bar{\boldsymbol{a}}_{L}, we deduce

 \displaystyle|\mathrm{Var}\left[\mathcal{Q}_{L,N}\right]\!|\lesssim N^{-1}\,% \mathbb{E}\left[\bigg{(}\sum_{b\in\mathcal{B}_{L}}\Big{(}L^{-\frac{d}{2}}\int_% {Q_{L}}\Delta_{b}\big{(}\boldsymbol{a}_{L}(\nabla\phi_{L}+\operatorname{Id})% \big{)}\Big{)}^{2}\bigg{)}^{2}\right].

Arguing as in the proof of Proposition 2.1 (with \varepsilon replaced by \frac{1}{L} and F_{\varepsilon} replaced by \operatorname{Id}), using the periodized version of Lemma 3.4 and the moment bounds of Lemma 6.1, the conclusion follows.

Step 2. Proof of (1.16) and of the systematic error estimate |\mathbb{E}\left[\mathcal{Q}_{L,N}\right]-\mathcal{Q}|\lesssim L^{-\frac{d}{2}% }\log^{\frac{d}{2}}L.

Since the realizations \bar{\boldsymbol{a}}_{L}^{(n)} are iid copies of \bar{\boldsymbol{a}}_{L}, the definition (1.17) of \mathcal{Q}_{L,N} yields after straightforward computations \mathbb{E}\left[\mathcal{Q}_{L,N}\right]=\mathrm{Var}\big{[}L^{\frac{d}{2}}% \bar{\boldsymbol{a}}_{L}^{*}\big{]}, that is,

 \displaystyle\mathbb{E}\left[(\mathcal{Q}_{L,N})_{ijkl}\right] \displaystyle= \displaystyle\operatorname{Cov}\left[{L^{\frac{d}{2}}\bar{\boldsymbol{a}}_{L,% ji}};{L^{\frac{d}{2}}\bar{\boldsymbol{a}}_{L,lk}}\right] (6.1) \displaystyle= \displaystyle L^{-d}\,\operatorname{Cov}\left[{\int_{Q_{L}}e_{j}\cdot% \boldsymbol{a}_{L}(\nabla\phi_{L,i}+e_{i})};{\int_{Q_{L}}e_{l}\cdot\boldsymbol% {a}_{L}(\nabla\phi_{L,k}+e_{k})}\right].

For b\in\mathcal{B}, we write b=(z_{b},z_{b}+e_{b}). Using the periodized corrector equation (2.10) and its vertical derivative, and recalling that \Delta_{b}\boldsymbol{a}_{L}(x)=(a(b)-a^{b}(b))\mathds{1}_{Q(z_{b})}(x)e_{b}% \otimes e_{b} for b\in\mathcal{B}_{L} and x\in Q_{L}, we find

 \displaystyle\Delta_{b}\int_{Q_{L}}e_{j}\cdot\boldsymbol{a}_{L}(\nabla\phi_{L,% i}+e_{i})=\int_{Q_{L}}e_{j}\cdot\Delta_{b}\boldsymbol{a}_{L}(\nabla\phi_{L,i}^% {b}+e_{i})+\int_{Q_{L}}e_{j}\cdot\boldsymbol{a}_{L}\nabla\Delta_{b}\phi_{L,i} (6.2) \displaystyle                 = \displaystyle\int_{Q_{L}}e_{j}\cdot\Delta_{b}\boldsymbol{a}_{L}(\nabla\phi_{L,% i}^{b}+e_{i})-\int_{Q_{L}}\nabla\phi_{L,j}^{*}\cdot\boldsymbol{a}_{L}\nabla% \Delta_{b}\phi_{L,i} \displaystyle                 = \displaystyle\int_{Q_{L}}(\nabla\phi_{L,j}^{*}+e_{j})\cdot\Delta_{b}% \boldsymbol{a}_{L}(\nabla\phi_{L,i}^{b}+e_{i}) \displaystyle                 = \displaystyle(a(b)-a^{b}(b))(e_{b}\cdot(\nabla\phi_{L,j}^{*}(z_{b})+e_{j}))(e_% {b}\cdot(\nabla\phi_{L,i}^{b}(z_{b})+e_{i})).

Applying the Helffer-Sjöstrand representation formula of Lemma 5.1 to the covariance in (6.1), we obtain by stationarity, as in the proof of Proposition 2.9(i),

 \displaystyle\mathbb{E}\left[(\mathcal{Q}_{L,N})_{ijkl}\right]\,=\,\sum_{n=1}^% {d}\mathbb{E}\left[M_{ij,L}^{n}\,\mathcal{T}\,M_{kl,L}^{n}\right],

where we have set

 M_{ij,L}^{n}\,:=\,\mathbb{E}_{a^{b_{n}}(b_{n})}\Big{[}(a(b_{n})-a^{b_{n}}(b_{n% }))(e_{n}\cdot(\nabla\phi_{L,j}^{*}(0)+e_{j}))(e_{n}\cdot(\nabla\phi_{L,i}^{b_% {n}}(0)+e_{i}))\Big{]}.

Noting that (6.2) implies \mathbb{E}\big{[}M_{ij,L}^{n}\big{]}=0, comparing the above identity for \mathbb{E}\left[(\mathcal{Q}_{L,N})_{ijkl}\right] with formula (5.1) for \mathcal{Q}, and using that the operator \mathcal{T} on \operatorname{L}^{2}(\Omega)/\mathbb{R} has operator norm bounded by 1, we deduce

 \displaystyle\big{|}\mathbb{E}\left[(\mathcal{Q}_{L,N})_{ijkl}\right]-\mathcal% {Q}_{ijkl}\big{|}\,\lesssim\,\mathbb{E}\left[|\nabla(\phi_{L}-\phi)(0)|^{4}% \right]^{\frac{1}{4}}\big{(}\,\mathbb{E}\left[|\nabla\phi_{L}|^{4}\right]+% \mathbb{E}\left[|\nabla\phi|^{4}\right]\big{)}^{\frac{3}{4}},

and the conclusion follows from Lemmas 3.3 and 6.1.

## Acknowledgements

The work of MD is supported by F.R.S.-FNRS (Belgian National Fund for Scientific Research) through a Research Fellowship. MD and AG acknowledge financial support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2014-2019 Grant Agreement QUANTHOM 335410). The authors acknowledge the hospitality of IHÉS, where this work was initiated in February 2015, and the support of the Chaire Schlumberger.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters