Green’s function for elliptic systems: existence and Delmotte–Deuschel bounds

This paper is divided into two parts: In the main deterministic part, we prove that for an open domain D⊂Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D \subset \mathbb {R}^d$$\end{document} with d≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \ge 2$$\end{document}, for every (measurable) uniformly elliptic tensor field a and for almost every point y∈D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y \in D$$\end{document}, there exists a unique Green’s function centred in y associated to the vectorial operator -∇·a∇\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\nabla \cdot a\nabla $$\end{document} in D. This result implies the existence of the fundamental solution for elliptic systems when d>2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d>2$$\end{document}, i.e. the Green function for -∇·a∇\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\nabla \cdot a\nabla $$\end{document} in Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^d$$\end{document}. In the second part, we introduce a shift-invariant ensemble ⟨·⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \cdot \rangle $$\end{document} over the set of uniformly elliptic tensor fields, and infer for the fundamental solution G some pointwise bounds for ⟨|G(·;x,y)|⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle |G(\cdot ; x,y)|\rangle $$\end{document}, ⟨|∇xG(·;x,y)|⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle |\nabla _x G(\cdot ; x,y)|\rangle $$\end{document} and ⟨|∇x∇yG(·;x,y)|⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle |\nabla _x\nabla _y G(\cdot ; x,y)|\rangle $$\end{document}. These estimates scale optimally in space and provide a generalisation to systems of the bounds obtained by Delmotte and Deuschel for the scalar case.

tensor field on R d : Our notion of ellipticity is slightly weaker than the more standard strong ellipticity, and this allows us to include the setting of linearised elasticity. We also need to impose an extra condition on the domain D in the case d = 2, namely that it has at least one bounded direction.
The paper is structured in two parts. In the first deterministic part we provide an existence and uniqueness result for Green functions. That is we prove for every a and almost every y ∈ D the function G D (a; ·, y) exists (in fact, in the case of systems it is a tensor field). In the case d > 2 this implies the existence of the full-space Green function, i.e. of G(a; ·, y) = G R d (a; ·, y). In the second stochastic part we introduce a shift-invariant probability measure on the coefficient fields a (or, equivalently, an ensemble · ), and consider when d > 2 the random field given by G(a; ·, y). In this random setting we establish for G optimal pointwise moment bounds: If · denotes expectation with respect to the ensemble and λ is the ellipticity ratio of a, we prove that with similar estimates for ∇ x G, and ∇ x ∇ y G.
In the scalar case it is a well-known result (see e.g. Grüter and Widman [20], Littman et al. [23]) that for any measurable and strongly uniformly elliptic a, the Green function exists and has optimal pointwise decay, e.g. as the Green function associated to the Laplacian (c.f. also the r.h.s. in (1)). This bound on the decay is a consequence of the De Giorgi-Nash-Moser theory, which does not hold in the case of systems. Moreover, when working with systems the existence of a Green's function is itself not ensured for this class of (possibly very rough) coefficient fields: A famous example of De Giorgi [10], disproving both the Liouville property and the C α -regularity theory for a-harmonic functions, also implies that there are strongly elliptic tensor fields and points y ∈ R d for which a Green's function centred in y does not exist.
Under additional regularity assumptions on the coefficient fields and/or on the domain D, the existence of the Green function has been proved also for systems: For any bounded and C 1 domain D, Fuchs [15] establishes existence of the Green function for strongly elliptic continuous coefficient fields a, and optimal pointwise bounds under the stronger assumption of Hölder continuity of a. Subsequently, Dolzmann and Müller [12] improve the previous result by obtaining for continuous a not only the existence of the Green function, but also optimal decay properties. In a series of works, Hoffman and Kim [21] and Kim and collaborators (see e.g. [22] and [7]) considerably weaken the assumptions on the domain D and on the regularity of a ( both in the elliptic and in the corresponding parabolic setting): In [21], they establish the existence of the Green function for an arbitrary open domain D ⊂ R d with d > 2, provided that the coefficient field is such that a-harmonic functions satisfy an interior Hölder continuity estimate (e.g. if a is VMO). In [22], Kang and Kim (see also Cho et al. [7] for the case d = 2) further develop the previous theory and in addition provide a necessary and sufficient condition on a in order to have for the Green function an optimal pointwise bound. We also mention that a result similar to [22] has been proved by Auscher and Tchamitchian [2] in the parabolic case via the introduction of a criterion [the Dirichlet Property (D)] for the parabolic Green function to have Gaussian bounds.
In this paper we adopt a different approach: Instead of restricting the class of coefficient fields a by further regularity assumptions, we show that the "bad" cases as the one of De Giorgi's example are exceptional for any coefficient field a. The main idea consists of extending our definition of the Green function to a two-variable object G D (a; ·, ·) which solves the equation almost surely in y: With this understanding, we may establish L 2 a-priori bounds in (x, y) on the gradients ∇ x G, ∇ y G and the mixed derivatives ∇ x ∇ y G. By an approximation procedure, we then obtain the existence statement.
The optimal stochastic estimates (c.f. (1)), obtained in the second part of this work, extend the inequalities established by Delmotte and Deuschel [11] for scalar equations to elliptic systems: Their methodology relies on the theory of De Giorgi, Nash and Moser for uniformly elliptic and parabolic PDEs in divergence form and therefore does not generalise to elliptic systems. Stochastic estimates on the (whole space) Green function have been largely used in the context of stochastic homogenization for elliptic PDE's in divergence form, in particular to quantify the decay of the homogenisation error, i.e. the difference between the solution of the heterogeneous equation and the solution of the homogenised equation. Inspired by the work of Naddaf and Spencer [26] on Gradient Gibbs Measures, the third author and Gloria and the third author, Gloria and Neukamm (see e.g. [17,19]) provide optimal estimates for the fluctuations (variance) of the corrector by controlling the decay of the moments of the gradients and mixed derivatives Green function, i.e. |∇ x G(·; x, y)| 2 p 1 p and |∇ x ∇ y G(·; x, y)| 2 p 1 p for p ≥ 1. There, an important role is played by the assumption on the ensemble of coefficient fields to satisfy a quantification of ergodicity in the form of a Spectral Gap. In [24], Marahrens and the third author rely on Logarithmic Sobolev's inequalities to upgrade the bounds obtained by Delmotte and Deuschel for any moment of G, ∇ x G and ∇ x ∇ y G and infer optimal estimates on the fluctuation of the homogenisation error. The result of this paper should therefore allow to obtain the optimal quantitative results cited above also in the case of stochastic homogenisation of elliptic systems. We also mention that in [4], Bella and the second author upgraded as well (1) to a bound for any moment in probability of G, its gradient and its mixed derivatives. Estimates (1) immediately imply optimal decay bounds for the averaged Green function G(a; ·, ·) . It is an interesting exercise to compare the methodology used in the present paper with the methodology used by the first author and Naddaf [8] to prove in the scalar case pointwise estimates on the averaged Green function. and their derivatives. While in this work we infer the stochastic bounds on the Green function directly from the deterministic existence result for G(a; ·, ·), in [8] a major part is played by the Fourier representation of the averaged Green function, which is a generalisation of the Fourier representation of the Green function for an elliptic PDE with constant coefficients. Naddaf and the first author then obtain Fourier space estimates strong enough to imply the pointwise estimates on averaged Green functions. We remark that the method in [8] does not make use of the scalar structure and therefore may be applied also to obtain a-priori estimates in the vectorial case. In the last section of this paper we indeed summarise how our main estimate in the proof of (1) can be obtained using this Fourier method.
We conclude this introduction by remarking that the issue of the regularity of averaged Green functions for PDEs with random coefficients plays an important role in statistical mechanics (in fact, [11] belongs itself to this context). In particular, it appears to have first come up in the study of the equilibrium statistical mechanics of the Coulomb dipole gas. Correlation inequalities for the dipole gas on the integer lattice were first obtained by Gawedzki and Kupiainen [16] by means of a complicated multi-scale induction argument known as the renormalisation group method [6]. A major drawback to the implementation of the renormalisation group method is that it always requires smallness in some parameters. In the case of the dipole gas this implies that the density of the gas must be extremely small, and with no reasonable estimate on how large the density is allowed to be. In [25] Naddaf and Spencer pioneered an alternative approach to estimating correlation functions for the dipole gas which was based on convexity theory. Their starting point was the observation that a correlation function closely related to the charge-charge correlation function for the dipole gas is equal to the integral over time of an averaged Green function for a parabolic PDE in divergence form with random coefficients. One of the main results of [16] follows immediately from this identity by application of a discrete version of the Aronson bounds [1] for Green functions of parabolic PDE. In addition, the argument applies for gases with density of order 1. The Aronson bounds make use of the main ingredient of Nash's argument for the De Giorgi-Nash-Moser theory, and are thus restricted to the scalar setting.
An important intuition in the study of the Coulomb gas is the notion of screening. That is the interaction between two particles of the gas is decreased by the presence of the other particles. In the case of two dipoles centred at x, y ∈ R d the interaction behaves like 1/|x − y| d at large distances. Hence one expects that for a dilute Coulomb gas of dipoles the charge-charge correlation between two dipoles at x, y also behaves like 1/|x − y| d at large distances. In the Naddaf-Spencer representation the charge-charge correlation function is approximately given by the averaged second derivative ∇ x ∇ y G of the Green function evaluated at (x, y). Hence one is motivated to expect pointwise bounds on averages of second derivatives of Green functions for certain parabolic PDE in divergence form with random coefficients, a conjecture formulated by Spencer [29] and proven in [11].

Notation and setting
In this section we introduce the elliptic systems' setting we are interested in and the notion of associated Green's function for a general open domain D ⊂ R d , d ≥ 2. In particular, we want to justify the scalar notation which is used throughout the paper. In order to do so, we first introduce our problem in a more general setting: In the rest of this section we consider a Hilbert space Y with dimY := m < ∞. We denote by zy and z · y respectively the inner product in Y and the natural one induced over Y d . In the same spirit, we write |z| = (zz) 1 2 and |y| = (y · y) 1 2 for z ∈ Y and y ∈ Y d .

Coefficient field
Let be the set of all symmetric and elliptic coefficient fields, i.e. all maps a such that We stress that (4) is a weaker condition than the uniform ellipticity assumption and it includes a larger class of elliptic systems to which belongs also the case of linearised elasticity. In this paper we mainly consider coefficient fields a ∈ , thus elliptic in the more general sense (4).

Generalised Green's function
For an open domain D ⊂ R d with d ≥ 2 and a given a ∈ , we refer to the map G D (a; ·, ·) : R d × R d → L(Y, Y ) as a Green's function, if there exists an α ∈ (0, d) 1 and R > 0 such that for every z ∈ Z d and for almost every y ∈ R d the application G(a, ·, y) satisfies in the sense that G D (a; ·, y) = 0 almost everywhere outside D or vanishes at infinity for D = R d , and that for any ξ ∈ Y , |ξ | = 1 it holds for every We note that if we obtain estimates over G(a; ·, y)ξ , independent of ξ , then we automatically deduce the same bounds for G(a; ·, y) itself. Therefore, as long as we estimate uniformly in ξ , instead of (9) we can adopt the formal notation Given G(a; ·, ·) as defined before, we may also consider ∇ y G D (a; ·, y) which, for almost every y ∈ R d , is formally a solution (with the same understanding of (8)) of where the distribution ∇ y δ(· − y) acts on any ζ ∈ C ∞ 0 (D) aŝ Throughout the rest of the paper we fix Y = R m with the canonical inner product and use the previous scalar-like notation. When no ambiguity occurs, we write ∇G D , ∇∇G D for the gradient ∇ x G(a; x, y) and the mixed derivatives ∇ x ∇ y G D (a; x, y). In the case D = R d , we also use the notation G(a; ·, ·) = G R d (a; ·, ·). In the estimates carried out in this paper, stands for ≤ C with a constant depending exclusively on the dimension d and the ellipticity ratio λ and thus in particular independent of the choice of the domain D; similarly, D stands for ≤ C with C depending on d, λ and also on the domain D: Unless stated otherwise, the dependence of C on the domain is merely through the size of the smallest bounded direction of D.
We remark that our definition of Green's function guarantees that for every coefficient field a ∈ and for every open domain D ⊂ R d with d ≥ 2, G D (a; ·, ·) is unique. More precisely, we have the following

Lemma 1 Let a ∈ and let D be an open domain in
. The proof of this result in the appendix is very similar to [18], Section A.3, Step 4.

Random coefficient fields
We restrict our definition of as (4) and (3) where the measurability is considered with respect to the coarsest σ -algebra F such that ∀ξ ∈ Y d , |ξ | = 1, the evaluation a →ˆ(a(x)ξ )χ(x)dx (12) is measurable for every χ ∈ C ∞ 0 (R d ) ( where R is equipped with the usual Borel σ -algebra). We define a random coefficient field by endowing the couple ( , F ) with a probability measure P, or equivalently by considering an ensemble · over symmetric, uniformly elliptic coefficient fields a. We assume the ensemble · to be stationary, namely that ∀z ∈ R d the coefficient fields {R d x → a(x)} and {R d x → a(x + z)} have the same distribution, and to be stochastically continuous, in the sense that for every x ∈ R d and δ > 0 This last condition ensures that the map ×R d (a, z) → a(·+z) ∈ is measurable with respect to the product topology of × R d . With this additional structure, we can consider the random map G : a → G(a; ·, ·) . We also remark that, by definition (12), F is countably generated and therefore for every p ∈ [1; +∞), the space

Main result and remarks
Throughout the paper, as a basic assumption, we consider the domain D ⊂ R d to be open and such that This condition basically ensures that if a function u is zero almost everywhere outside D, and admits weak derivatives up to order k, then the derivatives are as well almost surely zero outside the domain: This will prove to be useful when defining the approximate problem for (8), cf.(32), which calls for a higher order operator and thus Dirichlet boundary conditions also for the derivatives.
In this paper we mainly provide two existence results for the Green function in a domain D ⊂ R d with d ≥ 2. As introduced in Sect. 2, for a given a ∈ , we treat the Green function for a domain D as an object G D (a; ·, ·) in two space variables (x, y) ∈ R d × R d , which satisfies for almost every singularity point y ∈ R d the equation (8). It is with this generalised definition of Green function that we manage to prove its existence and uniqueness (cf. Lemma 1) and also to obtain optimal estimates for the L 2 -norm in both the space variables of G D (a; ·, ·), its gradient and its mixed derivatives.
In the first theorem, we show that if the domain D is open, and if d = 2 also bounded in at least one direction, then for every coefficient field a ∈ the Green function G D (a; ·, ·) exists; In particular, this result also includes the existence of the fundamental solution, i.e. the Green function for D = R d , with d > 2. In the case of open domains bounded in at least one direction and strongly elliptic (cf. (5) coefficient fields a ∈ , we also provide in Corollary 1 an improvement of the estimates obtained in Theorem 1, namely that the offdiagonal L 2 -norms of G D (a; ·, ·), G D (a; ·, ·) and ∇∇G D (a; ·, ·) decay exponentially in the unbounded directions. Corollary 2 deals with the introduction of a stationary ensemble on the space of coefficient fields and provides in the case of systems a generalisation, at the level of the first moments in probability, of the stochastic bounds obtained by Delmotte and Deuschel in [11] for the scalar case.
More precisely, we prove the following statements Theorem 1 Let d ≥ 2 and D ⊂ R d be a general open domain satisfying (13). Then (a) If d > 2, for every a ∈ there exists the Green function G D (a; ·, ·) and it satisfies for and for every z ∈ R d , (b) If d = 2 and D is bounded in at least one direction, then for every a ∈ the Green function G D (a; ·, ·) exists as well and satisfies the bounds (14)-(15), (18) and (16)- (17).
All the constants, with the exception of (18), depend also on the size of the smallest bounded direction of D and the bound (15) holds for radii R D 1.
be an open domain satisfying (13) and bounded in at least one direction. Then for every a ∈ satisfying also (5), there exists a constant C 1 D 1 such that for every z ∈ R d and R D 1 it holdŝ Corollary 2 Let d > 2 and let · be a stationary ensemble on . Then, the Green function G(a; ·, ·) for the whole space R d satisfies for almost every x, y ∈ R d the annealed pointwise bounds We recall that in all the previous inequalities and D stand for ≤ C with the constant C respectively depending on d and λ or on d,λ and the size of the smallest bounded component of D.
In the following remark we argue that the bounds (22) and (23) require the expectation · : Remark 1 (i) For d > 2, a coefficient field a ∈ and an associated Green's function G(a; ·, y) on R d , the bound implies that any finite energy a-harmonic application u is (locally) bounded. More precisely, the local boundedness of a-harmonic applications is also implied if assuming instead of (24) the weaker L 2 -bound |x−y|>R |∇G(a; x, y)| 2 R 2−d for every R > 0 and a.e. y ∈ R d .
While in the scalar case the bound (25) holds ( [17], Lemma 2.9), we cannot expect it to be true for every coefficient field a ∈ in the case of systems. The following example of De Giorgi [10] shows indeed that in d > 2, the unbounded function u : is locally of finite energy and a-harmonic for the symmetric and elliptic coefficient field (ii) Assuming that both (25) and |x−y|>R |∇∇G(a; x, y)| 2 R −d for every R > 0 and a.e. y ∈ R d (28) hold, implies that any a-harmonic application u is also locally Lipschitz. Hence, also in the scalar case both conditions (25) and (28) cannot be true for every coefficient field a ∈ ( [28], Example 3). (iii) For α ∈ (0, 1), even a suboptimal assumption on the decay of (28) aŝ cannot hold for every coefficient field a ∈ both in the scalar and in the systems' case. Indeed, (29) implies a Liouville property for a-harmonic functions, namely that for β ∈ ( α 2 , 1) for any a-harmonic u on R d where It can be shown that in the scalar case (30) doesn't even hold for uniformly elliptic and smooth coefficient fields ( [14], Proposition 21): For every ε > 0, there exists indeed a smooth a ∈ and an a-harmonic function u such that ffl |x|<R |u| 2 in the case of systems, De Giorgi's example (26) shows that a-harmonic functions can also (non trivially) vanish at infinity. We postpone the proofs of (i), (ii) and (iii) to the Appendix.
This paper is organised as follows: In Sect. 4 we give the argument for Theorem 1, part (a) and (b). The core of the proof for part (a), i.e. when d > 2, is an L 2 -off-diagonal bound for ∇∇G D and ∇G D , in both space variables x and y and depending only on the dimension and the ellipticity ratio. It is mainly obtained through a duality argument à la Avellaneda-Lin ( [3], Theorem 13) on standard energy estimates for solutions of −∇ · a∇u = ∇ · g, combined with an inner-regularity estimate for a-harmonic functions in the spirit of Lemma 4 of [5]. We stress here that this result is inspired by Lemma 2 of [5] and provides the new and pivotal ingredient for the first fundamental estimate for G D . This may be considered as the key ingredient for the whole argument of Theorem 1. Sobolev's inequality allows to extend the previous estimates also for ∇ y G D and G D . Finally, with the aid of rescaling and dyadic decomposition arguments, from the off-diagonal estimate on ∇∇G D we also infer bounds for ∇G D and G D close to the singularity x = y.
In dimension d = 2, it is well known that the fundamental solution, i.e. the Green function for D = R d , does not exist. We indeed restrict our attention to domains D which have at least one bounded direction: By substituting the scale-invariant Sobolev's inequality, which holds only for d > 2, with Poincaré-type inequalities, we may extend the arguments of part (a) to the two-dimensional case. We point out that the appeal to Poincaré-type inequalities introduces in the estimates for G D and its derivatives a dependence also on the minimal bounded direction of D.
We stress that our assumptions on include in this set also very rough coefficient fields for which the existence of G D is not a priori guaranteed. Therefore, we need to first approximate the problem (8), carry out and adapt to the approximate solution the aforementioned a priori bounds on G D , and then argue by standard weak-compactness of W 1,q loc -spaces. We approximate (8) through an ε-perturbation of the operator −∇ · a∇ with the hyper-elliptic term 2 and thus consider for ε > 0, a ∈ and y ∈ R d the problem The assumption (13) on the domain D and our understanding of the boundary conditions, i.e. that G ε,D vanishes almost surely outside D or vanishes at infinity for D = R d , imply that the same boundary conditions hold also for the higher-order derivatives up to index n − 1. For D = R d , the Dirichlet conditions on the boundary turn into the requirement for every ∂ α u with 0 ≤ |α| ≤ n − 1, to vanish at infinity. For n > d 2 , Riesz's representation theorem ensures the existence of a unique weak solution G ε,D for every ε > 0, a ∈ and singularity point y ∈ R d . Moreover, assuming n > d 2 + 1 also implies that there exists a unique ∇ y G ε,D (a; ·, y) which solves the approximate problem for (11) −∇ · a∇∇ y G ε,D (a; ·, y) + εL n ∇ y G D (a; ·, y) = ∇ y δ(· − y) in D ∇ y G ε,D (a; ·, y) = 0 o n ∂ D.
In Sect. 5 we provide the proof of Corollary 1 and Corollary 2; In the first corollary we show that in the case of domains with at least one bounded direction and strongly elliptic coefficient fields we improve the estimates of Theorem 1 again by a duality argument which this time relies on a refinement of the standard energy estimate for solutions of −∇ · a∇u = ∇ · g in domains which have a bounded direction. While the arguments of Theorem 1 and Corollary 1 are purely deterministic, in Corollary 2 we introduce a stationary ensemble · on and focus our attention on the fundamental solution G in d > 2 seen as a random map. The stationarity assumption on · provides an improvement of the estimates on G by upgrading the bounds of Theorem 1 from space-averaged in both variables x and y to annealed in a but pointwise in y. An a priori estimate for locally a-harmonic functions allows us to conclude the argument and obtain estimates averaged in a, but pointwise in x and y.
In the last section we present an alternative partial proof for Corollary 2 which makes use of the Fourier techniques developed in [8] and relies on a representation formula for the Fourier transform of the Green function. Finally, in the Appendix we give a self-contained proof of all the auxiliary results which are used in the arguments.

Proof of Theorem 1
Let a ∈ and D ⊂ R d be a generic open domain satisfying (13), with d ≥ 2. For a fixed y ∈ R d and ε > 0, we consider the approximate problem for (8) introduced in (32), i.e.
where L n is as in definition (31) and with n a fixed odd integer such that n > d 2 + 1.
Analogously we may consider solutions on outer domains by substituting in Definition 1 the domain {|x| < R} with {|x| > R}.
We start with two variants of Lemma 4 of [5]; while this last one is a statement for ensembles of locally a-harmonic functions, the following Lemma 2 takes into account the new perturbation term L n and the more general domain D. If d > 2, then Lemma 3 is a further generalisation to the case of functions solving −∇ · a∇u + εL n u = 0 on outer domains. We postpone the proofs of Lemma 2 and Lemma 3 to the Appendix. Lemma 2 For a radius R > 0 and a ∈ , we consider a σ -finite measure μ on functions u satisfying in the sense of Definition 1 with ε ≥ 0. Then we havêˆ| where the supremum runs over all linear functionals F bounded in the sense of with v satisfying (i) and (ii) of Definition 1.

Lemma 3 Let d > 2.
For a radius R > 0 and a ∈ , we consider a σ -finite measure μ on functions u with finite Dirichlet energy in {|x| > R} and satisfying in the sense of Definition 1 with ε ≥ 0. Then we havêˆ| where the supremum runs over all linear functionals F bounded in the sense of with v satisfying (i) and (ii) in the sense Definition 1 where the set {|x| < R} is substituted by {|x| > R}.
Analogously to Theorem 1, means ≤ C with a generic C = C(d, λ). We remark that the inequalitieŝˆ| hold trivially by conditions (36) and (39) and duality (in L 2 ). Roughly speaking, Lemma 2 and Lemma 3 state that the previous inequalities remain true also if we exchange in the r.h.s. the order of the integration in μ and the supremum over the functionals F. We may refer to the result of Lemma 4 of [5], which corresponds to Lemma 2 with D = R d and ε = 0, as a compactness statement for ensembles of locally a-harmonic functions. Indeed, as we show in the appendix, inequality (35) actually follows by an inner regularity estimate which allows to control the energy of an a-harmonic function u in an interior domain by the L 2 -norm on {|x| < 2R} of (− N ) − l 2 u for any even l ∈ N. Here, − N denotes the Laplacian with Neumann boundary conditions. This last estimate basically implies that in the space of locally a-harmonic functions, the local W 1,2 -norm (the strongest norm which is meaningful to consider for weak solutions of a variable-coefficient and uniformly elliptic second-order operator) is actually equivalent to much weaker norms, provided we consider a slightly bigger domain. Therefore, in this sense we may say that the space of locally a-harmonic functions is "close" to being finite-dimensional, in which case all the norms are equivalent. The previous lemmas state similar compactness results in the case of the approximate operator −∇ · a∇ + εL n .
Proof of Theorem 1 Throughout the whole proof we assume D ⊂ R d to be a generic open domain satisfying (13) which is also bounded in at least one direction if d = 2.
Given the Hilbert spaces is bounded thanks to (3) and coercive in the sense of thanks to (4). Let us first consider the case d > 2: Sobolev's embedding implies that i.e. B is non-degenerate. We now argue that B satisfies for every u ∈ X Thanks to the coercivity condition (40), inequality (42) is implied by which can be restated by passing to dual (i.e. Fourier) variables k aŝ As the conditions n > d 2 + 1 > d 2 and d > 2 imply that dk |k| 2 + |k| 2n < +∞, from (44) we obtain (43) and thus (42). Inequality (42) in particular yields that for the linear which implies by Riesz's representation theorem that there exists a unique G D,ε (a; ·, y) ∈ X weakly solving (32). As we have shown that for every a ∈ , the map D y → G ε,D (a, ·, y) is well defined, we now also show that for every a ∈ and y ∈ D there exists ∇ y G ε,D (a; ·, y), unique (weak) solution 4 of (33). We appeal again to Riesz's representation theorem: In this case, the linear functional that we need to bound with B is given byF(v) := ∇v(y). Once again, thanks to the coercivity condition (40) we conclude the argument if we show that or equivalently, by passing in Fourier variables, that Similarly to (44), we estimatê The fact that the solution of (33) is actually the y-gradient of G D (a; ·, ·) is rigorously proven by first showing that, on the one hand, by symmetry (47), the difference quotients 1 h G D,ε (a; ·, · + he i ) − G D,ε (a; ·, · + he i ) are uniformly bounded in X for h << 1, and thus weakly converge up to subsequences. On the other hand, letting h → 0 + in the equation solved by 1 h G D,ε (a; ·, y + he i ) − G D,ε (a; ·, y + he i ) , we recover (33) and conclude the argument by uniqueness of the solution.
Let us now assume that d = 2: Also in this case, by (3) and (4), B is respectively bounded and coercive in the sense of (40). Our assumption on D and the Dirichlet boundary conditions allow us to appeal to Poincaré's inequalitŷ and infer that B is non degenerate in the sense of (41). We note that once we prove that also in this case B satisfies inequalities (42) and (45), we may argue analogously to the case d > 2 and conclude that there exist unique G ε (a, ·, y) and ∇ y G ε (a, ·, y) solving in D respectively (32) and (33) . The argument used above for (45) is still valid; to show (42), we observe that by Poincaré's inequality and (40) it is enough to prove that We rewrite the previous inequality in Fourier variables aŝ Relying on our assumption n > d 2 + 1 > d 2 , we have that dk 1 + |k| 2n < +∞, and thus we infer (42).
In addition, we note that uniqueness and the symmetry of the operator −∇ · a∇ + εL n , cf. (2), yield for all a ∈ , z ∈ R d , R > 0, y ∈ R d and almost every x ∈ R d that Moreover, we claim that for every compactly supported f ∈ L 2 (R d ) and every g ∈ in the sense of Definition 1, then we have the representation formula By Hölder's inequality it is immediate to show that for the linear functional F : . Therefore, by Riesz's representation thoerem, there exists a unique solution in X of (50). It thus remains to show that the r.h.s. of (51) solves the equation: An application of Hölder's inequality, together with the assumptions on f , g and the fact that G ε,D (a; ·, y) ∈ X , implies that u in (51) is well-defined and belongs to X . By (32) and (33), u satisfies the boundary conditions and for every v ∈ X we havê We thus established (51).
Step 2: Uniform bounds for {G ε,D } ε>0 if d > 2. We presently argue that the family {G ε,D (a, ·, ·)} constructed in the previous step satisfies (14)- (15), and (16)-(17)- (18). By the properties (48)-(49), without loss of generality it is sufficient to fix z = 0 and R = 1, i.e. to prove that for all and that for every For a given L 2 -vector field g with support in {|x| > 2} ∩ D, the solution 5 of satisfies by (4) Thus, the previous estimate and the energy estimate respectively yield, thanks to the representation formula (51), that We now apply Lemma 3 to the families {G ε,D (a; ·, y)} |y|<1 and {∇ y G ε,D (a; ·, y)} |y|<1 , with functionals given by´·g and measure μ(dy) = dy |{|y|<1} . We specify that we are allowed to use Lemma 3 on both families since, by the previous step, they respectively solve the problems (32) and (33) and thus are a-harmonic in {|x| > 2} ∩ D for (μ-)almost every y such that |y| < 1. Therefore, from (57) and (58) we get which implies the bounds (52) and (56). We now turn to inequality (53): By the shift invariant property (48) and the invariance under scaling of the previous argument, it follows from (59) that for all w ∈ R d and r > 0 it holdŝ We appeal to the scale-invariant Sobolev Inequality in the exterior domain 6 {|x − w| > 8r } to obtain from (61) that Thus, Hölder's inequality in the x-variable yieldŝ We now cover the ball {|y| < 1} with the union of smaller balls of radius 0 Then, estimates (63) and (60) It follows from this that for any α > d − 2 and 0 < r < 1 and thus (53). We now claim that from (53) we obtain (55): If we smuggle in (55) the weight |x − y| α 2 q and apply Hölder's inequality first in x and then in y, we get |x|<1ˆ|y|<1 |∇ x G ε,D (a; x, y)| q + |∇ y G ε,D (a; x, y)| q 6 To show Sobolev's inequality in the outer domain {|x| > R} we argue as follows: By scale invariance, we may reduce ourselves to the domain {|x| > 1}; moreover, by standard approximation, we may assume u to be smooth and zero outside a ball big enough. We now extend u inside {|x| < 1} using the radial reflection x → x |x| 2 , apply Sobolev's inequality on the whole space and conclude by observing that, due to our choice of extension, the Dirichlet integral in {|x| < 1} can be controlled by the Dirichlet integral in {|x| > 1}.
and thus (55), as our assumption It only remains to establish (54): We first observe that if we prove an analogy of (63) for G ε,D , namely that for r > 0 |y|<rˆ8r <|x|<16r then by a scaling and covering argument similar to the one in (64) and (65) (14)- (15) and (18) satisfies by (4) the energy estimate λ´D |∇u| 2 ≤´|g| 2 and yields, thanks to the representation formula (51),ˆ| We now apply Lemma 2 to the family {∇G ε,D (a; x, ·)} {|x|>4}∩D , with functionals given bý ·g and measure μ(dx) = dx {|x|>4}∩D . We observe that we are allowed to use Lemma 2 on this family since, by (47), we can identify with ∇ x G ε,D (a; ·, x) constructed in Step 1 (with exchanged roles of the x and y variable). It follows from (33) that for x ∈ R d with |x| > 4, ∇ x G ε,D (a; ·, x) is solution of (37) in the domain {|y| < 2} ∩ D. Therefore, from (68) we get by Lemma 2 the desired bound (18). We remark that since the scale invariant Sobolev's inequality is no more available for d = 2 we cannot infer also (62). Appealing to our assumption on D to have at least one bounded direction, we may use as a replacement for Sobolev's inequality the following version of Poincarè-Sobolev's estimate: 8 Let D ⊂ R 2 be open and having at least one bounded direction. Then, for every 2 ≤ p < +∞, z ∈ R 2 and R > 0, it holds ˆ| for every u ∈ W 1,1 loc (R 2 ) and such that u = 0 almost everywhere outside D. Here, the constant depends on the size of the smallest bounded component of D.
With the same reasoning used in Step 2, once that we show that for every δ > 0, z ∈ R 2 and r > 0 we havê it follows by a covering argument, that and thus that for every α > d − 2 = 0 and for every α > d − 4 = −2 We may analogously argue for a general radius R and establish (72)-(73), and thus bound (14) for any R > 0. As shown in Step 2, these estimates also yield (16)-(17) Since the exponent p can be chosen arbitrarily large, we obtain (70) for ∇ y G ε,D . We now observe thatˆ| for every w and z such that {5R < |z − w| < 7R}. Indeed, this is implied by (74) and the inclusion For a fixed w ∈ R d , we choose n 1 balls of radius R which cover the annulus {5R < |y − w| < 7R} and whose centres {z i } n i=1 are contained in {5R ≤ |z − w| ≤ 7R}. Thus, from the previous inequality we infer By switching the labels x and y and using the symmetry property (47), this may be rewritten asˆ5 i.e. inequality (70) thanks to the arbitrariness of 2 ≤ p < +∞.
It thus remains to prove (71): By Poincaré's inequality in the x-variable we havê Therefore, thanks to (70), we conclude (71) once that we show that the second term on the r.h.s. of (76) satisfies for δ > 0 To do so, let us fix p and consider any function g ∈ L 2 (R d ) with supp(g) ⊂ {|x| < R} ∩ D. Let u be the solution of .
The arbitrariness of g allows to argue by duality that i.e. the desired bound (77) thanks to the arbitrariness of 2 ≤ p < +∞.
At last, we prove that the bound (18) implies (15): Modulo a change of coordinates, we may assume D ⊂ I × R, with I a bounded interval. Moreover, since by construction for almost every y / ∈ D, G ε,D (a; ·, y) = 0 almost surely in R 2 , we reduce ourselves to those z ∈ R 2 and R > 0 such that {|y − z| < R} ∩ D = ∅ and, without loss of generality we fix z = 0. Therefore, for every R D 1 the rectangle I × (−2R, 2R) is such that Since by (47) Step 4: Existence of G D (a, ·, ·). In this final step we do not distinguish between the cases d > 2 and d = 2. The uniform bounds (in ε) (18)- (14) and (16)-(17) for the family {G ε,D (a; ·, ·)} ε↓0 allow us to argue by weak-compactness that, modulo a subsequence, for Since G D (a; ·, ·) ∈ W 1,q loc (R d × R d ) with G D (a; ·, ·) ≡ 0 outside D × D, it follows respectively that for almost every y ∈ R d G D (a, ·, y) ∈ W 1,q loc (R d ), G D (a, ·, y) = 0 almost everywhere outside D.
We now show that for almost every y ∈ R d , the application G D (a, ·, y) solves (8): By construction of G ε,D , it holds indeed that for almost every y ∈ R d and every ζ ∈ C ∞ 0 (D) For every ρ ∈ C ∞ 0 (R d ), the previous identity yieldŝ so that for ε → 0, by weak convergence (80), we get The arbitrariness of the test function ρ ∈ C ∞ 0 (R d ) implies that for almost every y ∈ R d ∇ζ(x) · a(x)∇G D (a; x, y) = ζ(y).
We now appeal to the separability of C ∞ 0 (D) with respect to the C 1 topology to conclude that for almost every y ∈ R d and for every ζ ∈ C ∞ 0 (D) i.e. for almost every y ∈ R d a solution G D (a; ·, y) of (8) exists. Reasoning in the same way, from (33) and weak convergence (82) we also obtain that for every R > 0, z ∈ R d and almost every y ∈ {|y − z| > 2R}, the function ∇ y G D (a; ·, y) solves Furthermore, appealing to (80), (81) and (82) and the lower semiconinuity of the bounds (14)-(15), (18) and (16)- (17), we get that they hold also for G D (a; ·, ·); in particular, inequalities (14)- (15) imply that G D (a; ·, ·) satisfies bound (6) for any α ∈ (d − 2, d) as well as bound (7) for any R > 0 if d > 2 and for any R D 1 if d = 2. Thus, G D (a; ·, ·) is the Green function for the domain D. By uniqueness (cf. Lemma 1) and symmetry of the operator −∇ · a∇, cf.
(2), we also have that for all a ∈ , z ∈ R d , R > 0 and almost every x, y ∈ R d it holds Remark 2 We observe that also for G D holds a representation formula for weak solutions of with f ∈ L q (D), q > d and compactly supported. This may easily follow by uniqueness of the solution u (via Riesz's representation theorem) and the fact that the function is well defined and such that ∇û ∈ L 2 (D), thanks to the bounds (15), (16) and (17). Note thatû weakly solves (90): This may be shown as for (51) first for smooth test functions and then extended by standard approximation. We also have that for any f ∈ L 2 (D) and g ∈ [L 2 (D)] d with compact support, the weak solution of admits the representation whenever x is outside the support of both g and f . We first consider the family {u ε } ε>0 of solutions to the approximate problems (50) with the same r.h.s. : By standard weakcompactness arguments, (up to a subsequence) {u ε } ε>0 weakly converge in W 1,2 loc (D) to the solution u of (91). We thus conclude the identity above by using (51), together with (81)-(82), and the uniqueness of the (weak) limit.

Proof of Corollary 1 and 2
Proof of Corollary 1 Let D ⊂ R d with d ≥ 2 be as in the statement of Corollary 1. Modulo a change of coordinates, we can assume that there exists a bounded interval I ⊂ R such that D ⊂ I × R d−1 . In addition, without loss of generality we may suppose that |I | = 1: It will become clear along the proof that the estimates obtained depend on the size of I . For The main ingredient for the argument of Corollary 1 is the following elliptic regularity result ( [27], Lemma 2.2), adapted to elliptic systems with Dirichlet boundary conditions. We postpone its proof to the Appendix. (5). For g ∈ L 2 (D) d , let u solve (in the sense of Definition 1 with ε = 0 and R = +∞)

Lemma 4 Let D be as introduced above, and a ∈ such that it satisfies
Then, there exists a constant C 0 depending on d, λ (and the size of I ) such that for any We start by claiming that the previous lemma, together with an application of Lemma 2, yields (19). More precisely, we have for every z ∈ R d and R > 0 that We now apply to (94) Lemma 2, this time in the case ε = 0, with functionals given by´g, measure μ(dy) = exp |y −z | which implies inequality (93) since for R D 1 it holds the inclusion To obtain also (20) we argue similarly to Step 3. of Theorem 1; we first tackle the bound for the gradient of G D : Without loss of generality we may reduce ourselves to consider the case {|x −z| < R}∩ D = ∅ and fix z = 0. For every R D 1 the rectangle I ×(−2R, 2R) d−1 is such that Since by (87) This trivially yields (20) for ∇G D . The bound (20) for G D follows from (96) by an application of Poincaré's inequality, this time in the domain {|x − z| > R} ∩ D. 9 Proof of Corollary 2 Throughout this proof we assume d > 2 and recall that, for a ∈ ,we adopt the notation G(a, ·, ·) for the Green function for the whole space R d .
Step 1: · -almost sure solutions of (8) and (11). We show that with the additional structure · on , it holds that ∀ a.e. y ∈ R d and · -a.e. a ∈ , G(a; ·, y) solves (8), ∀ a.e. y ∈ R d and · -a.e. a ∈ , ∇ y G(a; ·, y) solves (11), and for every R > 0, z ∈ R d , almost every x, y ∈ R d and · -almost every a ∈ , In other words, we prove that the ensemble on , chosen to be such that the L 1 ( ) space is separable (cf. Sect. 2), allows to exchange in (8) and (11) as well as (87), (88) and (89) the order of the quantors a and x, y. This will be useful in the next steps, when we treat G(·; ·, y) and ∇ y G(·; ·, y) as almost sure solutions of respectively (8) and (11). From Theorem 1, we have that for every φ ∈ L 1 ( ) and ζ, ρ ∈ C ∞ 0 (R d ) it holds
We now appeal to the separability of C ∞ 0 (R d ) with respect to the C 1 topology, to conclude that for almost every y ∈ R d , · -almost every a ∈ and for every ζ ∈ C ∞ 0 (R d ),
In a similar way we obtain identities (99), (100) and (101): We show the argument only for (101) since the arguments for the other two are analogous. Identity (89) with a fixed R > 0 yields for any triple withˆ: → such thatâ(·) := a(R ·). By Fubini's theorem we may exchange the order of integration in the previous identity and obtain thatˆζ Therefore, separability of L 1 ( ) yields that for almost every x, y ∈ R d and · -almost every a ∈ , identity (101) holds.
Step 2: Spacially averaged annealed bounds. We argue that for almost every y ∈ R d and We claim that it is sufficient to prove (103),(104) and (105) for R = 1: Let us assume for instance that (103) holds for a R = 1, namely that ˆ1 <|x−y|<2 |G(a; x, y)| 2 1.
Since for almost every x, y ∈ R d , · -almost every a ∈ and every countable set of radii R identity (101) holds, we may infer from (106) that for almost every y ∈ R d , bound (103) is true for every R ∈ R. We now show that with an appropriate choice of R, we extend (103) to any R > 0: Picking The same reasoning holds for (104) and (105). Moreover, since the previous argument may be adapted to any fixed R 1, for convenience in the next estimates, we prove (103),(104) and (105) with R = 3.
We start with inequality (105): We claim that it is enough to prove that for almost every y ∈ R d and δ << 1, Indeed, using (18) we may send δ → 0 and conclude by Lebesgue Differentiation Theorem. We thus prove (107): We take the average · into inequality (18) with z = y, R = 2 and, after integrating in the x and y -variables, we obtain ˆ| and also ˆ| y −y|<1ˆ|x−y |>3 |∇∇G(a; x, y )| 2 1.
We now consider n ∼ δ −d disjoint balls of radius δ << 1 centred in {w i } n i=1 points and contained in the unitary ball centred at the origin: The previous inequality yields Moreover, thanks to (100) and stationarity, we rewrite the l.h.s. of the previous inequality as and, by the change of coordinates x = x − w i and y = y − w i , as Inserting this into the l.h.s. of (110) allows to conclude (107) and thus establish (105). The bound (104) for ∇G follows analogously from inequality (15 1. Since as we argue above we may assume that (101) holds for almost every x, y ∈ R d , ·almost every a ∈ and on a countable set of radii, we infer that from the above inequality we have also for every n ∈ N that ˆ2 so that summing over n ∈ N we conclude (104) also for ∇ y G. Inequality (103) follows from (104) for ∇G, again by Sobolev's and Hölder's inequality.
Step 3: Spatially pointwise estimates. We now post-process (103), (104) and (105) to obtain (21), (22) and (23). Reasoning as in Step 2, without loss of generality it suffices to prove (21), (22) and (23) Before proving (112), we show how to conclude the argument for (21), (22) and (23). We start with (21): By symmetry (99) and (97), for almost every x ∈ R d and · -almost every a ∈ the application u(a; ·) = G(a; x, ·) is a-harmonic in |x − y| > 2. Moreover, since we may select in this domain N 1 balls of radius 8, centred in {w i } N i=1 points such that their union covers the annulus 10 < |x − y| < 12, estimate (112) yields that for almost every y In addition, by the inclusion i.e. bound (21). In order to have also (22)-(23), we consider u(a; y) = ∇G(a; x, y) which, thanks to symmetry (99) and (98), for almost every x ∈ R d and · -almost every a ∈ is a-harmonic in {|y − x| > 2}. Therefore, reasoning as for bound (21), we may apply estimate (112) to this choice of u and get that for almost every y such that {10 < |x − y| < 12} We now argue that (112) Indeed, arguing again by separability of L 1 ( ), we also infer that the previous bounds hold for almost every y ∈ R d such that |y − w| < 1 and for · -almost every a ∈ . Therefore, we may take in (113)  .
By definition of η, v, f , and g together with conditions (4)

Fourier approach
Here we summarise how the Fourier method developed in [8] can be used to prove Corollary 2 provided the system is uniformly elliptic, so we shall assume that both (3) and (5) hold. The method is based on a representation of the Fourier transform of G in terms of a function Y ), which satisfies an elliptic PDE on .
To define the PDE for we introduce some notation. First observe that ξ ∈ R d can be regarded as being in the space L(Y, Y d ). In that case we denote its adjoint by ξ * ∈ L (Y d , Y ). Similarly the gradient operator D acts on functions F : → Y to yield a function D F : with e i denoting the standard ith-versor in R d . We denote by D * the corresponding divergence operator, which takes a function F : → Y d to a function D * F : → Y . Using this notation, the function is the solution to the equation where P is the projection operator on L 2 ( ) orthogonal to the constant. We can see using . To do this we apply the adjoint (a, ξ) * ∈ L(Y, Y d ) to (115) and take the expectation. This yields the inequality Next we define a function q : Then from (4), (5) it follows that q(ξ ) is Hermitian for ξ ∈ R d and From (118) we conclude that ξ * q(ξ )ξ ∈ L(Y, Y ) is invertible provided ξ = 0. Generalising the representation of [8] (see equation (2.4) of [8] or equation (8.1) of [9]) to the case of systems, we see that ∇ x G(a; x, y) is given by the Fourier inversion formula ∇ x G(a; x, y) Let H be a Hilbert space with norm · and consider functions f : The norm of f , which we denote by f w, p , is the smallest constant C for which (120) holds. The following lemma can be proved in the the same way as Lemma 3.5 of [8].
Comments on Lemma 5. We give below the main ideas to prove Lemma 5 as in [8] (Lemma 3.5). We use a scalar notation, but all the arguments and techniques hold also for systems. The proof mainly relies on the representation formula for [see Lemma 3.2 in [8]] where if ρ is a random variable ρ : → R d and b = I − a. This can be obtained from (115) by a standard perturbation argument applied to the operator (D * + iξ * )a(D − iξ). Being ||T b,ξ || ≤ ||b|| L ∞ ( ) < 1, (122) is well defined and its Neumann series converges. Moreover, derivatives with respect to ξ of (123) can be explicitly written as Note that in a rigorous formulation, to assure the convergence of the integrals in (123), one should first work with the massive Green function G T associated to the operator T −1 −∇ ·a∇ ( [18], Definition 2.4) and then pass to the limit T → +∞ and obtain G. To keep notation lean, we neglect this issue. For the same reason, we restrict our attention to the case d = 3. Due to (122) and (124), f (ξ ) in (121) is equal to a sum of terms containing derivatives as in (124). Hence, for d = 3 we get (|n| = 1 < 3 2 ) More precisely, each term on the r.h.s. of (125) may be rewritten an operator acting on Lemma 5 follows once it is proved that S 1 and S 2 are bounded from L p w (R 3 ) to L p w (R 3 , L 2 ( )) for every p ∈ (2, +∞). The most challenging operator is S 2 : The one associated to the second term on the r.h.s. of (125), where the derivative falls on (1 − P T b,ξ ) −1 .
To deal with it, it is convenient to first prove its boundedness from L p (R 3 ) to L p (R 3 , L 2 ( )) for p ∈ {2, +∞} and then use Hunt's interpolation theorem. The case p = +∞ follows by an application of Bochner's theorem (cf. [8], formula (3.14)), while for p = 2 the main idea relies on the fact that (1 − P T b,ξ ) −1 can be written in Neumann series and every term can be explicitly expressed. Once explicit, one can recognise that each term acts on a function by essentially taking multiple convolutions of its Fourier transform with the Hessian of the standard Green function G (I ; x, y). Such a convolution kernel does not increase the (Frobenius) norm of the function. A generalisation to higher dimensions is in the same spirit but has to deal with more involved operators S 1 , S 2 , . . . S N (for N = N (d)). The upper bound for the number of derivatives |n| is related to the strict condition p > 2 which ensures the boundedness of the operators between the L p -weak spaces (see Lemma 3.7 and Lemma 3.9 in [8]).
The following lemma implies (104) of Corollary 2 provided d is odd. In order to prove (104) when d is even we would need to extend Lemma 5 to include fractional derivatives, something that is also required in [8]. ≥ 3 and n = (n 1 . . . , n d )

Lemma 6 Let d
Proof We have from (119) on integration by parts that , it follows from (117), (118) and Lemma 5 that f : for some constant C d (λ) depending only on λ, d. Let φ be a cut-off function for {|x| < 1} in {|x| < 2}. Then from (129) we have that ˆR It follows from (130) that ˆR where g : We can estimate the RHS of (131) by using the inequalitŷ where the constant C q diverges as q → p. We consider for any A > 0 the integral Taking q = 1 in (132), we see that the first term on the RHS of (133) is bounded by where C is a constant depending only on d and g p,w . Taking p = d/|n|, we see that the sum in (134) converges provided |n| < d − 1. If this is the case then the first term on the RHS of (133) is bounded by C A 2(d−1−|n|) R 2(|n|+1)−d , where C depends only on λ, d. To estimate the second term in (133) we take q = 2 in (132). Thus it is bounded by where C is a constant depending only on d and g p,w . Taking p = d/|n| as before, we see that the sum in (135)

Appendix
Proof of Lemma 1 Let a ∈ and the domain D be fixed, and let us assume that G We first argue that for every y ∈ R d , By definition of u, it indeed holds that for almost every y ∈ R d and every ζ ∈ C ∞ 0 (D) D (x, y ) dx (8) = 0.
Multiplying this with ρ(y − y) and integrating in y , we thus have for every ζ ∈ C ∞ 0 (D) that i.e. for every y ∈ R d ,û(·, y) solves the Eq. (136). The boundary conditions are trivially satisfied since by definition of G (1) and G (2) we havê (8) = 0.
We now claim that for every y ∈ R d fixed it holdŝ We start by noting that sincê after smuggling into the r.h.s. the weight |x − y | α ∧ 1 1 2 and using Hölder's inequality in the y -variable combined with ρ ∈ C ∞ 0 ({|y| < 1}), this may be estimated bŷ For i ∈ {1, 2}, let α i < d be the exponent in the inequality (6) for G (i) . Without loss of generality, we assume that α 1 < α 2 < d so that for every z ∈ Z d we havê We now appeal to (6) and (7) to conclude (137). Thanks to (137), we may test the Eq. (136) withû(·, y) itself, and obtain by (4) and (8), that u(·, y) = 0 almost everywhere in R d . It thus follows by definition ofû that for every y ∈ R d and ζ ∈ C ∞ 0 (R d ) we havêˆρ Therefore, by the arbitrariness of y ∈ R d and of both ζ ∈ C ∞ 0 (R d ) and ρ ∈ C ∞ 0 ({|y| < 1}) we conclude that G (141) From this, it immediately follows that u and ∇u are bounded in {|z| < R}, provided u has finite energy in {|z| < 8R}.
To show (140) and (141), we use an argument similar to the one of Corollary 2, Step 3. Moreover, with a standard approximation argument we can assume that u and a are smooth. By scaling, without loss of generality we can reduce ourselves to consider the case R = 1. For a cut-off function η of {|x| < 2} in {|x| < 4}, we have 11 −∇ · a∇(ηu) = −∇ · a(∇ηu) − ∇η · a∇u − η∇ · a∇u.

∇u(x) =ˆ(u(z) − c)∇∇G(a; z, x) · a∇η(z)dz −ˆ∇ x G(a; z, x)∇η(z) · a∇u(z)dz,
for any constant c ∈ R. We now apply Hölder's inequality in the integrals on the r.h.s. of the previous inequalities and obtain by (4) and the choice of the test-function η and .
The two ingredients for (153) are the interpolation inequality for any function v of zero spatial average and even l ∈ N and the Caccioppoli estimate for m ≥ 2n. Before proving (155) and (156) we show how to obtain (153) from them: We insert the Caccioppoli estimate (156) with m = n 2 in the interpolation inequality (155) with 2 ) d u, η replaced by η n and l = n − 1. Appealing to Young's inequality we obtain 1 |k| 2l |Fu(k)| 2 and thus inequality (153) by another application of the Caccioppoli estimate (156) and the choice of the support of η.
To obtain the interpolation estimate (155), we rewrite it without Fourier transform, appealing to the representation of the Laplacian − N with Neumann boundary conditions through the Fourier cosine series by F (− N )w(k) = |k| 2 F w(k): For (155) it thus suffices to show that for an arbitrary function w It is easily seen that this family of interpolation estimates indexed by even l follows from the following two-tier family of interpolation inequalities indexed by m ∈ N and Indeed, plugging (159) in (158) yields and after an application of Young's inequality We apply once more Young's inequality to the first to last term on the r.h.s and get By iterating the previous estimates we conclude (157) from (158) and (159). Obviously, the two-tier family (158)-(159) reduces to the two estimates which by Young's inequality follow from Thanks to (154) and the choice of the support of η, these two last estimates immediately follow from integration by parts, the Cauchy-Schwarz and the triangle inequalities.
In proving (156) we also introduce the notation (·, ·) for the usual L 2 -inner product in (− π 2 , π 2 ) d . We first observe that by identity (147) in the proof of Remark 1 and properties (4) and (3), it holds for any test functionη in We now test (37) with η 2m (u − c); by the invariance of the equation under translations, we may fix without loss of generality c = 0. Thanks to the cut-off function η, we obtain and by (160) This last inequality implies (156) provided that we have that for every i = 1, . . . , d it holds To simplify the notation, we drop the index i. We claim that to obtain (162) it is sufficient to show that Indeed, the l.h.s of (162) can be estimated from below by and after an application of Cauchy-Schwarz's inequality followed by Young's inequality, by We now may plug (163) with α = l = i and k = n − i in the terms of the sum in (164) and apply Young's inequality to conclude (162). Inequality (163) is implied by the interpolation estimate for all l ≥ 1, combined with an iterated application of Young's inequality, similar to the one used to infer (157) from (158)-(159). We prove (165) by induction on l: For l = 1, we estimate and use Cauchy-Schwarz's inequality to get Young's inequality yields i.e. (165) with l = 1. Let us now assume that (165) holds for every i ≤ l and l > 1: First using (165) with i = 1, we have then, another application of (165) with i = l on ||η (α+1) ∂ 2 u|| yields ||η α ∂u|| 2 ||η (α−1) u|| 2 + ||η α ∂u||||η (α−1) u|| By Young's inequality we can absorb the term ||η α ∂u|| on the r.h.s and thus obtain (165) with l + 1 .

Proof of Lemma 3
Using translation and scaling invariance, up to a relabelling of the domain D, we may reduce ourselves to R 1. As in the proof of Lemma 2 it is convenient to work with boxes instead of balls and thus argue that for a function u satisfying (37) in the domain where the supremum runs over all linear functionals F bounded in the sense of with v satisfying (i) and (ii) in the sense Definition 1 in We observe that Lemma 2 also implies by Poincaré's inequality in (− π 4 , π 4 ) d that for all w satisfying (34) in where the linear functionals F satisfy (36) in (− π 2 , π 2 ) d . We want to argue that a similar estimate holds also for u solving (37) in the domain for all linear functionals satisfying (167). Before tackling (168) we show how to conclude the proof: Thanks to (168), if we define We now observe that since d > 2, we may apply first Hölder's inequality in (− π 4 , π 4 ) d \ (− π 6 , π 6 ) d and then Sobolev's inequality in the outer domain {|x| > π 6 } to infer that i.e. F 0 satisfies (167). Therefore, it follows from (169) thatˆ( whereη is any cut-off function for the set R d \ (− π 5 , π 5 ) d in R d \ (− π 6 , π 6 ) d satisfying (154) and m ∈ N , m > n. We remark that this last inequality is obtained in a similar way to (156), this time testing the Eq. (37) withη 2m u which is an admissible test function since u solves (37). This time, with a more careful look, we rewrite (161) as We now choose another functionη, satisfying (154) and such that it cuts off the set supp(∇η) in (− π 4 , π 4 ) d \ (− π 8 , π 8 ) d . We thus have for any α > 0 Therefore, we may bound ε ∂ n η 2m u , ∂ n u We now may apply the interpolation inequality (163) with η =ηη to the second and third term on the r.h.s. and appeal to (173) to estimate and conclude that ε ∂ n η 2m u , ∂ n u −ε||η (m−n) u1 (− π 4 , π 4 ) d \(− π 6 , π 6 ) d ||, so that the Caccioppoli's inequality (171) follows by the previous inequality and (172).
It thus remains only to prove inequality (168): We fix aη to be a cut-off function of R d \ (− π 6 , π 6 ) d in R d \ (− π 8 , π 8 ) d and set In the same spirit of Lemma 2, we define for k ∈ Z d − {0} with F (ũ)(k) the k-th coefficient of the cosine Fourier series of the functionũ defined in (152). We first note that the functionalF(u)(k) satisfies (167) for all k ∈ Z d \ {0}: Indeed, similarly to the proof of Lemma 2 we may write so that first Cauchy-Schwarz's inequality, then Poincaré's inequality and the definition ofũ yield Another application of Poincaré's inequality on the second term of the r.h.s implies (167). Therefore, analogously to Lemma 2, if we show that for any even l ∈ N we have inf c∈Rˆ( − π 4 , π we may conclude (168)  domain D = I × R d−1 with I ⊂ R bounded. Therefore since in this setting the implicit multiplicity constant in (92) depends on D through the size of I , we may substitute the notation D with I . We start by observing that it is sufficient to prove that there exists a constant C 0 such that for every i = 1, . . . , d − 1 where x = (x 1 , . . . , x d−1 ) and x 0 = (x 0,1 , . . . , x 0,d−1 ). Indeed, from (180)-(181) it follows thatˆd and thus also (92) thanks to the convexity of the function exp( |s| C 0 ) and the equivalence of the norms in R d−1 .
Without loss of generality, we fix x 0 = 0 in (180)-(181). We start with the argument for (180) for a fixed i, say i = 1. As D is unbounded, we are not allowed a priori to test the equation before (92) with exp ± where g M (x 1 , x ) = X (−(M−1),+∞) d−1 (x )g(x 1 , x ). We start by showing that if we prove (180) for every u M , namelŷ  Thus, {u M } M∈N is uniformly bounded in W 1,2 (R d ) and, up to subsequences, weakly converges to u, the solution of the equation before (92). In addition, since by construction D 1 ⊂ · · · ⊂ D N ⊂ D N +1 ⊂ · · · ⊂ D, from (183) it holds that for every N , If we send M → +∞ we obtain by weak lower semi-continuitŷ Taking also the limit N → +∞, we conclude (180).
By our choice of ρ, we thus obtain . Therefore, after an application of Young's inequality in the last two terms on the r.h.s., we may choose C 0 > C(λ, |I |), absorb the term´ρ|∇u| 2 and thus establish (183). Inequality (181) follows similarly by considering an approximate problem in D × (−∞; M) d−1 and ρ = exp x C 0 .