Higher-order pathwise theory of fluctuations in stochastic homogenization

TL;DR

This paper develops a higher-order, pathwise theory for Gaussian fluctuations in stochastic homogenization of elliptic equations, revealing the algebraic structure of correctors and providing new estimates for random coefficients.

Contribution

It extends the homogenization theory to higher-order fluctuations, offering a detailed algebraic framework and technical estimates for Gaussian random environments.

Findings

Fluctuations are Gaussian at order d/2 in dimension d.

Higher-order correctors and expansions clarify the structure of fluctuations.

Introduction of annealed Calderón-Zygmund estimates enhances technical analysis.

Abstract

We consider linear elliptic equations in divergence form with stationary random coefficients of integrable correlations. We characterize the fluctuations of a macroscopic observable of a solution to relative order $\frac{d}{2}$, where $d$ is the spatial dimension; the fluctuations turn out to be Gaussian. As for previous work on the leading order, this higher-order characterization relies on a pathwise proximity of the macroscopic fluctuations of a general solution to those of the (higher-order) correctors, via a (higher-order) two-scale expansion injected into the homogenization commutator, thus confirming the scope of this notion. This higher-order generalization sheds a clearer light on the algebraic structure of the higher-order versions of correctors, flux correctors, two-scale expansions, and homogenization commutators. It reveals that in the same way as this algebra provides a higher-order theory for microscopic spatial oscillations, it also provides a higher-order theory for macroscopic random fluctuations, although both phenomena are not directly related. We focus on the model framework of an underlying Gaussian ensemble, which allows for an efficient use of (second-order) Malliavin calculus for stochastic estimates. On the technical side, we introduce annealed Calder\'on-Zygmund estimates for the elliptic operator with random coefficients, which conveniently upgrade the known quenched large-scale estimates.

Equations

$$|\mathbf{a}(x)\xi |\le |\xi |,\xi \cdot \mathbf{a}(x)\xi \ge \lambda |\xi {|}^2,\text{for all x,ξ∈Rd},$$

$$-\nabla \cdot \mathbf{a}(\frac{\cdot }{\epsilon })\nabla u_{\epsilon }=\nabla \cdot f\text{in Rd},$$

$$-\nabla \cdot \bar{\mathbf{a}}\nabla \bar{u}=\nabla \cdot f\text{in Rd},$$

$$-\nabla \cdot \mathbf{a}(\nabla {\phi }_i+e_i)=0\text{in Rd,}$$

$$\big{\|}\nabla\big{(}u_{\varepsilon}-(1+\varepsilon\varphi_{i}(\tfrac{\cdot}{\varepsilon})\nabla_{i})\bar{u}\big{)}\big{\|}_{\operatorname{L}^{2}(\mathbb{R}^{d}\times\Omega)}=O(\varepsilon).$$

$$-\nabla \cdot \mathbf{a}\nabla (x_i+{\phi }_i)=0,$$

$$-\nabla \cdot (\bar{\mathbf{a}}+\epsilon {\bar{\mathbf{a}}}_i^2{\nabla }_i)\nabla {\bar{U}}_{\epsilon }^2=\nabla \cdot f\text{in Rd},$$

$$-\nabla \cdot (\bar{\mathbf{a}}\nabla {\tilde{u}}^2+{\bar{\mathbf{a}}}_i^2\nabla {\nabla }_i\bar{u})=0\text{in Rd}.$$

$${\Xi }_{\epsilon }^1[\nabla w_{\epsilon }]:=(\mathbf{a}(\frac{\cdot }{\epsilon })-\bar{\mathbf{a}})\nabla w_{\epsilon }.$$

$$\Xi_{\varepsilon}^{2}[\nabla w_{\varepsilon}]:=\big{(}\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})-\bar{\boldsymbol{a}}-\varepsilon\bar{\boldsymbol{a}}^{2}_{i}\nabla_{i}\big{)}\nabla w_{\varepsilon},$$

$${\Xi }_{\epsilon }^{\circ ,1}[\nabla \bar{w}]:={\nabla }_i\bar{w}(\mathbf{a}(\frac{\cdot }{\epsilon })-\bar{\mathbf{a}})(\nabla {\phi }_i(\frac{\cdot }{\epsilon })+e_i).$$

$$\Xi^{\circ,2}_{\varepsilon}[\nabla\bar{w}]:=\nabla_{i}\bar{w}\,\big{(}\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})-\bar{\boldsymbol{a}}-\varepsilon\bar{\boldsymbol{a}}^{2}_{j}\nabla_{j}\big{)}(\nabla\varphi_{i}(\tfrac{\cdot}{\varepsilon})+e_{i})\\ +\varepsilon\nabla^{2}_{ij}\bar{w}\,\Big{(}\big{(}\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})-\bar{\boldsymbol{a}}-\varepsilon\bar{\boldsymbol{a}}^{2}_{l}\nabla_{l}\big{)}\big{(}\nabla\varphi_{ij}^{2}(\tfrac{\cdot}{\varepsilon})+\varphi_{i}(\tfrac{\cdot}{\varepsilon})e_{j}\big{)}-\bar{\boldsymbol{a}}^{2}_{j}\big{(}\nabla\varphi_{i}(\tfrac{\cdot}{\varepsilon})+e_{i}\big{)}\Big{)}.$$

$${\Xi }_{\epsilon }^{\circ ,1}[\nabla \bar{w}](x)$$

$${\Xi }_{\epsilon }^{\circ ,2}[\nabla \bar{w}](x)$$

$$\mathbf{a}(x):=a_0(G(x)),$$

$$\int_{\mathbb{R}^{d}}\Big{(}\sup_{B(x)}|c_{0}|\Big{)}\,dx\,\leq\,1,$$

$$\qquad N_{p}\bigg{(}\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\nabla u_{\varepsilon}-\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}\nabla\bar{v}_{\varepsilon}^{n}\cdot\Xi_{\varepsilon}^{n}[\nabla u_{\varepsilon}]\bigg{)}\,\lesssim_{p,f,g}\,\varepsilon^{n},$$

$$\qquad N_{p}\bigg{(}\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\big{(}\Xi_{\varepsilon}^{n}[\nabla u_{\varepsilon}]-\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{u}_{\varepsilon}^{n}]\big{)}\bigg{)}\,\lesssim_{p,f,g}\,\varepsilon^{n}\mu_{d,n}(\tfrac{1}{\varepsilon}),$$

$$\qquad\varepsilon^{n}\mu_{d,n}(\tfrac{1}{\varepsilon})\,:=\,\left\{\begin{array}[]{lll}\varepsilon^{n}&:&\text{$n<\ell$},\\ \varepsilon^{\frac{d}{2}}|\!\log\varepsilon|^{\frac{1}{2}}&:&\text{$n=\ell$ and $d$ even},\\ \varepsilon^{\frac{d}{2}}&:&\text{$n=\ell$ and $d$ odd}.\end{array}\right.$$

$${\mathcal{Q}}_{\epsilon }^n[g,\nabla \bar{w}]:=\sum _{m=0}^{n-1}{\epsilon }^m\sum _{\genfrac{}{}{0ex}{}{l^{\prime },k^{\prime },l,k\ge 0}{l^{\prime }+k^{\prime }+l+k=m}}{\mathcal{Q}}_{l^{\prime }{k}^{\prime }lk}^m\odot {\int }_{{\mathbb{R}}^d}{\nabla }^{l^{\prime }}g\otimes {\nabla }^{k^{\prime }}\nabla \bar{w}\otimes {\nabla }^{l}g\otimes {\nabla }^k\nabla \bar{w},$$

$$\qquad\bigg{|}\mathrm{Var}\left[\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]\right]-\mathcal{Q}_{\varepsilon}^{n}[g,\nabla\bar{w}]\bigg{|}\lesssim_{g,\bar{w}}\varepsilon^{n}\mu_{d,n}(\tfrac{1}{\varepsilon}).$$

$$\qquad\operatorname{d}_{\mathcal{N}}\Big{(}{\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]}\Big{)}\,\lesssim_{g,\bar{w}}\,\frac{\varepsilon^{\frac{d}{2}}|\!\log\varepsilon|\big{(}|\!\log\varepsilon|\log^{2}\!|\!\log\varepsilon|\big{)}^{\frac{1}{2d}n(n-1)}}{\mathrm{Var}\big{[}\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]\big{]}},$$

$$d_{\mathcal{N}}(X):=W_2\left(\frac{X-\mathbb{E}\left[X\right]}{\mathrm{Var}{\left[X\right]}^{\frac{1}{2}}};\mathcal{N}\right)+d_{\mathrm{TV}}\left(\frac{X-\mathbb{E}\left[X\right]}{\mathrm{Var}{\left[X\right]}^{\frac{1}{2}}};\mathcal{N}\right),$$

$$X_{\varepsilon}^{\circ}(g,f):=\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}\nabla\bar{v}_{\varepsilon}^{\ell}\cdot\big{(}\Xi_{\varepsilon}^{\circ,\ell}[\nabla\bar{u}_{\varepsilon}^{\ell}]-\mathbb{E}\big{[}\Xi_{\varepsilon}^{\circ,\ell}[\nabla\bar{u}_{\varepsilon}^{\ell}]\big{]}\big{)},$$

$$\mathbb{E}\left[\Big{(}\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot(\nabla u_{\varepsilon}-\mathbb{E}\left[\nabla u_{\varepsilon}\right])-X_{\varepsilon}^{\circ}(g,f)\Big{)}^{p}\right]^{\frac{1}{p}}\,\lesssim_{p,f,g}\,\varepsilon^{\frac{d}{2}}\times\left\{\begin{array}[]{lll}|\!\log\varepsilon|^{\frac{1}{2}}&:&\text{$d$ even},\\ 1&:&\text{$d$ odd},\end{array}\right.$$

$$d\left(X_{\epsilon }^{\circ }(g,f)\text{leavevmode};\text{leavevmode}{\mathcal{Q}}_{\epsilon }^{\ell }[\nabla {\bar{v}}_{\epsilon }^{\ell },\nabla {\bar{u}}_{\epsilon }^{\ell }{]}^{\frac{1}{2}}\mathcal{N}\right){\lesssim }_{f,g}{\epsilon }^{\frac{d}{2}}|\log \epsilon {|}^{1+\rho }(\log |\log \epsilon |{)}^{2\rho },$$

$$\rho\,:=\,\frac{1}{8d}\times\left\{\begin{array}[]{lll}d(d-2)&:&\text{$d$ even},\\ (d+1)(d-1)&:&\text{$d$ odd}.\end{array}\right.\qed$$

$${\mathrm{Sym}}_{i_1\dots i_k}{T}_{i_1\dots i_k}:=\frac{1}{k!}\sum _{\sigma \in {\mathcal{S}}_k}{T}_{i_{\sigma (1)}\dots i_{\sigma (k)}},$$

$$\nabla\cdot\boldsymbol{a}\nabla F^{n}[\bar{q}]=\nabla\cdot\Big{(}\sum_{k=1}^{n-1}\bar{\boldsymbol{a}}^{k}_{i_{1}\ldots i_{k-1}}\nabla^{k-1}_{i_{1}\ldots i_{k-1}}\Big{)}\nabla\bar{q},$$

$$\displaystyle-\nabla\cdot\boldsymbol{a}\nabla\varphi^{n}_{i_{1}\ldots i_{n}}=\nabla\cdot\big{(}(\boldsymbol{a}\varphi^{n-1}_{i_{1}\ldots i_{n-1}}-\sigma^{n-1}_{i_{1}\ldots i_{n-1}})\,e_{i_{n}}\big{)},$$

$$\nabla {\phi }_{i_1\dots i_n}^n\text{leavevmode}\text{stationary},\mathbb{E}\left[|\nabla {\phi }_{i_1\dots i_n}^n{|}^2\right]<\infty ,\mathbb{E}\left[\nabla {\phi }_{i_1\dots i_n}^n\right]=0.$$

Full text

1. 1 Introduction
2. 2 Higher-order theory of oscillations
3. 3 Higher-order homogenization commutators
4. 4 Reminder on Malliavin calculus
5. 5 Representation formulas
6. 6 Annealed Calderón-Zygmund theory
7. 7 Proof of the higher-order pathwise result
8. 8 Proof of the convergence of the covariance structure
9. 9 Proof of the normal approximation result
10. A More Malliavin calculus

1. Introduction

1.1. General overview

Let $\boldsymbol{a}$ be a random coefficient field on $\mathbb{R}^{d}$ that is stationary and ergodic and satisfies the boundedness and ellipticity properties

[TABLE]

for some $\lambda>0$, and denote by $(\Omega,\mathcal{F},\mathbb{P})$ the underlying probability space. Given a deterministic vector field $f\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, we consider the random family $(\nabla u_{\varepsilon})_{\varepsilon>0}$ of unique Lax-Milgram solutions to the rescaled problems

[TABLE]

where the rescaled coefficients $\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})$ vary on the “microscopic” scale $\varepsilon$, and we study “macroscopic” observables of the form $\int_{\mathbb{R}^{d}}g\cdot\nabla u_{\varepsilon}$ with $g\in C_{c}^{\infty}(\mathbb{R}^{d})^{d}$ deterministic. Qualitative homogenization theory [45, 34] states that almost surely $\int_{\mathbb{R}^{d}}g\cdot\nabla u_{\varepsilon}\to\int_{\mathbb{R}^{d}}g\cdot\nabla\bar{u}$ as $\varepsilon\downarrow 0$, where $\nabla\bar{u}$ solves the (deterministic) homogenized equation

[TABLE]

and the homogenized coefficients $\bar{\boldsymbol{a}}\in\mathbb{R}^{d\times d}$ are given by $\bar{\boldsymbol{a}}e_{i}=\mathbb{E}\left[\boldsymbol{a}(\nabla\varphi_{i}+e_{i})\right]$ for $1\leq i\leq d$, in terms of the corrector $\varphi_{i}$, that is, the unique (up to a random additive constant) almost sure solution of the corrector equation,

[TABLE]

in the class of functions the gradient of which is stationary, centered (that is, $\mathbb{E}\left[\nabla\varphi_{i}\right]=0$), and has finite second moment. In other words, the field-flux constitutive relation $\nabla w_{\varepsilon}\mapsto\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})\nabla w_{\varepsilon}$ is replaced on the macroscopic scale by the effective relation $\nabla\bar{w}\mapsto\bar{\boldsymbol{a}}\nabla\bar{w}$.

A well-travelled question in homogenization is to further characterize the oscillations of the solution $u_{\varepsilon}$. As expected from formal asymptotics, the so-called two-scale expansion111We systematically use Einstein’s summation rule on repeated indices. $(1+\varepsilon\varphi_{i}(\tfrac{\cdot}{\varepsilon})\nabla_{i})\bar{u}$ precisely captures the oscillations of $u_{\varepsilon}$ to order $O(\varepsilon)$ (in dimension $d>2$), in the sense of

[TABLE]

This expansion has a natural geometric interpretation: by definition (1.3), the corrector $\varphi_{i}$ defines the correction of Euclidean coordinates $x\mapsto x_{i}$ into $\boldsymbol{a}$-harmonic coordinates $x\mapsto x_{i}+\varphi_{i}(x)$,

[TABLE]

so that the expansion $(1+\varepsilon\varphi_{i}(\tfrac{\cdot}{\varepsilon})\nabla_{i})\bar{u}$ amounts to locally correcting the homogenized solution $\bar{u}$ in terms of oscillating $\boldsymbol{a}$-harmonic coordinates. Such expansions are naturally pursued to higher order: while the (centered) first-order corrector $\varphi=\varphi^{1}$ is characterized by the property that $(1+\varphi_{i}\nabla_{i})\bar{\ell}$ is $\boldsymbol{a}$-harmonic for all affine functions $\bar{\ell}$, the (centered) second-order corrector $\varphi^{2}$ is characterized by the property that for all quadratic polynomials $\bar{q}$ the expansion $(1+\varphi_{i}\nabla_{i}+\varphi_{ij}^{2}\nabla_{ij}^{2})\bar{q}$ fully captures the oscillations imposed by the heterogeneous elliptic operator $-\nabla\cdot\boldsymbol{a}\nabla$ in the sense of $-\nabla\cdot\boldsymbol{a}\nabla(1+\varphi_{i}\nabla_{i}+\varphi_{ij}^{2}\nabla_{ij}^{2})\bar{q}$ being deterministic (see Section 2 for a precise statement). The second-order two-scale expansion $(1+\varepsilon\varphi_{i}(\tfrac{\cdot}{\varepsilon})\nabla_{i}+\varepsilon^{2}\varphi^{2}_{ij}(\tfrac{\cdot}{\varepsilon})\nabla^{2}_{ij})\bar{u}^{2}_{\varepsilon}$ then captures the oscillations of $u_{\varepsilon}$ to order $O(\varepsilon^{2})$ (in dimension $d>4$), where $\bar{u}^{2}_{\varepsilon}$ is a proxy for a solution to the generically ill-posed higher-order homogenized equation

[TABLE]

in terms of the second-order homogenized coefficients \bar{\boldsymbol{a}}^{2}_{i}e_{j}:=\mathbb{E}\big{[}\boldsymbol{a}(\nabla\varphi^{2}_{ij}+\varphi_{i}e_{j})\big{]}. In other words, the effective field-flux constitutive relation $\nabla\bar{w}\mapsto\bar{\boldsymbol{a}}\nabla\bar{w}$ is refined to its second-order version $\nabla\bar{w}\mapsto(\bar{\boldsymbol{a}}+\varepsilon\bar{\boldsymbol{a}}^{2}_{i}\nabla_{i})\nabla\bar{w}$. The proxy $\bar{u}_{\varepsilon}^{2}$ for a second-order homogenized solution can e.g. be chosen in form of $\bar{u}^{2}_{\varepsilon}:=\bar{u}+\varepsilon\tilde{u}^{2}$ with $\tilde{u}^{2}$ solving

[TABLE]

Such a description of oscillations of $u_{\varepsilon}$ via higher-order two-scale expansions is classical in the periodic setting [11]. It also holds in the random setting for large enough dimension and for decaying enough correlations [24, 9] (cf. [22, 20, 3] for the first-order level): under appropriate mixing conditions, higher-order expansions allow to capture oscillations up to the critical order $O(\varepsilon^{d/2})$ (with a logarithmic correction in even dimensions). The limitation is due to the fact that, in the random setting, not all higher-order correctors are well-behaved: $\varphi^{n}$ can be chosen stationary with finite moments only for $n<\frac{d}{2}$. Precise statements and short proofs are included in Section 2 below in the Gaussian setting.

While periodic homogenization boils down to the description of oscillations of $u_{\varepsilon}$, stochastic homogenization in addition means studying the random fluctuations of the macroscopic observables $\int_{\mathbb{R}^{d}}g\cdot\nabla u_{\varepsilon}$. It was recently shown [25] that the centered rescaled observables $\varepsilon^{-d/2}\int_{\mathbb{R}^{d}}g\cdot(\nabla u_{\varepsilon}-\mathbb{E}\left[\nabla u_{\varepsilon}\right])$ converge in law to a Gaussian. Note that this does not enter the above theory of oscillations since accuracy there is precisely limited to order $O(\varepsilon^{d/2})$. We may naturally look for a finer description of this convergence by means of a two-scale expansion. As observed in [25], however, the limiting variance of $\varepsilon^{-d/2}\int_{\mathbb{R}^{d}}g\cdot\nabla u_{\varepsilon}$ generically differs from that of $\varepsilon^{-d/2}\int_{\mathbb{R}^{d}}g\cdot\nabla(1+\varepsilon\varphi_{i}(\tfrac{\cdot}{\varepsilon})\nabla_{i})\bar{u}$, which shows that two-scale expansion cannot be applied naïvely when it comes to fluctuations. In [19], we unravelled the full mechanism behind this observation by means of the so-called homogenization commutator, which led to a new pathwise theory of fluctuations in stochastic homogenization (in agreement with the pathwise heuristics in [25]). More precisely, fluctuations of $\nabla u_{\varepsilon}$ were shown to coincide with fluctuations of a suitable deterministic (Helmholtz) projection of the corresponding homogenization commutator, which behaves like a Gaussian white noise on large scales and for which the two-scale expansion is accurate at leading order. In the present article we explain how this result naturally extends to higher order in parallel to the known theory of oscillations: we obtain a full characterization of fluctuations of $\varepsilon^{-d/2}\int_{\mathbb{R}^{d}}g\cdot\nabla u_{\varepsilon}$ to order $O(\varepsilon^{d/2})$ (with again a logarithmic correction in even dimensions). The main emphasis is on unravelling the suitable higher-order commutator structure.

The homogenization commutator can be viewed as a commutator between large-scale averaging and the field-flux relation; on first-order level it takes the form

[TABLE]

This expression first appeared in [5] (see also [7]) and is particularly natural since $H$-convergence for (1.2) (cf. [39]) is precisely equivalent to weak convergence of $\Xi_{\varepsilon}[\nabla u_{\varepsilon}]$ to [math] (cf. [19]), thus summarizing the whole qualitative homogenization theory. This convergence property is the mathematical formulation of the Hill-Mandel relation in mechanics [27, 28] and can further be made quantitative [19]. For a higher-order theory, we need a higher-order extension of the commutator. In view of the second-order effective field-flux relation $\nabla\bar{w}\mapsto(\bar{\boldsymbol{a}}+\varepsilon\bar{\boldsymbol{a}}^{2}_{i}\nabla_{i})\nabla\bar{w}$, the second-order commutator is naturally defined as

[TABLE]

and a higher-order Hill-Mandel relation can indeed be formulated in these terms (cf. Lemma 3.3 below). Note that $w_{\varepsilon}\mapsto\Xi_{\varepsilon}^{1}[\nabla w_{\varepsilon}]$ and $\Xi^{2}_{\varepsilon}[\nabla w_{\varepsilon}]$ are viewed as first- and second-order differential operators, respectively. Next, we define suitable two-scale expansions of these homogenization commutators. For the first order, we inject the first-order two-scale expansion of $w_{\varepsilon}$ into $\Xi_{\varepsilon}^{1}[\nabla\cdot]$ and truncate the obtained differential operator at first order, thus defining the following first-order standard commutator,

[TABLE]

Similarly, injecting the second-order two-scale expansion of $w_{\varepsilon}$ into $\Xi_{\varepsilon}^{2}[\nabla\cdot]$ and truncating the obtained differential operator at second order, we define the following second-order standard commutator,

[TABLE]

The standard commutators $\bar{w}\mapsto\Xi_{\varepsilon}^{\circ,1}[\nabla\bar{w}]$ and $\Xi_{\varepsilon}^{\circ,2}[\nabla\bar{w}]$ are first- and second-order differential operators, respectively, with $\varepsilon$-rescaled (distributional) stationary random coefficients. These stationary differential operators are viewed as intrinsic quantities, with coefficients given by suitable combinations of correctors and of their gradients. The above definitions are efficiently summarized in the form

[TABLE]

where $T_{x}^{1}\bar{w}$ and $T_{x}^{2}\bar{w}$ denote the first- and second-order Taylor polynomials of $\bar{w}$ at the basepoint $x$, respectively. Note that these quantities are centered: \mathbb{E}\big{[}\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]\big{]}=0 for $n=1,2$ and any smooth deterministic field $\bar{w}$.

Pursuing these constructions to higher order, our main result is threefold and is summarized as follows, under appropriate mixing conditions:

1. (i)

Up to order $O(\varepsilon^{d/2})$ in the fluctuation scaling, the fluctuations of $\nabla u_{\varepsilon}$ are determined by those of the corresponding higher-order commutators; cf. Theorem 1(i) below. 2. (ii)

Up to order $O(\varepsilon^{d/2})$ in the fluctuation scaling (with a logarithmic correction in even dimensions), the two-scale expansions of the commutators are accurate: the fluctuations of the commutators are equivalent to those of the corresponding standard commutators; cf. Theorem 1(ii). This fully extends the pathwise theory of [19] to higher order. 3. (iii)

The standard commutators are approximately local functions of the coefficients $\boldsymbol{a}$. This allows to infer a quantitative characterization of their limiting covariance structure, as well as a quantitative central limit theorem (CLT); cf. Theorem 1(iii)–(iv). Similarly as for the coefficient field $\boldsymbol{a}$ itself, while at first order the scaling limit is given by a Gaussian white noise, corrections converge to derivatives of white noise.

1.2. Main results

For simplicity, we focus on the model setting of a Gaussian coefficient field, in which case a powerful Malliavin calculus is available on the probability space, substantially simplifying the analysis. More precisely, we set

[TABLE]

where $G$ is some $\mathbb{R}^{\kappa}$-valued centered stationary Gaussian random field on $\mathbb{R}^{d}$ constructed on a probability space $(\Omega,\mathcal{F},\mathbb{P})$ and where $a_{0}\in C^{2}_{b}(\mathbb{R}^{\kappa})^{d\times d}$ is such that the ellipticity and boundedness assumptions (1.1) are satisfied. We use scalar notation $a_{0}^{\prime},a_{0}^{\prime\prime}$ for derivatives of $a_{0}$. In addition, we assume that the Gaussian field $G$ has integrable correlations: more precisely, it has the form $G=c_{0}\ast\xi$, where $\xi$ denotes $\kappa$-dimensional white noise on $\mathbb{R}^{d}$ and where the kernel $c_{0}:\mathbb{R}^{d}\to\mathbb{R}^{\kappa\times\kappa}$ is chosen even (that is, $c_{0}(x)=c_{0}(-x)$ for all $x$) and satisfies the integrability condition

[TABLE]

where $B(x)$ is the unit ball centered at $x$. The covariance function of $G$ is then given by $c:=c_{0}\ast c_{0}$ and satisfies $\int_{\mathbb{R}^{d}}(\sup_{B(x)}|c|)\,dx\,\leq\,1$. The case of non-integrable correlations could be treated similarly, but it would lead to different scalings, cf. [16]. Our main result is as follows.

Theorem 1**.**

Consider the above Gaussian setting with integrable correlations (1.5) and set $\ell:=\lceil\frac{d}{2}\rceil$, that is, the smallest integer $\geq\frac{d}{2}$. For all $1\leq n\leq\ell$, let $\bar{u}_{\varepsilon}^{n}$ and $\bar{v}_{\varepsilon}^{n}$ be the $n$th-order homogenized solutions of the primal and dual equations, respectively, as in Definition 2.7 below, and let the $n$th-order homogenization commutator $\Xi_{\varepsilon}^{n}$ and its standard version $\Xi_{\varepsilon}^{\circ,n}$ be as in Definitions 3.1 and 3.4. Then the following hold for all $0<\varepsilon\leq\frac{1}{2}$.

1. (i)

Reduction to commutators:* For all $f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, $1\leq n\leq\ell$, and $p<\infty$,*

[TABLE]

where $N_{p}(X):=\mathbb{E}\left[|X-\mathbb{E}\left[X\right]\!|^{p}\right]^{\frac{1}{p}}$. 2. (ii)

Two-scale expansion of commutators:* For all $f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, $1\leq n\leq\ell$, and $p<\infty$,*

[TABLE]

where

[TABLE] 3. (iii)

Local covariance structure:* Further assume that $\int_{\mathbb{R}^{d}}|y|^{\ell}|c(y)|dy\leq 1$. There exists an *$(m+4)$th-order tensor $\mathcal{Q}_{l^{\prime}k^{\prime}lk}^{m}$ for all $l^{\prime},k^{\prime},l,k\geq 0$ with $l^{\prime}+k^{\prime}+l+k=m\leq\ell-1$, such that for all $g,\nabla\bar{w}\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$ and $1\leq n\leq\ell$, setting

[TABLE]

where $\odot$ denotes the complete contraction of tensors of the same order, there holds

[TABLE] 4. (iv)

Normal approximation:* For all $g,\nabla\bar{w}\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$ and $1\leq n\leq\ell$,*

[TABLE]

where we have set

[TABLE]

where $\operatorname{W}_{2}\left({\cdot};{\mathcal{N}}\right)$ and $\operatorname{d}_{\operatorname{TV}}\left({\cdot};{\mathcal{N}}\right)$ denote the $2$-Wasserstein and the total variation distances to a standard Gaussian law, respectively.

In the estimates above, the notation $\lesssim_{\gamma}$ (with parameters $\gamma$) stands for $\leq$ up to a multiplicative constant $C_{\gamma}\geq 1$ that only depends on $d$, $\lambda$, $\|a_{0}\|_{W^{2,\infty}}$, and on the parameters $\gamma$, through an upper bound on suitable (weighted) Sobolev norms of $\gamma$ if $\gamma$ is a function. ∎

Choosing $n=\ell$, we are thus led to an intrinsic description of the fluctuations of macroscopic observables $\int_{\mathbb{R}^{d}}g\cdot(\nabla u_{\varepsilon}-\mathbb{E}\left[\nabla u_{\varepsilon}\right])$ with accuracy $O(\varepsilon^{d})$ (with a logarithmic correction in even dimensions), that is, up to the square of the CLT order. In particular, non-Gaussian corrections are only expected beyond that order.

Corollary 1**.**

With the assumptions and notation of Theorem 1, for all $f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, in terms of the intrinsic quantity

[TABLE]

there holds for all $p<\infty$,

[TABLE]

and, under the non-degeneracy condition $\mathcal{Q}_{0000}\odot\int_{\mathbb{R}^{d}}(\nabla\bar{v}\otimes\nabla\bar{u})^{\otimes 2}\neq 0$,

[TABLE]

where $\operatorname{d}\left({\cdot};{\cdot}\right)$ stands for the sum of the $2$-Wasserstein and the total variation distances, where $\mathcal{N}$ denotes a standard Gaussian random variable, and where

[TABLE]

Another question that one may formulate concerns the description of fluctuations of higher-order correctors — although this is somewhat artificial since only combinations of correctors in form of two-scale expansions are intrinsically relevant. As briefly described in Remark 3.6, the proof of the above main result can be adapted in a straightforward way to also describe the fluctuations of $\varepsilon^{n-\frac{d}{2}}\varphi^{n}(\tfrac{\cdot}{\varepsilon})$ and $\varepsilon^{n-\frac{d}{2}}\nabla\varphi^{n+1}(\tfrac{\cdot}{\varepsilon})$ with accuracy $O(\varepsilon^{\frac{d}{2}})$ for all $n\leq\ell-1$ (with a logarithmic correction in even dimensions).

Outline of the article

In Section 2, we start by recalling the theory of oscillations in stochastic homogenization in terms of higher-order correctors and two-scale expansions. Section 3 is devoted to the definition and elementary properties of higher-order homogenization commutators, including the (easy) proof of item (i) of Theorem 1 as well as the statement of quantitative higher-order Hill-Mandel relations. Next, we turn to the rest of the proof of Theorem 1, which is based on the following three main ingredients.

1. $\bullet$

Malliavin calculus: In the present Gaussian setting, given a random variable $X$, a Poincaré inequality holds in the form $\mathrm{Var}\left[X\right]\leq C\,\mathbb{E}\left[\|DX\|_{\mathfrak{H}}^{2}\right]$ (cf. [15, 29]), where the Malliavin derivative $DX$ encodes the infinitesimal variation of $X$ with respect to changes in the underlying Gaussian field $G$ and where $\|\cdot\|_{\mathfrak{H}}$ is the Hilbert norm associated with the covariance structure of $G$. Similarly, the approximate normality of $X$ can be estimated by the size of $D^{2}X$ with help of a so-called second-order Poincaré inequality [14, 43]. With these functional analytic tools at hand, we are thus reduced to estimating infinitesimal variations of the quantities of interest with respect to $G$, thus somehow linearizing the dependence on $G$, which is then particularly amenable to PDE methods. On a more technical level, rather than using ad hoc derivatives with respect to the coefficient field as in previous works in stochastic homogenization, we use here the full power of the Gaussian Malliavin calculus. Relevant definitions and results are recalled in Section 4. 2. $\bullet$

Representation formulas: In view of applying the above tools from Malliavin calculus, we consider the infinitesimal variations of the quantities of interest, for which we establish suitable representation formulas. This requires a good understanding of the algebra behind higher-order correctors and commutators, and is performed in Section 5. 3. $\bullet$

PDE ingredients: It remains to estimate the infinitesimal variations of interest from their representation formulas. Based on the large-scale Lipschitz regularity theory for $\boldsymbol{a}$-harmonic functions as developed in [7, 6, 5, 21], we introduce in Section 6 a new annealed Calderón-Zygmund theory for linear elliptic equations in divergence form with random coefficients. This constitutes an upgrade of the quenched large-scale Calderón-Zygmund theory of [4, 3, 21] (see also [31] for a direct approach); it happens to be particularly well suited for our purposes in this article and is of independent interest.

Finally, with these three ingredients at hand, we turn to the proof of items (ii), (iii), and (iv) of Theorem 1 in Sections 7, 8, and 9, respectively.

Notation

1. $\bullet$

We denote by $C\geq 1$ any constant that only depends on $d$, $\lambda$, and $\|a_{0}\|_{W^{2,\infty}}$. We use the notation $\lesssim$ (resp. $\gtrsim$) for $\leq C\times$ (resp. $\geq\frac{1}{C}\times$) up to such a multiplicative constant $C$. We write $\simeq$ when both $\lesssim$ and $\gtrsim$ hold. We add subscripts to $C$, $\lesssim$, $\gtrsim$, $\simeq$ in order to indicate dependence on other parameters. If the subscript is a function, then it is understood as dependence on an upper bound on suitable (weighted) Sobolev norms of the function. 2. $\bullet$

The ball centered at $x$ of radius $r$ in $\mathbb{R}^{d}$ is denoted by $B_{r}(x)$, and we simply write $B(x)$ if $r=1$. 3. $\bullet$

For a function $f$ and $1\leq p<\infty$, we write $[f]_{p}(x):=(\fint_{B(x)}|f|^{p})^{1/p}$ for the local moving $\operatorname{L}^{p}$ average, and similarly $[f]_{\infty}(x):=\sup_{B(x)}|f|$. 4. $\bullet$

For a $k$th-order tensor $T$ we define for any fixed indices $i_{1},\ldots,i_{k}$,

[TABLE]

with $\mathcal{S}_{k}$ denoting the set of permutations of the set $\{1,\ldots,k\}$. Although this notation may seem redundant at this point, it will prove very useful in the sequel. 5. $\bullet$

For $r\in\mathbb{R}$ we denote by $\lceil r\rceil$ the smallest integer $\geq r$ and by $\lfloor r\rfloor$ the largest integer $\leq r$. For $r,s\in\mathbb{R}$ we write $r\vee s:=\max\{r,s\}$ and $r\wedge s:=\min\{r,s\}$. We set $\langle x\rangle:=(1+|x|^{2})^{1/2}$ for $x\in\mathbb{R}^{d}$, and $\langle\nabla\rangle$ denotes the corresponding pseudo-differential operator.

2. Higher-order theory of oscillations

In this section, as a prelude to the study of higher-order fluctuations, we describe the classical higher-order theory of oscillations in stochastic homogenization in terms of two-scale expansions, including some new results and shorter arguments. We start by recalling the definition of higher-order correctors, homogenized coefficients, fluxes, and flux correctors as motivated by formal two-scale expansions (e.g. [11, 30]). A simple iterative way to view this definition is as follows: for $n\geq 1$, assuming that the previous correctors $\varphi^{1},\ldots,\varphi^{n-1}$ are well-defined as centered stationary objects, the next corrector gradient $\nabla\varphi^{n}$ is the unique centered stationary object that is characterized by the property that for all $n$th-order polynomials $\bar{q}$ the corrected polynomial $F^{n}[\bar{q}]:=\bar{q}+\sum_{k=1}^{n}\varphi^{k}_{i_{1}\ldots i_{k}}\nabla_{i_{1}\ldots i_{k}}^{k}\bar{q}$ captures oscillations of the heterogeneous operator in the sense that $\nabla\cdot\boldsymbol{a}\nabla F^{n}[\bar{q}]$ is deterministic. In that case, there must hold (cf. (2.9))

[TABLE]

which in addition characterizes the higher-order homogenized coefficients $\bar{\boldsymbol{a}}^{2},\ldots,\bar{\boldsymbol{a}}^{n-1}$ next to the first-order one $\bar{\boldsymbol{a}}^{1}:=\bar{\boldsymbol{a}}$. In other words, the higher-order corrector expansion $F^{n}[\cdot]$ precisely intertwines the heterogeneous elliptic operator with its higher-order homogenized version. Note that the above only characterizes the symmetric part $\operatorname{Sym}_{i_{1}\ldots i_{n}}\varphi_{i_{1}\ldots i_{n}}$ of correctors, which is indeed the only relevant quantity in view of two-scale expansions. It is however convenient to define correctors with a suitable non-symmetric part in such a way that the corresponding fluxes can be written in terms of skew-symmetric flux correctors, which is a well-known crucial tool in view of error estimates (e.g. [30, p.27]). As stated below, in contrast with the periodic case, only $\ell:=\lceil\frac{d}{2}\rceil$ correctors can be defined, with the $\ell-1$ first of them being stationary; this number would be further reduced in the case of coefficient fields with stronger correlations (e.g. [20, 9]).

Definition 2.1** (Higher-order correctors).**

Given $\ell:=\lceil\frac{d}{2}\rceil$, we inductively define the correctors $(\varphi^{n})_{0\leq n\leq\ell}$, the homogenized coefficients $(\bar{\boldsymbol{a}}^{n})_{1\leq n\leq\ell}$, the fluxes $(q^{n})_{1\leq n\leq\ell}$, and the flux correctors $(\sigma^{n})_{0\leq n\leq\ell}$ as follows:

1. $\bullet$

$\varphi^{0}:=1$ and for all $1\leq n\leq\ell$ we define $\varphi^{n}:=(\varphi^{n}_{i_{1}\ldots i_{n}})_{1\leq i_{1},\ldots,i_{n}\leq d}$ with $\varphi^{n}_{i_{1}\ldots i_{n}}$ a scalar field that satisfies

[TABLE]

For $n<\ell$ we can choose $\varphi^{n}$ itself stationary with $\mathbb{E}\left[\varphi^{n}\right]=0$, while for $n=\ell$ we choose the anchoring $\varphi^{\ell}(0)=0$. 2. $\bullet$

For all $1\leq n\leq\ell$ we define $\bar{\boldsymbol{a}}^{n}:=(\bar{\boldsymbol{a}}^{n}_{i_{1}\ldots i_{n-1}})_{1\leq i_{1},\ldots,i_{n-1}\leq d}$ with $\bar{\boldsymbol{a}}^{n}_{i_{1}\ldots i_{n-1}}$ the matrix given by

[TABLE] 3. $\bullet$

For all $1\leq n\leq\ell$ we define $q^{n}:=(q^{n}_{i_{1}\ldots i_{n}})_{1\leq i_{1},\ldots,i_{n}\leq d}$ with $q^{n}_{i_{1}\ldots i_{n}}$ the vector field given by

[TABLE]

where the definition of $\bar{\boldsymbol{a}}^{n}$ and the normalization of $\sigma^{n-1}$ below ensure $\mathbb{E}\left[q^{n}\right]=0$. 4. $\bullet$

$\sigma^{0}:=0$ and for all $1\leq n\leq\ell$ we define $\sigma^{n}:=(\sigma^{n}_{i_{1}\ldots i_{n}})_{1\leq i_{1},\ldots,i_{n}\leq d}$ with $\sigma^{n}_{i_{1}\ldots i_{n}}$ a skew-symmetric matrix field that satisfies

[TABLE]

with the notation $(\nabla\times X)_{ij}:=\nabla_{i}X_{j}-\nabla_{j}X_{i}$ for a vector field $X$ and with the notation $(\nabla\cdot Y)_{i}:=\nabla_{j}Y_{ij}$ for a matrix field $Y$. For $n<\ell$ we can choose $\sigma^{n}$ itself stationary with $\mathbb{E}\left[\sigma^{n}\right]=0$, while for $n=\ell$ we choose the anchoring $\sigma^{\ell}(0)=0$.∎

We first state that the above definition of higher-order correctors indeed makes sense. This easily follows from [9, proof of Proposition 9 and Lemma 12] together with Lemma 7.1 below in the considered Gaussian setting (see also [21, 20, 3]), and the proof is omitted. Note that the optimal stochastic integrability (that is, the optimal constants $c_{n}$’s below) is not required in our analysis, hence is not addressed here.

Proposition 2.2** (Higher-order correctors [21, 20, 9, 3]).**

All the quantities in Definition 2.1 exist and are uniquely defined. In addition, there exists a sequence $(c_{n})_{0\leq n\leq\ell}$ of positive numbers with $c_{0}=\frac{1}{2}$ such that for all $0\leq n\leq\ell$ and $1\leq p<\infty$,

[TABLE]

where in line with (1.6) we have set

[TABLE]

For later reference we also state the following result on the fluctuation scaling of correctors, showing that higher-order correctors have stronger correlations hence worse fluctuation scalings, as first emphasized in [24].

Lemma 2.3** (Fluctuation scaling of higher-order correctors).**

For all $g\in C^{\infty}_{c}(\mathbb{R}^{d})$, $0\leq n\leq\ell-1$, and $1\leq p<\infty$,

[TABLE]

If $\boldsymbol{a}$ is replaced by the pointwise transpose field $\boldsymbol{a}^{*}$, we write $(\varphi^{*,n})_{0\leq n\leq\ell}$, $(\bar{\boldsymbol{a}}^{*,n})_{1\leq n\leq\ell}$, $(q^{*,n})_{1\leq n\leq\ell}$, and $(\sigma^{*,n})_{0\leq n\leq\ell}$ for the corresponding objects. While it is well-known that $\bar{\boldsymbol{a}}^{*,1}=(\bar{\boldsymbol{a}}^{1})^{*}$, the following lemma extends this relation to higher order.

Lemma 2.4** (Symmetries of homogenized coefficients).**

For all $1\leq n\leq\ell$,222For $n=1$ the notation $\bar{\boldsymbol{a}}^{*,n}_{ji_{1}\ldots i_{n-2}}e_{i_{n-1}}$ stands for $\bar{\boldsymbol{a}}^{*,1}e_{j}$, and we use a similar unifying notation throughout in the sequel.

[TABLE]

In particular, if $\boldsymbol{a}$ is symmetric and if $n$ is even, there holds

[TABLE]

Remark 2.5**.**

The proof of Lemma 2.4 proceeds by letting the successive inverse operators $(\nabla\cdot\boldsymbol{a}\nabla)^{-1}\nabla\cdot$ defining $(\nabla\varphi^{n},\varphi^{n-1})$ “migrate” to the test function $\boldsymbol{a}$ in the definition (2.1) of $\bar{\boldsymbol{a}}^{n}$. Using the algebra of correctors, a full migration precisely leads to the dual expression $\bar{\boldsymbol{a}}^{*,n}$ and proves the relation (2.2). Rather stopping this migration argument at an intermediate step, we obtain for all $1\leq m<n\leq\ell$,

[TABLE]

Note that this turns into (2.2) for $m=n-1$. We deduce for $n=2m+1$ odd,

[TABLE]

and for $n=2m$ even, further using the algebra of correctors to symmetrize the expression,

[TABLE]

This shows that $\bar{\boldsymbol{a}}^{n}$ can be defined in terms of correctors $(\nabla\varphi^{s+1},\varphi^{s})$ with $s\leq\lfloor\frac{n}{2}\rfloor$ only. This simple observation was independently made in [1, Theorem 3.5] (see also [47]) and has the following striking consequence: Choosing a periodization in law $\boldsymbol{a}_{L}$ of the coefficient field $\boldsymbol{a}$ and recalling that the periodic Poincaré inequality ensures that all the correctors $\varphi_{L}^{n}$ are defined as stationary objects so that all the corresponding homogenized coefficients $\bar{\boldsymbol{a}}_{L}^{n}$ are well-defined, the above computation entails that \operatorname{Sym}_{i_{1}\ldots i_{n}}\big{(}\bar{\boldsymbol{a}}^{n}_{L;i_{1}\ldots i_{n-1}}e_{i_{n}}\big{)} has a limit as $L\uparrow\infty$ for all $n\leq d$ if $d$ is odd and for all $n<d$ if $d$ is even. This partially solves a weak version of the Bourgain-Spencer conjecture [12, 33] as stated in [17, Section 4.2]. Refinements in this direction are postponed to a future work. ∎

We turn to the accuracy of two-scale expansions in the homogenization regime. Given $f\in\operatorname{L}^{2}(\mathbb{R}^{d})^{d}$, we consider the unique Lax-Milgram solution $\nabla u_{\varepsilon}$ of the rescaled elliptic PDE (1.2),

[TABLE]

Standard two-scale expansion techniques [11] formally suggest

[TABLE]

where $\bar{U}_{\varepsilon}^{n}$ satisfies the $n$th-order homogenized equation

[TABLE]

The effective field-flux constitutive relation $\nabla\bar{w}\mapsto\bar{\boldsymbol{a}}\nabla\bar{w}$ is thus refined into the higher-order relation

[TABLE]

which includes dispersive corrections. However, equation (2.3) is ill-posed in general for $n$ even and a suitable proxy needs to be devised (cf. also [32]). We start by introducing a notation for two-scale expansions.

Definition 2.6** (Higher-order two-scale expansions).**

For $0\leq n\leq\ell$, given a deterministic smooth function $\bar{w}$, its $n$th-order two-scale expansion $F^{n}_{\varepsilon}[\bar{w}]$ is defined as

[TABLE]

At the level of gradients, we similarly define for $0\leq n\leq\ell$,

[TABLE]

In particular, note that the quantities $\nabla F^{n}_{\varepsilon}[\bar{w}]$ and $E^{n}_{\varepsilon}[\nabla\bar{w}]$ are related as follows: for all $0\leq n\leq\ell$,

[TABLE]

where all coefficients of $E^{n}_{\varepsilon}[\nabla\bar{w}]$ are stationary. We turn to the definition of a suitable proxy for the ill-posed higher-order homogenized equation (2.3). For later reference we simultaneously consider the dual equations.

Definition 2.7** (Higher-order homogenized equations).**

Given $f,g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, let $\nabla u_{\varepsilon}$ and $\nabla v_{\varepsilon}$ be the unique Lax-Milgram solutions of the rescaled elliptic PDEs

[TABLE]

For $1\leq n\leq\ell$, we define the corresponding $n$th-order homogenized solutions

[TABLE]

where $\nabla\tilde{u}^{1}=\nabla\bar{u}$ and $\nabla\tilde{v}^{1}=\nabla\bar{v}$ are the unique Lax-Milgram solutions in $\mathbb{R}^{d}$ of the first-order homogenized equations,

[TABLE]

and where for $2\leq n\leq\ell$ we inductively define the $n$th-order corrections $\nabla\tilde{u}^{n}$ and $\nabla\tilde{v}^{n}$ as the unique Lax-Milgram solutions in $\mathbb{R}^{d}$ of

[TABLE]

With the above notation, the main result from the higher-order corrector theory in stochastic homogenization takes the following guise. For $q=2$, it is a straightforward consequence of the definition of higher-order correctors and of the corrector estimates of Proposition 2.2. It extends [22, 20, 9, 3] to higher order and completes [24, Theorem 1.6] up to the optimal order with finer norms. For $n=\ell$, the obtained maximal accuracy is $O(\varepsilon^{d/2})$ (with a logarithmic correction in even dimensions), indicating that the study of fluctuations goes beyond the theory of oscillations.

Proposition 2.8** (Accuracy of higher-order two-scale expansions).**

For all $f\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, $1\leq n\leq\ell$, and $1<q\leq p<\infty$,

[TABLE]

2.1. Proof of Lemma 2.3

The result is easily obtained based on Malliavin calculus and on the annealed Calderón-Zygmund theory of Section 6. More precisely, for $0\leq n\leq\ell-1$, we start from Proposition 4.1(iii) together with (4.2) in the form

[TABLE]

and the conclusion follows after $\varepsilon$-rescaling from the estimates of Lemma 7.1 on the Malliavin derivatives of higher-order correctors. ∎

2.2. Proof of Lemma 2.4

As explained in Remark 2.5, the argument consists in letting the successive inverse operators $(-\nabla\cdot\boldsymbol{a}\nabla)^{-1}\nabla\cdot$ defining the correctors $(\nabla\varphi^{n},\varphi^{n-1})$ migrate to the test function $\boldsymbol{a}$ in the definition of the homogenized coefficient

[TABLE]

Taking into account the particular algebraic structure of the correctors, a full migration precisely leads to the dual expression $\bar{\boldsymbol{a}}^{*,n}$ up to additional terms involving higher-order flux correctors, which disappear when taking the symmetric part. We split the proof into two steps.

Step 1. Migration process: proof that for all $1\leq m<n\leq\ell$,

[TABLE]

We abundantly use the properties of higher-order correctors (cf. Definition 2.1). The equation for $\varphi^{n}$, in conjunction with the stationarity of all involved objects, and the skew-symmetry of $\sigma^{n-1}_{i_{1}\ldots i_{n-1}}$ yield

[TABLE]

hence, using the equation $\nabla\cdot\sigma^{n-1}=q^{n-1}$ and the choice $\mathbb{E}\left[\varphi^{*,m}\right]=0$,

[TABLE]

We focus on rewriting the second right-hand side term. Using the equation for $\varphi^{*,m+1}$ and the skew-symmetry of $\sigma^{*,m}_{j_{1}\ldots j_{m}}$, we find

[TABLE]

hence, using the equation for $\nabla\cdot\sigma^{*,m}$ and the choice $\mathbb{E}\left[\varphi^{n-1}\right]=0$,

[TABLE]

and the claim (2.8) follows.

Step 2. Conclusion.

For $1\leq n\leq\ell$, the definition of $\bar{\boldsymbol{a}}^{n}$ and the equation for $\varphi^{*,1}$ yield

[TABLE]

Iterating identity (2.8) then leads to

[TABLE]

By skew-symmetry of $\sigma^{*,l}_{ji_{n}\ldots i_{n-l+2}}$ and $\sigma^{n-l-2}_{i_{1}\ldots i_{n-l-2}}$ for all $l$, we deduce

[TABLE]

hence, using the equation for $\varphi^{1}$ and the definition of $\bar{\boldsymbol{a}}^{*,n}$,

[TABLE]

and the conclusion follows.∎

2.3. Proof of Proposition 2.8

We focus on the case $\varepsilon=1$ and drop it from all subscripts in the notation, while the final result is obtained after $\varepsilon$-rescaling. We split the proof into two steps.

Step 1. Equation for $F_{n}[\bar{w}]$: proof that for all $\bar{w}\in C^{\infty}_{c}(\mathbb{R}^{d})$ and $n\geq 0$,

[TABLE]

We argue by induction. The claim is obvious for $n=0$. Now, if it holds for some $n\geq 0$, we deduce

[TABLE]

The definition of $\sigma_{i_{1}\ldots i_{n+1}}^{n+1}$ (cf. Definition 2.1) yields

[TABLE]

Hence, using the skew-symmetry of $\sigma_{i_{1}\ldots i_{n+1}}^{n+1}$ and decomposing

[TABLE]

we obtain

[TABLE]

Injecting this into (2.10) leads to the claim (2.9) at level $n+1$.

Step 2. Conclusion.

Let $n\geq 1$ be fixed. Applying (2.9) with $\bar{w}:=\bar{u}^{n}$ and subtracting the equation (2.6) for $u$, we obtain

[TABLE]

We now examine the equation satisfied by $\bar{u}^{n}=\sum_{k=1}^{n}\tilde{u}^{k}$: summing the defining equations for $(\tilde{u}^{k})_{1\leq k\leq n}$ (cf. Defintion 2.7), we find

[TABLE]

or equivalently,

[TABLE]

(This is to be compared with the formal ill-posed equation (2.3).) Inserting this into (2.11) yields

[TABLE]

An energy estimate and the corrector estimates of Proposition 2.2 then imply for all $p<\infty$,

[TABLE]

where the last estimate follows from the weighted Calderón-Zygmund theory for the constant-coefficient equations defining $(\tilde{u}^{k})_{1\leq k\leq n}$. The conclusion (2.7) in $\operatorname{L}^{2}(\mathbb{R}^{d})$ follows after $\varepsilon$-rescaling. (Note that for the sake of shortness in the statement the norms of $\nabla^{k}f$ for $n\leq k\leq 2n-1$ are replaced by the norm of $\langle\nabla\rangle^{2n-1}f$, which is no longer scale invariant.) The corresponding result in $\operatorname{L}^{q}(\mathbb{R}^{d})$ for general $1<q<\infty$ requires to replace the energy estimate for (2.12) by some $\operatorname{L}^{q}$ theory, as provided e.g. by Theorem 6.1 below; details are omitted. ∎

3. Higher-order homogenization commutators

As discovered in [5, 19], the homogenization commutator $\Xi_{\varepsilon}^{1}[\nabla u_{\varepsilon}]:=(\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})-\bar{\boldsymbol{a}}^{1})\nabla u_{\varepsilon}$, which takes the form of a commutator between large-scale averaging and the field-flux constitutive relation, plays a key role in the study of fluctuations in stochastic homogenization. Indeed, while the two-scale expansion of the solution $\nabla u_{\varepsilon}$ is not accurate in the fluctuation scaling [25], the expansion of its commutator $\Xi_{\varepsilon}^{1}[\nabla u_{\varepsilon}]$ is accurate [19]. The leading-order fluctuations are then governed by the standard homogenization commutator, that is, the commutator of the corrector field, $\Xi_{\varepsilon}^{\circ,1}[\nabla\bar{u}^{1}]:=(\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})-\bar{\boldsymbol{a}}^{1})(\nabla\varphi_{i}^{1}(\tfrac{\cdot}{\varepsilon})+e_{i})\nabla_{i}\bar{u}^{1}$. For higher-order fluctuations, a suitable higher-order correction of the homogenization commutator $\Xi_{\varepsilon}^{1}[\nabla u_{\varepsilon}]$ is defined as follows.

Definition 3.1** (Higher-order homogenization commutators).**

For $0\leq n\leq\ell$, given a random function $w$, its $n$th-order homogenization commutator $\Xi^{n}_{\varepsilon}[\nabla w]$ is the (distributional) random vector field

[TABLE]

This definition is natural in view of the higher-order effective field-flux constitutive relation (2.4). Note that the symmetries of the homogenized coefficients (cf. Lemma 2.4) and an integration by parts lead to the following duality relation for all smooth and compactly supported random functions $w,w^{\prime}$,

[TABLE]

Similarly as for the first-order commutator [19, (1.8)], we state that the higher-order fluctuations of the field $\nabla u_{\varepsilon}$ are determined by those of the higher-order commutators. While at first order we have the exact relation $\int_{\mathbb{R}^{d}}g\cdot(\nabla u_{\varepsilon}-\nabla\bar{u})=\int_{\mathbb{R}^{d}}\nabla\bar{v}\cdot\Xi_{\varepsilon}^{1}[\nabla u_{\varepsilon}]$, this does no longer hold at higher orders due to the choice of a proxy for solutions of the generically ill-posed higher-order homogenized equations (2.3) (cf. (3.4) below). A corresponding formula can be deduced for the flux $\boldsymbol{a}(\tfrac{\cdot}{\varepsilon})\nabla u_{\varepsilon}$. This key principle corresponds to item (i) in Theorem 1.

Proposition 3.2** (Reduction to commutators).**

For all $1\leq n\leq\ell$ and $p<\infty$,

[TABLE]

where we recall the notation $N_{p}(X)=\mathbb{E}\left[|X-\mathbb{E}\left[X\right]|^{p}\right]^{\frac{1}{p}}$. ∎

Next, we wish to emphasize that the above definition of the higher-order commutator $\Xi_{\varepsilon}^{n}$ leads to a higher-order Hill-Mandel relation (e.g. [51] in the mechanics literature). First recall that the first-order commutator $\Xi^{1}_{\varepsilon}[\nabla u_{\varepsilon}]:=(\boldsymbol{a}-\bar{\boldsymbol{a}})(\tfrac{\cdot}{\varepsilon})\nabla u_{\varepsilon}$ was a natural quantity to consider in [5, 19] since the very definition of $H$-convergence [39] is equivalent to the weak convergence of $\Xi^{1}_{\varepsilon}[\nabla u_{\varepsilon}]$ to [math], which is the mathematical formulation of the classical Hill-Mandel relation [27, 28]. Interestingly, this can be made quantitative: a direct computation for the solution of (2.6) yields

[TABLE]

where the right-hand side is checked to be of order $O(\varepsilon)$ in a weak norm if $d>2$. This property extends to higher order as follows.

Lemma 3.3** (Higher-order Hill-Mandel relation).**

For all $f\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$ and $1\leq n\leq\ell$,

[TABLE]

hence, for all $g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$ and $p<\infty$,

[TABLE]

Finally, the relevant higher-order two-scale expansion of the higher-order commutator $\Xi^{n}_{\varepsilon}[\nabla u_{\varepsilon}]$ takes the form $\Xi^{\circ,n}_{\varepsilon}[\nabla\bar{u}^{n}_{\varepsilon}]$, where the so-called standard commutator $\Xi^{\circ,n}_{\varepsilon}[\nabla\cdot]$ is the $n$th-order differential operator with $\varepsilon$-rescaled (distributional) stationary random coefficients that is obtained by inserting the $n$th-order two-scale expansion $F_{\varepsilon}^{n}[\cdot]$ into the commutator $\Xi^{n}_{\varepsilon}[\nabla\cdot]$ and by truncating the obtained differential operator at order $n$. In other words, the standard commutator $\Xi^{\circ,n}_{\varepsilon}[\nabla\cdot]$ is the $n$th-order linear differential operator characterized by $\Xi^{\circ,n}_{\varepsilon}[\nabla\bar{q}]=\Xi^{n}_{\varepsilon}[\nabla F_{\varepsilon}^{n}[\bar{q}]]$ for all $n$th-order polynomials $\bar{q}$. This is efficiently expressed as follows.

Definition 3.4** (Standard higher-order homogenization commutators).**

For $0\leq n\leq\ell$, given a smooth deterministic function $\bar{w}$, we consider its $n$th-order Taylor polynomial with basepoint $x$,

[TABLE]

where we use standard multi-index notation, and the $n$th-order standard homogenization commutator $\Xi^{\circ,n}_{\varepsilon}[\nabla\bar{w}]$ is then defined as the (distributional) random vector field

[TABLE]

(Henceforth by a slight abuse of notation we similarly define $\Xi^{n}[H]$ as in Definition 3.1 even when $H$ is not a gradient field.) ∎

This definition ensures $\mathbb{E}\left[\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]\right]=0$ for all $0\leq n\leq\ell$ and $\bar{w}\in C^{\infty}_{c}(\mathbb{R}^{d})$. Note that $\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]$ in general differs from $\Xi_{\varepsilon}^{n}[E_{\varepsilon}^{n}[\nabla\bar{w}]]$ whenever $n>1$, even though $E^{n}_{\varepsilon}[\nabla T^{n}_{x}\bar{w}](x)=E^{n}_{\varepsilon}[\nabla\bar{w}](x)$. An explicit formula is as follows.

Lemma 3.5** (Explicit formula for $\Xi^{\circ,n}_{\varepsilon}$).**

For all $\bar{w}\in~C^{\infty}_{c}(\mathbb{R}^{d})$ and $1\leq n\leq\ell$,

[TABLE]

Remark 3.6**.**

While Proposition 3.2 reduces higher-order fluctuations of the field $\nabla u_{\varepsilon}$ to corresponding fluctuations of its higher-order homogenization commutator, we briefly argue that fluctuations of higher-order correctors can similarly be reduced to fluctuations of higher-order standard commutators. A simple adaptation of the proof of Proposition 3.2 yields for $0\leq n\leq\ell-1$,

[TABLE]

Inserting the formula of Lemma (3.5) and the definition of $E_{\varepsilon}^{n}$ (cf. Definition 2.6), the expression on the left-hand side can, after straightforward computations, be rewritten as follows,

[TABLE]

in terms of

[TABLE]

If we wish to describe only first-order fluctuations, we may replace $T_{\varepsilon}^{n,s}(f,g)$ by $g\nabla^{s+1}\bar{w}$ and $\nabla\bar{v}_{\varepsilon}^{n}$ by $\nabla\bar{v}$, and use Lemma 2.3 to estimate the remainder, to the effect of

[TABLE]

hence, recalling Definition 2.6 and taking differences, for all $0\leq n\leq\ell-1$,

[TABLE]

This characterizes the first-order fluctuations of $\operatorname{Sym}_{j_{1}\ldots j_{n}}(\nabla\varphi^{n}_{j_{1}\ldots j_{n}}+\varphi^{n-1}_{j_{1}\ldots j_{n-1}}e_{j_{n}})$ in terms of (differences of) standard commutators. Note that fluctuations of the latter are easily jointly characterized by repeating the analysis of Sections 8–9, and are seen to converge to suitable linear combinations of derivatives of Gaussian white noise. From here we may inductively infer first-order fluctuations of $\operatorname{Sym}_{j_{1}\ldots j_{n}}\varphi^{n}_{j_{1}\ldots j_{n}}$ for all $n\leq\ell-1$, while corresponding higher-order fluctuations are similarly extracted from (3.3). ∎

3.1. Proof of Proposition 3.2

We focus on the case $\varepsilon=1$ and drop it from all subscripts in the notation, while the final result is obtained after $\varepsilon$-rescaling. We split the proof into three steps.

Step 1. Proof of

[TABLE]

Testing equation $-\nabla\cdot\bar{\boldsymbol{a}}^{*,1}\nabla\bar{v}^{1}=\nabla\cdot g$ with $u$ and testing equation $-\nabla\cdot\boldsymbol{a}\nabla u=\nabla\cdot f$ with $\bar{v}^{n}$, we find

[TABLE]

Noting that the definition of $\Xi^{n}[\nabla u]$ (cf. Definition 3.1) together with (3.1) yields after integration by parts

[TABLE]

we deduce

[TABLE]

Decomposing $\bar{v}^{n}:=\sum_{j=1}^{n}\tilde{v}^{j}$ and using the defining equations for $(\tilde{v}^{j})_{2\leq j\leq n}$ (cf. Definition 2.7) in the resummed form

[TABLE]

the claim follows.

Step 2. CLT scaling: for all $g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$,

[TABLE]

This result is classical [37, 20, 18]. A short proof based on Malliavin calculus and on the annealed Calderón-Zygmund theory of Section 6 is easily obtained as in (7.9); details are omitted.

Step 3. Conclusion.

Since $\int_{\mathbb{R}^{d}}\nabla\bar{v}^{n}\cdot f$ is deterministic, the identity of Step 1 yields

[TABLE]

Applying the CLT scaling result of Step 2 and the Calderón-Zygmund theory for the constant-coefficient equations defining $(\tilde{v}^{k})_{1\leq k\leq n}$, this expression is estimated by

[TABLE]

and the conclusion follows after $\varepsilon$-rescaling. ∎

3.2. Proof of Lemma 3.3

We focus on the case $\varepsilon=1$ and drop it from all subscripts in the notation, while the final result is obtained after $\varepsilon$-rescaling. We split the proof into two steps.

Step 1. Proof that for all $1\leq n\leq\ell$,

[TABLE]

which coincides with (3.2) after $\varepsilon$-rescaling.

First, by Lemma 2.4, the $n$th-order commutator can alternatively be written as

[TABLE]

It remains to show that for all $0\leq n\leq\ell$,

[TABLE]

while the result (3.5) indeed follows after injecting the equation $-\nabla\cdot\boldsymbol{a}\nabla u=\nabla\cdot f$ and writing $\varphi^{*,k}_{ji_{1}\ldots i_{k-1}}\nabla\cdot f=\nabla\cdot(\varphi^{*,k}_{ji_{1}\ldots i_{k-1}}f)-\nabla\varphi^{*,k}_{ji_{1}\ldots i_{k-1}}\cdot f$. We prove (3.7) by induction. It is obvious for $n=0$ (recall that for $n=0$ the notation $\boldsymbol{a}^{*}\varphi_{ji_{1}\ldots i_{n-1}}^{*,n}e_{i_{n}}$ stands for $\boldsymbol{a}^{*}e_{j}$). Assume that it holds for some $n\geq 0$. In order to deduce it at level $n+1$, it suffices to prove that

[TABLE]

The equation for $\nabla\cdot\sigma^{*,n+1}$ (cf. Definition 2.1) yields

[TABLE]

and the claim (3.8) follows from the skew-symmetry of $\sigma^{*,n+1}_{ji_{1}\ldots i_{n}}$.

Step 2. Conclusion.

For $g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, combining (3.5) with the energy estimate $\|\nabla u\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}\lesssim\|f\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}$ and with the corrector estimates of Proposition 2.2, we find for all $p<\infty$,

[TABLE]

In order to optimally estimate the last contribution, we must exploit stochastic cancellations in form of Lemma 2.3, which yields after $\varepsilon$-rescaling, for all $1\leq n\leq\ell$,

[TABLE]

Since $\ell-1\leq\frac{1}{2}(d-1)$, using Hölder’s and Sobolev’s inequalities in the form

[TABLE]

the conclusion follows. ∎

3.3. Proof of Lemma 3.5

First note that, given a linear differential operator $O$ of order $m$, say $O[\bar{w}]=\sum_{k=0}^{m}P^{k}_{j_{1}\ldots j_{k}}(x)\nabla^{k}_{j_{1}\ldots j_{k}}\bar{w}$, the characterizing property of the Taylor polynomial yields

[TABLE]

which amounts to truncating the differential operator $O$ at order $n$. Now observe that $\Xi_{\varepsilon}^{n}[E_{\varepsilon}^{n}[\nabla\cdot]]$ is a linear differential operator of order $2n-1$ with (distributional) stationary random coefficients. Composing the explicit definitions of $E_{\varepsilon}^{n}$ and $\Xi_{\varepsilon}^{n}$ (cf. Definitions 2.6 and 3.1), we obtain

[TABLE]

Hence, relabelling $k-1$ as $k$ in the second right-hand side contribution, using the symmetry in $i_{1},\ldots,i_{k}$, and expanding the multiple derivative $\nabla^{k}_{i_{1}\ldots i_{k}}$ by using the general Leibniz rule,

[TABLE]

Now applying (3.9) to the definition of $\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]$ (cf. Definition 3.4) together with the above formula for $\Xi_{\varepsilon}^{n}[E_{\varepsilon}^{n}[\nabla\bar{w}]]$, we are led to

[TABLE]

Comparing this with (3.10) yields the conclusion.∎

4. Reminder on Malliavin calculus

In this section, we recall some basic definitions of the Malliavin calculus with respect to the Gaussian field $G$ (e.g. [36, 44, 42] for details). For simplicity, we focus on the explicit description of the Malliavin derivative on a simple dense subspace of random variables (cf. $\mathcal{R}$ hereafter), where it is seen to coincide with the formal $\operatorname{L}^{2}$-gradient with respect to $G$. Next, we give a precise statement of the different tools from Malliavin calculus that will be used in the sequel. We emphasize that the reader should not be intimidated by the somehow abstract appearance of the theory: for our application to stochastic homogenization all (averaged) quantities of interest happen to be Malliavin smooth, and their Malliavin derivatives are systematically computed by duality arguments so that in the end everything is expressed in terms of the (explicit) Malliavin derivative of the coefficient field $\boldsymbol{a}$ itself.

The Gaussian random field $G$ can be seen as a random Schwartz distribution, that is, as a random element in $\mathcal{S}^{\prime}(\mathbb{R}^{d})^{\kappa}$. Indeed, for all $\zeta_{1},\zeta_{2}\in C^{\infty}_{c}(\mathbb{R}^{d})^{\kappa}$, we define $G(\zeta_{1})$, $G(\zeta_{2})$ (or formally $\int_{\mathbb{R}^{d}}G\zeta_{1}$, $\int_{\mathbb{R}^{d}}G\zeta_{2}$) as centered Gaussian random variables with covariance

[TABLE]

We define $\mathfrak{H}$ as the closure of $C^{\infty}_{c}(\mathbb{R}^{d})^{\kappa}$ for the (semi)norm

[TABLE]

The space $\mathfrak{H}$ (up to taking the quotient with respect to the kernel of $\|\cdot\|_{\mathfrak{H}}$) is a separable Hilbert space. Under the integrability condition (1.5), it obviously contains $\operatorname{L}^{2}(\mathbb{R}^{d})^{\kappa}$ and there holds

[TABLE]

In view of the isometry relation $\operatorname{Cov}\left[{G(\zeta_{1})};{G(\zeta_{2})}\right]=\langle\zeta_{1},\zeta_{2}\rangle_{\mathfrak{H}}$, the random field $G$ is said to be an isonormal Gaussian process over $\mathfrak{H}$.

Without loss of generality, we work under the minimality assumption $\mathcal{F}=\sigma(G)$, which implies that the linear subspace

[TABLE]

is dense in $\operatorname{L}^{2}(\Omega)$. This allows to define operators and prove properties on the simpler subspace $\mathcal{R}$ before extending them to $\operatorname{L}^{2}(\Omega)$ by density. For $r\geq 1$ we similarly define

[TABLE]

which is dense in $\operatorname{L}^{2}(\Omega;\mathfrak{H}^{\otimes r})$. For a random variable $X\in\mathcal{R}$, say $X=g(G(\zeta_{1}),\ldots,G(\zeta_{n}))$, we define its Malliavin derivative $DX\in\operatorname{L}^{2}(\Omega;\mathfrak{H})$ as

[TABLE]

which coincides with the formal $\operatorname{L}^{2}$ gradient of $X$ with respect to $G$. For an element $X\in\mathcal{R}(\mathfrak{H}^{\otimes r})$ with $r\geq 1$, say $X=\sum_{i=1}^{n}\psi_{i}X_{i}$, the Malliavin derivative $DX\in\operatorname{L}^{2}(\Omega;\mathfrak{H}^{\otimes(r+1)})$ is defined as $DX=\sum_{i=1}^{n}\psi_{i}\otimes DX_{i}$. For $j\geq 1$, we iteratively define the $j$th-order Malliavin derivative $D^{j}:\mathcal{R}(\mathfrak{H}^{\otimes r})\to\operatorname{L}^{2}(\Omega;\mathfrak{H}^{\otimes(r+j)})$ for all $r\geq 0$. For all $r,m\geq 0$, we then set

[TABLE]

we define the Malliavin-Sobolev space $\mathbb{D}^{m,2}(\mathfrak{H}^{\otimes r})$ as the closure of $\mathcal{R}(\mathfrak{H}^{\otimes r})$ for this norm, and we extend the Malliavin derivatives $D^{j}$ by density to these spaces.

Next, we define a divergence operator $D^{*}$ as the adjoint of the Malliavin derivative $D$ and we start with its domain: for $r\geq 1$, we set

[TABLE]

and for all $X\in{\operatorname{dom}}\,(D^{*}|_{\operatorname{L}^{2}(\Omega;\mathfrak{H}^{\otimes r})})$ we define the divergence $D^{*}X\in\operatorname{L}^{2}(\Omega;\mathfrak{H}^{\otimes(r-1)})$ by duality via the relation

[TABLE]

One can show that ${\operatorname{dom}}\,(D^{*}|_{\operatorname{L}^{2}(\Omega;\mathfrak{H}^{\otimes r})})\supset\mathbb{D}^{1,2}(\mathfrak{H}^{\otimes r})$ (e.g. [44, Proposition 1.3.1]). Combining the Malliavin derivative and the associated divergence, we construct the so-called Ornstein-Uhlenbeck operator (or infinite-dimensional Laplacian)

[TABLE]

which is an essentially self-adjoint positive operator. The explicit actions of $D^{*}$ and $\mathcal{L}$ on $\mathcal{R}$ are easily computed (e.g. [42, p.34]): for $X\in\mathcal{R}(\mathfrak{H}^{\otimes r})$, say $X=\sum_{i=1}^{n}\psi_{i}X_{i}$, we have $\mathcal{L}X=\sum_{i=1}^{n}\psi_{i}\,\mathcal{L}X_{i}$, while for $X_{i}=g(G(\zeta_{1}),\ldots,G(\zeta_{n}))$,

[TABLE]

In particular, a direct computation (e.g. [42, p.35]) leads to the crucial commutator relation

[TABLE]

In addition, one checks that $\mathcal{L}$ is shift-invariant. While the positivity of $\mathcal{L}$ ensures that the inverse operator $(1+\mathcal{L})^{-1}$ is contractive on $\operatorname{L}^{2}(\Omega)$, it is easily checked to be also contractive on $\operatorname{L}^{p}(\Omega)$ for all $p>1$ (cf. e.g. [38, Proposition 3.2]). Although we focus here on the explicit description of operators on their core $\mathcal{R}$, the action of $D,D^{*},\mathcal{L}$ on general $\operatorname{L}^{2}$ random variables are naturally described in terms of their Wiener chaos expansion (e.g. [44, 42]).

Based on the above definitions, we state the following proposition collecting various useful results for the fine analysis of functionals of the Gaussian field $G$. Item (i) is classical [40, 15, 29], as well as item (iii) (e.g. [8] and references therein). Item (ii) is well-known in the discrete Gaussian setting [26] (cf. also [50, 38]). Item (iv) in total variation distance (and similarly in $1$-Wasserstein distance) is a consequence of Stein’s method: it was first obtained in the discrete setting by Chatterjee [14], while the present Malliavin analogue is due to [41, 43] (of which we give a slightly different formulation here). The corresponding result in $2$-Wasserstein distance is different and is due to [35]. Note the particular role of the so-called Stein kernel $\langle DX,(1+\mathcal{L})^{-1}DX\rangle_{\mathfrak{H}}$ in items (ii) and (iv). A short proof is included in Appendix A for the reader’s convenience.

Proposition 4.1** ([29, 26, 14, 41, 43, 35]).**

**

1. (i)

First-order Poincaré inequality:* For all $X\in\operatorname{L}^{2}(\Omega)$,*

[TABLE] 2. (ii)

Helffer-Sjöstrand identity:* For all $X,Y\in\mathbb{D}^{1,2}(\Omega)$,*

[TABLE] 3. (iii)

Logarithmic Sobolev inequality:* For all $X\in\operatorname{L}^{2}(\Omega)$,*

[TABLE]

hence in particular, for all $p\geq 1$,

[TABLE] 4. (iv)

Second-order Poincaré inequality:* For all $X\in\operatorname{L}^{2}(\Omega)$ with $\mathbb{E}\left[X\right]=0$ and $\mathrm{Var}\left[X\right]=1$,*

[TABLE]

where $\operatorname{W}_{2}\left({\cdot};{\mathcal{N}}\right)$ and $\operatorname{d}_{\operatorname{TV}}\left({\cdot};{\mathcal{N}}\right)$ denote the $2$-Wasserstein and the total variation distance to a standard Gaussian law, respectively, and where the operator norm of $D^{2}X$ is defined by

[TABLE]

Remark 4.2**.**

Malliavin-type calculus has also been successfully developed in some non-Gaussian settings, to which most of our conclusions could thus be extended. For instance, in the case of an i.i.d. discrete random field $G:=(G_{x})_{x\in\mathbb{Z}^{d}}$, a corresponding calculus is based on the so-called Glauber derivative $D_{x}X=\mathbb{E}\left[\left.X\,\right|\,G_{x}\right]-X$ (e.g. [19]). Another important example is that of a Poisson point process $\xi$ on $\mathbb{R}^{d}$, in which case the Malliavin derivative is chosen as the difference operator $D_{x}X(\xi)=X(\xi\cup\{x\})-X(\xi\setminus\{x\})$ (e.g. [46]). ∎

5. Representation formulas

As recalled in Proposition 4.1 above, tools from Malliavin calculus allow to linearize the dependence on the underlying Gaussian field $G$ and reduce various subtle questions about quantities of interest to the estimation of their Malliavin derivative, that is, of their local dependence with respect to $G$. In this spirit, the key ingredient for items (ii) and (iii) of Theorem 1 consists in suitable representation formulas for the Malliavin derivative of higher-order commutators.

1. $\bullet$

Locality of higher-order commutators:

The key property of homogenization commutators for fluctuations is their improved locality with respect to the coefficient field. While this was first discovered at first order in [19] (see also [25, 5, 23]), formula (5.1) below indicates that the higher-order commutator $\Xi^{n}_{\varepsilon}[\nabla u_{\varepsilon}]$ behaves like a local quantity at the level of its Malliavin derivative up to a higher-order error $O(\varepsilon^{n})$. More precisely, the formula takes the form

[TABLE]

where the main right-hand side terms are indeed local in view of $D_{z}\boldsymbol{a}=a_{0}^{\prime}(G(z))\,\delta(\cdot-z)$, and where the error can be checked to be of order $O(\varepsilon^{n})$ in $\operatorname{L}^{2}(\Omega;\mathfrak{H})$ when integrated with a smooth test function. While the full formula (5.1) below is given for simplicity for $\varepsilon=1$, the $\varepsilon$-scaling is transparent by counting the number of derivatives that fall on test functions. The proof consists in exploiting the algebraic structure of higher-order correctors and commutators in order to make this clean local structure appear.

Since the standard commutator $\Xi^{\circ,n}_{\varepsilon}[\nabla\bar{w}]$ is obtained by applying $\Xi_{\varepsilon}^{n}[\nabla\cdot]$ to the two-scale expansion $F_{\varepsilon}^{n}[\nabla\bar{w}]$ whenever $\bar{w}$ is an $n$th-order polynomial, we may deduce a corresponding formula (5.2) for its Malliavin derivative,

[TABLE]

In view of Proposition 4.1(ii), this improved locality property provides a simple access to the higher-order characterization of the asymptotic covariance structure of large-scale averages of commutators, cf. Section 8. 2. $\bullet$

Accuracy of the two-scale expansion of higher-order commutators:

Comparing the formulas for Malliavin derivatives of the commutator $\Xi^{n}_{\varepsilon}[\nabla u_{\varepsilon}]$ and of its standard version $\Xi^{\circ,n}_{\varepsilon}[\nabla\bar{w}]$, and choosing $\bar{w}=\bar{u}_{\varepsilon}^{n}$ the homogenized solution, we find for the difference,

[TABLE]

where we exactly recognize the higher-order two-scale expansion error $\nabla u_{\varepsilon}-E_{\varepsilon}^{n}[\nabla\bar{u}_{\varepsilon}^{n}]$, which is known to be of small order $O(\varepsilon^{n})$ in view of Proposition 2.8. This constitutes the core of the proof of the higher-order pathwise result, cf. Section 7.

We turn to the precise statement of the representation formulas. These are understood in a weak sense: when integrated in space with a smooth test function, the quantities in consideration are indeed Malliavin smooth and the formulas make sense after integrations by parts (see also Lemmas 7.1 and 9.1).

Proposition 5.1** (Formulas for first Malliavin derivatives).**

Let $\varepsilon=1$ and drop it from all subscripts in the notation. For all $1\leq n\leq\ell$ there holds for any random function $w$ with $\nabla\cdot\boldsymbol{a}\nabla w$ deterministic,

[TABLE]

and also, for any smooth deterministic function $\bar{w}$,

[TABLE]

In particular,

[TABLE]

∎

In view of the approximate normality result of Theorem 1(iv), the second-order Poincaré inequality of Proposition 4.1(iv) requires to further estimate the second Malliavin derivative of higher-order standard commutators, and we establish the following formula. Although some of the terms are proportional to $D^{2}\boldsymbol{a}$, hence display an exact locality, some others unavoidably take the form $D\boldsymbol{a}D\nabla\varphi^{1}$, hence are only approximately local. This is due to the nonlinear dependence of the solution operator on the coefficient field, and is the reason why some logarithms are lost in the convergence rate to normality in Theorem 1(iv), cf. Section 9.

Proposition 5.2** (Formulas for second Malliavin derivatives).**

Let $\varepsilon=1$ and drop it from all subscripts in the notation. For all $1\leq n\leq\ell$, there holds for any smooth deterministic function $\bar{w}$,

[TABLE]

5.1. Proof of Proposition 5.1

We split the proof into three steps.

Step 1. Proof that for any random function $w$ with $\nabla\cdot\boldsymbol{a}\nabla w$ deterministic there holds for all $0\leq n\leq\ell$,

[TABLE]

We argue by induction. The result is obvious for $n=0$ (recall $\Xi^{0}[\nabla w]=\boldsymbol{a}\nabla w$, cf. Definition 3.1). Assume that it holds for some $n\geq 0$. In order to deduce it at level $n+1$, we appeal to the alternative definition (3.6) of the commutator, so that it suffices to prove

[TABLE]

The definition of $\sigma^{*,n+1}$ (cf. Definition 2.1) yields

[TABLE]

Then using the skew-symmetry of $\sigma^{*,n+1}_{ji_{1}\ldots i_{n}}$ and the equation $D(\nabla\cdot\boldsymbol{a}\nabla w)=0$ in the form

[TABLE]

the claim (5.5) follows.

Step 2. Proof that for any random function $w$ with $\nabla\cdot\boldsymbol{a}\nabla w$ deterministic and for any smooth deterministic function $\bar{w}$ there holds for all $0\leq n,m\leq\ell$,

[TABLE]

where we recall that by a slight abuse of notation we similarly define $\Xi^{n}[H]$ as in Definition 3.1 even if $H$ is not a gradient field.

Let $n\geq 0$ be fixed. We argue by induction on $m$. For $m=0$, the result (5.6) coincides with (5.4) (recall $E^{0}[\nabla\bar{w}]=0$, cf. Definition 2.6). Now assume that (5.6) holds for some $m\geq 0$ and let us deduce it at the level $m+1$. Given a polynomial $\bar{q}$ of order $m+1$, noting that $\nabla^{m+2}\bar{q}$ vanishes, identity (2.9) ensures that $\nabla\cdot\boldsymbol{a}\nabla F^{m+1}[\bar{q}]$ is deterministic. We then apply (5.6) with $w=F^{m+1}[\bar{q}]$ and $\bar{w}=\bar{q}$. Abundantly using that $\nabla^{m+1}\bar{q}$ is constant and noting that identity (2.5) together with the definition of $E^{m}$ yields

[TABLE]

this leads to

[TABLE]

Note that both sides of this identity depend linearly on the (constant and deterministic) symmetric tensor $\nabla^{m+1}\bar{q}$. Hence, we may replace it by the (deterministic) symmetric tensor $\nabla^{m+1}\bar{w}$, to the effect of

[TABLE]

By definition of $E^{m}$ and $\Xi^{n}$ (cf. Definitions 2.6 and 3.1) together with Lemma 2.4, using the general Leibniz rule, we find

[TABLE]

Injecting identities (5.6) and (5.7) into this equality, using again the definition of $E^{m}$ (cf. Definition 2.6) and the general Leibniz rule in the forms

[TABLE]

and

[TABLE]

and recombining the terms, the result (5.6) follows at level $m+1$. For instance, the increment of the first triple sum over $(k,s,l)$ splits into three contributions,

[TABLE]

Step 3. Conclusion.

Applying (5.4) yields (5.1). We now turn to (5.3) and first note that by definition of $\Xi^{\circ,n}$ (cf. Definition 3.4) and by (2.5) we have

[TABLE]

so that by definition of $E^{n}$ (cf. Definition 2.6) and by (3.6) we obtain

[TABLE]

We now apply $D$ to this identity, add (5.6) with $m=n$, with $(w,\bar{w})$ replaced by $(u,\bar{u}^{n})$, and with $s$ replaced by $s+1$, and so obtain (5.3). Finally, (5.2) follows by subtracting (5.3) (with $(u,\bar{u}^{n})=(w,\bar{w})$) from (5.1) and inserting the formula for $DE^{n}[\nabla\bar{u}^{n}]$ already used in (5.1) for $m=n$. ∎

5.2. Proof of Proposition 5.2

The conclusion follows by repeating the proof of Proposition 5.1, taking advantage of the same algebraic identities, but now further keeping track of the mixed term in $D^{2}(\boldsymbol{a}\nabla w)=\boldsymbol{a}\nabla D^{2}w+2D\boldsymbol{a}\nabla Dw+D^{2}\boldsymbol{a}\,\nabla w$. ∎

6. Annealed Calderón-Zygmund theory

In this section, we establish the following new annealed Calderón-Zygmund estimate for linear elliptic equations in divergence form with random coefficients. This constitutes a useful upgrade of the quenched large-scale Calderón-Zygmund estimates of [4, 3, 21]. The merit of (6.2)–(6.4) below is that the stochastic $\operatorname{L}^{q}$ norm is inside the spatial $\operatorname{L}^{p}$ norm; while there is a tiny loss in stochastic integrability, the spatial integrability is as in the constant-coefficient case.

Theorem 6.1**.**

Consider the Helmholtz projection $T:=\nabla(\nabla\cdot\boldsymbol{a}\nabla)^{-1}\nabla\cdot\boldsymbol{a}$. There exists a $\frac{1}{8}$-Lipschitz stationary field $r_{*}\geq 1$ on $\mathbb{R}^{d}$ such that for all $h\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$ and $1<q\leq p<\infty$ there holds

[TABLE]

where $B_{*}(x):=B_{r_{*}(x)}(x)$. In addition, in the Gaussian setting with integrable correlations (1.5), the stationary random field $r_{*}$ satisfies $\mathbb{E}\left[\exp(\frac{1}{C}r_{*}^{d})\right]\leq 2$ for some $C\simeq 1$. In particular, for all $1<p,q<\infty$ and $0<\delta\leq\frac{1}{2}$,

[TABLE]

We shall also make use of the following consequence of this annealed estimate for the corresponding Riesz potentials. As opposed to the above, this can alternatively be deduced from the annealed Green’s function estimates in [37].

Corollary 6.2**.**

Consider the Riesz potentials $U_{e,\circ}:=\nabla\triangle^{-1}e\cdot$ and $U_{e}:=\nabla(\nabla\cdot\boldsymbol{a}\nabla)^{-1}e\cdot\boldsymbol{a}$ with $|e|=1$. For all $h\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$ and $1<p,q<\infty$ with $\frac{dp}{d+p}>1$,

[TABLE]

and for all $0<\delta\leq\frac{1}{2}$,

[TABLE]

6.1. Proof of Theorem 6.1

Our argument is based on a multiple use of the following version of the Calderón-Zygmund lemma due to Shen [48, Theorem 2.1] (cf. also [49, Theorem 2.4]), based on ideas by Caffarelli and Peral [13]. For a ball $B$, we set $B=B(x_{B},r_{B})$ and abusively write $\alpha B:=B(x_{B},\alpha r_{B})$ for $\alpha>0$.

Lemma 6.3** ([48, 13]).**

Let $1\leq p_{0}<p_{1}\leq\infty$, $C_{0}>0$, and $f,g\in\operatorname{L}^{p_{0}}\cap\operatorname{L}^{p_{1}}(\mathbb{R}^{d})$. Assume that for all balls $B\subset\mathbb{R}^{d}$ there exist measurable functions $f_{B,0}$ and $f_{B,1}$ such that $f=f_{B,0}+f_{B,1}$ on $B$ and

[TABLE]

Then, for all $p_{0}<q<p_{1}$,

[TABLE]

Combining this lemma with the quenched large-scale Lipschitz regularity theory due to [7, 6, 21] leads to the following quenched large-scale Calderón-Zygmund estimate. This slightly shortens the proof of [4] and [21, Corollary 4] (cf. also [3, Section 7]).

Proposition 6.4** (Quenched large-scale Calderón-Zygmund estimate).**

There exists a stationary random field $r_{*}$ as in the statement of Theorem 6.1 such that for all $h\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$ and $1<p<\infty$,

[TABLE]

where we recall the notation $B_{*}(x):=B_{r_{*}(x)}(x)$. ∎

Before proceeding to the proof, we state the following useful properties of averages on the family of balls $\{B_{*}(x)\}_{x\in\mathbb{R}^{d}}$. The point is that these balls have a radius that varies with their center $x$.

Lemma 6.5**.**

For any measurable function $f$ on $\mathbb{R}^{d}$, the following estimates hold.

1. (i)

For all balls $B\subset\mathbb{R}^{d}$ with $r_{B}\geq\frac{1}{4}r_{*}(x_{B})$, we have

[TABLE]

hence, in particular,

[TABLE] 2. (ii)

For all balls $B\subset\mathbb{R}^{d}$ and $x_{0}\in B$ with $r_{B}\leq\frac{1}{4}r_{*}(x_{0})$, we have

[TABLE]

hence, in particular,

[TABLE]

We may now proceed to the proof of Proposition 6.4, which amounts to combining the quenched large-scale Lipschitz regularity theory of [7, 6, 21] together with Lemma 6.3.

Proof of Proposition 6.4.

As in [21, Step 1 of the proof of Corollary 4], the random field $r_{*}$ is chosen as the largest $\frac{1}{8}$-Lipschitz lower bound on the so-called minimal radius defined in [21] (what is denoted here by $r_{*}$ thus corresponds to the notation $\underline{r_{*}}$ in [21]). We split the proof into two steps.

Step 1. For all balls $B\subset\mathbb{R}^{d}$, decomposing $Th=\nabla w_{B,0}+\nabla w_{B,1}$ with

[TABLE]

we prove for all $2\leq p\leq\infty$,

[TABLE]

We start with the proof of (6.10), and we first consider the case $r_{B}\geq\frac{1}{4}r_{*}(x_{B})$. An energy estimate on the equation for $w_{B,0}$ yields

[TABLE]

hence, combined with (6.5) and (6.6),

[TABLE]

which proves (6.10). Next, we consider the case $r_{B}\leq\frac{1}{4}r_{*}(x_{B})$. For $x_{0}\in B$, using the $\frac{1}{8}$-Lipschitz property of $r_{*}$ in form of $|B_{*}(x_{0})|\simeq|B_{*}(x_{B})|$, using the energy estimate (6.12), and applying (6.8) (with $x_{0}$ and $B$ replaced by $x_{B}$ and $\frac{1}{10}B$), we find

[TABLE]

(The choice of $\frac{1}{2}B$ in the right-hand side does not matter here and is only made for later reference.) Integrating this estimate over $x_{0}\in B$, the conclusion (6.10) follows.

We turn to the proof of (6.11), which is a consequence of quenched large-scale Lipschitz regularity theory, noting that $w_{B,1}$ is indeed $\boldsymbol{a}$-harmonic in $B$. It is obviously sufficient to prove for all $x_{0}\in\frac{1}{12}B$,

[TABLE]

If $r_{*}(x_{0})\geq\frac{1}{3}r_{B}$, this follows from (6.8) (with $B$ replaced by $\frac{1}{12}B$). Let now $x_{0}\in\frac{1}{12}B$ be fixed with $r_{*}(x_{0})\leq\frac{1}{3}r_{B}$. Since the ball $B(x_{0},r_{*}(x_{0})+\frac{1}{12}r_{B})$ is contained in $\frac{1}{2}B$, where $w_{B,1}$ is $\boldsymbol{a}$-harmonic, the quenched large-scale Lipschitz regularity theory (e.g. [21, Theorem 1]) implies

[TABLE]

Since for $r_{*}(x_{0})\leq\frac{1}{3}r_{B}$ the $\frac{1}{8}$-Lipschitz property of $r_{*}$ yields $r_{*}(x_{B})\leq r_{*}(x_{0})+\frac{1}{96}r_{B}\leq 2r_{B}$, the conclusion (6.14) follows from (6.5) (with $B$ replaced by $\frac{1}{2}B$).

Step 2. Conclusion.

Combining (6.10) and (6.11), and applying Lemma 6.3 with

[TABLE]

we conclude for all $2\leq p<\infty$,

[TABLE]

A standard duality argument allows to deduce the corresponding result for $1<p\leq 2$ (e.g. [21, Step 7 of the proof of Corollary 4]). ∎

We need two additional ingredients for the proof of Theorem 6.1. The first is a deterministic regularity result for $\boldsymbol{a}$-harmonic functions. It follows e.g. from [10, proof of Lemma 4], but a short proof is included below for the reader’s convenience.

Lemma 6.6**.**

If $w$ is $\boldsymbol{a}$-harmonic on a ball $B$, there holds

[TABLE]

The second ingredient is needed to pass from statements on averages at the scale $r_{*}$ to corresponding statements on unit scale. This can be done at no loss in the spatial, but a small controlled loss in the stochastic integrability.

Lemma 6.7**.**

Let the stationary random field $r_{*}$ satisfy $\mathbb{E}\left[\exp(\frac{1}{C}r_{*}^{d})\right]\leq 2$ for some $C\simeq 1$. For all $f\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))$ and $1\leq q\leq p<\infty$, there hold

1. (i)

for all $r>q$,

[TABLE] 2. (ii)

for all $r<q$,

[TABLE]

where we have set $L(t):=\log(2+\frac{1}{t})$. ∎

With these ingredients at hand, we turn the proof of Theorem 6.1. The argument is based on a second application of Lemma 6.3 starting with the quenched large-scale Calderón-Zygmund estimate of Proposition 6.4.

Proof of Theorem 6.1.

Let $r_{*}$ be chosen as in the proof of Proposition 6.4 above. We split the proof into two steps.

Step 1. For all balls $B\subset\mathbb{R}^{d}$, decomposing $Th=\nabla w_{B,0}+\nabla w_{B,1}$ as in (6.9), we prove for all $1<q<\infty$ and $1\leq p_{1}\leq\infty$,

[TABLE]

We start with the proof of (6.15). It suffices to show

[TABLE]

since the conclusion (6.15) then follows by taking the expectation. First consider the case $r_{B}\leq\frac{1}{4}r_{*}(x_{B})$. Integrating (6.13) over $x_{0}\in B$ yields

[TABLE]

The claim (6.17) follows from this estimate together with Jensen’s inequality if $q\geq 2$, but a different argument is needed if $q<2$. For $x_{0}\in\frac{1}{2}B$, since the condition $r_{B}\leq\frac{1}{4}r_{*}(x_{B})$ and the $\frac{1}{8}$-Lipschitz property of $r_{*}$ imply $r_{B}\leq\frac{1}{4}r_{*}(x_{0})+\frac{1}{64}r_{B}$, hence $r_{B}\leq\frac{1}{2}r_{*}(x_{0})$, it follows from (6.7) (with $B$ replaced by $\frac{1}{2}B$) that

[TABLE]

Combined with (6.18), this leads to (6.17) by Jensen’s inequality. It remains to consider the case $r_{B}\geq\frac{1}{4}r_{*}(x_{B})$. Applying Proposition 6.4 (with $q$ playing the role of $p$) to the equation for $w_{B,0}$ yields

[TABLE]

For $x\in B$, the balls $B_{*}(x)$ and $B$ are disjoint whenever $|x-x_{B}|\geq r_{*}(x)+r_{B}$. Since the condition $r_{B}\geq\frac{1}{4}r_{*}(x_{B})$ and the $\frac{1}{8}$-Lipschitz property of $r_{*}$ imply

[TABLE]

we deduce that $B_{*}(x)$ and $B$ are disjoint whenever $|x-x_{B}|\geq 5r_{B}+\frac{1}{8}|x-x_{B}|$, that is, whenever $|x-x_{B}|\geq 6r_{B}$. This implies

[TABLE]

so that the claim (6.17) follows from (6.19).

We turn to the proof of (6.16). For $2\leq q\leq\infty$, taking the $\frac{q}{2}$th power of (6.14), using Jensen’s inequality, and taking the expectation, we easily obtain (6.16). In order to conclude for all values $1\leq q\leq\infty$, it suffices to establish the following improved version of (6.14): for all $x_{0}\in\frac{1}{24}B$,

[TABLE]

If $r_{*}(x_{0})\geq\frac{1}{6}r_{B}$, this already follows from (6.7) (with $B$ replaced by $\frac{1}{24}B$). Let now $x_{0}\in\frac{1}{24}B$ be fixed with $r_{*}(x_{0})\leq\frac{1}{6}r_{B}$. Since $w_{B,1}$ is $\boldsymbol{a}$-harmonic on $B$ and since the ball $B(x_{0},r_{*}(x_{0})+\frac{1}{24}r_{B})$ is contained in $\frac{1}{4}B$, the quenched large-scale Lipschitz regularity theory (e.g. [21, Theorem 1]) implies

[TABLE]

and Lemma 6.6 (with $B$ replaced by $\frac{1}{2}B$) then leads to

[TABLE]

Since for $r_{*}(x_{0})\leq\frac{1}{6}r_{B}$ the $\frac{1}{8}$-Lipschitz property of $r_{*}$ yields $r_{*}(x_{B})\leq r_{*}(x_{0})+\frac{1}{192}r_{B}\leq 2r_{B}$, the conclusion (6.20) follows from (6.5) (with $B$ replaced by $\frac{1}{2}B$) and from Jensen’s inequality.

Step 2. Conclusion.

Combining (6.15) and (6.16), and applying Lemma 6.3 with

[TABLE]

we deduce for all $1<q\leq p<\infty$,

[TABLE]

which is the main stated result (6.1). It remains to post-process this result by rewriting both sides of the inequality: combining (6.1) with the different bounds stated in Lemma 6.7, and noting that for $q\leq p$ there holds $(\frac{1}{q}-\frac{1}{2})_{+}+(\frac{1}{2}-\frac{1}{p})_{+}=\frac{1}{q\wedge 2}-\frac{1}{p\vee 2}$, the conclusion (6.2) follows for $q\leq p$, while the corresponding result for $p\leq q$ is deduced by duality. ∎

6.2. Proof of Corollary 6.2

For $d>1$, appealing to the integral representation

[TABLE]

the triangle inequality yields

[TABLE]

and the result (6.3) follows from the Hardy-Littlewood-Sobolev inequality. Next, writing $U_{e}=T\boldsymbol{a}^{-1}U_{e,\circ}\boldsymbol{a}$, the estimate (6.4) follows from the combination of (6.2) and (6.3). ∎

6.3. Proof of Lemma 6.5

We start with the proof of (i), and thus let $r_{B}\geq\frac{1}{4}r_{*}(x_{B})$. For $|x-z|\leq\frac{1}{2}r_{*}(z)$, the $\frac{1}{8}$-Lipschitz property of $r_{*}$ implies $|x-z|\leq\frac{1}{2}r_{*}(x)+\frac{1}{16}|x-z|$, hence $|x-z|\leq r_{*}(x)$, and similarly $|B_{*}(x)|\simeq|B_{*}(z)|$. These observations lead to the lower bound

[TABLE]

Since $r_{*}(x_{B})\leq 4r_{B}$, we note that for all $z\in B$ the ball $\frac{1}{5}B_{*}(z)$ is included in

[TABLE]

and (6.5) follows. The upper bound in (6.6) is a consequence (6.5) with $r_{B}\uparrow\infty$, so that it remains to establish the lower bound. For that purpose, we write again

[TABLE]

Since for $|x-z|\leq r_{*}(x)$ the $\frac{1}{8}$-Lipschitz property of $r_{*}$ implies $r_{*}(x)\simeq r_{*}(z)$, the last integral in bracket is $\lesssim 1$, and the desired lower bound follows.

We turn to the proof of (ii). By Jensen’s inequality, it suffices to establish (6.7). Let $x_{0}\in B$ with $r_{B}\leq\frac{1}{4}r_{*}(x_{0})$. Noting that for all $x\in B$ the ball $B_{*}(x_{0})$ is contained in $B(x,r_{*}(x_{0})+2r_{B})$, we find

[TABLE]

The ball $B_{4r_{B}}$ at the origin can be covered by $N\simeq 1$ balls of radius $r_{B}$. Denote by $(z_{i})_{i=1}^{N}\subset B_{4r_{B}}$ their centers. The collection $(B(z_{i},r_{*}(x_{0})-r_{B}))_{i=1}^{N}$ then covers $B_{r_{*}(x_{0})+2r_{B}}$. By subadditivity of the square root, we may decompose

[TABLE]

hence, noting that $z_{i}+B\subset 5B$ for all $i$,

[TABLE]

For $x\in 5B$, the $\frac{1}{8}$-Lipschitz property of $r_{*}$ yields

[TABLE]

and (iii) follows.∎

6.4. Proof of Lemma 6.6

By scaling and translation invariance, we may assume that $B$ is the unit ball at the origin and that $\fint_{B}w=0$. Choose a smooth cut-off function $\chi$ with $\chi=1$ in $\frac{1}{2}B$ and $\chi=0$ outside $B$. Caccioppoli’s inequality yields

[TABLE]

Sobolev’s inequality combined with Poincaré’s inequality (with vanishing boundary condition) then implies for some suitable $p=p(d)>2$,

[TABLE]

Choosing $q:=2\frac{p-1}{p-2}$ and writing $\chi=\xi^{q}$, we deduce by Hölder’s inequality,

[TABLE]

hence, by Young’s inequality,

[TABLE]

Injecting this into (6.22) and using Poincaré’s inequality (with vanishing mean-value), we conclude

[TABLE]

as desired.∎

6.5. Proof of Lemma 6.7

Step 1. Proof of the upper bound (i).

Let $1\leq q\leq p<\infty$. We start with the following deterministic version of the claimed upper bound: for all $R\geq 1$ we have

[TABLE]

If $q\geq 2$, Jensen’s inequality on $t\mapsto t^{\frac{q}{2}}$ yields

[TABLE]

If $q\leq 2$, covering the ball $B_{R}$ by $N\simeq R^{d}$ balls of radius $\frac{1}{2}$ with centers $(z_{i})_{i=1}^{N}\subset B_{R-\frac{1}{4}}$ such that $\min_{i\neq j}|z_{i}-z_{j}|\geq\frac{1}{2}$, it follows from the discrete $\ell^{2}$–$\ell^{q}$ inequality that

[TABLE]

and hence, since the balls $(B_{\frac{1}{4}}(x+z_{i}))_{i=1}^{N}$ are disjoint and included in $B_{R}(x)$,

[TABLE]

These two estimates (6.25) and (6.26) lead to

[TABLE]

and the claim (6.24) follows from Jensen’s inequality on $t\mapsto t^{\frac{p}{q}}$.

We turn to the proof of (i). Conditioning with respect to the value of $r_{*}(x)$ on dyadic scale and using Jensen’s inequality on $t\mapsto t^{\frac{p}{q}}$ with $\sum_{n=0}^{\infty}2^{-nd\delta}\simeq\delta^{-1}$, we find for $0<\delta\leq 1$,

[TABLE]

Taking advantage of the $\frac{1}{8}$-Lipschitz property of $r_{*}$, we can rewrite

[TABLE]

and applying (6.24) with $R=2^{n+1}$ then leads to

[TABLE]

Hence, for $r>q$, using Hölder’s inequality and the moment bounds on the stationary $r_{*}$,

[TABLE]

where in the last estimate we used $\sum_{n}(2^{nd})^{\alpha}e^{-2^{nd}\theta}\lesssim_{\alpha,d}\theta^{-\alpha}$. Optimizing $\frac{1}{\delta}\theta^{-\delta}$ in $0<\delta\leq 1$ through $\delta=(\log\frac{1}{\theta})^{-1}$ yielding $\sim\log\frac{1}{\theta}$, the desired upper bound follows.

Step 2. Proof of the lower bound (ii).

Let $1\leq r\leq p\leq\infty$. We start with the following deterministic version of the claimed lower bound: for all $R\geq 1$ we have

[TABLE]

If $2\leq r\leq p$, using the discrete $\ell^{r}$–$\ell^{2}$ inequality and the discrete $\ell^{p}$–$\ell^{r}$ inequality, we indeed find

[TABLE]

If $r\leq 2\leq p$, the triangle inequality in $\operatorname{L}^{\frac{2}{r}}(B_{R}(x))$ and the discrete $\ell^{p}$–$\ell^{2}$ inequality lead to

[TABLE]

Finally, if $r\leq p\leq 2$, using again the triangle inequality in $\operatorname{L}^{\frac{2}{r}}(B_{R}(x))$ and Jensen’s inequality on $t\mapsto t^{\frac{2}{p}}$, we obtain

[TABLE]

This proves the claim (6.27) in all cases.

We turn to the proof of (ii). Conditioning with respect to the value of $r_{*}(x)$ on dyadic scale and using Jensen’s inequality on $t\mapsto t^{\frac{p}{r}}$ with $\sum_{n=0}^{\infty}2^{-nd\delta}\simeq\delta^{-1}$, we find for $0<\delta\leq 1$,

[TABLE]

Applying (6.27) with $R=(2^{n}-2)\vee 1$ and taking advantage of the $\frac{1}{8}$-Lipschitz property of $r_{*}$, we deduce

[TABLE]

Hence, for $q>r$, using Hölder’s inequality and the moment bounds on $r_{*}$, and proceeding as for (i),

[TABLE]

Optimizing in $0<\delta\leq 1$ yields the desired lower bound. ∎

7. Proof of the higher-order pathwise result

In order to establish the accuracy of the two-scale expansion of higher-order commutators (hence, the higher-order pathwise structure), we appeal to the representation formula (5.3) established in Proposition 5.1 and the conclusion easily follows from the following estimates on the Malliavin derivative of higher-order correctors.

Lemma 7.1**.**

Let $T$, $U_{e}$, and $U_{e,\circ}$ denote the operators defined in Theorem 6.1 and Corollary 6.2 with $\boldsymbol{a}$ replaced by its adjoint $\boldsymbol{a}^{*}$. For all $g\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$ and $0\leq s\leq\ell-1$,

[TABLE]

where $U_{\cdot}^{k}$ denotes any product of the form $U_{1}\ldots U_{k}$ with $U_{j}\in\{U_{e_{i}},U_{e_{i},\circ}\}_{1\leq i\leq d}$ for all $j$, and where we implicitly sum over all possible such products. Furthermore, for all $\delta>0$ and $2\leq p<\infty$,

[TABLE]

Proof.

Denote by $J_{s,1}(g)^{2}$, $J_{s,2}(g)^{2}$, and $J_{s,3}(g)^{2}$ the three left-hand side terms in (7.1), hence $J_{s}(g)^{2}=J_{s,1}(g)^{2}+J_{s,2}(g)^{2}+J_{s,3}(g)^{2}$. We split the proof into two steps.

Step 1. Proof that for all $g\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$ and $0\leq s\leq\ell-1$,

[TABLE]

We start with (7.3). Taking the Malliavin derivative of the equation for $\varphi^{s+1}$ (cf. Definition 2.1) and testing it with $(\nabla\cdot\boldsymbol{a}^{*}\nabla)^{-1}\nabla\cdot(\boldsymbol{a}^{*}g)$, we obtain by definition of $Tg=\nabla(\nabla\cdot\boldsymbol{a}^{*}\nabla)^{-1}\nabla\cdot(\boldsymbol{a}^{*}g)$,

[TABLE]

where for the first right-hand side term we used that by definition (1.4) the Malliavin derivative of $\boldsymbol{a}$ takes the form

[TABLE]

The claim (7.3) follows by applying $(\int_{\mathbb{R}^{d}}(\int_{B(x)}\cdot\,dz)^{2}dx)^{\frac{1}{2}}$ to both sides of the estimate. Next, by definition of $U_{e_{i_{s+1}}}g=\nabla(\nabla\cdot\boldsymbol{a}^{*}\nabla)^{-1}(e_{i_{s+1}}\cdot\boldsymbol{a}^{*}g)$, we find

[TABLE]

so that (7.4) follows directly. We turn to (7.5). Writing $g=\triangle\triangle^{-1}g$, integrating by parts, and using the Malliavin derivative of the equation for $\triangle\sigma^{s}$ (cf. Definition 2.1), we find

[TABLE]

and the claim (7.5) follows from taking the Malliavin derivative of the definition of $q^{s}$ and writing $\nabla\triangle^{-1}g_{i}=U_{e_{i},\circ}g$.

Step 2. Conclusion.

Using the notation introduced in the statement, Step 1 yields for all $g\in C_{c}^{\infty}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$,

[TABLE]

Since by definition $TU=U$ and $TU_{\circ}=U_{\circ}$, the decomposition (7.1) follows by induction. Next, combining (7.1) with the corrector estimates of Proposition 2.2, we find for all $g\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))^{d}$, $0\leq s\leq\ell-1$, $2\leq p<\infty$, and $\delta>0$,

[TABLE]

where the anchoring in $s=0$ follows from (7.3) because of $J_{0,2}=J_{0,3}=0$ in view of $\varphi^{0}=1$ and $\sigma^{0}=0$. Here we have used that for a stationary $\psi$ we have by the triangle inequality in $\operatorname{L}^{\frac{p}{2}}(\Omega)$ and Hölder’s inequality,

[TABLE]

Now, for $0\leq s\leq\ell-1$, $2\leq p<\infty$, $0\leq k\leq s$, $\beta\in\{0,1\}$, and $\delta>0$, we apply the bound of Corollary 6.2 $k$ times and the one of Theorem 6.1 once, increasing at every application the stochastic integrability by $\frac{\delta}{k+1}$, to the effect of

[TABLE]

given that the restriction $k\leq s\leq\ell-1$ indeed ensures $\frac{2d}{d+2k}>1$, and the conclusion follows. ∎

With these estimates at hand, we may now turn to the proof of Theorem 1(ii), that is, of the higher-order pathwise result.

Proof of Theorem 1(ii).

Let $1\leq p<\infty$. The starting point is the $\operatorname{L}^{p}$ inequality of Proposition 4.1(iii) together with (4.2), in the form

[TABLE]

Inserting identity (5.3) of Proposition 5.1, integrating by parts in the $z$-integral to put derivatives on $g$ and $\bar{u}^{n}$, using (7.6), and using the corrector estimates of Proposition 2.2, we deduce for all $1\leq n\leq\ell$,

[TABLE]

where we have set

[TABLE]

We separately treat the different right-hand side terms in (7.7). First, from equation (2.12), formula (2.5), the annealed Calderón-Zygmund estimate of Theorem 6.1, and the corrector estimates in Proposition 2.2, we obtain

[TABLE]

Second, recalling the notation $T=\nabla(\nabla\cdot\boldsymbol{a}^{*}\nabla)^{-1}\nabla\cdot\boldsymbol{a}^{*}$ from Theorem 6.1 and using the equation for $u$ in the form $-\nabla\cdot\boldsymbol{a}\nabla D_{z}u=\nabla\cdot D_{z}\boldsymbol{a}\nabla u$, we find

[TABLE]

hence, using Theorem 6.1 to bound the operator $T$ and to estimate $\nabla u$, and applying the corrector estimates of Proposition 2.2, we deduce

[TABLE]

Third, decomposing

[TABLE]

applying (7.2), and using the corrector estimates of Proposition 2.2,

[TABLE]

Fourth, applying (7.2) and using the corrector estimates of Proposition 2.2,

[TABLE]

Collecting the above estimates on the right-hand side terms of (7.7), we see that all of them are bilinear in $g$ and $f$. Crucial for us is their scaling in the underlying lengthscale, which has three origins: the (total) number of derivatives (to be counted negatively), the spatial integrability exponent $p$ (to be counted in form of $\frac{d}{p}$), and the weights $\mu_{d,n}$ (cf. Proposition 2.2). For instance, in view of (7.8), the first right-hand side term in (7.7) has scaling $\leq L^{-n+\frac{d}{4}+\frac{d}{4}}\mu_{d,n}(L)=L^{-n+\frac{d}{2}}\mu_{d,n}(L)$, where $L$ stands for the lengthscale, because $\nabla\tilde{u}^{l}$ has the same scaling as $\nabla^{l-1}f$ (cf. Definition 2.7). Let us take the last term as a second example: the scaling of its $(s,l)$-contribution is $L^{-(n+s)+\frac{d+2s}{2}}\mu_{d,n}(L)=L^{-n+\frac{d}{2}}\mu_{d,n}(L)$. In fact, all terms have a scaling at most $L^{-n+\frac{d}{2}}\mu_{d,n}(L)$, which after $\varepsilon$-rescaling ($L=\varepsilon^{-1}$) turns into the desired estimate. ∎

8. Proof of the convergence of the covariance structure

Combining the representation formula (5.2) of Proposition 5.1 with the Helffer-Sjöstrand identity of Proposition 4.1(ii), we establish Theorem 1(iii), that is, the convergence of the covariance structure of the higher-order standard homogenization commutator.

Proof of Theorem 1(iii).

Combining the Helffer-Sjöstrand identity of Proposition 4.1(ii) with the representation formula (5.2) tested with $g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, we obtain, using (7.6) and recalling that $(1+\mathcal{L})^{-1}$ has operator norm $\leq 1$,

[TABLE]

where we have set

[TABLE]

and

[TABLE]

Appealing to the definition of $E^{n}[\nabla\bar{w}]$ (cf. Definition 2.6), to (4.2), and to the corrector estimates of Proposition 2.2, we find

[TABLE]

and similarly,

[TABLE]

in terms of

[TABLE]

Then applying (7.2) and again using the corrector estimates of Proposition 2.2, we deduce for all $\delta>0$,

[TABLE]

As in the proof of Theorem 1(ii) (cf. Section 7), it is easily checked that these expressions in $(g,\nabla\bar{w})$ have the desired length scaling, hence lead to the claimed error estimate after $\varepsilon$-rescaling.

Now that the right-hand side of (8.1) is properly estimated, we turn to a suitable decomposition of $\mathbb{E}\left[\langle S,(1+\mathcal{L})^{-1}S\rangle_{\mathfrak{H}}\right]$. Inserting the definition of $E^{n}[\nabla\bar{w}]$ (cf. Definition 2.6) into $S$, using the definition (4.1) of the scalar product in $\mathfrak{H}$, expanding the product, and focusing on terms of order $\varepsilon^{k}$ with $k\leq n-1$, we write

[TABLE]

where we have set

[TABLE]

Combining this with (8.1) and with the above estimates on $\mathbb{E}\left[\|S\|_{\mathfrak{H}}^{2}\right]^{\frac{1}{2}}$ and $\mathbb{E}\left[\|R\|_{\mathfrak{H}}^{2}\right]^{\frac{1}{2}}$ and appealing to the stationarity of $T^{l,k}_{ji_{1}\ldots i_{l}j_{1}\ldots j_{k+1}}$ (recalling that $\mathcal{L}$ commutes with shifts), we deduce after $\varepsilon$-rescaling,

[TABLE]

Appealing to the boundedness of $(1+\mathcal{L})^{-1}$ and to the corrector estimates of Proposition 2.2 in the form

[TABLE]

the conclusion follows after Taylor-expanding $\nabla^{l^{\prime}}g_{j^{\prime}}(z+\varepsilon y)\nabla^{k^{\prime}+1}\bar{w}(z+\varepsilon y)$ to order $\varepsilon^{n-m}$ with the integrability assumption $\int_{\mathbb{R}^{d}}|y|^{n}|c(y)|dy\leq 1$. ∎

9. Proof of the normal approximation result

In order to establish the asymptotic normality of the higher-order standard homogenization commutators, we combine the representation formula of Proposition 5.2 with the second-order Poincaré inequality of Proposition 4.1(iv). We then need the following version of Lemma 7.1 for the estimation of second Malliavin derivatives of higher-order correctors.

Lemma 9.1**.**

For all $\zeta\in C^{\infty}_{c}(\mathbb{R}^{d})$, $g\in C^{\infty}_{c}(\mathbb{R}^{d})^{d}$, $0\leq s\leq\ell-1$, and $p\geq 1$,

[TABLE]

where the notation is taken from Lemma 7.1. ∎

Proof.

Denote by $K_{s,1}(g,\zeta)$, $K_{s,2}(g,\zeta)$, and $K_{s,3}(g,\zeta)$ the three left-hand side terms in the claimed estimate, hence $K_{s}(g,\zeta)=K_{s,1}(g,\zeta)+K_{s,2}(g,\zeta)+K_{s,3}(g,\zeta)$. We prove that for all $g\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\mathbb{R}^{d}))^{d}$, $\zeta\in C^{\infty}_{c}(\mathbb{R}^{d})$, $0\leq s\leq\ell-1$, and $p\geq 1$,

[TABLE]

These estimates combine to

[TABLE]

Noting as in the proof of Lemma 7.1 that $TU=U$ and $TU_{\circ}=U_{\circ}$, the conclusion follows by induction.

Let $0\leq s\leq\ell-1$. We start with the proof of (9.1). Taking the second Malliavin derivative of the equation for $\varphi^{s+1}$ (cf. Definition 2.1), we find

[TABLE]

Testing with $(\nabla\cdot\boldsymbol{a}^{*}\nabla)^{-1}\nabla\cdot(\boldsymbol{a}^{*}g)$, recalling the definition of $T$, and using (7.6) and

[TABLE]

we find after integration against $\zeta\otimes\zeta$,

[TABLE]

and the claim (9.1) follows from Hölder’s inequality. Next, (9.2) follows from the definition of $U$. We turn to (9.3). Writing $g=\triangle\triangle^{-1}g$, integrating by parts, and using the equation for $\triangle\sigma^{s}$ (cf. Definition 2.1), we find

[TABLE]

Hence, by the definition of $q^{s}$ to which we apply $D^{2}_{xy}$, we obtain an analogous structure as in (9.5), and the claim (9.3) follows. ∎

With these estimates at hand, we now turn to the proof of Theorem 1(iv), that is, of the normal approximation result for the higher-order standard homogenization commutators.

Proof of Theorem 1(iv).

We split the proof into two steps.

Step 1. Proof that for all $8\leq p,q<\infty$ and $1\leq n\leq\ell$,

[TABLE]

In order to not overestimate the logarithmic correction in the final result, we have to be specific here on the $p$-dependence; however, this is only necessary in the leading-order terms (that is, in the first two terms). Let $1\leq n\leq\ell$. Applying the representation formula of Proposition 5.2, using (7.6) and (9.4), inserting the definition of $E_{\varepsilon}^{n}[\nabla\bar{w}]$ (cf. Definition 2.6), and ordering terms by the number of Malliavin derivatives on $\boldsymbol{a}$, we find

[TABLE]

in terms of

[TABLE]

where we identify the operators $U_{i}$ with their kernels and we have set

[TABLE]

with the implicit understanding that the contribution of $\sum_{l=0}^{k+s-n}$ vanishes if $k<n-s$. We analyze the three right-hand side terms in (9.6) separately and we start with $U_{1}$. Decomposing the covariance function as $c=c_{0}\ast c_{0}$ so that $\|\zeta\|_{\mathfrak{H}}=\|c_{0}\ast\zeta\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}$, the definition (4.7) of the norm $\|\cdot\|_{\operatorname{op}}$ is equivalent to

[TABLE]

For $p\geq 4$, inserting the special form of $U_{1}(x,y)=\delta(x-y)\,S^{1}(x)$ and noting that the discrete $\ell^{\frac{2p}{p-2}}$–$\ell^{2}$ inequality and the assumption (1.5) imply

[TABLE]

we get

[TABLE]

hence,

[TABLE]

By definition of $S^{1}$, appealing to the corrector estimates of Proposition 2.2, we deduce

[TABLE]

We turn to $U_{2}$. By symmetry, it suffices to estimate the norm of $\tilde{U}^{2}$. From (9.7) and $\|[c_{0}\ast\zeta]_{\infty}\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}\leq\|\zeta\|_{\operatorname{L}^{2}(\mathbb{R}^{d})}$, we see

[TABLE]

Taking the expectation and arguing by duality, we find

[TABLE]

where

[TABLE]

Let $\zeta\in C^{\infty}_{c}(\mathbb{R}^{d};\operatorname{L}^{\infty}(\Omega))$ be fixed with $\|[\zeta]_{\infty}\|_{\operatorname{L}^{4}(\Omega;\operatorname{L}^{2}(\mathbb{R}^{d}))}=1$. Noting that Jensen’s inequality yields

[TABLE]

and that the discrete $\ell^{4}$–$\ell^{2}$ inequality gives

[TABLE]

we deduce by interpolation, for all $2\leq r\leq 4$,

[TABLE]

By definition of $\tilde{U}_{2}$ and by Lemma 7.1 in the form (7.1), we get

[TABLE]

and the corrector estimates of Proposition 2.2 yield for all $p\geq 1$,

[TABLE]

For $p\geq 2$, applying $m$ times the bound of Corollary 6.2 and once the one of Theorem 6.1 with $\delta=\frac{1}{m+1}\frac{2p}{(p-1)(p-2)}\sim\frac{1}{p}$ in order to pass from the stochastic integrability $\frac{2p}{p-1}$ to $\frac{2p}{p-2}$, we obtain the following,

[TABLE]

where the exponent $\frac{m(m+1)}{d}$ of $p\log^{2}p$ comes from the large prefactors in Corollary 6.2 and Theorem 6.1 and is due to the successive deviations of spatial integrability from $2$. Hölder’s inequality together with (9.10) yields for $p\geq 8$,

[TABLE]

By definition of $S^{2;s}$, using the corrector estimates of Proposition 2.2, we conclude

[TABLE]

It remains to analyze $U_{3}$. For $p\geq 2$, we first appeal to the definition of $U_{3}$ and to Lemma 9.1 (with $p$ replaced by $\frac{p}{2}$) in the form

[TABLE]

From (9.7) and (9.8), we deduce

[TABLE]

Taking the $\operatorname{L}^{4}(\Omega)$ norm and using the corrector estimates of Proposition 2.2 in the first right-hand side term, and arguing by duality as in (9.9) and repeating the proof of (9.12) for the second right-hand side term, we obtain for all $8\leq p<\infty$,

[TABLE]

A multiple use of the annealed bounds of Theorem 6.1 and Corollary 6.2 then leads to

[TABLE]

By definition of $S^{3;s}$, using the corrector estimates of Proposition 2.2, we conclude

[TABLE]

Gathering the above estimates and replacing $p$ by $q$ in the estimate for $U^{3}$, the claim follows.

Step 2. Conclusion.

By scaling, the conclusion of Step 1 yields for all $8\leq p,q<\infty$ and $1\leq n\leq\ell$,

[TABLE]

hence, choosing $p=|\!\log\varepsilon|$ and $q=8\vee 4d$, and recalling that $c_{0}=\frac{1}{2}$,

[TABLE]

In terms of $X_{\varepsilon}^{n}:=\varepsilon^{-\frac{d}{2}}\int_{\mathbb{R}^{d}}g\cdot\Xi_{\varepsilon}^{\circ,n}[\nabla\bar{w}]$, we now apply Proposition 4.1(iv) in the form

[TABLE]

Inserting the above estimate on $\|D^{2}X_{\varepsilon}^{n}\|_{{\operatorname{op}}}$ and noting that the proof of Theorem 1(iii) yields $\mathbb{E}\left[\|DX_{\varepsilon}^{n}\|_{\mathfrak{H}}^{4}\right]\,\lesssim_{f,g}\,1$, the conclusion follows. ∎

Appendix A More Malliavin calculus

In this appendix, we include for completeness a short, essentially self-contained proof of Proposition 4.1.

Proof of Proposition 4.1.

We split the proof into four steps.

Step 1. Proof of (i).

In terms of the Ornstein-Uhlenbeck semigroup $e^{-t\mathcal{L}}$ (e.g. [44, Section 1.4]), we may write

[TABLE]

where Mehler’s formula (e.g. [44, (1.67)]) indeed ensures that $\mathbb{E}\left[(e^{-t\mathcal{L}}X)^{2}\right]\to\mathbb{E}\left[X\right]^{2}$ as $t\uparrow\infty$. Computing the derivative in $t$ yields

[TABLE]

hence, appealing to the commutation relation (4.5) in the form $De^{-t\mathcal{L}}=e^{-t}e^{-t\mathcal{L}}D$,

[TABLE]

and the positivity of $\mathcal{L}$ leads to the conclusion,

[TABLE]

Step 2. Proof of (ii).

Let $X,Y\in\mathbb{D}^{1,2}$ with $\mathbb{E}\left[X\right]=\mathbb{E}\left[Y\right]=0$. First note that the Poincaré inequality (i) implies that the restriction of $\mathcal{L}$ to $\operatorname{L}^{2}(\Omega)/\mathbb{C}:=\{U\in\operatorname{L}^{2}(\Omega):\mathbb{E}\left[U\right]=0\}$ is invertible. In particular, there exists $Z\in\operatorname{L}^{2}(\Omega)$ such that $Y=\mathcal{L}Z$. By definition of the adjoint $D^{*}$ (cf. (4.4)), we then find

[TABLE]

Appealing to the commutator relation (4.5) in the form $DY=D\mathcal{L}Z=(1+\mathcal{L})DZ$, we conclude

[TABLE]

Step 3. Proof of (iii).

Similarly as in Step 1, in terms of the Ornstein-Uhlenbeck semigroup $e^{-t\mathcal{L}}$, we may write

[TABLE]

hence, by the commutation relation (4.5) in the form $De^{-t\mathcal{L}}=e^{-t}e^{-t\mathcal{L}}D$,

[TABLE]

Appealing to Mehler’s formula for the Ornstein-Uhlenbeck semigroup $e^{-t\mathcal{L}}$ (e.g. [44, (1.67)]), the Cauchy-Schwarz inequality leads to

[TABLE]

so that the above becomes

[TABLE]

which proves the logarithmic Sobolev inequality. Next, higher integrability is deduced by integration as e.g. in [2, Theorem 3.4]: we write

[TABLE]

so that the logarithmic Sobolev inequality implies

[TABLE]

and the conclusion follows from the Poincaré inequality (i).

Step 4. Proof of (iv).

Let $X\in\operatorname{L}^{2}(\Omega)$ with $\mathbb{E}\left[X\right]=0$ and $\mathrm{Var}\left[X\right]=1$. For $h\in\operatorname{L}^{\infty}(\mathbb{R})$, we define its Stein transform $S_{h}$ as the solution of Stein’s equation

[TABLE]

As in Step 2, there exists $Z\in\operatorname{L}^{2}(\Omega)$ such that $Y=\mathcal{L}Z$. We then compute

[TABLE]

and hence, integrating by parts and using (A.1) in the form $\mathbb{E}\left[\langle DZ,DX\rangle_{\mathfrak{H}}\right]=\mathrm{Var}\left[X\right]=1$,

[TABLE]

Noting that $\|S_{h}^{\prime}\|_{\operatorname{L}^{\infty}}\leq 2\|h\|_{\operatorname{L}^{\infty}}$ (e.g. [42, Theorem 3.3.1]) and taking the supremum over all $h\in\operatorname{L}^{\infty}(\mathbb{R})$, we deduce

[TABLE]

Noting that for $h$ Lipschitz-continuous there holds $\|S_{h}^{\prime}\|_{\operatorname{L}^{\infty}}\leq\sqrt{\frac{2}{\pi}}\|h^{\prime}\|_{\operatorname{L}^{\infty}}$ (e.g. [42, Proposition 3.5.1]), we can deduce a similar bound on the $1$-Wasserstein distance. The corresponding bound on the $2$-Wasserstein distance takes the form

[TABLE]

its proof is of a different nature and is based on an optimal transport argument in density space (cf. [35, Proposition 3.1]).

It remains to estimate the variance $\mathrm{Var}\left[\langle DZ,DX\rangle_{\mathfrak{H}}\right]$. For that purpose, we apply the first-order Poincaré inequality (i) in the form

[TABLE]

Noting that the commutation relation (4.5) leads to $DX=(1+\mathcal{L})DZ$ and $D^{2}X=(2+\mathcal{L})D^{2}Z$, we deduce

[TABLE]

Noting as in [38, Proposition 3.2] that $(1+\mathcal{L})^{-1}$ and $(2+\mathcal{L})^{-1}$ have operator norms on $\operatorname{L}^{4}(\Omega)$ bounded by $1$ and $\frac{1}{2}$, respectively, the conclusion follows. ∎

Acknowledgements

The work of MD was supported by F.R.S.-FNRS and by the CNRS-Momentum program.