# Boundary Layer Estimates in Stochastic Homogenization

Peter Bella  TU Dortmund, DE-44227, Dortmund, Germany  Julian Fischer  IST
Austria, AT-3400, Klosterneuburg, Austria  Marc Josien  CEA, DES, IRESNE, DEC,
Cadarache, F-13108, Saint-Paul-Lez-Durance, France  Claudia Raithel  TU
Dresden, DE-01217, Dresden, Germany

###### Abstract

We prove quantitative decay estimates for the boundary layer corrector in
stochastic homogenization in the case of a half-space boundary. Our estimates
are of optimal order and show that the gradient of the boundary layer
corrector features nearly fluctuation-order decay; its expected value decays
even one order faster. As a corollary, we deduce estimates on the accuracy of
the representative volume element method for the computation of effective
coefficients: our understanding of the decay of boundary layers enables us to
improve the order of convergence of the RVE method for $d\geq 3$. ††This
project has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No 948819).

###### Contents

  1. 1 Introduction
  2. 2 Discussion of main results and strategy
    1. 2.1 Setting of our results
    2. 2.2 Theorem 1: Optimal decay estimate for the boundary layer
    3. 2.3 Theorem 2: Improved error estimate for representative volume element method.
  3. 3 Overview of our arguments
    1. 3.1 Strategy for Theorem 1
    2. 3.2 Argument for Proposition 2: Iterative improvement of decay estimate for the boundary layer
    3. 3.3 Proof of Theorem 2
  4. 4 Argument for Proposition 1: A suboptimal decay estimate for the boundary layer
    1. 4.1 Proof of Proposition 1
    2. 4.2 Proof of auxiliary lemmas for Proposition 1
  5. 5 Argument for Proposition 5: Fluctuation estimates for the boundary layer corrector
    1. 5.1 Identification of the Malliavin / Fréchet derivative
    2. 5.2 Estimates on the auxiliary functions: A weighted Meyers estimate
    3. 5.3 Proof of Proposition 5
    4. 5.4 Proofs of auxiliary lemmas for Proposition 5
  6. 6 Argument for Proposition 6: Bound for expectation of weighted average of the boundary layer
    1. 6.1 Proof of Proposition 6
    2. 6.2 Proof of auxiliary lemma for Proposition 6
  7. 7 Proof of Theorem 1
  8. A Regularity of random elliptic operators on $\mathbb{R}^{d}$ and corrector estimates
  9. B Regularity estimates for random elliptic operators on a domain
  10. C Proof of the weighted Meyers estimate in Lemma 14
  11. D Exponentially localized boundary layer

##  1 Introduction

Both in periodic homogenization and in stochastic homogenization, the presence
of a domain boundary should be expected to lead to intricate effects and
substantial new challenges: In the periodic setting, the domain boundary
breaks the periodic structure, while in stochastic homogenization, the domain
boundary breaks the stationarity (invariance with respect to spatial shifts)
of the problem.

Indeed, in particular in quantitative homogenization theory on bounded
domains, it is important to understand the impact of boundary layers on the
solution: For instance, for periodic coefficient fields
$a(\cdot)=a(\cdot+k),k\in\mathbb{Z}^{d}$, with $a=a^{*}$, solutions to the
linear elliptic PDE on the whole space

$\displaystyle-\nabla\cdot\big{(}a(\tfrac{\cdot}{\varepsilon})\nabla
u_{\varepsilon}\big{)}=\nabla\cdot f\quad\quad\text{on }\mathbb{R}^{d}$ (1.1)

may be approximated by the solution to a constant-coefficient effective
equation $-\nabla\cdot(\bar{a}\nabla\bar{u})=\nabla\cdot f$ up to an error of
second order111Provided that the right-hand side $f$ is sufficiently regular.

$\displaystyle||u_{\varepsilon}-\bar{u}||_{H^{-1}(B_{1})}\lesssim\varepsilon^{2}.$
(1.2)

In contrast, the analogous estimate _fails_ for the Dirichlet problem on
bounded domains, even for the homogeneous Dirichlet problem and even strictly
in the interior of the domain: On bounded domains the estimate
$||u_{\varepsilon}-\bar{u}||_{L^{p}}\lesssim\varepsilon$ cannot be improved by
passing to weak norms. The underlying reason for this failure is boundary
layer effects that provide an $\varepsilon$-order contribution to the
_effective boundary data_ for $\bar{u}$.

On a mathematical level, the failure of the derivation of estimates like (1.2)
in the case of the Dirichlet problem on domains may be seen as follows: In
order to prove an estimate like (1.2), one attempts to approximate
$u_{\varepsilon}$ by the so-called two-scale expansion; written up to second
order, it reads for periodic and symmetric $a=a^{*}$

$\displaystyle\bar{u}(x)+\varepsilon\phi_{i}(\tfrac{x}{\varepsilon})\partial_{i}\bar{u}(x)+\varepsilon^{2}\psi_{ij}(\tfrac{x}{\varepsilon})\partial_{i}\partial_{j}\bar{u}(x),$
(1.3)

where $\phi_{i}$ and $\psi_{ij}$ are the so-called first-order and second-
order homogenization correctors. The homogenization corrector $\phi_{i}$ is
defined as the $1$-periodic solution to the PDE
$-\nabla\cdot(a(e_{i}+\nabla\phi_{i}))=0$; note that it is unique up to
additive constants. If we attempt to apply the two-scale expansion ansatz
(1.3) on a bounded domain $\Omega\subset\mathbb{R}^{d}$, we immediately
observe that it fails to attain the correct boundary data $\bar{u}$ on the
domain boundary. While the second-order term $O(\varepsilon^{2})$ does not
affect the accuracy, the first-order contribution
$\varepsilon\phi_{i}(\tfrac{x}{\varepsilon})\partial_{i}\bar{u}(x)$ must be
corrected for. In other words, we would need to subtract a term of the form
$\varepsilon v(x)$ with $v$ solving

$\displaystyle-\nabla\cdot\big{(}a\big{(}\tfrac{\cdot}{\varepsilon}\big{)}\nabla
v\big{)}$ $\displaystyle=0$ $\displaystyle\text{in }\Omega,$ $\displaystyle
v(x)$ $\displaystyle=\phi_{i}(\tfrac{x}{\varepsilon})\partial_{i}\bar{u}(x)$
$\displaystyle\text{on }\partial\Omega.$

This is an example of the so-called _oscillating Dirichlet problem_ : The
boundary data oscillates on the same scale as the coefficient field. As a
consequence, the interactions of the oscillations in the boundary data and in
the coefficient field give rise to quite nontrivial nonlinear interactions; it
turns out that under suitable assumptions, $v$ may be described in the
interior in terms of a suitable $\bar{a}$-harmonic function with certain
_effective boundary conditions_ on $\partial\Omega$. Note that in particular,
these effective boundary conditions are not given in terms of a simple average
of the boundary data.

In the context of regularity theory for linear elliptic operators with
oscillating coefficients, a particularly important instance of an oscillating
Dirichlet problem is given by the _boundary layer corrector_
$\theta_{i}^{\varepsilon}$, defined as the solution to the problem

$\displaystyle-\nabla\cdot\big{(}a\big{(}\tfrac{\cdot}{\varepsilon}\big{)}\nabla\theta_{i}^{\varepsilon}\big{)}$
$\displaystyle=0$ $\displaystyle\text{in }\Omega,$
$\displaystyle\theta_{i}^{\varepsilon}(x)$
$\displaystyle=\phi_{i}(\tfrac{x}{\varepsilon})$ $\displaystyle\text{on
}\partial\Omega.$

By its definition, subtracting the boundary layer corrector
$\theta_{i}^{\varepsilon}$ from the whole space corrector
$\phi_{i}^{\varepsilon}(x):=\phi_{i}(\tfrac{x}{\varepsilon})$ corrects the
boundary conditions of $\phi_{i}^{\varepsilon}$ while maintaining the defining
corrector equation
$-\nabla\cdot(a(e_{i}+\nabla(\phi_{i}^{\varepsilon}-\theta_{i}^{\varepsilon})))=0$
– in other words, $\phi_{i}^{\varepsilon}-\theta_{i}^{\varepsilon}$ yields a
first-order homogenization corrector in $\Omega$ with homogeneous Dirichlet
boundary conditions on $\partial\Omega$.

In the case of periodic homogenization, the study of the oscillating Dirichlet
problem has a long history and by now a relatively complete picture has
emerged in the literature. On half-spaces with rationally aligned normals
(with respect to the underlying periodic structure), Allaire and Amar [1] have
shown exponential decay of the gradient of the boundary layer corrector away
from the boundary; they subsequently used this result to homogenize the
oscillating Dirichlet problem on polygonal domains with rationally aligned
boundaries 222However, this homogenization result only holds along
subsequences.. For half-spaces with Diophantine normals, Gérard-Varet and
Masmoudi [14, 15] established superpolynomial decay of the boundary layer
corrector gradients; as a consequence, they were able to formulate
quantitative homogenization error estimates for the oscillating boundary data
problem on uniformly convex domains, albeit with suboptimal rate. Recently,
Armstrong, Kuusi, Mourrat, and Prange [2] have obtained optimal convergence
estimates in $L^{p}$ for $p\geq 2$ for the oscillating Dirichlet problem on
convex domains, employing also a subsequent regularity result for the
effective boundary data by Shen and Zhuge [28]. The case of the Neumann
problem has recently been treated by Shen and Zhuge [27].

In principle, similar issues arise in the homogenization of random media: For
stationary and suitably decorrelating random coefficient fields $a$ with
$a=a^{*}$, an analogue of the higher-order error estimate (1.2) has been
derived in [6]. More precisely, in the case of such random coefficient fields,
solutions to the linear elliptic PDE (1.1) on $\mathbb{R}^{d}$ with smooth
right-hand side $f$ may be approximated by a constant-coefficient effective
PDE $-\nabla\cdot(\bar{a}\nabla\bar{u})=\nabla\cdot f$ up to an error of
order333Note that for $d=3$ and $d=4$ this error estimate is of fluctuation
order and therefore optimal.

$\displaystyle||u_{\varepsilon}-\bar{u}||_{H^{-1}(B_{1})}\leq\mathcal{C}\begin{cases}\varepsilon^{3/2}&\text{
for }d=3,\\\ \varepsilon^{2}|\log\varepsilon|&\text{ for }d=4,\\\
\varepsilon^{2}&\text{ for }d\geq 5.\end{cases}$

Again, in the case of the Dirichlet problem on bounded domains no analogue of
this higher-order convergence result holds and the best possible
homogenization error estimate for the Dirichlet problem is
$||u_{\varepsilon}-\bar{u}||_{L^{p}}\leq\mathcal{C}\varepsilon$ (for $d\geq
3$).

While the periodic homogenization problem with fluctuating boundary data is by
now well-studied, in the case of stochastic homogenization only few results on
the boundary layer corrector have been available in the literature. For
instance, in [13] boundary layer correctors have been constructed for half-
spaces to develop a large-scale regularity theory; however, the available
estimates on the boundary layer corrector are far from optimal for $d\geq 3$.
In [23], correctors for harmonic functions around lower-dimensional features
like cusps have been studied. For correctors for interface problems, we refer
to [22]. For interface problems in periodic homogenization we refer to [9, 21,
29].

Quantitative estimates on boundary layers in stochastic homogenization. In the
present work, we establish quantitative estimates on the decay of the boundary
layer corrector $\theta_{i}$ on a half-space $\mathbb{H}^{d}_{+}$, given as
the solution to the PDE

$\displaystyle-\nabla\cdot(a\nabla\theta_{i})$ $\displaystyle=0$
$\displaystyle\text{in }\quad\,\mathbb{H}^{d}_{+},$ (1.4)
$\displaystyle\theta_{i}$ $\displaystyle=\phi_{i}$ $\displaystyle\text{on
}\quad\partial\mathbb{H}^{d}_{+}.$

More precisely, we consider an ensemble $\langle\cdot\rangle$ of coefficient
fields $a$ on $\mathbb{R}^{d}$ for $d\geq 3$ (that is, a probability measure
on the space of coefficient fields on $\mathbb{R}^{d}$) that is subject to a
standard set of conditions of quantitative stochastic homogenization: We
require uniform ellipticity and boundedness, stationarity (invariance of the
law under spatial shifts), and quantitative decorrelation on scales larger
than a microscale $\varepsilon>0$ in form of a spectral gap inequality. Under
an additional small-scale regularity condition, we show that the boundary
layer corrector decays as

$\displaystyle\big{|}\nabla\theta_{i}(x)\big{|}\leq\mathcal{C}\bigg{(}\frac{\varepsilon}{\operatorname{dist}(x,\partial\mathbb{H}^{d}_{+})}\bigg{)}^{d/2-\delta}$
(1.5)

with a random prefactor $\mathcal{C}$ subject to uniform (stretched
exponential) moment bounds (and with any $\delta>0$). One should expect this
estimate to be of optimal order (up to the loss of $\delta$), as the average
of correctors on the boundary $\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{\partial\mathbb{H}^{d}_{+}\cap B_{r}}\phi_{i}\,dS$
displays fluctuations of order $r^{1-d/2}$ and $\theta_{i}$ is the
$a$-harmonic extension of $\phi_{i}|_{\partial\mathbb{H}^{d}_{+}}$ to
$\mathbb{H}^{d}_{+}$. For the expected value $\mathbb{E}[\nabla\theta_{i}]$,
our arguments morally establish a higher-order decay

$\displaystyle\big{|}\mathbb{E}[\nabla\theta_{i}(x)]\big{|}\leq
C\bigg{(}\frac{\varepsilon}{\operatorname{dist}(x,\partial\mathbb{H}^{d}_{+})}\bigg{)}^{d/2+1-\delta},$

although we make this rigorous only in a spatially averaged sense (see
Proposition 6). To the best of our knowledge, the estimate (1.5) is the first
optimal-order estimate on the decay of boundary layers in stochastic
homogenization.

Application to the representative volume approximation. An interesting
application of these estimates for boundary layers arises in the analysis of
the representative volume element (RVE) method for the approximation of the
effective coefficient $\bar{a}$. The RVE method proceeds by taking a sample of
the random coefficient field – say, on a box $[0,L]^{d}$ with side length
$L\gg\varepsilon$ – and approximates the effective coefficient $\bar{a}$ as
$\bar{a}\approx a^{\operatorname{RVE}}$ with

$\displaystyle
a^{\operatorname{RVE}}e_{i}:=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[0,L]^{d}}a(e_{i}+\nabla\phi_{i}^{L})\,\mathrm{d}x.$

Here, $\phi_{i}^{L}$ solves the corrector equation
$-\nabla\cdot(a(e_{i}+\nabla\phi_{i}^{L}))=0$ on the representative volume
$[0,L]^{d}$, for instance with homogeneous Dirichlet boundary conditions. It
turns out that this RVE approximation with Dirichlet boundary data features
random fluctuations of order
$|a^{\operatorname{RVE}}-\mathbb{E}[a^{\operatorname{RVE}}]|\sim(\tfrac{\varepsilon}{L})^{d/2}$
and a systematic error of order
$|\mathbb{E}[a^{\operatorname{RVE}}]-\bar{a}|\sim\tfrac{\varepsilon}{L}$. The
fluctuations are of central limit theorem scaling and therefore of optimal
order. Note, however, that for $d\geq 3$ the systematic error dominates. It
turns out that the reason for this rather large systematic error contribution
is in fact a boundary layer phenomenon. Naively one might expect that by
averaging only in the interior (“oversampling”), i. e. by setting

$\displaystyle{\tilde{a}}^{\operatorname{RVE}}e_{i}:=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i}^{L})\,\mathrm{d}x$ (1.6)

for some constant $0<\kappa\leq\tfrac{1}{4}$, one may eliminate this issue of
boundary layers. Interestingly, this turns out _not_ to be the case: while
this procedure of oversampling decreases the systematic error by a constant
factor, the order of convergence remains the same. The underlying reason is
that the Dirichlet corrector $\phi_{i}^{L}$ differs from the whole space
corrector $\phi_{i}$ (the latter of which defines the effective coefficient as
$\bar{a}e_{i}:=\mathbb{E}[a(e_{i}+\nabla\phi_{i})]$) in terms of a boundary
contribution: Indeed, $\phi_{i}-\phi_{i}^{L}$ is an $a$-harmonic function in
$[0,L]^{d}$ with boundary data $\phi_{i}$, i. e. we need to understand a
homogenization problem with fluctuating boundary data. The failure of higher-
order convergence of (1.6) is due to _effective boundary conditions_ imposed
on $\phi_{i}-\phi_{i}^{L}$ due to the interaction of the _fluctuating boundary
data_ $\phi_{i}$ with the fluctuations of the coefficient field $a$. In
Theorem 2, we show that the slight modification of the RVE approximation given
by

$\displaystyle
a^{{\operatorname{RVE}},new}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{i}+\nabla\phi_{i}^{L}\,\mathrm{d}x:=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i}^{L})\,\mathrm{d}x$ (1.7)

with overwhelming probability improves the overall error estimate to
$(\tfrac{\varepsilon}{L})^{d/2-\delta}$ for $d=3$ and $d=4$. The underlying
reason is that the new formula takes into account the slight change of the
field gradient in the interior that is caused by the effective boundary
conditions of $\phi_{i}-\phi_{i}^{L}$. The new formula (1.7) is also rather
natural, as the effective coefficient should provide the relation between a
given macroscopic field gradient and the corresponding macroscopically
averaged flux.

Notation. Throughout the paper, to simplify notation and without loss of
generality we mostly consider the half-space of the form
$\mathbb{H}^{d}_{+}:=\\{x\in\mathbb{R}^{d}:x_{1}>0\\}$. Thus, the natural
(interior) unit normal vector is $e_{1}:=(1,0,\ldots,0)$. For points
$x\in\mathbb{R}^{d}$ we decompose their coordinate vector $x$ as
$x:=(x^{\perp},x^{\parallel})$ with $x^{\perp}\in\mathbb{R}$ and
$x^{\parallel}\in\mathbb{R}^{d-1}$. Occasionally, we denote a fixed point as
$x_{0}\in\mathbb{H}^{d}_{+}$ –this should not be confused with the notation
for the coordinates $x_{i}:=x\cdot e_{i}$ for $i=1,\ldots,d$.

We frequently use the convention that $``\lesssim"$ means $``\leq
C(d,\lambda),"$ where $d$ denotes the dimension and where $\lambda>0$ is the
ellipticity ratio of the coefficient field $a$.

By $c$ and $C$, we denote generic constants whose value may change from
appearance to appearance; typically, $C$ is used to denote a large constant,
while $c$ denotes a small constant.

Similarly, we use the notation $\mathcal{C}$ to denote a generic random
constant (respectively, $\mathcal{C}(x)$ to denote a generic random field)
whose value may change from appearance to appearance. Unless otherwise
specified, the constant $\mathcal{C}$ (respectively the random field
$\mathcal{C}(x)$) will be subject to a uniform stretched exponential moment
bound.

By $B_{r}(x)$ we denote the ball of radius $r$ centered at $x$, while by
$\mathcal{O}$ we denote a (typically $C^{1,1}$, but not necessarily bounded)
domain.

Throughout the paper, we consider an ensemble $\langle\cdot\rangle$ (a
probability distribution) of elliptic and bounded coefficient fields $a$; we
denote the corresponding probability space as $\Omega$ and the associated
expected value of a random variable $G$ by $\langle G\rangle$. Regarding
terminology, a _quenched_ estimate is a bound that holds almost surely, while
an _annealed_ estimate is a bound on suitable stochastic moments.

We also use standard notation for Sobolev spaces; in particular, by
$\dot{H}^{1}(\mathbb{H}^{d}_{+})$ we denote the space of decaying functions
with finite Dirichlet energy equipped with the norm
$||v||_{\dot{H}^{1}(\mathbb{H}^{d}_{+})}^{2}:=\int_{\mathbb{H}^{d}_{+}}|\nabla
v|^{2}\,\mathrm{d}x$.

##  2 Discussion of main results and strategy

###  2.1 Setting of our results

To begin this section, we first list the assumptions on $\langle\cdot\rangle$,
the ensemble of coefficient fields $a$ on $\mathbb{R}^{d}$, under which we
obtain our results. Towards further illustrating our setting, we then give a
few examples of ensembles on $\mathbb{R}^{d}$ that satisfy these assumptions.
Having clarified the setting of our results we then discuss, in turn, Theorem
1 on the decay of the boundary layer and Theorem 2 on the improved convergence
rate for the representative volume approximation of effective coefficients.

#### Assumptions on the ensemble.

Throughout this article we will assume that:

  * (A1)

The ensemble $\langle\cdot\rangle$ is _uniformly elliptic and bounded_ : There
exists some $\lambda>0$ such that $\langle\cdot\rangle$-almost surely we have

$\displaystyle|a(x)v|$ $\displaystyle\leq|v|,$ $\displaystyle a(x)v\cdot v$
$\displaystyle\geq\lambda|v|^{2}$

for all $v\in\mathbb{R}^{d}$ and a. e. $x\in\mathbb{R}^{d}$.

  * (A2)

The ensemble $\langle\cdot\rangle$ is _stationary_ , that is, the law of
$a(\cdot+x)$ coincides with the law of $a(\cdot)$ for any
$x\in\mathbb{R}^{d}$. In other words, the statistics of the coefficient field
do not depend on the region in space.

  * (A3)

The ensemble is subject to decorrelation on scales larger than the microscale
in the sense that it satisfies a spectral gap inequality: For any random
variable $\xi(a)$ it holds that

$\displaystyle\left<(\xi-\left<\xi\right>)^{2}\right>\leq\varepsilon^{d}\left<\int_{\mathbb{R}^{d}}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\varepsilon}(x)}\big{|}\frac{\partial\xi}{\partial
a}\big{|}{\textrm{d}}y\bigg{)}^{2}\,\mathrm{d}x\right>,$ (2.1)

where $\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\varepsilon}(x)}\big{|}\frac{\partial\xi}{\partial
a}\big{|}{\textrm{d}}y$ is used to denote

$\displaystyle\sup_{\delta a}\limsup_{h\rightarrow 0}\frac{\xi(a+h\delta
a)-\xi(a)}{h}$

with the supremum taken over perturbations $\delta
a:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ that are supported in
$B_{\varepsilon}(x)$ and satisfy $\|\delta
a\|_{L^{\infty}(\mathbb{R}^{d})}\leq 1$. The parameter $\varepsilon>0$ is the
correlation length of the ensemble.

  * (A4)

The coefficient field $a$ is Hölder continuous with stretched exponential
moments for some Hölder exponent $\alpha\in(0,1]$. In particular, there exists
a stationary random field $\mathcal{C}(a,x)$ with
$\mathbb{E}[\exp(\nu\mathcal{C}^{\nu})]\leq 2$ for some $\nu>0$ such that the
Hölder estimate

$\displaystyle\sup_{y,z\in
B_{\varepsilon}(x)}\frac{|a(y)-a(z)|}{|y-z|^{\alpha}}\leq\mathcal{C}(a,x)\varepsilon^{-\alpha}$
(2.2)

holds.

Note that (A1)–(A3) are a standard set of assumptions in quantitative
stochastic homogenization; in particular, assumption (A3) is a standard
quantification of the qualitative ergodicity assumption that is used in
classical qualitative stochastic homogenization results.

Assumption (A4) is a small-scale regularity assumption that we require mostly
for convenience; it enables us to state pointwise gradient bounds (which may
otherwise fail due to a possible failure of regularity on scales smaller than
the microscale).

Examples of ensembles satisfying (A1) - (A4). We give two examples of random
coefficient fields, both of which are standard examples used to illustrate the
applicability of quantitative stochastic homogenization results.

Example 1: We consider a Gaussian ensemble that satisfies a certain
decorrelation estimate. In particular, let $\tilde{a}(x)$ be a scalar Gaussian
field that is centered and stationary and, furthermore, satisfies a
decorrelation estimate of the form

$\displaystyle\int_{\mathbb{R}^{d}}\sup_{|x|=|y|}\big{|}\left<\tilde{a}(x)\tilde{a}(0)\right>\big{|}\,{\textrm{d}}y\leq
1.$

Defining $a:=g(\tilde{a}(x))$, where $g$ is a bounded Lipschitz map from
$\mathbb{R}$ into the uniformly elliptic matrices $\mathbb{R}^{d\times d}$, it
is shown that (A1)–(A3) are satisfied e. g. in [11]. For sufficiently regular
covariances $\langle\tilde{a}(x)\tilde{a}(0)\rangle$, the ensemble also
satisfies (A4).

Example 2: The second example involves random spherical inclusions. More
precisely, we consider a Poisson point process on
$\mathbb{R}^{d}\times[0,\infty)$. Let $(X_{i},h_{i})_{i\in\mathbb{N}}$ be a
realization; then consider the points $X_{i}\in\mathbb{R}^{d}$ in the order of
their coordinate $h_{i}$ and for each point $X_{i}$ place a ball of radius $1$
around $X_{i}$ if the ball does not intersect with any of the previously
placed balls. Denote the union of the balls obtained by this procedure by $U$.
Then define a coefficient field as $a(x):=\lambda\chi_{\mathbb{R}^{d}\setminus
U}+\chi_{U}$. While such ensembles of coefficient fields consisting of random
high-conductivity inclusions do not satisfy the spectral gap inequality (A3),
they have been shown to satisfy a slightly different version of a spectral gap
inequality in [11]. We expect that a relatively straightforward adaption of
our arguments would make our results applicable also to ensembles satisfying
such an “oscillation spectral gap” estimate, though it would involve some
additional work to obtain the necessary small-scale (microscopic) regularity
theory (as derived from assumption (A4) below in our present work).

###  2.2 Theorem 1: Optimal decay estimate for the boundary layer

In our first main result we show that if the ensemble $\langle\cdot\rangle$ of
coefficient fields on $\mathbb{R}^{d}$ satisfies (A1)-(A4), then we are able
to obtain an almost-optimal decay estimate for the gradient of the boundary
layer corrector $|\nabla\theta_{i}|$ away from the boundary
$\partial\mathbb{H}^{d}_{+}$. We remark that this estimate is quenched, in the
sense that it is stated for a. e. coefficient field $a$ and includes a random
constant. Here is our first main result:

######  Theorem 1.

Let $d\geq 3$ and $\langle\cdot\rangle$ be an ensemble of coefficient fields
on $\mathbb{R}^{d}$ that satisfies the assumptions (A1)-(A4). Denote by
$\theta_{i}\in H^{1}_{\rm{loc}}(\mathbb{H}^{d}_{+};\mathbb{R}^{d})$ the half-
space corrector, that is, the unique444Uniqueness follows from the Liouville
principle for random elliptic operators on half-spaces proven in [13]. weak
solution to (1.4) with the property that $\nabla\theta_{i}$ decays away from
$\partial\mathbb{H}^{d}_{+}$ –in the sense that

$\displaystyle\lim_{r\rightarrow\infty}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}\cap\mathbb{H}^{d}_{+}}|\nabla\theta_{i}|^{2}{\textrm{d}}x=0.$

Then there exists a random field $\mathcal{C}(a,x)$, defined on
$\Omega\times\mathbb{H}^{d}_{+}$, that is stationary with respect to
tangential shifts and has stretched exponential moments –in the sense that
there exists $C=C(d,\lambda,\nu,\alpha,\delta)$ with

$\displaystyle\left<\exp\bigg{(}\frac{\mathcal{C}(a,x)^{1/C}}{C}\bigg{)}\right>\leq
2,$ (2.3)

for any $x\in\mathbb{H}^{d}_{+}$– such that the gradient of the boundary layer
corrector satisfies the decay estimate

$\displaystyle|\nabla\theta_{i}(x)|\leq\mathcal{C}(a,x)\bigg{(}1+\frac{x^{\perp}}{\varepsilon}\bigg{)}^{-\frac{d}{2}+\delta},$
(2.4)

for any $x\in\mathbb{H}^{d}_{+}$ and $0<\delta\ll 1$.

We remark that the decay of $|\nabla\theta_{i}|$ away from
$\partial\mathbb{H}^{d}_{+}$ that is shown in (2.4) is, up to the
$\delta$-loss, the optimal scaling that may be expected. This can be seen in
the following heuristic argument: Now fixing $x_{0}\in\mathbb{H}^{d}_{+}$, let
$h$ solve the divergence form equation

$\displaystyle-\nabla\cdot(a^{*}\nabla h)$
$\displaystyle=\nabla\cdot(\delta_{x_{0}})$ $\displaystyle\text{in
}\mathbb{H}^{d}_{+},$ (2.5) $\displaystyle h$ $\displaystyle=0$
$\displaystyle\text{on }\partial\mathbb{H}^{d}_{+},$

and rewrite $\nabla\theta_{i}(x_{0})$ as

$\displaystyle\nabla\theta_{i}(x_{0})=\int_{\mathbb{H}^{d}_{+}}-\nabla\cdot(\delta_{x_{0}})\theta_{i}\,\,\mathrm{d}x=\int_{\partial\mathbb{H}^{d}_{+}}\phi_{i}{\mathbf{n}}\cdot
a^{*}\nabla h\,{\textrm{d}}S,$ (2.6)

where ${\mathbf{n}}$ is the outward unit normal vector on
$\partial\mathbb{H}^{d}_{+}$. Notice that to obtain the above relation we have
used the equation (1.4), that $h=0$ on $\partial\mathbb{H}^{d}_{+}$, and the
identity $\theta_{i}=\phi_{i}$ on $\partial\mathbb{H}^{d}_{+}$. Now, we expect
the right-hand side of (2.6) to have the same asymptotic scaling behaviour as
the corresponding quantity with the Green’s function $h$ replaced by its
homogenized counterpart $\bar{h}$ solving
$-\nabla\cdot(\bar{a}^{*}\nabla\bar{h})=\nabla\cdot(\delta_{x_{0}})$ in
$\mathbb{H}^{d}_{+}$. Using that $|\nabla\bar{h}(\cdot)|\leq C|\cdot-
x_{0}|^{-d}$, we obtain

$\displaystyle\int_{\partial\mathbb{H}^{d}_{+}}\phi_{i}{\mathbf{n}}\cdot\bar{a}\nabla\bar{h}\,{\textrm{d}}S\sim(x_{0}^{\perp})^{-1}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{x_{0}^{\perp}}(x_{0}^{\parallel})}\phi_{i}\,{\textrm{d}}S,$
(2.7)

where we have used $B_{x_{0}^{\perp}}(x_{0}^{\parallel})$ to denote the $d-1$
dimensional ball with radius $x_{0}^{\perp}$ centered around
$x_{0}^{\parallel}$ in $\partial\mathbb{H}^{d}_{+}$. Since under the
assumptions $(A1)-(A3)$ volume averages of the gradient of the whole space
corrector display central limit theorem scaling [19], one expects that

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{x_{0}^{\perp}}(x_{0}^{\parallel})}\phi\,{\textrm{d}}S\sim|x_{0}^{\perp}|\bigg{(}\frac{|x_{0}^{\perp}|}{\varepsilon}\bigg{)}^{-\frac{d}{2}}.$

Combining this observation with (2.6) and (2.7) suggests that (2.4) cannot be
improved –aside from the $\delta$-loss.

###  2.3 Theorem 2: Improved error estimate for representative volume element
method.

In the second main result of our paper we use Theorem 1 to obtain an improved
error estimate for the representative volume element method, which is the
standard method used to approximate the homogenized coefficients $\bar{a}$.
The error in the representative volume element method can be split into two
contributions, the random error and the systematic error. While the former is
well-understood to be of central limit theorem scaling, without additional
modifications such as “screening” terms in the corrector equation [8, 17, 20]
the systematic error would dominate (for $d\geq 3$). As discussed above, our
understanding of boundary layers enables us to propose a minor modification of
the RVE method that improves the bound on the systematic error, thereby
improving the overall rate of convergence for $d\geq 3$.

######  Theorem 2.

Let $d\geq 3$ and let $\langle\cdot\rangle$ be an ensemble of coefficient
fields on $\mathbb{R}^{d}$ that satisfies the assumptions (A1)-(A4) with
$\varepsilon=1$. Let $\bar{a}\in\mathbb{R}^{d\times d}$ be the corresponding
effective coefficient. Denote by $\phi_{i}^{L}$ the homogenization corrector
with homogeneous Dirichlet boundary conditions on the box $[0,L]^{d}$, i. e.
let $\phi_{i}^{L}$ satisfy

$\displaystyle-\nabla\cdot(a(e_{i}+\nabla\phi_{i}^{L}))$ $\displaystyle=0$
$\displaystyle\text{in }(0,L)^{d},$ $\displaystyle\phi_{i}^{L}$
$\displaystyle=0$ $\displaystyle\text{on }\partial[0,L]^{d}.$

Given any $\kappa\in(0,\tfrac{1}{4}]$, define the RVE approximation
$\bar{a}^{L}\in\mathbb{R}^{d\times d}$ for the effective coefficient as the
solution to the system of linear equations

$\displaystyle\bar{a}^{L}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{i}+\nabla\phi_{i}^{L}\,\mathrm{d}x=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i}^{L})\,\mathrm{d}x\quad\text{for all
}1\leq i\leq d.$ (2.8)

Then for any $\delta>0$ and any $0<s\leq
c(d,\lambda,\nu)(\tfrac{L}{\varepsilon})^{\min\\{d/2,2\\}-\delta}$, we have
the error estimate

$\displaystyle\mathbb{P}\big{[}\big{|}\bar{a}^{L}-\bar{a}\big{|}>s\big{(}\tfrac{\varepsilon}{L}\big{)}^{\min\\{d/2,2\\}-\delta}\big{]}\leq
2\exp(-s^{1/C}/C).$

In particular, with at least this probability $\bar{a}^{L}$ is well-defined,
i. e. the defining linear system has a unique solution.

Overall, in dimensions $3$ and $4$ we see that with overwhelming probability,
the error $|\bar{a}^{L}-\bar{a}|$ is at most of order
$(\tfrac{\varepsilon}{L})^{d/2-\delta}$, improving the error estimate of the
standard RVE approximation to (nearly) the fluctuation scaling.

Note that as a slight downside of the formula (2.8), in principle we cannot
rule out that there may be a small probability that the linear system for
$\bar{a}^{L}$ given by (2.8) may fail to be invertible; thus, we cannot make a
statement on the error of the expected value $\mathbb{E}[\bar{a}^{L}]$.
However, the probability for this to happen is exceedingly low, as shown by
our stretched exponential estimate.

##  3 Overview of our arguments

###  3.1 Strategy for Theorem 1

Recall that a random field $\mathcal{C}(a,x)$, defined on
$\Omega\times\mathbb{H}^{d}_{+}$, having stretched exponential moments is
equivalent to there existing constants $m,C>0$ such that
$\left<|\mathcal{C}(a,x)|^{q}\right>^{1/q}\leq Cq^{m}$ for any $q\geq 1$ and
$x\in\mathbb{H}^{d}_{+}$. (By Jensen’s inequality it is, in fact, sufficient
to obtain the previous moment bound on $\mathcal{C}(a,x)$ only for $q\geq
q_{0}$ for some $q_{0}>1$.) For (2.4) it is, therefore, sufficient to show
that there exist constants $m,C>0$ and $q_{0}\geq 1$ such that

$\displaystyle\left<|\nabla\theta_{i}(x)|^{q}\right>^{\frac{1}{q}}\leq
Cq^{m}(1+x^{\perp})^{-\frac{d}{2}}\quad\text{ for any
}x\in\mathbb{H}^{d}_{+}\text{ and }q\geq q_{0}.$ (3.1)

Notice that, while Theorem 1 is a quenched result, the estimate (3.1) is
annealed, which is more convenient for using the spectral gap inequality (2.1)
from (A3). We also remark that, since in this contribution Theorem 1 is only
stated for $\mathbb{H}^{d}_{+}$ (rather than for a general domain), we may
w.l.o.g. set $\varepsilon=1$ –since the result for a general correlation
length follows from rescaling $a$.

Our argument for (3.1) comes in three main steps:

Step 1: We show a quenched suboptimal decay of $|\nabla\theta_{i}(x)|$ away
from $\partial\mathbb{H}^{d}_{+}$ –this is contained in Proposition 1. The
argument is deterministic, relying on already available estimates for the
whole space corrector $\phi_{i}$ and a Campanato-iteration type argument. The
result of Proposition 1 is slightly post-processed to obtain a corresponding
suboptimal annealed decay of a spatial average of $|\nabla\theta_{i}|^{2}$
(around $x$).

Step 2: We show that a decay estimate on moments of a spatial average of
$|\nabla\theta_{i}|^{2}$ (around $x$) is self-improving (but saturates at the
desired decay $-d/2+\delta$ in Theorem 1) –this is contained in Proposition 2.

Proposition 2 is the workhorse of our strategy: For notational convenience we
now fix $x_{0}\in\mathbb{H}^{d}_{+}$. To prove Proposition 2 we control the
moments of $\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{x_{0}^{\perp}/8}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x$
in terms of moments of suitable weighted averages
$F:=\int_{\mathbb{H}^{d}_{+}}g\cdot\nabla\theta_{i}\,\mathrm{d}x$. To control
the moments of the $F$, we control the moments of the fluctuations $F-\langle
F\rangle$ (Proposition 5) and the expected value $\langle F\rangle$
(Proposition 6). Proposition 5 relies on the spectral gap inequality (A3) and
sensitivity estimates –the arguments are contained in Section 5. The argument
for Proposition 6 is contained in Section 6; it relies on a careful two-scale
expansion argument, which exploits the stationarity of $\theta_{i}$ and the
standard correctors with respect to shifts tangential to
$\partial\mathbb{H}^{d}_{+}$.

Initializing the use of Proposition 2 with the result of Step 1, we
iteratively apply Proposition 2 –this yields Corollary 3.

Step 3: To finish, we post-process the result of Corollary 3 using the whole
space large-scale regularity theory from random linear elliptic operators
contained in [18] and the assumption (A4) on the regularity of the coefficient
field $a$. This post-processing is performed in Section 7 –and yields (3.1).

We now discuss Steps 1 and 2 in a bit more detail:

#### Step 1: A suboptimal quenched decay for $|\theta_{i}(x)|$ away from
$\partial\mathbb{H}^{d}_{+}$.

As already mentioned above, the result of this step is used to initialize the
iterative application of Proposition 2 in Step 2. For this purpose, we show
the following:

######  Proposition 1.

Under the assumptions of Theorem 1, there exists a random field
$\mathcal{C}(a,x)$, defined on $\Omega\times\mathbb{H}^{d}_{+}$, that is
stationary with respect to shifts tangential to $\partial\mathbb{H}^{d}_{+}$
and with stretched exponential moments in the sense of (2.3) such that

$\displaystyle|\nabla\theta_{i}(x)|\lesssim_{d,\lambda,\alpha}\mathcal{C}(a,x)\bigg{(}1+\frac{x^{\perp}}{\varepsilon}\bigg{)}^{-\frac{1}{3}},$
(3.2)

for any $x\in\mathbb{H}^{d}_{+}$.

The proof of Proposition 1 is given in Section 4 and, while relying on
stochastic ingredients from [18] –which are summarized in Appendix A–, is
entirely deterministic here. The argument, in particular, relies on the
comparison of $\theta$ to the solution of the homogenized problem with the
same boundary data (on a half-ball) via a homogenization-inspired Campanato
iteration. This strategy is quite standard, having been introduced by
Avellaneda and Lin in the setting of periodic homogenization in the 80s (see
e. g. [5]). More recently, this strategy has also been fruitfully exploited in
stochastic homogenization (see e. g. [3, 18, 4]). The version of these
techniques that we use in this contribution is closest to those in [18,
Theorem 1].

#### Step 2: Iterative improvement of the decay estimate.

The workhorse of our iterative strategy for proving Theorem 1 is the following
proposition:

######  Proposition 2.

Let assumptions (A1)-(A4) hold and assume that there exist constants $m,n>0$
such that

$\displaystyle\left<\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\frac{x_{0}^{\perp}}{2}}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\bigg{)}^{\frac{q}{2}}\right>^{\frac{1}{q}}\lesssim_{d,\lambda,\alpha}q^{m}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-n},$
(3.3)

for any $x_{0}\in\mathbb{H}^{d}_{+}$ and $q\geq 1$. Then there exists a
constant $\tilde{m}>0$ such that

$\displaystyle\left<\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\frac{x_{0}^{\perp}}{8}}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\bigg{)}^{\frac{q}{2}}\right>^{\frac{1}{q}}\lesssim_{d,\lambda,\alpha}q^{\tilde{m}}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-\min\\{n+1-\delta,\frac{d}{2}-\delta\\}},$
(3.4)

for any $x_{0}\in\mathbb{H}^{d}_{+}$, $q\geq 1$, and $0<\delta\ll 1$.

Slightly postponing a thorough outline of the proof of Proposition 2 to
Section 3.2, we first discuss how we apply this result. To initiate the
iterative application of Proposition 2, we require access to an estimate like
(3.3). As already mentioned above, for this we post-process (3.2) of
Proposition 1 to obtain

$\displaystyle\begin{split}\left<\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\frac{x_{0}^{\perp}}{8}}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\bigg{)}^{\frac{q}{2}}\right>^{\frac{1}{q}}&\leq\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\frac{x_{0}^{\perp}}{8}}(x_{0})}\left<|\nabla\theta_{i}|^{q}\right>\,\mathrm{d}x\bigg{)}^{\frac{1}{q}}\\\
&\stackrel{{\scriptstyle\eqref{intermediate}}}{{\lesssim}}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-\frac{1}{3}}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\frac{x_{0}^{\perp}}{8}}(x)}\left<\mathcal{C}(a,x_{0})^{q}\right>{\textrm{d}}y\bigg{)}^{\frac{1}{q}}\\\
&\lesssim
q^{m}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-\frac{1}{3}},\end{split}$
(3.5)

for any $x_{0}\in\mathbb{H}^{d}_{+}$ and $q\geq 1$. Using (3.5), we can
initialize the use of Proposition 2 with $n=\frac{1}{3}$. Proposition 2 then
yields:

######  Corollary 3.

Under the assumptions of Theorem 1, there exists a constant $m>0$ such that

$\displaystyle\left<\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\frac{x_{0}^{\perp}}{8}}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\bigg{)}^{\frac{q}{2}}\right>^{\frac{1}{q}}\lesssim
q^{m}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-\frac{d}{2}+\delta},$
(3.6)

for any $x_{0}\in\mathbb{H}^{d}_{+},\,q\geq 1$, and $0<\delta\ll 1$.

Since the argument for Corollary 3 is obvious given the discussion proceeding
its statement, we forgo a formal proof.

###  3.2 Argument for Proposition 2: Iterative improvement of decay estimate
for the boundary layer

As elaborated on in the previous section, the core of our argument for Theorem
1 is the iterative application of Proposition 2. In this section we state the
other results required for our proof of Proposition 2 and, assuming that these
hold, give the proof of the proposition.

The first ingredient that we require is an estimate that controls the gradient
of an $a$-harmonic function in terms of linear functionals of the gradient,
plus a small remainder term.

######  Lemma 4.

Let $a$ be a uniformly elliptic and bounded coefficient field in the sense of
Assumption (A1). Let $u$ be $a$-harmonic in a ball $B_{R}$, $R\geq 1$. Then
for any $\delta>0$ and any $A<\infty$ there exists a set of functions
$g_{m,k}$, $1\leq m\leq M=M(\lambda,\delta,p)$, $1\leq
k\leq\lfloor\log_{2}\frac{R}{R^{1-\delta}}\rfloor$, subject to the bounds

$\displaystyle|g_{k,m}|$ $\displaystyle\leq(2^{k}R^{1-\delta})^{-d},$
$\displaystyle|\nabla g_{k,m}|$ $\displaystyle\leq(2^{k}R^{1-\delta})^{-1-d},$
$\displaystyle\operatorname{supp}g_{k,m}$ $\displaystyle\subset
B_{2^{k}R^{1-\delta}},$

and such that the estimate

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R^{1-\delta}}}|\nabla u|^{2}{\textrm{d}}x\leq
R^{-A}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R}}|\nabla
u|^{2}{\textrm{d}}x+C\sum_{k=1}^{\lfloor\log_{2}\frac{R}{R^{1-\delta}}\rfloor}(2^{k})^{-A/\delta}\sum_{m=1}^{M}\bigg{|}\int_{B_{R}}g_{k,m}\cdot\nabla
u{\textrm{d}}x\bigg{|}^{2}$ (3.7)

holds.

As the proof is rather short, we directly provide it here.

###### Proof.

We first notice that for any $\tilde{\delta}>0$ there exists a finite set of
functions $g_{m}$, $1\leq m\leq M$, supported in $B_{R}$ and satisfying
$|g_{m}|\leq R^{-d}$, $|\nabla g_{m}|\leq R^{1-d}$, such that any $a$-harmonic
function $u$ on $B_{R}$ satisfies

$\displaystyle\int_{B_{R/2}}|\nabla
u|^{2}{\textrm{d}}x\leq\tilde{\delta}\int_{B_{R}}|\nabla
u|^{2}{\textrm{d}}x+C\sum_{m=1}^{M}\bigg{|}\int_{B_{R}}g_{m}\cdot\nabla
u{\textrm{d}}x\bigg{|}^{2}.$ (3.8)

In fact, this may be obtained by a suitable variant of e. g. [6, Lemma 12].

Iterating this bound $K:=\lfloor\log_{2}\frac{R}{R^{1-\delta}}\rfloor$ many
times, we arrive at

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2^{-K}R}}|\nabla
u|^{2}{\textrm{d}}x\leq(2^{d}\tilde{\delta})^{I}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R}}|\nabla
u|^{2}{\textrm{d}}x+C\sum_{k=1}^{K}\sum_{m=1}^{M}(2^{d}\tilde{\delta})^{k-1}\bigg{|}\int_{B_{R}}g_{k,m}\cdot\nabla
u{\textrm{d}}x\bigg{|}^{2}.$

Choosing $\tilde{\delta}>0$ small enough, this implies the desired estimate. ∎

Lemma 4 reduces the proof of Proposition 2 to estimating linear functionals of
$\nabla\theta_{i}$: The main challenge is to control the $2q$-th stochastic
moments of the
$F_{k,m}:=\int_{B_{x_{0}^{\perp}}(x_{0})}g_{k,m}\cdot\nabla\theta_{i}\,\mathrm{d}x$.
Towards this goal, we begin by applying the triangle inequality as

$\displaystyle\left<|F_{k,m}|^{2q}\right>^{\frac{1}{2q}}\leq\left<|F_{k,m}-\left<F_{k,m}\right>|^{2q}\right>^{\frac{1}{2q}}+\big{|}\left<F_{k,m}\right>\big{|}.$
(3.9)

As already alluded to above, if we were deriving estimates for the whole space
corrector $\phi_{i}$ instead of $\theta_{i}$, that is, if we were considering
$F_{k,m}=\int g_{k,m}\cdot\nabla\phi_{i}\,\mathrm{d}x$ instead of
$F_{k,m}=\int g_{k,m}\cdot\nabla\theta_{i}\,\mathrm{d}x$, the expected value
$\left<F_{k,m}\right>$ would vanish due to the stationarity of
$\nabla\phi_{i}$ and the sublinear growth of $\phi_{i}$ (see e. g. [26, 24]).
Thus, the problem would entirely reduce to that of fluctuation bounds for the
$F_{k,m}$. In our situation, however, since $\nabla\theta_{i}$ is not
stationary in the direction perpendicular to $\partial\mathbb{H}^{d}_{+}$, we
need to control both terms on the right-hand side of (3.9): fluctuation bounds
for $F_{k,m}$ are obtained in Proposition 5 and the expected value $\langle
F_{k,m}\rangle$ is treated in Proposition 6.

The statements of Propositions 5 and 6 are given below, where we remark that,
for use in future work, we have proven Proposition 5 not only on
$\mathbb{H}^{d}_{+}$, but also on a generic bounded $C^{1,1}$-domain. Since
the our treatment of the expectations $\left<F_{k,m}\right>$ in Proposition 6
depends quite heavily on our estimate being on the half-space, we do not
generalize this statement to a generic domain in this contribution.

######  Proposition 5.

For $d\geq 3$, let $\langle\cdot\rangle$ be an ensemble of coefficient fields
on $\mathbb{R}^{d}$ that satisfies the assumptions (A1)-(A4). Furthermore, let
${\mathcal{O}}$ denote either a bounded $C^{1,1}$-domain or
$\mathbb{H}^{d}_{+}$; let $f\in W^{1,\infty}(\mathbb{R}^{d})$ be smooth such
that, in the case that ${\mathcal{O}}=\mathbb{H}^{d}_{+}$, $f\equiv 1$ and, if
${\mathcal{O}}$ is bounded, $f$ is supported in a bounded neighborhood of
${\mathcal{O}}$. Let $p$ satisfy $0<p-1\ll 1$, $r\geq\varepsilon$, and
$x_{0}\in{\mathcal{O}}$. Let $g\in L^{2p}(\mathbb{R}^{d};\mathbb{R}^{d})$
satisfy $\operatorname{supp}(g)\subset B_{r}(x_{0})$ with
$\operatorname{dist}(B_{r}(x_{0}),\partial{\mathcal{O}})\geq\frac{r}{4}$ and

$\displaystyle\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}(x_{0})}|g|^{2p}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\leq
r^{-d}.$ (3.10)

Consider $u$ solving

$\displaystyle-\nabla\cdot(a\nabla u)$ $\displaystyle=0$
$\displaystyle\textrm{ in }{\mathcal{O}}$ (3.11) $\displaystyle u$
$\displaystyle=f\phi_{i}$ $\displaystyle\textrm{ on }\partial{\mathcal{O}}$

and the random variable $F$ defined as

$\displaystyle F:=\int_{\mathcal{O}}g\cdot\nabla u\,\mathrm{d}x.$ (3.12)

Then there exists a constant $m=\kappa(d,\lambda,{\mathcal{O}},f)<\infty$ such
that

$\displaystyle\left<\big{|}F-\left<F\right>\big{|}^{2q}\right>^{\frac{1}{2q}}\lesssim_{d,\lambda,{\mathcal{O}},f}q^{\kappa}\bigg{(}\Big{(}\frac{\varepsilon}{r}\Big{)}^{\beta\frac{d}{2}}\bigg{(}\frac{\varepsilon}{\operatorname{dist}(x_{0},\partial{\mathcal{O}})+\varepsilon}\bigg{)}^{(1-\beta)\frac{d}{2}}+\frac{\varepsilon}{\operatorname{dist}(x_{0},\partial{\mathcal{O}})+\varepsilon}\bigg{(}\frac{\varepsilon}{r}\bigg{)}^{\frac{d}{2}}\bigg{)},$
(3.13)

for any $q\geq 1$ and $0<\beta\ll 1$.

As already mentioned, to prove Proposition 5 we rely on the spectral gap
inequality (2.1) –in particular, for $\xi=F$, we express the Malliavin
derivative on the right-hand side of (2.1) in terms of auxiliary functions
$h$, $\hat{h}$, $v_{1}$, and $v_{2}$ solving (5.2), (5.3), (5.4), and (5.5)
respectively. We then proceed to estimate these auxiliary functions via
elliptic regularity estimates. The argument is contained in Section 5 and is
similar to e. g. [19, 12, 6, 7].

Here is the statement of Proposition 6, which we use to treat the expected
values $\left<F_{k,m}\right>$ on the right-hand side of (3.9):

######  Proposition 6.

For $d\geq 3$, let $\langle\cdot\rangle$ be an ensemble of coefficient fields
that satisfies the assumptions (A1)-(A4). Suppose that there exists $n>0$ such
that

$\left<\left|\nabla\theta_{i}(x)\right|^{2}\right>^{\frac{1}{2}}\lesssim_{d,\lambda}\bigg{(}1+\frac{x^{\perp}}{\varepsilon}\bigg{)}^{-n},$
(3.14)

for any $x\in\mathbb{H}^{d}_{+}$. Let $r\geq\varepsilon$,
$x_{0}\in\mathbb{H}^{d}_{+}$, and $g\in
L^{2p}(\mathbb{H}^{d}_{+};\mathbb{R}^{d})$ satisfy
$\operatorname{supp}(g)\subset B_{r}(x_{0})\cap\mathbb{H}^{d}_{+}$ with
$\operatorname{dist}(B_{r}(x_{0}),\partial\mathbb{H}^{d}_{+})\geq\frac{r}{4}$
and

$\displaystyle\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}(x_{0})}|g|^{2p}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\leq
r^{-d}$ (3.15)

as well as

$\displaystyle\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}(x_{0})}|\nabla
g|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}\leq r^{-(d+1)}.$ (3.16)

Then the expected value $\langle F\rangle$ of the linear functional
$F:=\int_{\mathbb{H}^{d}_{+}}g\cdot\nabla\theta\,\mathrm{d}x$ satisfies the
improved decay estimate

$\displaystyle\left|\left\langle\int_{\mathbb{H}^{d}_{+}}\nabla\theta_{i}\cdot
g{\textrm{d}}x\right\rangle\right|\lesssim_{d,\lambda}\frac{\varepsilon}{r}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-n}.$
(3.17)

We argument for Proposition 6 is given in Section 6. The underlying idea is to
use a two-scale expansion approximation for the “Greens function analogue”
$-\nabla\cdot(a^{*}\nabla h)=\nabla\cdot g$ to reduce – up to the higher-order
error term resulting in the right-hand side of (3.17) – the expression for $F$
to a corresponding expression involving only the correctors and the
homogenized analogue $-\nabla\cdot(\bar{a}^{*}\nabla\bar{h})=\nabla\cdot g$.
One then exploits stationarity in the tangential direction to conclude that
the remaining expressions in fact vanish.

With all of the ingredients for the proof of Proposition 2 ready, we now
proceed to its proof.

###### Proof of Proposition 2.

Set $R:=\frac{x_{0}^{\perp}}{8}$. Inserting our assumption (3.3) into Lemma 4
(applied for $A=1$), we obtain

$\displaystyle\bigg{\langle}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R^{1-\delta}}(x_{0})}|\nabla\theta_{i}|^{2}{\textrm{d}}x\bigg{)}^{q}\bigg{\rangle}^{\frac{1}{q}}$
$\displaystyle\lesssim_{d,\lambda,\alpha}q^{m}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-n-1}+\sum_{k=1}^{\lfloor\log_{2}\frac{R}{R^{1-\delta}}\rfloor}(2^{k})^{-1/\delta}\bigg{\langle}\bigg{(}\sum_{m=1}^{M}\bigg{|}\int_{B_{R}}g_{k,m}\cdot\nabla\theta_{i}{\textrm{d}}x\bigg{|}^{2}\bigg{)}^{q}\bigg{\rangle}^{\frac{1}{q}}.$

To bound the fluctuations of the functionals on the right-hand side, we apply
Proposition 5 with $r=2^{k}R^{1-\delta}$ and $g=g_{k,m}$. Furthermore, to
bound the expected values of these functionals we apply Proposition 6 with
$g=g_{k,m}$ and $r=2^{k}R^{1-\delta}$; note that we obtain the pointwise bound
(3.14) required for its application from our hypothesis (3.3) and a down-
propagation to a pointwise bound via (B.3). Inserting the resulting estimates
into the previous bound, we obtain

$\displaystyle\bigg{\langle}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R^{1-\delta}}(x_{0})}|\nabla\theta_{i}|^{2}{\textrm{d}}x\bigg{)}^{q}\bigg{\rangle}^{\frac{1}{q}}$
$\displaystyle\lesssim_{d,\lambda,\alpha}q^{m}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-n-1}$
$\displaystyle\,\,+q^{\kappa}\sum_{k=1}^{\lfloor\log_{2}\frac{R}{R^{1-\delta}}\rfloor}(2^{k})^{-1/\delta}\bigg{(}\Big{(}\frac{\varepsilon}{2^{k}R^{1-\delta}}\Big{)}^{\beta\frac{d}{2}}\bigg{(}\frac{\varepsilon}{x_{0}^{\perp}+\varepsilon}\bigg{)}^{(1-\beta)\frac{d}{2}}+\frac{\varepsilon}{x_{0}^{\perp}+\varepsilon}\bigg{(}\frac{\varepsilon}{r}\bigg{)}^{\frac{d}{2}}\bigg{)}$
$\displaystyle\,\,+\sum_{k=1}^{\lfloor\log_{2}\frac{R}{R^{1-\delta}}\rfloor}(2^{k})^{-1/\delta}\frac{\varepsilon}{2^{k}R^{1-\delta}}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-n}.$

Estimating the sums, inserting the definition of $R$, and choosing $\beta$
sufficiently small (depending on $\delta$), we obtain

$\displaystyle\bigg{\langle}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R^{1-\delta}}(x_{0})}|\nabla\theta_{i}|^{2}{\textrm{d}}x\bigg{)}^{q}\bigg{\rangle}^{\frac{1}{q}}\lesssim_{d,\lambda,\alpha}q^{m}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-n-1+\delta}+q^{\kappa}\bigg{(}1+\frac{x_{0}^{\perp}}{\varepsilon}\bigg{)}^{-d/2+\delta}.$

Letting $x_{0}$ vary over a ball of radius $x_{0}^{\perp}/8$ around $x_{0}$
and integrating proves Proposition 2. ∎

###  3.3 Proof of Theorem 2

As the proof of our second main result Theorem 2 on the accuracy of the
representative volume approximation is rather short, we directly state it
here.

###### Proof of Theorem 2.

Without loss of generality, in the proof we assume $\varepsilon=1$.
_Step 1: Representation formula for $\phi_{i}^{L}$._ Denote by
$\mathbb{H}_{m}$, $1\leq m\leq 2d$, the half-spaces
$\\{x\in\mathbb{R}:x_{m}>0\\}$, $1\leq m\leq d$, and
$\\{x\in\mathbb{R}:x_{m-d}<L\\}$, $d+1\leq m\leq 2d$. Denote by slight abuse
of notation by $\theta_{i}^{\mathbb{H}_{m}}(\infty)$ the constants that the
boundary layer correctors $\theta_{i}^{\mathbb{H}_{m}}$ converge to (almost
surely) as the distance to the half-space boundary becomes infinite.
Furthermore, denote by
$E:=(\\{0,L\\}\times\\{0,L\\}\times[0,L]^{d-2})\cup\ldots\cup([0,L]^{d-2}\times\\{0,L\\}\times\\{0,L\\})$
the set of all edges in case $d=3$ respectively for $d\geq 4$ the set of all
$d-2$-dimensional hyperedges of the cube.

Note that by the defining equations, we may write

$\displaystyle\phi_{i}^{L}=\phi_{i}-\sum_{m=1}^{2d}(\theta_{i}^{\mathbb{H}_{m}}-\gamma_{i,m}^{L}-\tau_{i,m}^{L}),$
(3.18)

where $\gamma_{i,m}^{L}$ is defined as the solution to the boundary value
problem

$\displaystyle-\nabla\cdot(a\nabla\gamma_{i,m}^{L})$ $\displaystyle=0$
$\displaystyle\text{in }(0,L)^{d},$ $\displaystyle\gamma_{i,m}^{L}$
$\displaystyle=\theta_{i}^{\mathbb{H}_{m}}(\infty)$ $\displaystyle\text{on
}\partial[0,L]^{d}\setminus\partial\mathbb{H}_{m},$
$\displaystyle\gamma_{i,m}^{L}$ $\displaystyle=0$ $\displaystyle\text{on
}\partial[0,L]^{d}\cap\partial\mathbb{H}_{m},$

and where $\tau_{i,m}^{L}$ solves $-\nabla\cdot(a\nabla\tau_{i,m}^{L})=0$ in
$(0,L)^{d}$ with

$\displaystyle\tau_{i,m}^{L}$
$\displaystyle=\theta_{i}^{\mathbb{H}_{m}}-\theta_{i}^{\mathbb{H}_{m}}(\infty)$
$\displaystyle\text{on }\partial[0,L]^{d}\setminus\partial\mathbb{H}_{m},$
$\displaystyle\tau_{i,m}^{L}$ $\displaystyle=0$ $\displaystyle\text{on
}\partial[0,L]^{d}\cap\partial\mathbb{H}_{m}.$

Observe in particular that
$\theta_{i}^{\mathbb{H}_{m}}-\gamma_{i,m}^{L}-\tau_{i,m}^{L}\equiv 0$ on
$\partial[0,L]^{d}\setminus\partial\mathbb{H}_{m}$ and
$\theta_{i}^{\mathbb{H}_{m}}-\gamma_{i,m}^{L}-\tau_{i,m}^{L}=\phi_{i}$ on
$\partial[0,L]^{d}\cap\partial\mathbb{H}_{m}$. Notice that the boundary
conditions for $\tau_{i}^{L}$ in particular implies in view of (2.4) from
Theorem 1 that

$\displaystyle|\tau_{i,m}^{L}|(x)\leq\mathcal{C}(1+\operatorname{dist}(x,E))^{1-d/2+\delta}$
$\displaystyle\text{on }\partial[0,L]^{d}.$ (3.19)

Note also that some mathematical issues may arise due to the discontinuity of
the boundary data for $\gamma_{i,m}^{L}$ and $\tau_{i,m}^{L}$ at the edges; to
maintain strict mathematical rigor, a smoothing argument on a small scale
would be required. However, none of the below arguments will be affected by
the smoothing procedure, hence we omit details of this technicality.

The representation (3.18) enables us to split

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i}^{L}){\textrm{d}}x$
$\displaystyle=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i}){\textrm{d}}x-\sum_{m=1}^{2d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a\nabla\theta_{i}^{\mathbb{H}_{m}}{\textrm{d}}x$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\sum_{m=1}^{2d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a\nabla\gamma_{i,m}^{L}{\textrm{d}}x+\sum_{m=1}^{2d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a\nabla\tau_{i,m}^{L}{\textrm{d}}x.$

While the first term on the right-hand side is the desired one and the second
term may be estimated by (2.4) from Theorem 1, it in particular remains to
account for the contributions of the last two terms.

_Step 2: Estimate on $\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int\nabla\tau_{i,m}^{L}{\textrm{d}}x$._ Our first goal
is to bound the contributions

$\displaystyle\bigg{|}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}\nabla\tau_{i,m}^{L}{\textrm{d}}x\bigg{|}+\bigg{|}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a\nabla\tau_{i,m}^{L}{\textrm{d}}x\bigg{|}\leq\mathcal{C}L^{-3/2}.$
(3.20)

Heuristically, by the Poisson kernel formula, we expect (except at the
corners)

$\displaystyle|\nabla\tau_{i,m}^{L}|(x)$
$\displaystyle\leq\mathcal{C}(x)\int_{\partial[0,L]^{d}}\frac{\mathcal{C}(\tilde{x})}{|x-\tilde{x}|^{d}}|\tau_{i,m}^{L}(\tilde{x})|\,dS(\tilde{x})$
$\displaystyle\leq\mathcal{C}(x)\operatorname{dist}(x,E)^{1-d/2+\delta}\operatorname{dist}(x,\partial[0,L]^{d})^{-1}+\mathcal{C}(x)\operatorname{dist}(x,E)^{-2}.$

However, this is in fact a more precise estimate than actually required for
our conclusion. To make the required estimates rigorous without too much
effort, we define $h$ as the solution to $-\nabla\cdot(a^{*}\nabla
h)=\nabla\cdot(\chi_{[\kappa L,(1-\kappa)L]^{d}}e_{k})$ in $(0,L)^{d}$ with
$h=0$ on $\partial[0,L]^{d}$. Using the equations for $h$ and (in the second
step) for $\tau_{i,m}^{L}$, we obtain

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{k}\cdot\nabla\tau_{i,m}^{L}{\textrm{d}}x$
$\displaystyle=(1-2\kappa)^{-d}L^{-d}\int_{[0,L]^{d}}\tau_{i,m}^{L}{\mathbf{n}}\cdot
a^{*}\nabla h{\textrm{d}}S$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
-(1-2\kappa)^{-d}L^{-d}\int_{[0,L]^{d}}a\nabla\tau_{i,m}^{L}\cdot\nabla
h{\textrm{d}}x$
$\displaystyle=(1-2\kappa)^{-d}L^{-d}\int_{[0,L]^{d}}\tau_{i,m}^{L}{\mathbf{n}}\cdot
a^{*}\nabla h{\textrm{d}}S.$

A regularity theory estimate for $\nabla h$ (being $a$-harmonic in a
neighborhood of the boundary) yields (see for instance [23] for the two-
dimensional case, the higher-dimensional case being analogous)

$\displaystyle|\nabla
h|(x)\leq\mathcal{C}(x)\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[0,L]^{d}}|\nabla
h|^{2}{\textrm{d}}x\bigg{)}^{1/2}\leq\mathcal{C}(x).$

Plugging this into the previous bound and using the estimate (3.19) for the
boundary data of $\tau_{i,m}^{L}$, we arrive at

$\displaystyle\bigg{|}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{k}\cdot\nabla\tau_{i,m}^{L}{\textrm{d}}x\bigg{|}\leq\mathcal{C}(L^{-d/2+\delta}+L^{-2}).$

Proceeding analogously for the flux (simply inserting a factor of $a^{*}$ in
the right-hand side of the equation for $h$), we obtain (3.20).

_Step 3: Estimate on $\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int\nabla\gamma_{i,m}^{L}{\textrm{d}}x$._ Concerning
$\gamma_{i,m}^{L}$, we employ a two-scale expansion ansatz

$\displaystyle\gamma_{i,m}^{L}=\bar{\gamma}_{i,m}^{L}+\phi_{j}^{L}\partial_{j}\bar{\gamma}_{i,m}^{L}+w_{i,m}^{L}$

with $\bar{\gamma}_{i,m}^{L}$ solving

$\displaystyle-\nabla\cdot(\bar{a}\nabla\bar{\gamma}_{i,m}^{L})$
$\displaystyle=0$ $\displaystyle\text{in }(0,L)^{d},$
$\displaystyle\bar{\gamma}_{i,m}^{L}$
$\displaystyle=\theta_{i}^{\mathbb{H}_{m}}(\infty)$ $\displaystyle\text{on
}\partial[0,L]^{d}\setminus\partial\mathbb{H}_{m},$
$\displaystyle\bar{\gamma}_{i,m}^{L}$ $\displaystyle=0$ $\displaystyle\text{on
}\partial[0,L]^{d}\cap\partial\mathbb{H}_{m},$

and with $w_{i}^{L}$ solving

$\displaystyle-\nabla\cdot(a\nabla
w_{i,m}^{L})=\nabla\cdot((a\phi_{j}^{L}-\sigma_{j})\partial_{j}\nabla\bar{\gamma}_{i,m}^{L})+a(\nabla\phi_{j}^{L}-\nabla\phi_{j})\cdot\nabla\partial_{j}\bar{\gamma}_{i,m}^{L}$

with homogeneous Dirichlet boundary data. Note that the equation satisfied by
$w_{i,m}^{L}$ is obtained by the usual property of the two-scale expansion

$\displaystyle-\nabla\cdot\bigg{(}a\nabla\bigg{(}\bar{\gamma}_{i,m}^{L}+\sum_{j}\phi_{j}\partial_{j}\bar{\gamma}_{i,m}^{L}\bigg{)}\bigg{)}=-\nabla\cdot((a\phi_{j}-\sigma_{j})\partial_{j}\nabla\bar{\gamma}_{i,m}^{L}).$

Due to the piecewise constant boundary data we have
$|D^{2}\bar{\gamma}_{j}^{L}|(x)\leq C\operatorname{dist}(x,E)^{-2}$.

Defining the function $h$ as the solution to the equation
$-\nabla\cdot(a^{*}\nabla h)=\nabla\cdot(\chi_{[\kappa
L,(1-\kappa)L]^{d}}e_{k})$ in $(0,L)^{d}$ with homogeneous Dirichlet boundary
data, we deduce using first this equation and then the equation for
$w_{i,m}^{L}$

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa L,(1-\kappa)L]^{d}}e_{k}\cdot\nabla
w_{i,m}^{L}{\textrm{d}}x=-(1-2\kappa)^{-d}L^{-d}\int_{(0,L)^{d}}a^{*}\nabla
h\cdot\nabla w_{i,m}^{L}{\textrm{d}}x$
$\displaystyle=(1-2\kappa)^{-d}L^{-d}\int_{(0,L)^{d}}(a\phi_{j}^{L}-\sigma_{j})\partial_{j}\nabla\bar{\gamma}_{i,m}^{L}\cdot\nabla
h{\textrm{d}}x$ $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
-(1-2\kappa)^{-d}L^{-d}\int_{(0,L)^{d}}ha(\nabla\phi_{j}^{L}-\nabla\phi_{j})\cdot\nabla\partial_{j}\bar{\gamma}_{i,m}^{L}{\textrm{d}}x.$

Again, by regularity theory we have $|\nabla h|(x)\leq\mathcal{C}(x)$ and thus
$|h(x)|\leq\mathcal{C}(x)\operatorname{dist}(x,\partial[0,L]^{d})$.
Furthermore, applying Lemma 25 to $\zeta(\phi_{j}^{L}-\phi_{j})$ for a
suitable cutoff $\zeta$ and inserting the bound from Theorem 4, we see that
$|\nabla\phi_{j}^{L}-\nabla\phi_{j}|\leq\mathcal{C}(x)L^{-1+\delta}$ on
$[\kappa L,(1-\kappa)L]^{d}$. We thus obtain (using again an analogous
argument for the flux contribution)

$\displaystyle\bigg{|}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}\nabla\gamma_{i,m}^{L}{\textrm{d}}x-\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}(e_{j}+\nabla\phi_{j}^{L})\partial_{j}{\bar{\gamma}}_{i,m}^{L}{\textrm{d}}x\bigg{|}$
$\displaystyle+\bigg{|}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a\nabla\gamma_{i,m}^{L}{\textrm{d}}x-\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{j}+\nabla\phi_{j}^{L})\partial_{j}{\bar{\gamma}}_{i,m}^{L}{\textrm{d}}x\bigg{|}$
(3.21) $\displaystyle\leq\mathcal{C}L^{-2+\delta}.$

_Step 4: Conclusion._ Using the splitting (3.18), the decay estimate for the
boundary layer corrector gradient (2.4), the bounds (3.20) and (3.21), and the
classical fluctuation estimate $|\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}\nabla\phi_{i}^{L}\,\mathrm{d}x|\leq\mathcal{C}L^{-d/2}$,
we deduce

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{i}+\nabla\phi_{i}^{L}\,\mathrm{d}x=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{i}+\sum_{m=1}^{2d}\nabla\bar{\gamma}_{i,m}^{L}\,\mathrm{d}x+\mathcal{R}_{1}$
(3.22a) and $\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i}^{L})\,\mathrm{d}x$ (3.22b)
$\displaystyle=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{i}+\nabla\phi_{i})\,\mathrm{d}x+\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{j}+\nabla\phi_{j}^{L})\sum_{m=1}^{2d}\partial_{j}\bar{\gamma}_{i,m}^{L}\,\mathrm{d}x+\mathcal{R}_{2}$
$\displaystyle=\bar{a}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}e_{i}+\sum_{m=1}^{2d}\nabla\bar{\gamma}_{i,m}^{L}\,\mathrm{d}x+\mathcal{R}_{3}$
(3.22c)

with random variables $\mathcal{R}_{k}$ estimated as
$|\mathcal{R}_{1}|+|\mathcal{R}_{2}|+|\mathcal{R}_{3}|\leq\mathcal{C}(L^{-d/2+\delta}+L^{-2+\delta})$
for some random variable $\mathcal{C}$ with uniformly bounded stretched
exponential moments. Here in the last step we have used the (classical) result
that

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}a(e_{j}+\nabla\phi_{j})\,\mathrm{d}x=\bar{a}e_{j}+\mathcal{R}$

with $|\mathcal{R}|\leq\mathcal{C}L^{-d/2}$, see e. g. [19].

Note that if all the $\mathcal{R}_{j}$ were zero, the equations (3.22) and the
definition (2.8) together would imply $\bar{a}^{L}=\bar{a}$. In general, we
see that the quantity $\bar{a}^{L}$ is obtained by inverting a
$d^{2}$-dimensional system of linear equations, and in fact this system is
(with overwhelming probability) a perturbation of order
$L^{-\min\\{d/2,2\\}+\delta}$ of well-conditioned deterministic system with
solution $\bar{a}$ (note that the deterministic system is well-conditioned
since $\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{[\kappa
L,(1-\kappa)L]^{d}}|\nabla\bar{\gamma}_{i,m}^{L}|\,\mathrm{d}x\lesssim
L^{-1}$). The perturbation series for linear systems of equations now yields
the desired assertion. ∎

##  4 Argument for Proposition 1: A suboptimal decay estimate for the boundary
layer

The main step towards proving Proposition 1 is the following lemma:

######  Lemma 7.

Under the assumptions of Theorem 1, there exist random fields
$\mathcal{C}(a,x)$ and $\tilde{C}(a,x)$, defined on
$\Omega\times\mathbb{H}^{d}_{+}$, that are stationary with respect with
respect to tangential shifts and with stretched exponential moments –in the
sense of (2.3)– such that

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\leq\mathcal{C}(a,x_{0})\bigg{(}1+\frac{r}{\varepsilon}\bigg{)}^{-\frac{2}{3}}$
(4.1)

for any $x_{0}\in\partial\mathbb{H}^{d}_{+}$ and
$r\geq\tilde{C}(a,x_{0})\varepsilon$.

To obtain (4.1) we morally “transfer the regularity” from a corresponding
homogenized solution (with the same boundary data as $\theta_{i}$ on a ball)
onto $\theta_{i}$. We do this via a homogenization-inspired Campanato
iteration.

The argument for Lemma 7 uses four auxiliary lemmas, which are proven in
Section 4.2. In the first auxiliary lemma, we state an estimate for the whole
space corrector that is obtained by combining Theorem 4 of Appendix A with the
maximal ergodic theorem.

######  Lemma 8.

Under the assumptions of Theorem 1, there exists a stationary random field
$\mathcal{C}(a,x)$, defined on $\Omega\times\partial\mathbb{H}^{d}_{+}$, with
stretched exponential moments –in the sense of (2.3)–, such that

$\displaystyle\sup_{r\geq 1}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\phi|^{2}+|\sigma|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}\lesssim_{d,\lambda}\mathcal{C}(x_{0},a)\varepsilon,$
(4.2)

for any $x_{0}\in\partial\mathbb{H}^{d}_{+}$.

The second auxiliary lemma is a weighted energy estimate that is a consequence
of the classical Hardy inequality.

######  Lemma 9.

Let the assumption (A1) hold and $x_{0}\in\partial\mathbb{H}^{d}_{+}$. Then
there exists $\kappa=\kappa(d,\lambda)>0$ such that, for any $r>0$ and $w\in
H^{1}_{\rm{loc}}(\mathbb{H}^{d}_{+})$ a solution of

$\displaystyle-\nabla\cdot(a\nabla u)$ $\displaystyle=\nabla\cdot g+f$
$\displaystyle\text{in }B^{+}_{r}(x_{0}),$ $\displaystyle u$ $\displaystyle=0$
$\displaystyle\text{on }\partial\mathbb{H}^{d}_{+}\cap B^{+}_{r}(x_{0}),$

there holds

$\displaystyle\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa}|\nabla
u|^{2}\,\mathrm{d}x\lesssim_{d,\lambda,\kappa}\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa}(|g|^{2}+r^{2}f^{2})\,\mathrm{d}x$
(4.3)

The last auxiliary lemma gives the constant-coefficient regularity results
that we “transfer” onto the heterogenous solution $\theta_{i}$. Both Lemmas 9
and 10 are also used in the proof of Proposition 24.

######  Lemma 10.

Let $d\geq 2$ and $x_{0}\in\overline{\mathbb{H}}^{d}_{+}$. Furthermore, let
$\bar{a}\in\mathbb{R}^{d\times d}$ be uniformly elliptic and bounded with
ellipticity constant $\lambda>0$ and $||\bar{a}||\leq 1$. Let $v$ be
$\bar{a}$-harmonic in $B^{+}_{R}(x_{0})$ with homogeneous Dirichlet boundary
data on $\partial\mathbb{H}^{d}_{+}\cap B_{R}^{+}(x_{0})$. Then, for any
$\rho\in(0,\frac{R}{2}]$, we have that

$\displaystyle\sup_{B^{+}_{R-\rho}(x_{0})}(\rho^{2}|\nabla^{2}v|^{2}+|\nabla
v|^{2})\lesssim_{d,\lambda}\bigg{(}\frac{R}{\rho}\bigg{)}^{d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\nabla v|^{2}\,\mathrm{d}x.$
(4.4)

This lemma follows from using the Sobolev embedding and iterating the
Caccoppoli inequality in directions tangential to
$\partial\mathbb{H}^{d}_{+}$. To recover the desired estimate in the
$e_{1}$-direction, one uses the equation solved by $v$ –this strategy is also
used in Step 3 of the proof of Lemma 22. The full proof can be found in [13,
Lemma 4.2], we do not repeat it here.

Once we have access to Lemma 7, it remains to post-process (4.1) to obtain the
pointwise estimate (3.2). For this we use the following lemma, which is also
proved in Section 4.2 and is a simple consequence of standard Schauder theory.

######  Lemma 11.

Let the assumption (A4) be satisfied and $u$ be a weak solution of

$\displaystyle-\nabla\cdot(a\nabla u)$ $\displaystyle=\nabla\cdot g+f$
$\displaystyle\text{in }\mathbb{H}^{d}_{+},$ $\displaystyle u$
$\displaystyle=0$ $\displaystyle\text{on }\partial\mathbb{H}^{d}_{+}.$

Then, for any $x_{0}\in\mathbb{H}^{d}_{+}$ and $\alpha\in(0,1)$, we have that

$\displaystyle\begin{split}&\|\nabla
u\|_{C^{0,\alpha}(B_{\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+})}\\\
&\leq\mathcal{C}(a,x_{0})\bigg{[}\varepsilon^{-\alpha}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla
u|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}+\|g\|_{C^{0,\alpha}(B_{2\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+})}+\|f\|_{C^{0,\alpha}(B_{2\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+})}\bigg{]},\end{split}$
(4.5)

where $\mathcal{C}(a,x)$ has stretched exponential moments –in the sense of
(2.3)– which depend only on $d,\lambda,\alpha,\nu$.

###  4.1 Proof of Proposition 1

Here is the Proof of Lemma 7:

###### Proof of Lemma 7.

Let $1\leq i\leq d$ and assume w.l.o.g. that $\varepsilon=1$. Our argument
comes in three steps: First, in Step 1, we initialize the Campanato iteration.
Then, in Step 2, we show the requisite excess decay. In Step 3, we conclude
our argument by iterating the decay estimate obtained in the previous step.

Throughout this argument we use $\mathcal{C}(a,x)$ to denote a generic random
field on $\Omega\times\mathbb{H}^{d}_{+}$ that is stationary with respect to
tangential shifts and has stretched exponential moments –by tracking these
random fields in the proof below, one finds that all of the random-fields that
appear are polynomial in the random-field from Theorem 4.

Step 1: Initiation of the iteration. We show that $\langle\cdot\rangle$-a.s.
there exists a large enough radius $r>0$ such that

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\leq\frac{1}{r^{\frac{2}{3}}}.$
(4.6)

To show this, we introduce the event that (4.6) fails for all $r>0$ and denote
it by $\mathcal{B}$. We then have that

$\displaystyle\mathbb{P}\big{[}\mathcal{B}\big{]}\leq\mathbb{P}\Big{[}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\geq\frac{1}{r^{\frac{2}{3}}}\Big{]}$
(4.7)

for any fixed $r>0$. To estimate the right-hand probability, we notice that by
the stationarity of $\nabla\theta_{i}$ tangentially to
$\partial\mathbb{H}^{d}_{+}$ we have

$\displaystyle\left<\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\right>\lesssim\frac{1}{r}\left<\int_{0}^{\infty}|\nabla\theta_{i}(x_{0}+z)|^{2}\,{\textrm{d}}z\right>\lesssim\frac{1}{r}.$
(4.8)

We then apply Chebyshev’s inequality to obtain

$\displaystyle\mathbb{P}\Big{[}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\geq\frac{1}{r^{\frac{2}{3}}}\Big{]}\lesssim\frac{1}{r^{\frac{2}{3}}}.$
(4.9)

Since this is true for any $r>0$, we find that $\mathbb{P}[\mathcal{B}]=0$.

Step 2: Comparison to constant coefficient problem. We will now show that
there exist constants $\nu>0$ and $C_{q}>0$ such that the bound

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{r}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x\leq\frac{1}{2}\nu^{-1}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x+\mathcal{C}(a,x_{0})\frac{1}{r}+C_{q}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x\bigg{)}^{2}$
(4.10)

holds for any $R\geq r\geq\mathcal{C}(a,x_{0})$ such that $\frac{r}{R}=\nu$
–notice that $C_{q}$ and also the random fields may depend on $\nu$.

We first set up our argument: notice that w.l.o.g. $1\leq r\leq\frac{R}{8}$,
as otherwise (4.10) is trivial. Let $\zeta$ be a vertical cutoff along
$\partial\mathbb{H}^{d}_{+}$ such that $\zeta\equiv 1$ in $\\{x^{\perp}\leq
1\\}$, $\zeta\equiv 0$ in $\\{x^{\perp}\geq 2\\}$, and $|\nabla\zeta|\lesssim
1$; and, for some $\rho\leq R/8$ to be determined, we take $\eta$ to be a
smooth cutoff for $B_{R/2-2\rho}(x_{0})$ in $B_{R/2-\rho}(x_{0})$ such that
$|\nabla\eta|\leq 2/\rho$. We then define the “homogenization error” as

$\displaystyle
w:=\theta_{i}-\zeta\phi_{i}-(v+\eta(\phi_{j}-\theta_{j})\partial_{j}v),$
(4.11)

where $v$ is the solution to the constant coefficient problem

$\displaystyle-\nabla\cdot(\bar{a}\nabla v)$ $\displaystyle=0$
$\displaystyle\text{in }B^{+}_{R/2}(x_{0}),$ (4.12) $\displaystyle v$
$\displaystyle=\theta_{i}-\zeta\phi_{i}$ $\displaystyle\text{on }\partial
B_{R/2}^{+}(x_{0}).$

Notice that $w\equiv 0$ on $\partial B_{R/2}^{+}(x_{0})$ and satisfies the
relation

$\displaystyle-\nabla\cdot(a\nabla w)$
$\displaystyle=-\underbrace{\nabla\cdot(a\nabla\theta_{i})}_{=0}+\nabla\cdot(a\nabla(\phi_{i}\zeta))+\nabla\cdot(a(\phi_{j}-\theta_{j})\nabla(\eta\partial_{j}v))$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\nabla\cdot(a(e_{j}+\nabla\phi_{j}-\nabla\theta_{j})\eta\partial_{j}v)+\nabla\cdot((1-\eta)a\nabla
v)$
$\displaystyle=\nabla\cdot(a\nabla(\phi_{i}\zeta))+\nabla\cdot(a(\phi_{j}-\theta_{j})\nabla(\eta\partial_{j}v))$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+\underbrace{a(e_{j}+\nabla\phi_{j}-\nabla\theta_{j})\cdot\nabla(\eta\partial_{j}v)}_{=\nabla\cdot(\bar{a}\eta\nabla
v)+(\nabla\cdot\sigma_{j}-a\nabla\theta_{j})\cdot\nabla(\eta\partial_{j}v)}+\nabla\cdot((1-\eta)a\nabla
v)$

in $B_{R/2}^{+}(x_{0})$. In particular, we find that $w$ solves

$\displaystyle-\nabla\cdot(a\nabla w)$
$\displaystyle=\nabla\cdot\big{(}a\nabla(\zeta\phi_{i})+(a(\phi_{j}-\theta_{j})-\sigma_{j})\nabla(\eta\partial_{j}v)+(1-\eta)(a-\bar{a})\nabla
v\big{)}$ (4.13) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \,-a\nabla\theta_{j}\cdot\nabla(\eta\partial_{j}v)$
$\displaystyle\text{in }B^{+}_{R/2}(x_{0}),$ $\displaystyle w$
$\displaystyle=0$ $\displaystyle\text{on }\partial B^{+}_{R/2}(x_{0}),$ (4.14)

where $\sigma$ is the whole space flux-corrector.

We now begin with our argument. Notice that

$\displaystyle\theta_{i}-\zeta\phi_{i}-v(x_{0}^{\parallel})-\partial_{j}v(x_{0}^{\parallel})(x_{j}+\phi_{j}-\theta_{j})=\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1}),$

thanks to the boundary data of $v$. Furthermore,
$\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})$
solves the equation

$\displaystyle-\nabla\cdot(a\nabla(\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})))$
$\displaystyle=\nabla\cdot(a\nabla(\zeta\phi_{i}))$ $\displaystyle\text{in
}B^{+}_{R/2}(x_{0}),$
$\displaystyle\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})$
$\displaystyle=0$ $\displaystyle\text{on
}B_{R/2}^{+}(x_{0})\cap\partial\mathbb{H}^{d}_{+}.$

Emulating the classical proof for the Caccioppoli estimate, we find that

$\displaystyle\begin{split}&\int_{B_{r}^{+}(x_{0})}\big{|}\nabla\theta_{i}-\partial_{1}v(x_{0}^{\parallel})(e_{1}+\nabla\phi_{1}-\nabla\theta_{1})\big{|}^{2}\,\mathrm{d}x\\\
&\lesssim\int_{B_{2r}^{+}(x_{0})}|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x+\frac{1}{r^{2}}\int_{B_{2r}^{+}(x_{0})}\big{|}\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})-c\big{|}^{2}\,\mathrm{d}x\end{split}$
(4.15)

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

Towards showing (4.10), we use the definition of $w$ (see (4.11)) and the
Poincaré-Wirtinger inequality to obtain using also that $R\geq 8r$

$\displaystyle\frac{1}{r^{2}}\int_{B_{2r}^{+}(x_{0})}\big{|}\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})-c\big{|}^{2}\,\mathrm{d}x$
$\displaystyle\lesssim\int_{B^{+}_{R/2}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{R}\bigg{)}^{\kappa}|\nabla
w|^{2}\,\mathrm{d}x+r^{d-2}\sup_{x\in
B^{+}_{2r}(x_{0})}\big{|}v(x)-(v(x_{0}^{\parallel})+\partial_{j}v(x_{0}^{\parallel})x_{j})\big{|}^{2}$
$\displaystyle\,\,\,\,\,+r^{-2}\sup_{x\in B^{+}_{2r}(x_{0})}\big{|}\nabla
v(x)-\nabla
v(x_{0}^{\parallel})\big{|}^{2}\int_{B_{2r}^{+}(x_{0})}|\phi|^{2}+|\theta|^{2}\,\mathrm{d}x,$

for the choice $c=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2r}^{+}(x_{0})}w\,\mathrm{d}x$; note that in
this computation we have also used the fact that the weight in the second line
satisfies $(1-|x-x_{0}|/R)\geq\frac{1}{2}$ in $B_{R/2}^{+}(x_{0})$. This
implies

$\displaystyle\frac{1}{r^{2}}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2r}^{+}(x_{0})}\big{|}\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})-c\big{|}^{2}\,\mathrm{d}x$
$\displaystyle\lesssim
r^{-d}\int_{B^{+}_{R/2}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{R}\bigg{)}^{\kappa}|\nabla
w|^{2}\,\mathrm{d}x+\sup_{B_{2r}^{+}(x_{0})}|\nabla^{2}v|^{2}\bigg{(}r^{2}+\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2r}^{+}(x_{0})}|\phi|^{2}+|\theta|^{2}\,\mathrm{d}x\bigg{)}.$

Applying (4.3) from Lemma 9 to the function $w$ in $B^{+}_{R/2}(x_{0})$ and
using the equation (4.13), we obtain

$\displaystyle\begin{split}&\frac{1}{r^{2}}\int_{B_{2r}^{+}(x_{0})}\big{|}\theta_{i}-\zeta\phi_{i}-\partial_{1}v(x_{0}^{\parallel})(x_{1}+\phi_{1}-\theta_{1})-c\big{|}^{2}\,\mathrm{d}x\\\
&\lesssim\int_{B_{R/2}^{+}(x_{0})}|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\sup_{B^{+}_{R/2-\rho}(x_{0})}(\rho^{-2}|\nabla
v|^{2}+|\nabla^{2}v|^{2})\int_{B_{R/2}^{+}(x_{0})}(|\phi|^{2}+|\theta|^{2}+|\sigma|^{2}+R^{2}|\nabla\theta|^{2})\,\mathrm{d}x\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+r^{d}\sup_{B_{2r}^{+}(x_{0})}|\nabla^{2}v|^{2}\bigg{(}r^{2}+\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2r}^{+}(x_{0})}|\phi|^{2}+|\theta|^{2}\,\mathrm{d}x\bigg{)}+\bigg{(}\frac{\rho}{R}\bigg{)}^{\kappa}\int_{B_{R/2}^{+}(x_{0})}|\nabla
v|^{2}\,\mathrm{d}x.\end{split}$ (4.16)

Plugging the bound (4.16) into (4.15), we obtain:

$\displaystyle\begin{split}&\int_{B_{r}^{+}(x_{0})}\big{|}\nabla\theta_{i}-\partial_{1}v(x_{0}^{\parallel})(e_{1}+\nabla\phi_{1}-\nabla\theta_{1})\big{|}^{2}\,\mathrm{d}x\\\
&\lesssim\int_{B_{2r}^{+}(x_{0})}|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x+\int_{B_{R/2}^{+}(x_{0})}|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+R^{d}\sup_{B^{+}_{R/2-\rho}(x_{0})}(\rho^{-2}|\nabla
v|^{2}+|\nabla^{2}v|^{2})\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R/2}(x_{0})}(|\phi|^{2}+|\theta|^{2}+|\sigma|^{2}+R^{2}|\nabla\theta|^{2})\,\mathrm{d}x\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+R^{d}\sup_{B_{2r}^{+}(x_{0})}|\nabla^{2}v|^{2}\bigg{(}r^{2}+\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2r}^{+}(x_{0})}(|\phi|^{2}+|\theta|^{2})\,\mathrm{d}x\bigg{)}+\bigg{(}\frac{\rho}{R}\bigg{)}^{\kappa}\int_{B_{R/2}^{+}(x_{0})}|\nabla
v|^{2}\,\mathrm{d}x.\end{split}$ (4.17)

We now collect some ingredients: First, recall that

$\displaystyle\sup_{r\geq 1}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\phi_{i}|^{2}\,\mathrm{d}x\leq\mathcal{C}(a,x_{0}),$
(4.18)

which comes from (A.2) of Theorem 4. Using this in conjunction with the
Caccioppoli estimate applied to $\phi_{i}+x_{i}$, we obtain

$\displaystyle\sup_{r\geq 1}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}(x_{0})}|\nabla\phi_{i}|^{2}\,\mathrm{d}x+\sup_{r\geq
1}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\nabla\phi_{i}|^{2}\,\mathrm{d}x\leq\mathcal{C}(a,x_{0}).$
(4.19)

An estimate analogous to (4.18) also holds for $\sigma$. We also use the
constant-coefficient regularity estimate (4.4) of Lemma 10 applied to $v$.
Notice that when $r\leq R/4$, the estimate

$\displaystyle\sup_{B^{+}_{2r}(x_{0})}|\nabla^{2}v|^{2}\lesssim\frac{1}{r^{2}}\bigg{(}\frac{R}{r}\bigg{)}^{d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\nabla v|^{2}\,\mathrm{d}x$
(4.20)

follows from (4.4) by taking $\rho=2r$.

Using these tools and taking averages on both side and choosing $0<\rho\ll r$,
we rewrite (4.17) as

$\displaystyle\begin{split}&\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}\big{|}\nabla\theta_{i}-\partial_{1}v(x_{0}^{\parallel})(e_{1}+\nabla\phi_{1}-\nabla\theta_{1})\big{|}^{2}\,\mathrm{d}x\\\
&\lesssim\bigg{(}\frac{R}{r}\bigg{)}^{d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R/2}^{+}(x_{0})}|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x+r^{-1}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\bigg{(}\frac{R}{r}\bigg{)}^{d}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R}^{+}(x_{0})}|\nabla
v|^{2}\,\mathrm{d}x\bigg{)}\bigg{[}\bigg{(}\frac{R}{\rho}\bigg{)}^{2}\bigg{(}\mathcal{C}(a,x_{0})r^{-2}+R^{-2}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R}^{+}(x_{0})}(|\theta|^{2}+R^{2}|\nabla\theta|^{2})\,\mathrm{d}x\bigg{)}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ +r^{-2}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2r}^{+}(x_{0})}|\theta|^{2}\,\mathrm{d}x+\bigg{(}\frac{\rho}{R}\bigg{)}^{\kappa}\bigg{]}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\bigg{(}\frac{r}{R}\bigg{)}^{2}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R}^{+}(x_{0})}|\nabla
v|^{2}\,\mathrm{d}x.\end{split}$ (4.21)

We continue to process the above expression by using the bound

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R/2}(x_{0})}|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x\lesssim\mathcal{C}(a,x_{0})\frac{1}{R},$
(4.22)

which follows from (4.18), (4.19), and the definition of the cutoff $\zeta$.
Furthermore, since $\theta_{i}-\zeta\phi_{i}$ has homogeneous Dirichlet
boundary data on $\partial\mathbb{H}^{d}_{+}$, the Poincaré inequality and
(4.22) yield that

$\displaystyle\begin{split}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\theta_{i}|^{2}\,\mathrm{d}x&\lesssim\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\theta_{i}-\zeta\phi_{i}|^{2}\,\mathrm{d}x+\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\zeta\phi_{i}|^{2}\,\mathrm{d}x\\\
&\lesssim R^{2}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\nabla(\theta_{i}-\zeta\phi_{i})|^{2}\,\mathrm{d}x+\mathcal{C}(a,x_{0})\frac{1}{R}\\\
&\lesssim R^{2}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{R}(x_{0})}|\nabla\theta_{i}|^{2}\,\mathrm{d}x+\mathcal{C}(a,x_{0})R.\end{split}$
(4.23)

This bound is also valid for $R$ replaced by $2r$.

We then finally notice that $v-(\theta_{i}-\zeta\phi_{i})$ solves the equation

$\displaystyle-\nabla\cdot(\bar{a}\nabla(v-(\theta_{i}-\zeta\phi_{i})))$
$\displaystyle=\nabla\cdot(\bar{a}\nabla(\theta_{i}-\zeta\phi_{i}))$
$\displaystyle\text{in }B^{+}_{R/2}(x_{0}),$ (4.24) $\displaystyle
v-(\theta_{i}-\zeta\phi_{i})$ $\displaystyle=0$ $\displaystyle\text{on
}\partial B_{R/2}^{+}(x_{0}),$

whereby the energy estimate yields

$\displaystyle\int_{B_{R}^{+}(x_{0})}|\nabla
v|^{2}\,\mathrm{d}x\lesssim\int_{B_{R}^{+}(x_{0})}|\nabla\theta_{i}|^{2}+|\nabla(\zeta\phi_{i})|^{2}\,\mathrm{d}x.$
(4.25)

Observe also that by the Caccioppoli inequality for $x_{1}+\phi_{1}$, the
estimate (4.18), and (4.20) (the latter for $r=R/4$) we have for
$r\geq\mathcal{C}(a,x_{0})$

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{r}^{+}(x_{0})}|\partial_{1}v(x_{0})||e_{1}+\nabla\phi_{1}-\nabla\theta_{1}|^{2}\,\mathrm{d}x\lesssim\int_{B_{R}^{+}(x_{0})}|\nabla
v|^{2}\,\mathrm{d}x\Bigg{[}1+\bigg{(}\frac{R}{r}\bigg{)}^{d}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{R}^{+}(x_{0})}|\nabla\theta_{1}|^{2}\,\mathrm{d}x\Bigg{]}.$
(4.26)

Combining (4.21) with (4.22), (4.23), (4.25), and (4.26) yields (4.10) after
choosing $\nu:=r/R$ and $\rho/R\ll\nu$ small enough.

Step 3: Iteration to smaller scales. The result (4.10) of the previous step
implies that

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{\nu^{k+1}R_{0}}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x\leq\frac{1}{2}\nu^{-1}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{\nu^{k}R_{0}}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x+\mathcal{C}(a,x_{0})\frac{1}{r}+C_{q}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{\nu^{k}R_{0}}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x\bigg{)}^{2}$

for any $R_{0}>0$ as long as $\nu^{k+1}R_{0}\geq\mathcal{C}(a,x_{0})$. This
enables us to inductively propagate an estimate of the form

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B^{+}_{\nu^{k}R_{0}}(x_{0})}|\nabla\theta|^{2}\,\mathrm{d}x\leq\frac{4\mathcal{C}(a,x_{0})}{(\nu^{k}R_{0})^{\frac{2}{3}}}$

as long as

$\displaystyle
C_{q}\bigg{(}\frac{4\mathcal{C}(a,x_{0})}{(\nu^{k}R_{0})^{\frac{2}{3}}}\bigg{)}^{2}\leq\frac{\mathcal{C}(a,x_{0})}{(\nu^{k+1}R_{0})^{\frac{2}{3}}},$

or equivalently $\nu^{k}R_{0}\geq(16\mathcal{C}(a,x_{0})C_{q})^{3/2}\nu$. The
start of the induction is provided by (4.6). This directly entails our desired
result. ∎

With the result of Lemma 7 in hand, the argument for Proposition 1 is now a
simple matter of applying the regularity statement of Proposition 24.

###### Proof of Proposition 1.

Without loss of generality we may assume $\varepsilon=1$. Let
$x_{0}\in\mathbb{H}^{d}_{+}$ and $i=1,\ldots,d$. For $x_{0}^{\perp}\geq 1$ the
statement of Proposition 1 is an immediate consequence of (4.1) applied for
$r=2x_{0}^{\perp}$ and the regularity estimate (B.3) applied with
$r:=\frac{1}{2}x_{0}^{\perp}$.

To see the statement in case $x_{0}^{\perp}\leq 1$, let $\zeta$ be a cutoff
with $\zeta=1$ on $\partial\mathbb{H}^{d}_{+}$ and $\zeta\equiv 0$ outside of
$\\{x:\operatorname{dist}(x,\partial\mathbb{H}^{d}_{+})<\tfrac{1}{2}\\}$.
Applying (B.3) to $\theta_{i}-\zeta\phi_{i}$ solving
$-\nabla\cdot(a\nabla(\theta_{i}-\zeta\phi_{i}))=\nabla\cdot(a\nabla(\zeta\phi_{i}))$,
we find by Lemma 23 that

$\displaystyle\begin{split}\sup_{x\in
B_{1/2}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla\theta_{i}(x)|&\leq\|\nabla\theta_{i}\|_{C^{0,\alpha}(B_{1/2}(x_{0})\cap\mathbb{H}^{d}_{+})}\\\
&\lesssim_{d,\lambda,\alpha,\nu}\mathcal{C}(a,x_{0})\Bigg{[}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{1}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla(\theta_{i}-\zeta\phi_{i})|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}+\mathbbm{1}_{x^{\perp}\leq
2}\|a\nabla\phi_{i}\|_{C^{0,\alpha}(B_{1}(x_{0})\cap\mathbb{H}^{d}_{+})}\Bigg{]}.\end{split}$
(4.27)

To treat the first term on the right-hand side of (4.27) we notice that by
Lemma 23 and Theorem 4, we have

$\displaystyle\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{1}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla(\theta_{i}-\zeta\phi_{i})|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}\lesssim_{d,\lambda}\bar{r}(a,x_{0})^{\frac{d}{2}}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\bar{r}(a,x_{0})}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}+\mathcal{C}(a,x_{0})$
(4.28)

where we take
$\bar{r}(a,x_{0}):=\max(r^{*}(a,x_{0}),\tilde{C}(a,x^{\parallel}_{0}),1)$ with
$r^{*}_{\mathcal{O}}(a,x_{0})$ taken from Proposition 24 and with
$\tilde{C}(a,x^{\parallel}_{0})$ defined in the statement of Lemma 7 –both
have stretched exponential moments, which guarantees that $\bar{r}(a,x_{0})$
does as well.

To finish the argument we consider two cases: First, let
$x_{0}^{\perp}\leq\bar{r}(a,x_{0})$. Using (4.1) of Lemma 7, we obtain that

$\displaystyle\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{1}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla(\theta_{i}-\zeta\phi_{i})|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}$
$\displaystyle\leq\bar{r}(a,x_{0})^{\frac{d}{2}}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{2\bar{r}(a,x_{0})}(x^{\parallel}_{0})\cap\mathbb{H}^{d}_{+}}|\nabla\theta_{i}|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}$
$\displaystyle\leq\mathcal{C}(a,x_{0}^{\parallel})\bar{r}(a,x_{0})^{\frac{d}{2}}(1+2\bar{r}(a,x_{0}))^{-\frac{1}{3}}\leq\mathcal{C}(a,x_{0}^{\parallel})\bar{r}(a,x_{0})^{\frac{d}{2}}(1+x_{0}^{\perp})^{-\frac{1}{3}}.$

To control the second term on the right-hand side of (4.27), we use assumption
(A4) to write

$\displaystyle\|a\nabla\phi_{i}\|_{C^{0,\alpha}(B_{2}(x_{0})\cap\mathbb{H}^{d}_{+})}\lesssim_{\alpha,\nu}\mathcal{C}(a,x_{0})\|\nabla\phi_{i}\|_{C^{0,\alpha}(B_{2}(x_{0}))}\lesssim_{d,\lambda,\alpha,\nu}\mathcal{C}(a,x_{0}).$

This concludes the proof. ∎

###  4.2 Proof of auxiliary lemmas for Proposition 1

###### Proof of Lemma 9.

Letting $\eta(x)=\big{(}1-\frac{|x-x_{0}|}{r}\big{)}^{\frac{\kappa}{2}}$, we
test the equation for $u$ with $\eta^{2}u$. Using the uniform ellipticity of
$a$, we obtain that

$\displaystyle\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa}|\nabla
u|^{2}\,\mathrm{d}x$
$\displaystyle\lesssim\frac{\kappa^{2}}{r^{2}}\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa-2}u^{2}\,\mathrm{d}x+\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa}(|g|^{2}+r^{2}f^{2})\,\mathrm{d}x.$

When $x_{0}^{\perp}\geq r$, then the first term on the right-hand side may be
absorbed thanks to the Hardy inequality in [25, Theorem 1.6]. We use that
$x_{0}\in\partial\mathbb{H}^{d}_{+}$ to write the domain $B_{r}^{+}(x_{0})$ in
terms of $1$-dimensional filaments that originate at $x_{0}$ and terminate on
the round part of $\partial B_{r}^{+}$. On each of the rays we apply the
$1$-dimensional Hardy’s inequality:

$\int_{0}^{r^{\prime}}\rho^{\kappa-2}h^{2}(\rho)\,\mathrm{d}\rho\lesssim\frac{1}{(1-\kappa)^{2}}\int_{0}^{r^{\prime}}\rho^{\kappa}(h^{\prime}(\rho))^{2}\,\mathrm{d}\rho,$

for $r^{\prime}>0$, which holds under the assumption $h(0)=0$, to the effect
of

$\int_{B^{+}_{r}(x_{0})\setminus
B_{r/3}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa-2}u^{2}\,\mathrm{d}x\lesssim\frac{r^{2}}{(1-\kappa)^{2}}\int_{B^{+}_{r}(x_{0})\setminus
B_{r/3}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa}|\nabla
u|^{2}\,\mathrm{d}x.$

Notice that we have removed the ball $B_{r/3}$ from the domain of integration
on the left-hand side. To finish the argument we use the Poincaré-Wirtinger
type inequality

$\displaystyle\int_{B^{+}_{2r/3}(x_{0})}u^{2}\,\mathrm{d}x\lesssim
r^{2}\int_{B^{+}_{r}(x_{0})}|\nabla
u|^{2}\,\mathrm{d}x+\int_{B^{+}_{2r/3}(x_{0})\setminus
B_{r/3}(x_{0})}u^{2}\,\mathrm{d}x,$

combined with the fact that for $|x-x_{0}|<2r/3$ the weight
$(1-|x-x_{0}|/r)^{\kappa}$ has values between $1/2$ and $1$ (for $\kappa$
small enough), to obtain

$\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa-2}u^{2}\,\mathrm{d}x\lesssim
r^{2}\int_{B^{+}_{r}(x_{0})}\bigg{(}1-\frac{|x-x_{0}|}{r}\bigg{)}^{\kappa}|\nabla
u|^{2}\,\mathrm{d}x.$

∎

###### Proof of Lemma 11.

By standard Schauder theory, see e. g. [16, Theorem 5.19], we have that

$\displaystyle\|\nabla
u\|_{C^{0,\alpha}(B_{\frac{1}{2}}(x_{0})\cap\mathbb{H}^{d}_{+})}$
$\displaystyle\leq
C\bigg{(}d,\lambda,\|a\|_{C^{0,\alpha}(B_{1}(x_{0})\cap\mathbb{H}^{d}_{+})}\bigg{)}\Bigg{[}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{1}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla
u|^{2}\,\bigg{)}^{\frac{1}{2}}+\|g\|_{C^{0,\alpha}(B_{1}(x_{0})\cap\mathbb{H}^{d}_{+})}+\|f\|_{C^{0,\alpha}(B_{1}(x_{0})\cap\mathbb{H}^{d}_{+})}\Bigg{]},$

which, via a simple scaling argument, becomes

$\displaystyle\begin{split}&\|\nabla
u\|_{C^{0,\alpha}(B_{\frac{\delta}{2}}(x_{0})\cap\mathbb{H}^{d}_{+})}\leq
C\bigg{(}d,\lambda,\delta^{\alpha}\|a\|_{C^{0,\alpha}(B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+})}\bigg{)}\\\
&\hskip
99.58464pt\times\Bigg{[}\delta^{-\alpha}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla
u|^{2}\,\bigg{)}^{\frac{1}{2}}+\|g\|_{C^{0,\alpha}(B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+})}+\|f\|_{C^{0,\alpha}(B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+})}\Bigg{]},\end{split}$
(4.29)

for any $\delta>0$. To ensure that the constant in (4.29) is uniform in
$x_{0}$, we use the assumption (A4) and set
$\delta=\mathcal{C}(a,x_{0})^{-\frac{1}{\alpha}}\varepsilon$ so that

$\displaystyle\delta^{\alpha}\|a\|_{C^{0,\alpha}(B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+})}\leq
1.$

We then use a covering argument and that $\mathcal{C}(a,x)$ in (A4) is locally
uniform in $x$ to write

$\displaystyle\|\nabla
u\|_{C^{0,\alpha}(B_{\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+})}$
$\displaystyle\leq\sup_{y\in B_{\varepsilon}(x_{0})}\|\nabla
u\|_{C^{0,\alpha}(B_{\frac{\delta}{2}}(y)\cap\mathbb{H}^{d}_{+})}$
$\displaystyle\leq C(d,\lambda)\sup_{y\in
B_{\varepsilon}(x_{0})}\Bigg{[}\delta^{-\alpha}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla
u|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}+\|g\|_{C^{0,\alpha}(B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+})}+\|f\|_{C^{0,\alpha}(B_{\delta}(x_{0})\cap\mathbb{H}^{d}_{+})}\Bigg{]}$
$\displaystyle\leq
C(d,\lambda)\Bigg{[}\varepsilon^{-\alpha}\mathcal{C}^{d/\alpha}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+}}|\nabla
u|^{2}\,\mathrm{d}x\bigg{)}^{\frac{1}{2}}+\|g\|_{C^{0,\alpha}(B_{2\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+})}+\|f\|_{C^{0,\alpha}(B_{2\varepsilon}(x_{0})\cap\mathbb{H}^{d}_{+})}\Bigg{]},$

which completes our argument. ∎

##  5 Argument for Proposition 5: Fluctuation estimates for the boundary layer
corrector

In our proof of Proposition 5 we use that the spectral gap inequality (A3)
entails a corresponding estimate on arbitrary polynomial moments.

######  Lemma 12.

Let $\langle\cdot\rangle$ be an ensemble satisfying the assumption (A3). Then
for any random variable $F$ and $q\geq 1$ the estimate

$\displaystyle\left<(F-\left<F\right>)^{2q}\right>^{\frac{1}{2q}}\leq
Cq^{\frac{1}{2}}\varepsilon^{\frac{d}{2}}\left<\bigg{(}\int_{\mathbb{R}^{d}}\bigg{(}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.25003pt}}\\!\int_{B_{\varepsilon}(x)}\Big{|}\frac{\partial
F}{\partial
a}\Big{|}\,{\textrm{d}}y\bigg{)}^{2}\,\mathrm{d}x\bigg{)}^{q}\right>^{\frac{1}{2q}}$
(5.1)

holds.

This standard lemma has been shown e. g. in [10, Prop. 3.1]. Applying (5.1) to
the random variable $F$ defined in (3.12) and deriving an expression for the
“vertical derivatives” on the right-hand side –we find that obtaining (3.13)
reduces to estimating certain auxiliary functions that are used to describe
$\frac{\partial F}{\partial a}$.

We organize the argument into three parts: First, in Section 5.1, we derive an
expression for $\frac{\partial F}{\partial a}$ in terms of auxiliary functions
$h$, $v_{1}$, and $v_{2}$ –each solving an elliptic equation with coefficient
field $a^{*}$ on either ${\mathcal{O}}$ with homogeneous Dirichlet boundary
data or on $\mathbb{R}^{d}$ with a decay condition. In Section 5.2, we prove a
weighted Meyers estimate (Lemma 14) –this follows from the large-scale
regularity theory for $a$-harmonic functions with homogeneous Dirichlet
boundary data on ${\mathcal{O}}$, which is obtained in Proposition 24. The
result of Lemma 14, which has become a standard tool in stochastic
homogenization, yields estimates for the auxiliary functions $h$, $v_{1}$, and
$v_{2}$. In Section 5.3, we finish our argument by combining these estimates
with Lemma 12 and the expression for $\frac{\partial F}{\partial a}$ that has
been derived in Section 5.1.

###  5.1 Identification of the Malliavin / Fréchet derivative

In this section we calculate $\frac{\partial F}{\partial a}$ in terms of
auxiliary functions $h$, $v_{1}$, and $v_{2}$.

######  Lemma 13.

We adopt the assumptions of Proposition 5. Let $h\in
H^{1}_{\rm{loc}}({\mathcal{O}})$ solve

$\displaystyle-\nabla\cdot(a^{*}\nabla h)$ $\displaystyle=\nabla\cdot g$
$\displaystyle\text{in }{\mathcal{O}},$ (5.2) $\displaystyle h$
$\displaystyle=0$ $\displaystyle\text{on }\partial{\mathcal{O}},$

and $\hat{h}\in H_{\rm{loc}}^{1}(\mathbb{R}^{d})$ be the decaying solution of

$\displaystyle-\nabla\cdot(a^{*}\nabla\hat{h})=\nabla\cdot g\quad\text{on
}\mathbb{R}^{d}.$ (5.3)

When ${\mathcal{O}}=\mathbb{H}^{d}_{+}$, $h$ is the decaying solution of
(5.2); when ${\mathcal{O}}$ is a bounded $C^{1,1}$ domain then $h\in
H^{1}({\mathcal{O}})$. We, furthermore, define $v_{1}$ and $v_{2}\in
H^{1}(\mathbb{R}^{d})$ as the decaying solutions of

$\displaystyle-\nabla\cdot a^{*}\nabla v_{1}$ $\displaystyle=\nabla f\cdot
a^{*}(\chi_{\mathcal{O}}\nabla h-\nabla\hat{h})\quad\textrm{in
}\mathbb{R}^{d},$ (5.4)

and

$\displaystyle-\nabla\cdot a^{*}\nabla v_{2}$
$\displaystyle=-\nabla\cdot(fa^{*}(\chi_{\mathcal{O}}\nabla
h-\nabla\hat{h}))\quad\textrm{in }\mathbb{R}^{d},$ (5.5)

respectively. Then, for $F$ defined in (3.12), the vertical derivative
$\frac{\partial F}{\partial a}$ is given by

$\displaystyle\frac{\partial F}{\partial a}=\chi_{\mathcal{O}}\nabla
u\otimes\nabla h-(e_{i}+\nabla\phi_{i})\otimes(\nabla v_{1}+\nabla v_{2}).$
(5.6)

###### Proof.

Calculating the sensitivity of $u$ solving (3.11) to compactly supported
perturbations of $a$, we find that

$\displaystyle-\nabla\cdot(a\nabla\delta u)$ $\displaystyle=\nabla\cdot(\delta
a\nabla u)$ $\displaystyle\textrm{ in }{\mathcal{O}},$ (5.7)
$\displaystyle\delta u$ $\displaystyle=f\delta\phi_{i}$ $\displaystyle\textrm{
on }\partial{\mathcal{O}}.$

Similarly, for $\phi_{i}$ we obtain

$\displaystyle-\nabla\cdot(a\nabla\delta\phi_{i})=\nabla\cdot(\delta
a(e_{i}+\nabla\phi_{i}))\quad\text{on }\mathbb{R}^{d}.$ (5.8)

The relations (5.2) and (5.7) then allow us to write

$\displaystyle\delta F$ $\displaystyle=\int_{\mathcal{O}}g\cdot\nabla\delta
u\,\mathrm{d}x$ $\displaystyle=-\int_{\mathcal{O}}\nabla\delta u\cdot
a^{*}\nabla h\,\mathrm{d}x+\int_{\partial{\mathcal{O}}}\delta
u{\mathbf{n}}\cdot a^{*}\nabla h\,\mathrm{d}S$
$\displaystyle=\int_{\mathcal{O}}\delta a\nabla u\cdot\nabla
h\,\mathrm{d}x+\int_{\partial{\mathcal{O}}}f\delta\phi_{i}{\mathbf{n}}\cdot
a^{*}\nabla h\,\mathrm{d}S.$

Again using the equation (5.2) and now also that $g$ is compactly supported in
${\mathcal{O}}$, the second term on the right-hand side may be rewritten as

$\displaystyle\int_{\partial{\mathcal{O}}}f\delta\phi_{i}{\mathbf{n}}\cdot
a^{*}\nabla h\,\mathrm{d}S$
$\displaystyle=\int_{\mathcal{O}}\nabla(f\delta\phi_{i})\cdot(a^{*}\nabla
h+g)\,\mathrm{d}x.$

By (5.3) we see that

$\displaystyle\int_{\mathcal{O}}\nabla(f\delta\phi_{i})\cdot
g\,\mathrm{d}x=-\int_{\mathbb{R}^{d}}\nabla(f\delta\phi_{i})\cdot
a^{*}\nabla\hat{h}\,\mathrm{d}x,$

which then leads to

$\displaystyle\int_{\partial{\mathcal{O}}}f\delta\phi_{i}{\mathbf{n}}\cdot
a^{*}\nabla h\,\mathrm{d}S=$
$\displaystyle\int_{\mathcal{O}}\nabla(f\delta\phi_{i})\cdot a^{*}(\nabla
h-\nabla\hat{h})\,\mathrm{d}x-\int_{\mathbb{R}^{d}\setminus{\mathcal{O}}}\nabla(f\delta\phi_{i})\cdot
a^{*}\nabla\hat{h}\,\mathrm{d}x.$

Using the equations (5.4), (5.5) and (5.8), we rewrite the right-hand side as

$\displaystyle\int_{\mathcal{O}}\nabla(f\delta\phi_{i})\cdot a^{*}(\nabla
h-\nabla\hat{h})\,\mathrm{d}x-\int_{\mathbb{R}^{d}\setminus{\mathcal{O}}}\nabla(f\delta\phi_{i})\cdot
a^{*}\nabla\hat{h}\,\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\delta\phi_{i}\nabla f\cdot
a^{*}(\chi_{\mathcal{O}}\nabla
h-\nabla\hat{h})\,\mathrm{d}x+\int_{\mathbb{R}^{d}}f\nabla\delta\phi_{i}\cdot
a^{*}(\chi_{\mathcal{O}}\nabla h-\nabla\hat{h})\,\mathrm{d}x$
$\displaystyle=\int_{\mathbb{R}^{d}}\nabla\delta\phi_{i}\cdot a^{*}(\nabla
v_{1}+\nabla v_{2})\,\mathrm{d}x$
$\displaystyle=-\int_{\mathbb{R}^{d}}\bigg{(}\delta
a(\nabla\phi_{i}+e_{i})\cdot(\nabla v_{1}+\nabla v_{2})\bigg{)}\,\mathrm{d}x.$

Combining the above identities, we obtain

$\displaystyle\delta F=\int_{\mathcal{O}}\delta a\nabla u\cdot\nabla
h\,\mathrm{d}x-\int_{\mathbb{R}^{d}}\delta a(\nabla\phi_{i}+e_{i})\cdot(\nabla
v_{1}+\nabla v_{2})\,\mathrm{d}x,$

which concludes our argument. ∎

###  5.2 Estimates on the auxiliary functions: A weighted Meyers estimate

Having calculated the derivative of $F$ with respect to changes in the
coefficient field $a$ in the previous step, we now seek to control the terms
on the right-hand side of (5.6) by obtaining suitable estimates for
$v_{1},v_{2}$, and $h$. For this, the main tool is the following weighted
Meyers estimate:

######  Lemma 14 (Weighted Meyers Estimate).

Let $\langle\cdot\rangle$ be an ensemble of coefficient fields on
$\mathbb{R}^{d}$ that satisfies the assumptions (A1)-(A4). Furthermore, let
$x_{0}\in{\mathcal{O}}$ –where ${\mathcal{O}}$ is a bounded $C^{1,1}$-domain,
$\mathbb{H}^{d}_{+}$, or $\mathbb{R}^{d}$– and $R\geq
r^{{\mathcal{O}}}_{*}(x_{0})$ (defined in Proposition 24). For $f_{0}\in
L^{2}({\mathcal{O}})$ and $F_{0}\in L^{2}({\mathcal{O}},\mathbb{R}^{d})$, let
$v$ denote the unique solution of

$\displaystyle-\nabla\cdot(a^{*}\nabla v)$ $\displaystyle=\nabla\cdot
F_{0}+f_{0}\quad$ $\displaystyle\text{in }{\mathcal{O}},$ (5.9) $\displaystyle
v$ $\displaystyle=0\quad$ $\displaystyle\text{on }\partial{\mathcal{O}}$

and introduce the weight

$\omega_{\alpha,R}(x):=\bigg{(}\frac{|x-x_{0}|}{R}+1\bigg{)}^{\alpha}.$ (5.10)

Then, in the case that $f_{0}=0$, there exists
$\bar{p}(\lambda,d,{\mathcal{O}})>1$ such that, for $1\leq p<\bar{p}$ and
$0\leq\alpha_{0}<\alpha_{1}<d(2p-1)$, the estimate

$\displaystyle\begin{split}&\bigg{(}\int_{{\mathcal{O}}}|\nabla
v|^{2p}\omega_{\alpha_{0},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\leq
C(d,\lambda,p,{\mathcal{O}},\alpha_{0},\alpha_{1})\bigg{(}\int_{{\mathcal{O}}}|F_{0}|^{2p}\omega_{\alpha_{1},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\end{split}$
(5.11)

holds. Likewise, in the case that $F_{0}=0$, there exists
$\bar{p}(\lambda,d,{\mathcal{O}})>1$ such that, for $1\leq p<\bar{p}$ and
$0\leq\alpha_{0}<\alpha_{1}-2p<\alpha_{1}<d(2p-1)$, the estimate

$\displaystyle\begin{split}&\bigg{(}\int_{{\mathcal{O}}}|\nabla
v|^{2p}\omega_{\alpha_{0},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\leq
C(d,\lambda,p,{\mathcal{O}},\alpha_{0},\alpha_{1})R\bigg{(}\int_{{\mathcal{O}}}|f_{0}|^{2p}\omega_{\alpha_{1},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\end{split}$
(5.12)

holds. We emphasize that when ${\mathcal{O}}=\mathbb{R}^{d}$ or
$\mathbb{H}^{d}_{+}$, the boundary condition in (5.9) is augmented by a decay
condition in the far-field of the infinite part of the domain.

In the case that ${\mathcal{O}}=\mathbb{R}^{d}$, the statement has already
been shown in [6, Lemma 7]. For completeness, we give a proof of Lemma 14 in
Appendix C –we remark that it depends on the use of the classical Meyers
estimate on the contributions of a dyadic decomposition of ${\mathcal{O}}$.
Since $v$ solving (5.9) does not have Dirichlet boundary conditions on the
various members of the decomposition, when we apply the classical Meyers
estimate on each contribution, we pick up an additional term on the right-hand
side (see e. g. (C.8)), the sum of which we may handle by employing the large-
scale regularity theory for $a$-harmonic functions on ${\mathcal{O}}$.

An immediate consequence of Lemma 14 is the following estimate:

######  Corollary 15.

Adopt the assumptions of Proposition 5. Then, for $h$ and $\hat{h}$ solving
(5.2) and (5.3) respectively, there exists $\bar{p}(d,\lambda)>1$ such that,
for $1\leq p<\bar{p}$ and $0\leq\alpha<d(2p-1)$, the estimates

$\displaystyle\bigg{(}\int_{{\mathcal{O}}}|\nabla
h|^{2p}\omega_{\alpha,r+r^{{\mathcal{O}}}_{*}(x_{0})}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}$
$\displaystyle\leq
C(d,\lambda,p,{\mathcal{O}},\alpha)\bigg{(}\int_{B_{r}(x_{0})}|g|^{2p}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}},$
(5.13a)
$\displaystyle\bigg{(}\int_{\mathbb{R}^{d}}|\nabla\hat{h}|^{2p}\omega_{\alpha,r+r^{{\mathcal{O}}}_{*}(x_{0})}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}$
$\displaystyle\leq
C(d,\lambda,p,{\mathcal{O}},\alpha)\bigg{(}\int_{B_{r}(x_{0})}|g|^{2p}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}$
(5.13b)

hold.

Since this corollary follows easily from Lemma 14, we forego a formal proof
–notice that to obtain the result one simply applies Lemma 14 to $h$ and
$\hat{h}$ and uses the definition of the weight (5.10), along with the
property that $g$ is supported in $B_{r}(x_{0})$.

To obtain satisfactory control of $\nabla v_{1}$ and $\nabla v_{2}$ (solving
(5.4) and (5.5) respectively), we require a better estimate for
$\nabla(h-\hat{h})$ than that following from the triangle inequality and
Corollary 15. Towards this end, we prove the following lemma:

######  Lemma 16.

Adopt the assumptions of Proposition 5, and let
$R:=\operatorname{dist}(x,\partial{\mathcal{O}})+r+r^{{\mathcal{O}}}_{*}(x_{0})$
and $h$ and $\hat{h}$ solve (5.2) and (5.3) respectively. Then there exists
$\bar{p}>1$ such that, for all $1\leq p<\bar{p}$ and
$0\leq\alpha_{0}<\alpha_{1}<d(2p-1)$ additionally satisfying $\alpha_{1}>2p$,
the estimate

$\displaystyle\begin{split}&\bigg{(}\int_{{\mathcal{O}}}|\nabla(h-\hat{h})|^{2p}\omega_{\alpha_{0},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\leq
C(d,\lambda,p,{\mathcal{O}},\alpha_{0},\alpha_{1})\bigg{(}\frac{\operatorname{dist}(x_{0},\partial{\mathcal{O}})}{r+r^{{\mathcal{O}}}_{*}(x_{0})}+1\bigg{)}^{-\frac{\alpha_{1}}{2p}}\bigg{(}\int_{B_{r}(x_{0})}|g|^{2p}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\end{split}$
(5.14)

holds.

We make use of Lemma 16, which is proven at the end of this section, in the
form of the next corollary. In particular, combining the estimate (5.14) on
$\nabla(h-\hat{h})$ inside ${\mathcal{O}}$ with (5.13b) for $\nabla\hat{h}$
outside of ${\mathcal{O}}$, we obtain:

######  Corollary 17.

Adopt the assumptions of Proposition 5, and let
$R:=\operatorname{dist}(x,\partial{\mathcal{O}})+r+r^{{\mathcal{O}}}_{*}(x_{0})$
and $v_{1}$ and $v_{2}$ solve (5.4) and (5.5) respectively. Then, for $1\leq
p<\bar{p}$ and $0\leq\alpha_{0}<\alpha_{1}-2p<\alpha_{1}<d(2p-1)$, the
estimate

$\displaystyle\bigg{(}\int_{\mathbb{R}^{d}}|\nabla
v_{i}|^{2p}\omega_{\alpha_{0},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\leq
C(d,\lambda,p,{\mathcal{O}},\alpha_{0},\alpha_{1},f)\bigg{(}\frac{\operatorname{dist}(x_{0},\partial{\mathcal{O}})}{r+r^{{\mathcal{O}}}_{*}(x_{0})}+1\bigg{)}^{\frac{-\alpha_{1}}{2p}}\bigg{(}\int_{B_{r}(x_{0})}|g|^{2p}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}$
(5.15)

holds for $i=1,2$. In fact, in the case that $i=2$, (5.15) holds for
$0\leq\alpha_{0}<\alpha_{1}<d(2p-1)$.

###### Proof.

We may assume that $R\leq\operatorname{diam}({\mathcal{O}})$, since otherwise
(5.15) follows from the energy estimate and the Meyers estimate applied to
(5.4) and (5.5), for $i=1,2$ respectively.

We then first consider the case that $i=1$. To obtain (5.15), we apply the
whole space version of the weighted Meyers estimate (5.12) from Lemma 14 to
$v_{1}$ solving (5.4). This yields that

$\displaystyle\bigg{(}\int_{\mathbb{R}^{d}}|\nabla
v_{1}|^{2p}\omega_{\alpha_{0},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}\lesssim
R\bigg{(}\int_{\mathbb{R}^{d}}|\nabla f\cdot a^{*}(\chi_{\mathcal{O}}\nabla
h-\nabla\hat{h}))|^{2p}\omega_{\alpha_{1},R}\,\mathrm{d}x\bigg{)}^{\frac{1}{2p}}$

for $1<p\leq\bar{p}$ and $0<\alpha_{0}<\alpha_{1}-2p<\alpha_{1}<d(2p-1)$. In
the case that ${\mathcal{O}}=\mathbb{H}^{d}_{+}$, the right-hand side vanishes
since $f\equiv 1$. When alternatively ${\mathcal{O}}$ is a bounded $C^{1,1}$
-domain, (5.15) follows from (5.14), (5.13b), and that $|\nabla f|\leq C$.

The estimate for $i=2$ is obtained in essentially the same way, with the
application of (5.11) to $v_{2}$ solving (5.5). Since we use (5.11) instead of
(5.12), the result (5.15) holds for $0\leq\alpha_{0}<\alpha_{1}<d(2p-1)$. ∎

To complete this section, we now give the proof of Lemma 16. For this we use
Lemma 14 applied to $h^{T}-\eta\hat{h}^{T}$, where the cut-off $\eta$ ensures
homogeneous Dirichlet boundary data on ${\mathcal{O}}$. Here comes the
argument:

###### Proof of Lemma 16.

Let $L=\operatorname{dist}(x_{0},\partial{\mathcal{O}})$. Since (5.14) follows
from (5.13a) and (5.13b) by the triangle inequality when $L\leq
2(r^{{\mathcal{O}}}_{*}(x_{0})+r)$, we may assume that $L\geq
2(r^{{\mathcal{O}}}_{*}(x_{0})+r)$.

The main idea of our argument is to split $h-\hat{h}$ into a function which
vanishes on $\partial{\mathcal{O}}$ and a remainder, and show that (5.14)
holds for both of these functions. To this end, let
$\hat{\eta}:[0,\infty)\to[0,1]$ be a smooth function satisfying
$\hat{\eta}(t)=1$ for $t\in[0,\frac{1}{2}]$, $\hat{\eta}(t)=0$ for
$t\geq\frac{5}{8}$, and $|\hat{\eta}^{\prime}|\leq 9$, and use it to define a
cutoff function $\eta(x):=\hat{\eta}(\frac{|x-x_{0}|}{L})$. Since $L\geq
2(r^{{\mathcal{O}}}_{*}(x_{0})+r)$ and $g$ is supported in $B_{r}(x_{0})$, we
see that $\eta\equiv 1$ on the support of $g$; furthermore, $\eta=0$ in a
neighborhood of the boundary $\partial{\mathcal{O}}$ and
$\operatorname{supp}(\nabla\eta)\subseteq B_{\frac{3L}{4}}(x_{0})\setminus
B_{\frac{L}{2}}(x_{0})$.

Defining $w:=h-\eta\hat{h}$ and using the equations (5.2) and (5.3), we find
that

$\displaystyle-\nabla\cdot(a^{*}\nabla w)$
$\displaystyle=-\nabla\cdot(a^{*}\nabla