Green's function for elliptic systems: Moment bounds

1.   Introduction

This paper is a contribution to the recently very active area of quantitative stochastic homogenization of second order uniformly elliptic operators, the main goal of which is to quantify how close is the large scale behavior of the heterogeneous operator $(-\nabla \cdot A(x) \nabla)^{-1}$ to the behavior of the constant-coefficient solution operator $(-\nabla \cdot A_{\rm hom} \nabla)^{-1}$. Here $A(x)$ stands for the non-constant (random) coefficient field defined on $\mathbb{R}^d$ and $A_{\rm hom}$ is called the matrix of homogenized coefficients.

As originally realized in the seminal papers by Papanicolaou and Varadhan [23] and, independently, by Kozlov [20], the central object in the homogenization of elliptic operators with random coefficients is the corrector $\phi_\xi$, defined for each direction $\xi \in \mathbb{R}^d$ as a solution of the following elliptic problem

in the whole space $\mathbb{R}^d$. The function $\phi$ is called corrector since it corrects the linear function $x \cdot \xi$, which is clearly solution to the constant-coefficient equation, to be a solution of the equation with heterogeneous coefficients. Since $\phi$ serves as a correction of a linear function, it should naturally be smaller, i.e., sublinear. Assuming the distribution of random coefficient fields $A$ is stationary (meaning it is shift-invariant in the sense that $A$ and $A(\cdot+x)$ have the same law for any $x \in \mathbb{R}^d$) and ergodic (meaning any shift-invariant random variable is almost surely constant, a property encoding decorrelation of coefficient fields over large scales), they showed that correctors are almost surely sublinear and can be used to define the homogenized coefficient

Since the problem is linear, it clearly suffices to study the $d$ correctors $\phi_i : = \phi_{e_i}$ for $i = 1,\ldots,d$. Secondly, borrowing notation from the statistical physics, $\left\langle {\cdot} \right\rangle $ stands for the ensemble average (expected value) with respect to a probability distribution on the space of coefficient fields $A$. Here and also later, we will often drop the argument $A$ in random quantities like the corrector as well as the argument $x$ related to the spatial dependence of quantities like the coefficient field $A$ or the corrector $\phi$.

Both mentioned works [20,23] were purely qualitative in the sense that they showed the sublinearity of the corrector in the limit of large scales without any rate. Assuming that the correlation of the coefficient fields decays with a specific rate (either encoded by some functional inequality like the Spectral Gap estimate or the Logarithmic Sobolev Inequality, or by some mixing conditions or even assuming finite range of dependence), one goal of quantitative theory is to quantify the sublinearity (smallness) of the corrector and consequences thereof.

Though the present result is purely deterministic in the sense that it translates the fact that the energy of any $A$-harmonic function satisfies a "mean value property" from some scale on (a fact that follows from the sublinearity of the corrector) into estimates on Green's function and its derivatives, we will first mention some recent results related to sublinearity of the corrector without going too much into details.

A central assumption in our result involves a minimal radius, a notion introduced by Gloria, Neukamm, and Otto [16]: for given fixed $\delta = \delta(d,\lambda) > 0$ (here $\lambda$ denotes the ellipticity contrast), the random variable $r_* = r_*(A)$ is defined as

Here $(\phi,\sigma)$ stands for the augmented corrector, where the new element $\sigma$ (called vector potential) can be used for obtaining good error estimates and was originally introduced for the periodic homogenization (see, e.g., [5]). Here and in what follows $B_R$ stands for a ball of radius $R$ centered at the origin and $\rlap{-} \smallint $ denotes the average integral.

The introduction of the minimal radius $r_*$ allows to split the arguments nicely into a deterministic and a stochastic part. For example, Gloria, Neukamm, and Otto [16] showed that

the sublinearity of the corrector, as encoded in the definition of the random variable $r_*$, implies the mean value property (large-scale $C^{1,0}$-regularity) for $A$-harmonic functions. This means that for $R \ge r \ge r_*$ and for any $A$-harmonic function $u$ on $B_R$ (i.e., a solution of $-\nabla \cdot A \nabla u = 0$ in $B_R$) one has

The idea that the (large-scale) regularity theory of $A$-harmonic functions is closely related to the error estimates in stochastic homogenization was first realized by Armstrong and Smart [4]. Building on the work of Avellaneda and Lin [5] in the periodic setting, Armstrong and Smart studied homogenization of uniformly convex integral functionals under an independence assumption. For a quantity similar to $r_*$ (which they denoted $\mathcal{Y}$), assuming finite range of dependence of the coefficient fields they obtained almost Gaussian bounds of the form $\left\langle { \exp(\mathcal{Y})^s } \right\rangle < \infty$ for any $s<d$ (see also [1,2,3] for further developments). While $\mathcal{Y}$ and $r_*$ are quite related, the definition of $r_*$ is more intuitive and explicit, and as such it is easier to work with.

Assuming that the ensemble on the coefficient fields satisfies a coarsened version of the Logarithmic Sobolev Inequality, Gloria, Neukamm, and Otto [16] showed that the minimal radius $r_*$ has stretched exponential moments

where $0 \leq \beta < 1$ appearing in the exponent is related to the coarsening rate in the Logarithmic Sobolev Inequality. Observe that in the case $\beta = 0$, i.e., the case when we consider the Logarithmic Sobolev Inequality without coarsening, this is the optimal Gaussian bound.

Recently, reviving the parabolic approach used in the discrete setting [17], which has the benefit of conveniently disintegrating contributions to the corrector from different scales, Gloria and Otto [18] obtained a similar results assuming the coefficient fields have finite range of dependence. As a by-product, assuming finite range of dependence Gloria and Otto got the estimates for the minimal radius $r_*$ with optimal stochastic integrability of the form

As already said, using completely different methods, such almost Gaussian bounds for a related quantity $\mathcal{Y}$ in the case of scalar but possibly nonlinear equations were obtained before by Armstrong and Smart [4].

Finally, on the other side of the spectrum, Fischer and Otto [14] combined Meyer's estimate together with sensitivity analysis to show that for strongly correlated coefficient fields (more precisely, they consider coefficient fields which are $1$-Lipschitz images of a stationary Gaussian field with correlations bounded by $|x|^{-\beta}$, where $\beta > 0$ has to be much smaller that $1$, so that Meyer's estimate holds for $p = 2+2\beta$) it holds

In the present paper we will obtain deterministic estimates for the Green's function based on the minimal radii $r_*$ at different points. More precisely, fixing two points $x_0, y_0 \in \mathbb{R}^d$, we take as the input the coefficient field $A$ and the corresponding minimal radii $r_*(x_0), r_*(y_0)$ (here and in what follows $r_*(x_0)$ stays for the minimal radius $r_*$ of the shifted coefficient field $A(\cdot - x_0)$), and produce estimates on the Green's function $G$ and its derivatives $\nabla_xG, \nabla_yG, \nabla_x\nabla_yG$, averaged over small scale around the points $x_0$ and $y_0$. This averaging is necessary since we do not assume any smoothness of the coefficient fields.

An obvious advantage of the present approach is that it clearly separates the random effects, described by $r_*$, from the average large-scale behavior of the Green's function. Indeed, our result is modular in the sense that it can be combined with estimates on stochastic moments of $r_*$, as for example obtained under different ergodicity assumptions on the ensemble in [16,18], to yield estimates on the Green's function and its derivatives for different models of random coefficients.

Our only goal in this paper is to obtain bounds, and not to show existence (or other properties) of the Green's function. In fact, a well known counterexample of De Giorgi [11] shows that there are uniformly elliptic coefficient fields for which the Green's function does not exist. Nevertheless, as recently shown in [10] by Conlon, Otto, and the second author, this is not a generic behavior. More precisely, in [10] they show that for any uniformly elliptic coefficient field $A$ the Green's function $G = G(A;x,y)$ exists at almost every point $y \in \mathbb{R}^d$, provided the dimension $d \ge 3$. Therefore, in the case $d \ge 3$, we will assume that the Green's function $G(A;\cdot,y) \in L^1_{\textrm{loc}}(\mathbb{R}^d)$ exists, at least in the almost everywhere sense (i.e., for a.e. $y \in \mathbb{R}^d$), and focus solely on the estimates. Since in $\mathbb{R}^2$ the Green's function does not have to exist, but its "gradient" can possibly exists, using a reduction from $3D$ (where the Green's function exists) in Section 4 we construct and estimate $\nabla G$.

There are several works studying estimates on the Green's function in the context of uniformly elliptic equations with random coefficients. Using De Giorgi-Nash-Moser approach for a parabolic equation (which is naturally restricted to the scalar case), Delmotte and Deuschel [12] obtained annealed estimates on the first and second gradient of the Green's function, in $L^2$ and $L^1$ in probability respectively, under mere assumption of stationarity of the ensemble (see also [22] for a different approach). Using different methods, Conlon, Otto, and the second author [10] recently obtained similar estimates, together with other properties of the Green's function, but without the restriction to the scalar case. Compared to that work, our more quantitative assumption in terms of $r_*$ implies local and quenched estimates, and not only estimates in average. In the discrete setting, assuming that the spatial correlation of the coefficient fields decays sufficiently fast to the effect that the Logarithmic Sobolev Inequality is satisfied, Marahrens and Otto [21] upgraded the Delmotte-Deuschel bounds to any stochastic moments. Recently, again only in the scalar case this work was extended by Gloria and Marahrens [15] into the continuum setting.

Both works [15,21] used De Giorgi-Nash-Moser-type argument, and as such were restricted to a single equation. In contrast, our result is not restricted to the scalar case, a reason why we had to develop different techniques to obtain the estimates.

Before we state the main result, let us mention other works relating the smallness of the corrector and the properties of solutions to the heterogeneous equation. Together with Otto [8], we compare the finite energy solution $u$ of

with $g \in L^2(\mathbb{R}^d;\mathbb{R}^d)$ being supported in a unit ball around the origin, with twice corrected solution $u_\textrm{hom}$ of the homogenized equation

Here by twice corrected we mean that first the right-hand side $g$ from the heterogeneous equation is replaced by $\tilde g = g(\textrm{Id} + \nabla \phi)$ in the constant-coefficient equation, and second, we compare $u$ with $(\textrm{Id} + \phi_i \partial_i)v$ at the level of gradients. Using duality argument together with a compactness lemma, this gives an estimate of the difference between $\nabla_x \nabla_y G(x,y)$ and $\partial_i \partial_j G_\textrm{hom} (e_i + \nabla \phi_i(x)) \otimes (e_j + \nabla \phi_j(y))$ (averaged over small balls both in $x$ and $y$). In order to get such estimates, it is not enough to assume that the corrector is at most linear with small slope (as in (1)), but rather we need to assume that for some $\beta \in (0,1)$, it grows in the $L^2$-sense at most like $|x|^{1-\beta}$:

Compared to the condition (1) which we use in the present paper, the above condition (2) is obviously stronger. Indeed, while for example $r_*$ is almost surely finite under pure assumptions of stationarity and ergodicity of the ensemble, this is not true for $r_{*,\beta}$. Moreover, while one would believe that the stochastic estimates on $r_*$ (e.g. those in [16]) can be with some farther work modified to estimates on $r_{*,\beta}$, we are not aware of any such results.

Hence, in comparison with the present work, in [8] we get a stronger statement (since we estimate the difference between the heterogeneous Green's function and corrected constant-coefficient Green's function while in the present paper we only control the heterogeneous Green's function alone), at the expense of stronger assumptions on the smallness of the corrector and a more involved proof. More precisely, here we show that the second mixed derivative of the Green's function $\nabla_x \nabla_y G$ behaves like $C|x-y|^{-d}$ (clearly this estimate is sharp in scaling since it agrees with the behavior of the constant-coefficient Green's function), while in [8] we show that the homogenization error, i.e., the difference between $\nabla_x \nabla_y G$ and twice corrected mixed second derivative of the constant-coefficient Green's function, is estimated by $C|x-y|^{-(d+\beta)}$ -that means we gain a factor of $|x-y|^{-\beta}$, where the exponent $\beta \in (0,1)$ is the one appearing in (2).

Since we are dealing with linear equations, we make use of a duality argument, first introduced by Avellaneda and Lin in [5]. This allows us to obtain estimates on the $L^2$-norm (in one variable) of the averages (in the other variable) for $G$, $\nabla G$ and $\nabla_x\nabla_y G$. In particular, the main pivotal estimate is the one on the $L^2$-norm in both space variables for the mixed derivatives $\nabla_x \nabla_y G$. Since we do not assume any smoothness on the coefficient field $A$, and therefore we may not appeal to Schauder estimates, we obtain an off-diagonal $L^2$-estimate in both variables out of estimates on averaged quantities by what we call a compactness lemma for $A$-harmonic functions. This ingredient was first introduced in [8] and later exploited and extended also in [10]. In this work, we make use of the same result and provide a different proof which is based on a refined iteration of Caccioppoli's inequalities. In our opinion, this new argument admits an easier extension outside of the $L^2$-framework, as it is required in the case of non-uniformly elliptic coefficients as considered in [6]. We plan to investigate this question in the future. Besides, we remark that a simplified version of the previous compactness argument appears already in the works by Otto and the first author [9] and in [19]: There, the quantities involved are stationary (in fact they are related to the gradient of the corrector) and, under this assumption, the proof of such compactness result turns out to be much simpler. In our case, as well as in [8] and [10], one deals with Green's functions and it is no more possible to appeal to stationarity in one space variable, since a translation in space for these objects involves also a translation of the singularity.

Last, let us mention the work of Otto and the authors [7], where we push farther the results of [8] using higher order correctors. The second and higher order correctors were introduced into the stochastic homogenization setup by Fischer and Otto [13], in order to extend the $C^{1,\alpha}$ regularity estimates on large scales [16] to $C^{2,\alpha}$ estimates and $C^{k,\alpha}$ estimates respectively. In [7], under the assumption of smallness of the corrector, we obtain two results about $A$-harmonic function in exterior domains: first, for any integer $k$ we construct a finite dimensional space of functions such that the distance between any $A$-harmonic function in the exterior domain and this space is bounded by $C|x-y|^{-(d+k)}$ (this statement can be seen as an analogue of Liouville statements for finite energy solutions in the exterior domain). Second, assuming smallness of both the first and second order augmented correctors (i.e., also including the second order vector potential, which we had to introduce), compared with [8] we improve by $1$ the exponent in the estimate between the solution of the heterogeneous equation in the exterior domain and some corrected solution of the constant-coefficient equation.

The paper is organized as follows: In the next section we will state our assumptions together with the main result, Theorem 1, and its corollaries, Corollary 1, Corollary 2, and Corollary 3. In Section 3 we prove Theorem 1 and in Section 4 we give the argument for Corollary 1, which is the only corollary which does not immediately follow from the theorem.

Notation. Throughout the article, we denote by $C$ a positive generic constant which is allowed to depend on the dimension $d$ and the ellipticity contrast $\lambda$, and which may be different from line to line of the same estimate. By $\lesssim$ we will mean $\le C$. Finally, the integrals without specified domain of integration are meant as integrals over the whole space $\mathbb{R}^d$.

2.   The main result

We fix a coefficient field $A \in L^\infty(\mathbb{R}^d;\mathbb{R}^{d\times d})$, which we assume to be uniformly elliptic in the sense that

where $\lambda \in (0,1)$ is fixed throughout the paper. Then we have the following result:

Theorem 1. Let $d \ge 3$, let $A$ be a uniformly elliptic coefficient field on $\mathbb{R}^d$ in the sense of (3), and let $x_0,y_0 \in \mathbb{R}^d$ with $|x_0 - y_0| \ge 10$. For a point $x \in \mathbb{R}^d$, let $r_*(x) = r_*(A, x) \ge 1$ denote a radius such that for $R \ge r \ge r_*(x)$ and any $A$-harmonic function $u$ in $B_R(x)$ we have

Let $G = G(A;x,y)$ be the Green's function defined through

assuming it exists for a.e. $y \in \mathbb{R}^d$. Then we have

where $r_*'(y) = r_*(A^t, y) \ge 1$ denotes the minimal radius for the adjoint coefficient field $A^t$ at a point $y$.

Though the Green's function does not have to exist in $2D$, with the help of the Green's function in $3D$ we can at least define and estimate "its gradient & second mixed derivatives":

Corollary 1. Let $d = 2$. Let $A$ be a uniformly elliptic coefficient field on $\mathbb{R}^2$ in the sense of (3), such that for its extension into $\mathbb{R}^3$ of the form

there exists two points $\bar X, \bar Y \in \mathbb{R}^3$ so that the minimal radii for $\bar A^t$ and $\bar A$ at those points are finite, respectively (i.e., $r_*(\bar A^t,\bar X) < \infty$, $r_*(\bar A,\bar Y) < \infty$).

Then for a.e. $y \in \mathbb{R}^2$ there exists a function on $\mathbb{R}^2$, which we denote $\nabla G(A;\cdot,y)$, so that it satisfies in a weak sense

Moreover, given $x_0,y_0 \in \mathbb{R}^2$ with $|x_0 - y_0| \ge 10$, we have estimates for $\nabla G$ as well as for $\nabla_y \nabla G$:

Remark 1. Assuming that the coefficient field $A$ in the statement of Corollary 1 is chosen at random with respect to a stationary and ergodic ensemble, by the standard ergodic argument (see, e.g., [16]), applied in $3D$ to the ensemble obtained as a push-forward of the $2D$ ensemble through (9), the assumption on the finiteness of the minimal radii is almost surely satisfied (say with $\bar X = \bar Y = (0,0,0)$).

Remark 2. It is clear from the proof of Theorem 1 that all the above estimates, i.e. (5)-(11), are true also if the domains of integration $B_1(x_0)$ and $B_1(y_0)$ are replaced by larger balls with the corresponding radii $r_*$. Moreover, the radii of the balls could be even larger than the minimal radii (as long as these new radii are not larger than one third of a distance between centers of those balls), in which case we need to replace the minimal radii on the right-hand sides of those estimates with the actual radii of the balls.

Remark 3. The appearance of different minimal radii in (10) and (11) (in (10) the minimal radii are related to the equation in $2D$, while in (11) they are the minimal radii for the equation in $3D$) is not a typo. The reason is that while (10) is proved directly in $2D$, the proof of (11) passes through $3D$ -hence the need to consider the minimal radii in $3D$. In view of the relation $r_*(A,x) \le r_*(\bar A,\bar x)$, which easily follows from the fact that any $A$-harmonic function in $B_R \subset \mathbb{R}^2$ can be trivially extended to an $\bar A$-harmonic function in $\bar B_R \subset \mathbb{R}^3$, the estimate (11) seems to be less optimal.

For notational convenience we state the result for a single equation. Since in the proof of Theorem 1 we do not use any scalar methods (like for example De Giorgi-Nash-Moser iteration), the result holds also in the case of elliptic systems -for that one just considers that $u$ has values in some finite-dimensional Hilbert space. Naturally, in that case all the constants will depend on the dimension of this Hilbert space.

Using the Gaussian bounds on $r_*$ for the case of coefficient fields with finite range of dependence, which were obtained recently in [18], Theorem 1 implies the following bounds:

Corollary 2. Suppose $\left\langle {\cdot} \right\rangle $ is an ensemble of $\lambda$-uniformly elliptic coefficient fields which is stationary and of unit range of dependence, and let $d \ge 2$. Then there exist $C(d,\lambda)$ such that for every two points $x_0,y_0 \in \mathbb{R}^d$, $|x_0 - y_0| \ge 10$, and every $\epsilon > 0$ we have

and in $d \ge 3$ also

In the case of coefficient fields with stronger correlations we can use the result from [16]:

Corollary 3. Suppose $d \ge 2$, and that the ensemble $\left\langle {\cdot} \right\rangle $ is stationary and satisfies a logarithmic Sobolev inequality of the following type: There exists a partition $\{ D \}$ of $\mathbb{R}^d$ not too coarse in the sense that for some $0 \le \beta < 1$ it satisfies

Moreover, let us assume that there is $0 < \rho \le 1$ such that for all random variables $F$

where the carré-du-champ of the Malliavin derivative is defined as

Then there exists a constant $0 < C < \infty$, depending only on $d, \lambda, \rho, \beta$, such that

and in $d \ge 3$ also

3.   Proof of Theorem 1

The proof is inspired by a duality argument of Avellaneda and Lin [5,Theorem 13], which they used to obtain Green's function estimates in the periodic homogenization. After stating and proving two auxiliary lemmas, we first prove the estimate on the second mixed derivative (5). Then, (6) will follow from (5) using Poincaré inequality and one additional estimate. Next we observe that (7) can be obtained from (6) by replacing the role of $x$ and $y$, which can be done by considering the adjoint $A^t$ instead of $A$. Finally, (8) will follow from (7) in a similar way as (6) follows from (5).

We thus start with the following two auxiliary lemmas. The first one is very standard:

Lemma 1 (Caccioppoli inequality). Let $\rho > 0$, $\delta > 0$, and let $u$ be a solution of a uniformly elliptic equation $-\nabla \cdot A\nabla u = 0$ in $B_{(1+\delta)\rho}$. Then

for any $c \in \mathbb{R}$.

Proof. By considering $u-c$ instead of $u$, it is enough to show estimate (12) with $c = 0$. We test the equation for $u$ with $\eta^2 u$, where $\eta$ is a smooth cut-off function for $B_\rho$ in $B_{(1+\delta)\rho}$ with $|\nabla \eta| \lesssim (\delta\rho)^{-1}$, use (3) and Young's inequality to get

Since $|\nabla \eta| \le \frac{C}{\rho \delta}$, (12) immediately follows.

Lemma 2. Let $R_0 \ge r_*(0)$, and let $u$ be an $A$-harmonic function in $B_{R_0}$. Then we have

Proof. Throughout the proof we write $r_*$ instead of $r_*(0)$. We assume that $2r_* < R_0$, since otherwise (13) is trivial. For $r \in [r_*,R_0]$ we denote $u_r : = \rlap{-} \smallint _{B_r} u$. We have

Hence, to prove (4) it is enough to show

To prove it, we use the following estimate

which in fact holds for any function $u \in W^{1,2}(B_{2r})$.

We first argue how to obtain (14) thanks to estimate (15): Let $n \ge 0$ be the largest integer that satisfies $2^n r_* \le R_0/2$; using (15) multiple times we get

Using Jensen's inequality and the fact that $R_0 \le 2^{n+2}r_*$ we get

Combination of the two previous estimates then gives (14).

It remains to prove (15). Using Jensen's and Poincaré's inequalities we get

3.1.   Proof of (5)

We denote $R_0 : = |x_0 - y_0|/3$. We split the proof of (5) into 4 steps. In the first step we show that

for any $\rho \in [1,2]$ and any functional $F_\rho$ on $L^2(B_\rho(x_0))$ which satisfies

for any $v \in W^{1,2}(B_\rho(x_0))$. In the second step, using Neumann eigenfunctions on a ball, we define a family of functionals $F_k$ satisfying (17), which will play the role of Fourier coefficients. Using these $F_k$ we then estimate $\int |v|^2$ with a sum of $N$ terms $|F_k(v)|^2$ plus a residuum in the form $\frac{1}{\lambda_N} \int |\nabla v|^2$. Here, $\lambda_N$ denotes the $N$th Neumann eigenvalue of Laplacian on a ball. Using (17) together with the second step we get an estimate on $\int | \nabla_x \nabla_y G|^2$ in terms of a (good) term and a small prefactor times $\int | \nabla_x \nabla_y G|^2$, integrated over a slightly larger ball. In the last step, we apply iteratively this estimate.

Step 1. Proof of (16) (inspired by the duality argument of Avellaneda and Lin [5]).

Let $f \in L^2(B_{R_0}(y_0);\mathbb{R}^d)$, and let $u$ be the finite energy solution of

in $\mathbb{R}^d$, for which holds the energy estimate

Then on the one hand, the Green's function representation formula yields

If $r_*(x_0) \le R_0$ (w. l. o. g. we assume $r_*(x_0) \ge \rho$), we use this in (4) to get

If $r_*(x_0) \ge R_0$, we simply have

Since $F_\rho$ is linear, using (19) we have

where the dot means that $F_\rho$ acts on the first variable. The previous relations then give

Using definition of the norm $L^2(B_{R_0}(y_0))$ by duality we get

Let $r_*'(y_0)$ play the same role as $r_*(x_0)$ but for the adjoint equation. In the case $r_*'(y_0) \le R_0$, since the mapping $y \mapsto F_\rho(\nabla_x G(\cdot,y))$ solves the adjoint equation $-\nabla \cdot A^t \nabla F_\rho(\nabla_x G(\cdot,y)) = 0$ in $B_{R_0}(y_0)$, arguing similarly as we did after (19), an analogue of (4) solutions of the adjoint equation implies

If $r_*'(y_0) \ge R_0$, we get the same conclusion for free:

Step 2. Let $\rho \in [1,2)$ and $\delta > 0$ be fixed such that $(1+\delta)\rho \le 2$. For $n \ge 0$, let $F_n$ denote the functional on $L^2(B_{(1+\delta)\rho})$, defined as an inner product of the $n$th eigenfunction of Neumann Laplacian with $v$, and let $\lambda_n$ be the associated eigenvalue. By writing any $v \in W^{1,2}(B_{(1+\delta)\rho})$ with zero mean in terms of Neumann eigenfunctions (which can be done since they form an orthonormal basis in $L^2$) we see that for any $N \ge 1$

where the last inequality follows from the monotonicity of $\lambda_k$. Hence we see

where we used that $\lambda_k \ge 1$.

Step 3. Combination of Step 1 and Step 2 (applied to $\nabla_y G(\cdot,y)$) and use of (21) yields

where we defined $\mathcal{G}_{\textrm{avg}}(y) = \rlap{-} \smallint _{B_{(1+\delta)\rho}(x_0)} \nabla_y G(x,y) {\rm d}x$.

Step 4. For a given sequence $\delta_k > 0$ such that $\rho \Pi_{k = 1}^\infty (1+\delta_k) \le 2$ we consider the following iteration procedure. Let $\rho_0 : = 1$, and for $k\ge 1$ set $\rho_k : = (1+\delta_k)\rho_{k-1}$. We denote

For any $N_k \ge 1$, estimate (23) in Step 3 yields

where the values of $\delta_k$ and $N_k$ are at our disposal. We choose $\delta_k : = (2k)^{-2}$ and $N_k : = \alpha k^{2d} 2^d$. Since $\Pi_{k = 1}^\infty (1+\frac{1}{4k^2}) \sim 1.46 \le 2$, for this choice of $\delta_k$ for all $k\ge 1$ we have $\rho_k \in [1,2]$. Using lower bound on the Neumann eigenvalues for the ball in the form $\lambda_k \ge Ck^{\frac{2}{d}}$ (in the case of a cube one can use trigonometric functions to explicitly write down the formula for eigenfunctions and eigenvalues; for the ball, one uses the monotonicity of the eigenvalues with respect to the domain, which follows from the variational formulation of the eigenvalues), we can find large enough $\alpha$ such that the prefactor in front of $M_{k+1}$ above satisfies

For this choice (24) turns into

Iterating this we get

Assuming we have $\sup_k M_k < \infty$, we send $K \to \infty$ to get

Since the sum on the right-hand side is summable, we get that $M_1 \lesssim 1$.

It remains to justify the assumption $\sup_k M_k < \infty$. For any $\Lambda \ge 1$, let $\chi_\Lambda(y)$ be the characteristic function of the set $\left\{ y \in B_1(y_0) : \int_{B_2(x_0)} |\nabla_x \nabla_y G(x,y)|^2 \le \Lambda \right\}$. Using the previous arguments, applied to $\nabla_x\nabla_y G(x,y) \chi_\Lambda(y)$, we get that

where the right-hand side does not depend on $\Lambda$. Now we send $\Lambda \to \infty$, and get

by the Monotone Convergence Theorem. This completes the proof of (5).

3.2.   Proof of (6)

We first observe that using Poincaré's inequality we can control the difference between $\nabla_y G$ and its averages over $B_1(x_0)$ by the $L^2$-norm of $\nabla_x \nabla_y G$, which we already control by (5). Hence, to obtain (6) it is enough to estimate averages $\rlap{-} \smallint _{B_1(x_0)} \nabla_y G(x,y) {\rm d}x$. Such estimate will follow from an analogue of (16) applied to one particular functional $F$. Since in this setting we need to work with $\int |u|^2$ and not with previously used $\int |\nabla u|^2$, we will need to use Lemma 2.

Step 1. By Poincaré inequality in the $x$-variable we have

By the triangle inequality we have

and so (6) follows from (25) provided we show

Step 2. Proof of (26). Similarly as for (5), consider arbitrary $f \in L^2(B_{R_0}(y_0);\mathbb{R}^d)$ and the finite energy solution $u$ of

in $\mathbb{R}^d$, which satisfies the energy estimate

Let $F$ be a linear functional on $L^2(B_1(x_0))$ such that $|F(v)|^2 \le \int_{B_1(x_0)} |v|^2$ for any $v \in L^2(B_1(x_0))$. Then, if $r_*(x_0) \leq R_0$

If otherwise $r_*(x_0) > R_0$, we do not need anymore to appeal to Lemma 2 and may directly bound

and proceed as in the previous inequality. As before, we use linearity of $F$ and write

Since $f \in L^2(B_{R_0}(y_0);\mathbb{R}^d)$ was arbitrary, combination of the two previous estimates yields

As before, it remains to argue that by going from $\int_{B_{R_0}(y_0)}$ to $\int_{B_1(y_0)}$ we gain a factor $R_0^{-d}$. We define $v(y) : = F(G(\cdot,y))$, and observe that $-\nabla\cdot A^t \nabla v = 0 $ in $B_{R_0}(y_0)$, where $A^t$ denotes the adjoint coefficient field. Then by definition of $v$ estimate (28) implies

For the choice $F(v) = \rlap{-} \smallint _{B_1(x_0)} v$ (29) is exactly (26).

3.3.   Proof of (8)

Similarly to the proof of (6), we use Poincaré's inequality (Step 1) to show that (8) follows from (7) provided we control averages of $G$ (Step 2).

Step 1. By Poincaré's inequality in the $x$-variable we have

Then by the triangle inequality we have

and so (8) follows provided we show

Step 2. Proof of (30). Similarly as for (6), consider arbitrary $f \in L^2(B_{R_0}(y_0))$, but this time $u$ being a finite energy solution of

in $\mathbb{R}^d$. In order to get the energy estimate, we test the equation with $u$ to obtain:

and so

We point out that compared to the proof of (5) or (6), we got additional $R_0^2$ due to the right-hand side of the equation being $f$ and not $\nabla \cdot f$.

Let $F$ be a linear functional on $L^2(B_1(x_0))$ such that $|F(v)|^2 \le \int_{B_1(x_0)} |v|^2$. If $r_*(x_0) \leq R_0$, then

If otherwise $r_*(x_0) > R_0$, then we directly bound

and proceed analogously to the other case. Using the Green's function representation formula we have $u(x) = \int_{B_{R_0}(y_0)} G(x,y) f(y) {\rm d}y$, and thus the linearity of $F$ yields

Since $f \in L^2(B_{R_0}(y_0))$ was arbitrary, we may combine the two previous estimates and conclude

As before, it remains to argue that by going from $\int_{B_{R_0}(y_0)}$ to $\int_{B_1(y_0)}$ we gain a factor $R_0^{-d}$. We define $v(y) : = F(G(\cdot,y))$, and observe that $-\nabla A^t \nabla v = 0$ in $B_{R_0}(y_0)$.

Now we use Lemma 2 with $v$ to get

For the choice $F(v) = \rlap{-} \smallint _{B_1(x_0)} v$, relation (33) is exactly (30).

4.   Proof of Corollary 1

We provide a generalization of (6)-(7) in the two-dimensional case. When $d = 2$, the Green's function for the whole space $\mathbb{R}^2$ does not have to exist; nevertheless, we may give a definition for $\nabla G$ via the Green's function on $\mathbb{R}^3$. To this purpose we introduce the following notation: If $\bar x \in \mathbb{R}^3$, we write $\bar x = (x, x_3) \in \mathbb{R}^2 \times \mathbb{R}$ and denote by $\bar B_r \subset \mathbb{R}^3$ and $B_r \subset \mathbb{R}^2$ the balls of radius $r$ and centered at the origin. For a given bounded and uniformly elliptic coefficient field $A$ in $\mathbb{R}^2$, recall that its trivial extension $\bar A$ to $\mathbb{R}^3$ was defined in (9) by

and the three-dimensional Green's function $\bar{G} = \bar{G}(\bar{A};\bar{x},\bar{y})$ is defined as a solution of

It will become clear below that the argument for the representation formula for $\nabla G$ through $\nabla_x \bar G$ calls for the notion of pointwise existence in $\bar y \in \mathbb{R}^3$ of the Green's function $\bar G (\bar A; \cdot , \bar y)$. As mentioned in Section 1, in the case of systems we may only rely on a definition of the Green's function for almost every singularity point $\bar y$. Therefore, differently from the previous sections, we need to bear in mind this weaker notion of existence of $\bar G$.

Step 1. We argue that for almost every $y \in \mathbb{R}^2$ the function $\nabla G$ (since $G$ does not exist, $\nabla G$ should be understood as a symbol for a function and not as a gradient of some function $G$), defined through

satisfies for every $\zeta \in C^\infty_0(\mathbb{R}^2)$

i.e., in a weak sense it solves $-\nabla_x \cdot A \nabla G(A; \cdot, y) = \delta(\cdot -y)$.

By definition of $\bar G(\bar A; \cdot, \cdot)$, we have for almost every $\bar y\in \mathbb{R}^3$ and every $\bar \zeta \in C^\infty_0( \mathbb{R}^3)$

Thus, for any $\bar \rho \in C^\infty_0(\mathbb{R}^3)$ this yields

We now choose a sequence $\{\bar \zeta_n \}_{n \in \mathbb{N}}$ of test functions $\bar \zeta_n = \eta_n \zeta$, with $\zeta = \zeta(x) \in C^\infty_0( \mathbb{R}^2)$ and $\eta_n = \eta_n(x_3)$ smooth cut-off function for $\{|x_3|< n\}$ in $\{|x_3| < n+1 \}$: From the previous identity and definition (9) it follows

We now want to send $n \rightarrow +\infty$ in the previous identity : By our assumptions on $\bar \rho$ and $\bar \zeta_n$, if we show that

then by the Dominated Convergence Theorem we may conclude that

and thus (35) by the arbitrariness of the test function $\bar \rho$ and the separability of $C^\infty_0(\mathbb{R}^2)$.

To argue inequality (36) we proceed as follows: We define a finite radius $M$ such that

and observe that inequality (36) is implied by

Since $\bar A$ is translational invariant in $x_3$, the minimal radius $r_*(\bar A^t,\cdot)$ is independent of $x_3$. Then, by the definition of $M$ and Remark 2 we have

provided $|X_3 - Y_3| \ge 3M$.

We now cover the cylinder $B_{M/2}(X) \times \mathbb{R}$ with countably many balls of radius $M$ centered at the points $(X, \pm Mn) \in \mathbb{R}^3$. By translational invariance we can w. l. o. g. assume that $Y_3 = 0$. We thus bound the integral in (37) by

We claim that $\nabla_{\bar x} \bar G( \bar A; \cdot, \cdot) \in L^1_{loc}( \mathbb{R}^3 \times \mathbb{R}^3)$, and so the first integral on the r. h. s. of the previous identity is finite.

Here we only sketch the idea why $\bar \nabla_{\bar x} \bar G \in L^1_{loc}(\mathbb{R}^3 \times \mathbb{R}^3)$; for the proof with all the details we refer to the proof of [10,Theorem 1]. To show that $\bar \nabla_{\bar x} \bar G \in L^1_{loc}$ it suffices to show that $\int_{\bar B_R(0)} \int_{\bar B_R(0)} |\bar \nabla_{\bar x}\bar G| < \infty$. In order to do that we observe that for given two distinct points $\bar x, \bar y \in \mathbb{R}^3$, the proof of Theorem 1 (without the use of $r_*$ to go to smaller scales; see also Remark 2) implies in $3D$

where $r = |\bar x-\bar y|/3$, which by Hölder's inequality turns into

Using a simple covering argument, the above estimate holds also in the case when the balls are replaced by cubes. Since $\bar B_R(0) \times \bar B_R(0)$ can be written as a null-set plus a countable union of pairs of open cubes $\bar Q_{r_n}(\bar x_n) \times \bar Q_{r_n}(\bar y_n)$, each with size $r_n : = |\bar x_n - \bar y_n|/3$ and such that each pair of points $(\bar x,\bar y) \in \bar B_R(0) \times \bar B_R(0)$ belongs to at most one such pair of cubes, we conclude

where we used that for $(\bar x,\bar y) \in \bar Q_{r_n}(\bar x_n) \times \bar Q_{r_n}(\bar y_n)$ we have $|\bar x - \bar y| \sim r_n$.

Going back to the second term on the right-hand side of (39), an application of Hölder's inequality in both variables $\bar x $ and $\bar y$ yields for the the sum over $n$

We now may apply to the r.h.s. the bound (38) and thus obtain

We have established (36).

Before concluding Step 1, we show that the representation formula (34) does not depend on the choice of the coordinate $y_3 \in \mathbb{R}$, namely that for almost every two values $y_{0,3},y_{1,3} \in \mathbb{R}$, for almost every $y_0, x_0 \in \mathbb{R}^2$

Without loss of generality we assume $y_{0,3} = 0$: Since by the uniqueness of $\bar G(\bar A; \cdot , \cdot)$, for every $\bar z \in \mathbb{R}^3$ and almost every $\bar x, \bar y \in \mathbb{R}^3$

by choosing $\bar z = (0, z_3)$ and using definition (9) for $\bar A$, we get

Let $x_0, y_0\in \mathbb{R}^2$ and $y_{1,3} \in \mathbb{R}^3$ be fixed: For every $\delta > 0$ we may write

and use (40) with $\bar z = (0, y_{1,3})$ to get

We now appeal to Lebesgue's theorem and conclude (38).

Step 2. Proof of (11). For this part we denote $r_x : = r_*(\bar A^t,(x_0,0))$ and $r_y : = r_*(\bar A, (y_0,0))$. By translational invariance of $\bar A$ and $\bar A^t$ we have $r_x = r_*(\bar A^t,(x_0,x_3))$ and $r_y = r_*(\bar A,(y_0,y_3))$ for any $x_3, y_3 \in \mathbb{R}$. Denoting $\mathcal{B} : = B_1(y_0) \times (-r_y/2,r_y/2)$, the independence of (34) from $y_3$ yields

Since $\mathcal{B} \subset \bar B_{r_y}((y_0,0))$, the previous identity implies

We define a sequence

and observe that

Since

where for simplicity we assumed $|x_0-y_0| \ge r_x$, we combine the three above relations to infer

which is exactly (11).

Concerning (10), there are two possible ways how to proceed. For the first we observe that (35) implies for every test function $\phi \in C^\infty_c(\mathbb{R}^2)$

where $f \in L^2(\mathbb{R}^2;\mathbb{R}^2)$ and $u$ is a solution of $-\nabla \cdot A \nabla u = -\nabla \cdot f$. Therefore we have that

and the proof of (5) applies verbatim. A different way would be to mimic the argument for (11), i.e., to define $\nabla_y \nabla G$ as an integral of the second mixed derivative of the Green's function in three dimension. Unfortunately, this way we would obtain the estimate where the minimal radii in $2D$ appearing on the right-hand side of (10) would need to be replaced (with possibly larger) minimal radii for $3D$.

Acknowledgment

We warmly thank Felix Otto for introducing us into the world of stochastic homogenization and also for valuable discussions of this particular problem. This work was begun while both authors were at the Max Planck Institute for Mathematics in the Sciences in Leipzig.