# Convergence to the thermodynamic limit for random-field random surfaces

Paul Dario  Université Claude Bernard Lyon 1, Institut Camille Jordan, 69622
Villeurbanne, France<EMAIL_ADDRESS>
###### Abstract.

We study random surfaces with a uniformly convex gradient interaction in the
presence of quenched disorder taking the form of a random independent external
field. Previous work on the model has focused on proving existence and
uniqueness of infinite-volume gradient Gibbs measures with a given tilt and on
studying the fluctuations of the surface and its discrete gradient.

In this work we focus on the convergence of the thermodynamic limit,
establishing convergence of the finite-volume distributions with Dirichlet
boundary conditions to translation-covariant (gradient) Gibbs measures.
Specifically, it is shown that, when the law of the random field has finite
second moment and is symmetric, the distribution of the gradient of the
surface converges in dimensions $d\geq 4$ while the distribution of the
surface itself converges in dimensions $d\geq 5$. Moreover, a power-law upper
bound on the rate of convergence in Wasserstein distance is obtained. The
results partially answer a question discussed by Cotar and Külske [19].

###### Contents

  1. 1 Introduction
  2. 2 Notation and preliminary results
  3. 3 Translation-covariant gradient Gibbs states in dimensions $d\geq 4$
  4. 4 Translation-covariant Gibbs states in dimensions $d\geq 5$
  5. A Gradient of the Neumann Green’s function
  6. B De Giorgi-Nash-Moser regularity

##  1\. Introduction

###  1.1. Description of the model and main results

Random surfaces in statistical mechanics are used to model the interface
separating two pure thermodynamic phases. In classical effective interface
models of this type, the interface is represented by a real-valued field
$\left\\{\phi(x)\,:\,x\in\mathbb{Z}^{d}\right\\}$ to which one associates an
energy determined by an interaction potential applied to the discrete gradient
of the surface. In this article, we consider a model of random interfaces in
the presence of a random disorder taking the form of a random, independent
field. Let us define this model precisely. Given finite set
$\Lambda\subseteq\mathbb{Z}^{d}$, a map $\eta:\Lambda\to\mathbb{R}$ and an
intensity parameter $\lambda>0$, we assign to each surface
$\phi:\Lambda\to\mathbb{R}$ satisfying the Dirichlet boundary condition
$\phi=0$ on $\partial\Lambda$ the energy

(1.1)
$H_{\Lambda}^{\eta}\left(\phi\right):=\sum_{\begin{subarray}{c}x,y\in\Lambda^{+}\\\
|x-y|=1\end{subarray}}V(\phi(x)-\phi(y))-\lambda\sum_{x\in\Lambda}\eta(x)\phi(x).$

where $\partial\Lambda$ is the external vertex boundary of $\Lambda$,
$\Lambda^{+}$ is the set $\Lambda\cup\partial\Lambda$, $|\cdot|$ denotes the
Euclidean norm and $V:\mathbb{R}\to\mathbb{R}$ is an interaction potential
satisfying suitable assumptions. Typically, one considers even $V$ which are
twice-continuously differentiable and uniformly convex, i.e., there exist two
constants $c_{-}$ and $c_{+}$ satisfying

(1.2) $0<c_{-}\leq V^{\prime\prime}(t)\leq c_{+}<\infty.$

The law of random surface is then given by

(1.3)
$\mu_{\Lambda}^{\eta}(\mathrm{d}\phi):=\frac{1}{Z_{\Lambda}^{\eta}}\exp\left(-H_{\Lambda}^{\eta}(\phi)\right)\prod_{x\in\Lambda}d\phi(x),$

where $d\phi(x)$ denotes the Lebesgue measure on $\mathbb{R}$ and
$Z_{\Lambda}^{\eta}$ is the constant ensuring that $\mu_{\Lambda}^{\eta}$ is a
probability measure. We assume that the external field $\eta$ is random, with
the random variables $\left(\eta(x)\right)_{x\in\mathbb{Z}^{d}}$ independent
with symmetric distribution and unit variance, and wish to study the
properties of the random surface distributed according to
$\mu_{\Lambda}^{\eta}$ for a typical realization of the random field $\eta$.
We denote by $\mathbb{P}$ the law of the random field, and by $\mathbb{E}$ the
corresponding expectation.

A natural question pertaining to the model is whether the distribution (1.3)
of the random surface or the one of its discrete gradient converges upon
considering an increasing sequence of domains
$\Lambda_{n}\subseteq\mathbb{Z}^{d}$ tending to $\mathbb{Z}^{d}$ as $n$ tends
to infinity. In the case of the pure system (i.e., when $\eta\equiv 0$),
Funaki and Spohn [29] proved the existence of this limit for the discrete
gradient of the random surface in any dimension $d\geq 1$, as well as the
uniqueness of the infinite-volume ergodic gradient Gibbs measures with a
specific tilt. Regarding the distribution of the surface, it is known that the
limit does not exist in dimensions $d\leq 2$ [13] and, by using the Brascamp-
Lieb concentration inequality [11, 12] (see e.g., [28, Section 9], [24,
Section 6]) that the limit exists in dimensions $d\geq 3$.

These qualitative properties are modified upon the addition of a quenched
random external field to the system. This phenomenon can be observed by
considering the specific case of the random-field Gaussian free field, that is
the model (1.3) with the interaction potential $V(x)=x^{2}/2$. In this
setting, for any fixed realization of the external field $\eta$, the
disordered random surface has a multivariate Gaussian distribution whose mean
vector and covariance matrix are explicit. One can then rely on explicit
computations to prove that, for almost every realization of the disorder, the
law of the discrete gradient of the random surface converges weakly to an
infinite-volume gradient Gibbs measure in dimensions $d\geq 3$ and does not
converge in dimensions $d\leq 2$. Similarly, one can show that, for almost
every realization of the disorder, the law of the random surface converges
weakly to an infinite-volume Gibbs measure in dimensions $d\geq 5$ and does
not converge in dimensions $d\leq 4$ (see e.g., [19, Appendix A.1] where the
non-existence of Gaussian Gibbs measures with finite first moment is
established in dimensions $d=3,4$).

It is natural to expect that these properties persist for the random-field
random surface model (1.3) with a uniformly-convex interaction potential $V$
(for which the distribution of the random surface is not explicitly known).
Indeed, Cotar and Külske [19, 20] proved the existence and uniqueness of
infinite-volume translation-covariant gradient Gibbs measures in $d\geq 3$
with a fixed tilt (the uniqueness is established under the necessary condition
that the corresponding annealed measure is ergodic with respect to spatial
shifts). However, their work left open the question of the convergence of the
finite-volume distributions to the infinite-volume Gibbs measures; see the
discussion in [19, Appendix A]. In this article, we work under the assumption
that the law of the random field is symmetric and establish that, for almost
every realization of the disorder, the distribution of the gradient of the
random surface converges weakly to an infinite-volume gradient Gibbs measure
in dimensions $d\geq 4$, and that the distribution of the random surface
converges to an infinite-volume Gibbs measure in dimensions $d\geq 5$. This
resolves the convergence question in the case of a random field with symmetric
distribution and with uniformly elliptic interaction potential, apart from the
gradient convergence in dimension $d=3$, which is expected but does not follow
from our proof; see the discussion below Theorem 1. Further discussion of the
relation of these results to the previous literature can be found in Section
1.2.

In order to state the results, we introduce a few additional notations and
definitions. Write $\Lambda_{L}:=\left\\{-L,\ldots,L\right\\}^{d}$, let
$\vec{E}(\mathbb{Z}^{d})$ denote the set of oriented edges of the lattice
$\mathbb{Z}^{d}$ (see Section 2.1). For each function
$\phi:\mathbb{Z}^{d}\to\mathbb{R}$, we define its discrete gradient
$\nabla\phi:\vec{E}(\mathbb{Z}^{d})\to\mathbb{R}$ by, for any
$e=(x,y)\in\vec{E}(\mathbb{Z}^{d})$, $\nabla\phi(e):=\phi(y)-\phi(x)$. We
denote by $\Omega:=\mathbb{R}^{\mathbb{Z}^{d}}$ the set of functions defined
on $\mathbb{Z}^{d}$ and valued in $\mathbb{R}$, by
$\mathcal{P}\left(\Omega\right)$ the set of probability measures on the space
$\Omega$ and by $C_{b}\left(\Omega\right)$ the set of functions
$f:\Omega\to\mathbb{R}$ which depend on finitely many coordinates, are
continuous and bounded. For a measure $\mu\in\mathcal{P}\left(\Omega\right)$,
we denote by $\left\langle\cdot\right\rangle_{\mu}$ the expectation with
respect to $\mu$. Given a pair of measures
$\mu_{0},\mu_{1}\in\mathcal{P}\left(\Omega\right)$, we denote by
$\Gamma(\mu_{0},\mu_{1})$ the set of couplings between $\mu_{0}$ and
$\mu_{1}$, that is, the set of probability measures on the product space
$\Omega\times\Omega$ whose first and second marginals are equal to $\mu_{0}$
and $\mu_{1}$ respectively.

We denote by $\mathcal{X}$ the set of _gradient fields_ , that is

(1.4)
$\mathcal{X}:=\left\\{\chi:\vec{E}(\mathbb{Z}^{d})\to\mathbb{R}\,:\,\exists\phi:\mathbb{Z}^{d}\to\mathbb{R},\,\chi=\nabla\phi\right\\}.$

We let $\mathcal{P}\left(\mathcal{X}\right)$ be the set of probability
measures on the space $\mathcal{X}$ and $C_{b}\left(\mathcal{X}\right)$ be the
set of functions $f:\mathcal{X}\to\mathbb{R}$ which depend on finitely many
coordinates, are continuous and bounded. For a measure
$\mu\in\mathcal{P}\left(\mathcal{X}\right)$, we denote by
$\left\langle\cdot\right\rangle_{\mu}$ the expectation with respect to $\mu$.
Given a pair of measures
$\mu_{0},\mu_{1}\in\mathcal{P}\left(\mathcal{X}\right)$, we denote by
$\Gamma_{\nabla}(\mu_{0},\mu_{1})$ the set of couplings between $\mu_{0}$ and
$\mu_{1}$.

For $y\in\mathbb{Z}^{d}$, we denote by $\tau_{y}:\Omega\to\Omega$ the shift
operator, i.e., $\tau_{y}\phi(x)=\phi(x-y)$. We extend the notation to the set
of gradient fields by writing $\tau_{y}\chi(e)=\chi(e-y)$ and to the set of
external fields by writing $\tau_{y}\eta(x)=\eta(x-y)$. For any measure
$\mu\in\mathcal{P}\left(\Omega\right)$ (resp.
$\mathcal{P}\left(\mathcal{X}\right)$), we denote by $(\tau_{y})_{*}\mu$ the
pushforward of the measure $\mu$ by the shift $\tau_{y}$. A sequence of
probability measures $(\mu_{n})_{n\geq 0}\in\mathcal{P}\left(\Omega\right)$
(resp. $\mathcal{P}\left(\mathcal{X}\right)$) converges weakly in
$\mathcal{P}(\Omega)$ (resp. $\mathcal{P}(\mathcal{X})$) to a measure
$\mu\in\mathcal{P}\left(\Omega\right)$ (resp.
$\mathcal{P}\left(\mathcal{X}\right)$) if, for any function $f\in
C_{b}\left(\Omega\right)$ (resp. $f\in C_{b}\left(\mathcal{X}\right)$),

$\left\langle
f\right\rangle_{\mu_{n}}\underset{n\to\infty}{\longrightarrow}\left\langle
f\right\rangle_{\mu}.$

We remark that, for any finite subset $\Lambda\subseteq\mathbb{Z}^{d}$ and any
realization of the random field, one can see the measure
$\mu_{\Lambda}^{\eta}$ as an element of the set
$\mathcal{P}\left(\Omega\right)$ by extending the field $\phi$ by the value
$0$ outside the set $\Lambda$. Similarly, if one wishes to study the discrete
gradient of the random surface, one can see the measure $\mu_{\Lambda}^{\eta}$
as an element of the set $\mathcal{P}\left(\mathcal{X}\right)$ by considering
the discrete gradient of the random surface and extending it to be equal to
$0$ outside the set $\Lambda$. We may specify a boundary condition on the law
of the disordered random surface as follows. For any
$\psi:\mathbb{Z}^{d}\to\mathbb{R}$, $\eta:\mathbb{Z}^{d}\to\mathbb{R}$ and
$\Lambda\subseteq\mathbb{Z}^{d}$ finite, let
$\mu_{\Lambda}^{\eta,\psi}\in\mathcal{P}(\Omega)$ be defined by

$\mu_{\Lambda}^{\eta,\psi}(\mathrm{d}\phi):=\frac{1}{Z_{\Lambda}^{\eta}}\exp\left(-H_{\Lambda}^{\eta}(\phi)\right)\prod_{v\in\Lambda}d\phi(v)\prod_{v\in\mathbb{Z}^{d}\setminus\Lambda}\delta_{\psi(v)}(d\phi(v)),$

where the symbol $\delta_{x}$ refers to the Dirac measure at the point
$x\in\mathbb{R}$.

The first theorem of this article establishes that, for almost every
realization of the random field, the sequence of measures
$(\mu_{\Lambda_{L}}^{\eta})$ converges weakly in
$\mathcal{P}\left(\mathcal{X}\right)$ as $L$ tends to infinity in dimensions
$d\geq 4$.

######  Theorem 1 (Infinite-volume translation-covariant gradient Gibbs states
in $d\geq 4$).

Suppose $d\geq 4$, let $\lambda>0$ and assume that
$(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent with
$\mathbb{E}\left[\eta(x)^{2}\right]=1$ and a symmetric distribution. Then, for
almost every realization of the disorder $\eta$, the sequence of measures
$(\mu_{\Lambda_{L}}^{\eta})_{L\geq 0}$ converges weakly in
$\mathcal{P}\left(\mathcal{X}\right)$ to a measure
$\mu^{\eta}_{\nabla}\in\mathcal{P}\left(\mathcal{X}\right)$ which satisfies
the following properties:

  * •

Quantitative convergence: There exist an exponent $\alpha>0$ depending on the
dimension $d$ and the ellipticity constants $c_{-},c_{+}$ and a constant $C$
depending on the parameters $d,c_{+},c_{-}$ and $\lambda$ such that, for any
side length $L\geq 1$,

(1.5)
$\mathbb{E}\left[\inf_{\kappa^{\eta}\in\Gamma_{\nabla}\left(\mu^{\eta}_{\Lambda_{L}},\mu^{\eta}_{\nabla}\right)}\int_{\mathcal{X}\times\mathcal{X}}\sup_{e\in
E(\Lambda_{L/2})}\left|\chi_{0}(e)-\chi_{1}(e)\right|^{2}\kappa^{\eta}\left(d\chi_{0},d\chi_{1}\right)\right]\leq
CL^{-\alpha}.$

  * •

Translation covariance: For any vertex $y\in\mathbb{Z}^{d}$ and almost every
realization of the random field, we have

(1.6)
$\mu^{\tau_{y}\eta}_{\nabla}=\left(\tau_{y}\right)_{*}\mu^{\eta}_{\nabla}.$

  * •

Dobrushin–Lanford–Ruelle (DLR) equation: For any finite set
$\Lambda\subseteq\mathbb{Z}^{d}$, any function $f\in C_{b}(\mathcal{X})$ and
almost every realization of the disorder, one has the identity

(1.7)
$\int_{\mathcal{X}}\mu^{\eta}_{\nabla}(d\chi)\int_{\mathcal{X}}f(\nabla\phi)\mu_{\Lambda}^{\eta,\psi(\chi)}(d\phi)=\int_{\mathcal{X}}f(\chi)\mu^{\eta}_{\nabla}(d\chi),$

where $\psi(\chi)$ is any function whose discrete gradient is equal to $\chi$.

As mentioned above, the result of Theorem 1 should hold in any dimension
$d\geq 3$, and the three dimensional case is not covered by our argument
(while it is for instance included in the results of [19, 20]). The technical
reason behind is that, as established in [21, Theorem 2], the fluctuations of
the disordered random surface distributed according to the measure
$\mu_{\Lambda_{L}}^{\eta}$ are of order $\sqrt{L}$ in three dimensions, while
they are of order $\sqrt{\ln L}$ in four dimensions and are bounded in
dimensions $5$ and higher. In dimension $d=3$, these fluctuations are too
large to be compensated by the $\alpha$-Hölder regularity provided by the De
Giorgi-Nash-Moser regularity theory (see Proposition 2.2), which only holds
for a small exponent $\alpha>0$.

The second theorem of this article states that, for almost every realization
of the random field, the sequence of measures $(\mu_{\Lambda_{L}}^{\eta})$
converges weakly in $\mathcal{P}\left(\Omega\right)$ as $L$ tends to infinity,
and establishes the existence of infinite-volume translation-covariant Gibbs
states in dimensions $d\geq 5$.

######  Theorem 2 (Infinite-volume translation-covariant Gibbs states in
$d\geq 5$).

Suppose $d\geq 5$, let $\lambda>0$ and assume that
$(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent with
$\mathbb{E}\left[\eta(x)^{2}\right]=1$ and a symmetric distribution. Then, for
almost every realization of the disorder $\eta$, the sequence of measures
$(\mu_{\Lambda_{L}}^{\eta})_{L\geq 0}$ converges weakly in
$\mathcal{P}\left(\Omega\right)$ to a measure
$\mu^{\eta}\in\mathcal{P}\left(\Omega\right)$ which satisfies the following
properties:

  * •

Quantitative convergence: There exists an exponent $\alpha>0$ depending on the
dimension $d$ and the ellipticity constants $c_{-},c_{+}$ and a constant $C$
depending on the parameters $d,c_{+},c_{-}$ and $\lambda$ such that, for any
$L\geq 1$,

(1.8)
$\sup_{x\in\Lambda_{L/4}}\mathbb{E}\left[\inf_{\kappa^{\eta}\in\Gamma\left(\mu^{\eta}_{\Lambda_{L}},\mu^{\eta}\right)}\int_{\Omega\times\Omega}\left|\phi_{0}(x)-\phi_{1}(x)\right|^{2}\kappa^{\eta}\left(d\phi_{0},d\phi_{1}\right)\right]\leq
CL^{-\alpha}.$

  * •

Translation covariance: For any vertex $y\in\mathbb{Z}^{d}$ and almost every
realization of the random field, we have

(1.9) $\mu^{\tau_{y}\eta}=\left(\tau_{y}\right)_{*}\mu^{\eta}.$

  * •

Dobrushin–Lanford–Ruelle (DLR) equation: For any finite set
$\Lambda\subseteq\mathbb{Z}^{d}$, any function $f\in C_{b}(\Omega)$ and almost
every realization of the random field, one has the identity

(1.10)
$\int_{\Omega}\mu^{\eta}(d\psi)\int_{\Omega}f(\phi)\mu_{\Lambda}^{\eta,\psi}(d\phi)=\int_{\Omega}f(\phi)\mu^{\eta}(d\phi).$

We conclude this section by mentioning that a follow-up question to Theorem 1
and Theorem 2 is the convergence to a scaling limit for the random-field
random surface model. It was pointed out in [21, Section 7] that the ground
state of random-field Gaussian free field coincides with the lattice _membrane
model_ , and that it is natural to expect that the scaling limit is the
continuum membrane model.

###  1.2. Related results

####  1.2.1. Random surfaces without disorder

The study of random surfaces (without disorder) was initiated in the 1970s by
Brascamp, Lieb and Lebowitz [13] who obtained sharp localization and
delocalization estimates for the random surface under the uniform convexity
assumption (1.2). In 1998, Funaki and Spohn [29] established that the large-
scale limit of the Langevin dynamics associated with the model is a
deterministic quantity governed by a parabolic equation taking the form of a
mean curvature flow. The scaling limit of the random surface was identified by
Naddaf and Spencer [41], and reworked probabilistically by Giacomin, Olla and
Spohn [30] in order to establish that the fluctuations of the space-time
dynamics associated with the model behave like an infinite-dimensional
Ornstein-Uhlenbeck process over large scales. We mention that, in the earlier
work [15], Brydges and Yau established, among other results, a statement
similar to the one of [41] in a perturbative setting by using a
renormalization group argument. In [25], Deuschel, Giacomin and Ioffe obtained
large deviation estimates and concentration inequalities for the model. Sharp
decorrelation estimates for the surface and its discrete gradient were
established by Delmotte and Deuschel [24]. More recently, the techniques of
[41] were generalized by Miller [38] to finite-volume measures with general
tilt. They were also combined with techniques of quantitative stochastic
homogenization by Armstrong and Wu [6] to obtain the $C^{2}$ regularity of the
surface tension of the $\nabla\phi$-model.

The case when the potential $V$ is not (uniformly) convex may exhibit a
different phenomenology, and a first-order phase transition from uniqueness to
non-uniqueness of infinite-volume (gradient) Gibbs measure (with zero tilt)
may be observed. In [7], Biskup and Kotecký studied a class of models with
potentials of the form

(1.11)
$e^{-V_{e}(\chi)}=pe^{-\kappa_{e}\chi^{2}}+(1-p)e^{-\kappa_{e}^{\prime}\chi^{2}}~{}\mbox{with}~{}\kappa_{e},\kappa_{e}^{\prime}>0,\,p\in[0,1]$

and established the coexistence of infinite-volume, shift-ergodic gradient
Gibbs measures when the parameters $\kappa_{e},\kappa_{e}^{\prime},p$ are
chosen suitably. We also refer to [47, 48] where a related phenomenon was
observed for spin systems with continuous symmetry and sufficiently nonlinear
interaction potential. In [8], Biskup and Spohn proved that, for a general
category of non-convex potentials including (1.11), the scaling limit of the
model is a Gaussian free field. Cotar and Deuschel [17] then proved the
uniqueness of ergodic Gibbs measures, obtained sharp estimates on the decay of
covariance and identified the scaling limit for a class of random surfaces
with non-convex potentials. In [18], Cotar, Deuschel and Müller established
the strict convexity of the surface tension of the $\nabla\phi$-model with a
non-convex potential which is a perturbation of a uniformly-convex potential.
The uniform convexity of the surface tension for a class of random surfaces
with non-convex potentials was obtained by Adams, Kotecky and Müller [1]
through a renormalization group argument. Further details on random surfaces
without external field can be found in the surveys of Funaki [28], Sheffield
[44] and Velenik [49].

####  1.2.2. Disordered random surfaces

In their seminal work [35], Imry and Ma predicted that the first-order phase
transitions in low-dimensional spin systems are rounded upon the addition of a
quenched random external field. Their predictions were confirmed for a large
class of spin systems by Aizenman and Wehr [4, 5], and recently quantified in
a series of works in the case of the two-dimensional random-field Ising model
[16, 3, 27, 2, 26]. While the phenomenon has mostly been investigated in the
context of compact spin spaces, a related effect occurs for the
$\nabla\phi$-model when a quenched random field is suitably added to the
system.

Motivated by disordered spin systems, Bovier and Külske [9, 10] studied the
effect of the addition of a $(d+1)$-dimensional disorder to the SOS-model for
integer-valued random surfaces. Specifically, the disordered Hamiltonian they
considered takes the following form:

$H^{\eta}(\phi):=\sum_{x\sim
y}|\phi(x)-\phi(y)|+\lambda\sum_{x,k}\eta_{k}(x)\mathds{1}_{\\{\phi(x)=k\\}},$

for integer-valued surfaces $\phi:\Lambda\to\mathbb{Z}$ and where the sum is
considered over the vertices $x\in\mathbb{Z}^{d}$ and the integers
$k\in\mathbb{Z}$, and the random variables $(\eta_{k}(x))_{k,x}$ are i.i.d.
with mean zero and unit variance. Based on the renormalization group argument
for the three dimensional random-field Ising model [14], they established in
[9] the existence of infinite-volume Gibbs measures in dimensions $d\geq 3$
when the strength of the external field and the temperature are small enough,
and, by adapting the strategy of Aizenman and Wehr [5, 4], they proved that
infinite-volume, translation-covariant Gibbs states cannot exist for the model
in dimensions $d\leq 2$ for any temperature and any strength of the random
field.

Another class of random surfaces for which the phenomenon is observed is the
model (1.3) which was first studied by van Enter and Külske [46]. They
established that, contrary to the pure model (with $\eta\equiv 0$), infinite-
volume gradient Gibbs measures cannot exist in two dimensions when the
potential $V$ is even, continuously differentiable and has super-linear
growth. They additionally obtained a quantitative estimate on the decay of
correlation for infinite-volume translation-covariant gradient Gibbs measures
in three dimensions.

In the article [19], Cotar and Külske studied two models of disordered random
surfaces. The first one is the model (1.3). In the second one, the disorder is
added to the model through the interacting potential as follows: they consider
a collection of independent and identically distributed functions
$\left(V_{e}\right)_{e\in E(\mathbb{Z}^{d})}$ satisfying suitable growth
condition and study the random surface whose distribution is formally given by

(1.12)
$\mu:=\exp\left(-\sum_{\begin{subarray}{c}e=(x,y)\end{subarray}}V_{e}(\phi(x)-\phi(y))\right)/Z.$

In this setting, they investigate the question of the existence of infinite-
volume, translation-covariant gradient Gibbs states for the model. They prove
that such measures exist when the external field satisfies
$\mathbb{E}\left[\eta(x)\right]=0$ in dimensions $d\geq 3$ for (1.3) and in
dimensions $d\geq 1$ for (1.12), as well as the existence of an infinite-
volume surface tension under the same assumptions. They additionally establish
that the infinite-volume surface tension does not exist for (1.3) in
dimensions $d\leq 2$, and that none of these quantities exist if the condition
$\mathbb{E}\left[\eta(x)\right]=0$ is not satisfied.

Their result is closely related to Theorem 1, but the techniques implemented
in the two articles rely on different strategies. The article [19] uses
surface tension bounds to obtain tightness of a suitable sequence of spatially
averaged finite-volume gradient Gibbs measures. The techniques developed in
the proof of Theorem 1 are based on a coupling argument for the Langevin
dynamics associated with the model (initially introduced by Funaki and Spohn
[29]), together regularity estimates and upper bounds on the fluctuations of
the random surface.

In the article [20], Cotar and Külske proved that infinite-volume translation-
covariant and ergodic gradient Gibbs states are unique when the tilt of the
interface is specified in any dimension $d\geq 1$ for the model (1.12), and in
any dimension $d\geq 3$ when the distribution of the disorder is symmetric for
the model (1.3). They additionally obtained sharp decorrelation estimates for
general functionals of the discrete gradient of the field for both models. We
mention that some of the tools needed to prove the result are used in the
proof of Theorem 2 (see Lemma 4.2 and Lemma 4.3 below).

Regarding the height of the random surface, Külske and Orlandi [36] considered
the model (1.3) in two dimensions when the domain $\Lambda$ is a box. They
establish that the typical height of the random surface at the center of the
box is at least logarithmic in the side length of the box when the potential
$V$ grows super-linearly for any realization of the disorder. In [37], Külske
and Orlandi established that, when the potential $V$ has super linear growth
and satisfies $V^{\prime\prime}\leq C$, the random surface is delocalized
(quantitatively) in two dimensions even when the interface is subjected to a
strong $\delta$-pinning. They additionally obtained a lower bound on the
height of the interface in dimensions $d\geq 3$ (exhibiting localization this
time) as well as a lower bound on the fraction of pinned sites.

In the article [21], the author, Harel and Peled obtained lower and upper
bounds for the typical height of the random surface (1.3) under the convexity
assumption (1.2), exhibiting delocalization in dimensions $d\leq 4$ and
localization in dimensions $d\geq 5.$

###  1.3. Strategy of the proofs

In this section, we present of the main arguments developed in this article.

####  1.3.1. Translation-covariant gradient Gibbs states

To highlight the main ideas of the proof of Theorem 1, we only describe the
strategy in the case of the ground state of the model defined to be the
minimizer of the variational problem

$\inf_{\begin{subarray}{c}v:\Lambda_{L}^{+}\to\mathbb{R}\\\
v=0\,\mathrm{on}\,\partial\Lambda_{L}\end{subarray}}\sum_{\begin{subarray}{c}x,y\in\Lambda_{L}^{+}\\\
x\sim
y\end{subarray}}V\left(v(y)-v(x)\right)-\lambda\sum_{x\in\Lambda_{L}}\eta(x)v(x).$

We denote by $v_{L,\eta}$ the unique minimizer; it can be equivalently
characterized as the solution of the discrete nonlinear elliptic equation

(1.13) $\left\\{\begin{aligned} -\sum_{y\sim
x}V^{\prime}\left(v_{L,\eta}(y)-v_{L,\eta}(x)\right)&=\lambda\eta(x)&~{}\mbox{for}~{}x\in\Lambda_{L},\\\
v_{L,\eta}(x)&=0&~{}\mbox{for}~{}x\in\partial\Lambda_{L}.\end{aligned}\right.$

In this setting, the question of the weak convergence of the disordered Gibbs
measure $\mu^{\eta}_{\Lambda_{L}}$ (Theorem 1) can be reformulated as follows:
for almost every realization of the disorder, the sequence of functions
$(\nabla v_{L,\eta})_{L\geq 0}$ converges pointwise as $L$ tends to infinity.
The argument relies on the observation that, for any fixed realization of the
disorder $\eta$, the difference $w_{L}:=v_{2L,\eta}-v_{L,\eta}$ solves _the
uniformly elliptic linear equation_ (see (2.1))

$-\nabla\cdot\mathbf{a}\nabla w_{L}=0~{}\mbox{in}~{}\Lambda_{L},$

where the environment $\mathbf{a}:E(\Lambda_{L})\to\mathbb{R}$ is explicit,
depends on the functions $v_{\eta,L}$ and $V$, and satisfies the uniform
ellipticity estimates $c_{-}\leq\mathbf{a}\leq c_{+}$ (see (3.5)). The De
Giorgi-Nash-Moser regularity estimate in the discrete setting (the version
below is an immediate consequence of [23, Proposition 5.3 and Proposition
6.2]) implies that there exists an exponent $\alpha>0$ such that

(1.14) $\sup_{\begin{subarray}{c}x,y\in\Lambda_{L/2}\\\ x\neq
y\end{subarray}}\frac{\left|w_{L,\eta}(x)-w_{L,\eta}(x)\right|}{|x-y|^{\alpha}}\leq\frac{C}{L^{\alpha}}\left(\frac{1}{\left|\Lambda_{L}\right|}\sum_{x\in\Lambda_{L}}w_{L,\eta}(x)^{2}\right)^{\frac{1}{2}},$

and thus, by restricting the supremum to the pairs of neighboring vertices,

(1.15) $\sup_{e\in E\left(\Lambda_{L/2}\right)}\left|\nabla
w_{L,\eta}(e)\right|\leq\frac{C}{L^{\alpha}}\left(\frac{1}{\left|\Lambda_{L}\right|}\sum_{x\in\Lambda_{L}}w_{L,\eta}(x)^{2}\right)^{\frac{1}{2}}.$

Using the results on the fluctuations of the random surface in the presence of
a random field established in [21, Theorem 2] and stated in Proposition 2.3
below (the result stated below only covers the case of the $\nabla\phi$-model,
but its extension to the ground state follows from [21, Theorem 2]), we can
estimate the term in the right-hand side of (1.15) and obtain

(1.16) $\mathbb{E}\left[\sup_{e\in E\left(\Lambda_{L/2}\right)}\left|\nabla
v_{2L,\eta}(e)-\nabla
v_{L,\eta}(e)\right|\right]\leq\frac{C}{L^{\alpha}}\left(\frac{1}{\left|\Lambda_{L}\right|}\sum_{x\in\Lambda_{L}}\mathbb{E}\left[\left|w_{L,\eta}(x)\right|^{2}\right]\right)^{\frac{1}{2}}\leq\frac{C}{L^{\alpha}}.$

The inequality (1.16) implies that the sequence $\nabla v_{2^{l},\eta}(e)$
converges almost surely by summing over the scales.

The extension of the result from the ground state to the disordered Gibbs
measure $\mu^{\eta}_{\Lambda_{L}}$ (as stated in Theorem 1) is achieved by
extending the argument from the elliptic setting to the parabolic setting,
using the Langevin dynamics associated with the model (see Section 2.2.1) and
the notion of Wasserstein metric (see Section 2.1.3) to estimate the distance
between two probability distributions.

####  1.3.2. Translation-covariant Gibbs states

As in Section 1.3.1, we only outline the argument in the case of the ground
state of the model. In this setting, Theorem 2 can be restated as follows: for
almost every realization of the disorder, the sequence of functions
$(v_{L,\eta})_{L\geq 0}$ converges pointwise as the side length $L$ tends to
infinity. To establish the claim, we use standard results on the gradient of
the Neumann Green’s function on $\mathbb{Z}^{d}$ (see Appendix A) which imply
that, for any $x\in\Lambda_{L/2}$ and any $\ell\in\mathbb{N}$, there exists a
map
$\mathbf{f}_{x+\Lambda_{\ell}}:\vec{E}\left(x+\Lambda_{\ell}\right)\to\mathbb{R}$
which is bounded by a constant depending only on the dimension and satisfies

$\displaystyle v_{2L,\eta}(x)-v_{L,\eta}(x)$
$\displaystyle=\sum_{e\in\vec{E}\left(x+\Lambda_{\ell}\right)}\mathbf{f}_{x+\Lambda_{\ell}}(e)\left(\nabla
v_{2L,\eta}(e)-\nabla v_{L,\eta}(e)\right)$
$\displaystyle\quad+\frac{1}{\left|\Lambda_{\ell}\right|}\sum_{y\in(x+\Lambda_{\ell})}\left(v_{2L,\eta}(y)-v_{L,\eta}(y)\right).$

By choosing the side length $\ell$ to be a suitable power of $L$, we may use
(1.16) to estimate the first term in the right-hand side, and a decorrelation
estimate (see Section 4.1) to establish that the spatial averages of the
functions $v_{2L,\eta}$ and $v_{L,\eta}$ over the cube $x+\Lambda_{\ell}$ are
small in dimensions $d\geq 5$. Combining the different arguments, we obtain
that there exists an exponent $\alpha>0$ such that, for any vertex
$x\in\Lambda_{L/4}$,

$\mathbb{E}\left[\left|v_{2L,\eta}(x)-v_{L,\eta}(x)\right|\right]\leq\frac{C}{L^{\alpha}}.$

The rest of the argument is then similar to the one of Theorem 1 presented in
Section 1.3.1.

###  1.4. Organisation of the article

The rest of the article is organized as follows. Section 2 introduces some
additional notations and preliminary results used in the proofs of Theorem 1
and Theorem 2. Section 3 and Section 4 contain the proofs of Theorem 1 and
Theorem 2 respectively. Appendix A contains the proof of a technical (and
standard) result pertaining to the gradient of the Neumann Green’s function
used in the proof of Theorem 1. Appendix B builds upon standard results from
the theory of elliptic regularity and contains the proof of the version of the
De Giorgi-Nash-Moser regularity estimate used in the article.

###  1.5. Convention for constants

Throughout this article, the symbols $C$ and $c$ denote positive constants
which may vary from line to line, with $C$ increasing and $c$ decreasing. The
symbol $\alpha$ denotes a positive exponent which may vary from line to line.
These constants and exponents may depend only on the dimension $d$, the
ellipticity constants $c_{+}$ and $c_{-}$ and the intensity $\lambda$ of the
random field. We specify the dependency of the constants and exponents by
writing, for instance, $C:=C(d,c_{+})$ to mean that the constant $C$ depends
on the parameters $d$ and $c_{+}$. We simply write $C$ to mean
$C:=C(d,c_{+},c_{-},\lambda)$.

#### Acknowledgements

The research of the author was supported in part by Israel Science Foundation
grants 861/15 and 1971/19 and by the European Research Council starting grant
678520 (LocalOrder). The author would like to thank the anonymous referees who
provided useful and detailed comments on an earlier version of the manuscript,
as well as for suggesting a simplification in the proofs of the estimates
(1.5) and (1.8).

##  2\. Notation and preliminary results

###  2.1. Notation

In this section, we introduce some general notations pertaining to lattices,
discrete functions and probability distributions.

####  2.1.1. General notation

We consider the hypercubic lattice $\mathbb{Z}^{d}$ in dimension $d\geq 4$ and
denote by $\vec{E}(\mathbb{Z}^{d})$ (resp. $E(\mathbb{Z}^{d})$) its set of
directed (resp. undirected) edges. Given a directed edge
$e\in\vec{E}\left(\mathbb{Z}^{d}\right)$, we denote by $x_{e}$ and $y_{e}$ the
first and second endpoint of $e$. We denote by $\left|\cdot\right|$ the
Euclidean norm on $\mathbb{Z}^{d}$ and set
$\left|\cdot\right|_{+}:=\max(\left|\cdot\right|,1)$. We write $x\sim y$ to
denote that $\\{x,y\\}\in E(\mathbb{Z}^{d})$. We identify a subset
$\Lambda\subseteq\mathbb{Z}^{d}$ with the induced subgraph of $\mathbb{Z}^{d}$
on $\Lambda$ and, in particular, denote by $\vec{E}\left(\Lambda\right)$ and
$E\left(\Lambda\right)$ the sets of directed and undirected edges of
$\Lambda$. We denote by $\partial\Lambda$ the _external vertex boundary_ of
$\Lambda$, that is,

$\partial\Lambda:=\left\\{x\in\mathbb{Z}^{d}\setminus\Lambda\,\colon\,\exists
y\in\Lambda,y\sim x\right\\}.$

We define $\Lambda^{+}:=\Lambda\cup\partial\Lambda$ and let
$\left|\Lambda\right|$ be the cardinality of $\Lambda$. We write
$\Lambda_{L}:=\\{-L,\ldots,L\\}^{d}\subseteq\mathbb{Z}^{d}$ for any integer
$L\geq 0$. This is extended to all $L\in[0,\infty)$ by setting
$\Lambda_{L}:=\Lambda_{\lfloor L\rfloor}$, where $\lfloor L\rfloor$ denotes
the floor of $L$. Given a point $x\in\mathbb{Z}^{d}$, we denote by
$x+\Lambda_{L}:=\left\\{y\in\mathbb{Z}^{d}\,:\,y-x\in\Lambda_{L}\right\\}$. A
box $\Lambda$ is then a subset of $\mathbb{Z}^{d}$ of the form
$\Lambda:=x+\Lambda_{L}$, for $x\in\mathbb{Z}^{d}$ and $L\in\mathbb{N}$.

Given a subset $\Lambda\subseteq\mathbb{Z}^{d}$ and a function
$\phi:\Lambda\to\mathbb{R}$, we define the discrete gradient
$\nabla\phi(e):=\phi(y)-\phi(x)$ for each directed edge
$e=(x,y)\in\vec{E}(\Lambda)$. We extend the definition to _time-dependent_
functions $\phi:[t_{-},t_{+}]\times\Lambda\to\mathbb{R}$ by setting
$\nabla\phi(t,e):=\phi(t,y)-\phi(t,x)$ for directed edges
$e=(x,y)\in\vec{E}(G)$ and times $t\in[t_{-},t_{+}]$. In expressions which do
not depend on the orientation of the edge, such as $|\nabla\phi(e)|^{2}$ or
$V(\nabla\phi(e))$, we allow the edge $e$ to be undirected. Given a vertex
$x\in\mathbb{Z}^{d}$, we use the notation

$\sum_{e\ni x}:=\sum_{e\in\vec{E}(\mathbb{Z}^{d}),x_{e}=x}.$

Given a time-dependent function $\phi:[s,t]\times\Lambda\to\mathbb{R}$, we
denote by $\partial_{t}\phi$ the derivative of $\phi$ with respect to the time
variable. Given a vertex $x\in\mathbb{Z}^{d}$ and a function $f$ depending on
the random field $\eta$, we denote by $\frac{\partial f}{\partial\eta(x)}$ the
derivative of the function $f$ with respect to the variable $\eta(x)$.

####  2.1.2. Environments and elliptic operators

A time-dependent environment is a map $\mathbf{a}:[0,\infty)\times
E(\Lambda)\to[0,\infty)$. We say that an environment is uniformly elliptic if
it satisfies the inequalities $0<c_{-}\leq\mathbf{a}\leq c_{+}<\infty$ for any
$(t,e)\in[0,\infty]\times E\left(\mathbb{Z}^{d}\right)$.

Given an environment $\mathbf{a}$, we define the time-dependent elliptic
operator $-\nabla\cdot\mathbf{a}\nabla$ by the formula: for a function
$\phi:[t_{-},t_{+}]\times\Lambda\to\mathbb{R}$, a vertex $x\in\Lambda$ and a
time $t\in[t_{-},t_{+}]$,

(2.1) $\nabla\cdot\mathbf{a}\nabla\phi(t,x)=\sum_{y\sim
x}\mathbf{a}(t,\\{x,y\\})(\phi(t,y)-\phi(t,x)).$

####  2.1.3. Coupling and Wasserstein metric

As in the introduction, we let $\Omega=\mathbb{R}^{\mathbb{Z}^{d}}$ be the
space of infinite-volume surfaces. Following the presentation of Funaki–Spohn
[29], we introduce the norm

$\left|\phi\right|_{1}:=\left(\sum_{x\in\mathbb{Z}^{d}}|\phi(x)|^{2}e^{-|x|}\right)^{\frac{1}{2}}.$

We denote by
$\Omega_{1}:=\left\\{\phi\in\Omega\,:\,\left|\phi\right|_{1}<\infty\right\\}$
and let $\mathcal{P}_{2}\left(\Omega_{1}\right)$ be the set of probability
measures on $\Omega_{1}$ with finite second moment. For any pair
$\mu_{0},\mu_{1}\in\mathcal{P}_{2}\left(\Omega_{1}\right)$, we define the
Wasserstein distance between $\mu_{0}$ and $\mu_{1}$ by the formula

$W_{\Omega}(\mu_{0},\mu_{1}):=\inf_{\kappa\in\Gamma(\mu_{0},\mu_{1})}\left(\int_{\Omega_{1}\times\Omega_{1}}\left|\phi_{0}-\phi_{1}\right|_{1}^{2}d\kappa(\phi_{0},\phi_{1})\right)^{\frac{1}{2}}.$

where $\Gamma(\mu_{0},\mu_{1})$ denotes the set of probability measures on the
space $\Omega_{1}\times\Omega_{1}$ whose first and second marginals are the
measures $\mu_{0}$ and $\mu_{1}$ respectively.

By [51, Theorem 6.18], the space
$(\mathcal{P}_{2}\left(\Omega_{1}\right),W_{\Omega})$ is metric, complete and
separable. Additionally, if a sequence of measures $(\mu_{n})_{n\geq 0}$
converges to a measure $\mu$ in the metric space
$(\mathcal{P}_{2}\left(\Omega_{1}\right),W_{\Omega})$ then it converges weakly
in the sense that, for any $f\in C_{b}(\Omega)$, $\left\langle
f\right\rangle_{\mu_{n}}\underset{n\to\infty}{\longrightarrow}\left\langle
f\right\rangle_{\mu}$ (see [51, Theorem 6.9]).

We denote by $\mathcal{X}$ the set of gradient fields, that is

$\mathcal{X}:=\left\\{\chi:\vec{E}(\mathbb{Z}^{d})\to\mathbb{R}\,:\,\exists\phi:\mathbb{Z}^{d}\to\mathbb{R}\,:\,\chi=\nabla\phi\right\\}.$

As it was the case in Section 2.1.1, in expressions which do not depend on the
orientation of the edge, we allow the edge $e$ to be undirected. We introduce
the norms

$\left|\chi\right|_{1}:=\left(\sum_{e\in\vec{E}(\mathbb{Z}^{d})}|\chi(e)|^{2}e^{-|x_{e}|}\right)^{\frac{1}{2}}$

We denote by
$\mathcal{X}_{1}:=\left\\{\chi\in\mathcal{X}\,:\,\left|\chi\right|_{1}<\infty\right\\}$,
and by $\mathcal{P}_{2}\left(\mathcal{X}_{1}\right)$ the set of probability
measures on $\mathcal{X}_{1}$ with finite second moment. For any pair
$\mu_{0},\mu_{1}\in\mathcal{P}_{2}\left(\mathcal{X}_{1}\right)$, we define the
Wasserstein distance between $\mu_{0}$ and $\mu_{1}$ according to the formula

(2.2)
$W_{\mathcal{X}}(\mu_{0},\mu_{1}):=\inf_{\kappa\in\Gamma_{\nabla}(\mu_{0},\mu_{1})}\left(\int_{\mathcal{X}_{1}\times\mathcal{X}_{1}}\left|\chi_{0}-\chi_{1}\right|_{1}^{2}d\kappa(\chi_{0},\chi_{1})\right)^{\frac{1}{2}}$

It is also clear from the definition of the Wasserstein metric
$W_{\mathcal{X}}$ and the norm $\left|\cdot\right|_{1}$ that, for any
$y\in\mathbb{Z}^{d}$ and any pair
$\mu_{0},\mu_{1}\in\mathcal{P}_{2}\left(\mathcal{X}_{1}\right)$,

(2.3)
$W_{\mathcal{X}}\left(\left(\tau_{y}\right)_{*}\mu_{0},\left(\tau_{y}\right)_{*}\mu_{1}\right)\leq
Ce^{C|y|}W_{\mathcal{X}}\left(\mu_{0},\mu_{1}\right),$

where the notation $\left(\tau_{y}\right)_{*}\mu_{0}$ denotes the pushforward
of the measure $\mu_{0}$ by the shift operator $\tau_{y}$ as introduced in
Section 1.1.

By [51, Theorem 6.18], the space
$(\mathcal{P}_{2}\left(\mathcal{X}_{1}\right),W_{\mathcal{X}})$ is metric,
complete and separable. Additionally, if a sequence of measures
$(\mu_{n})_{n\geq 0}$ converges to a measure $\mu$ in the metric space
$(\mathcal{P}_{2}\left(\mathcal{X}_{1}\right),W_{\mathcal{X}})$ then it
converges weakly in the sense that, for any $f\in C_{b}(\mathcal{X})$,
$\left\langle
f\right\rangle_{\mu_{n}}\underset{n\to\infty}{\longrightarrow}\left\langle
f\right\rangle_{\mu}$.

###  2.2. Preliminary results

In this section, we collect some tools pertaining to random surfaces and the
regularity theory for solutions of parabolic equations which are used in the
proofs of Theorem 1 and Theorem 2.

####  2.2.1. Langevin dynamics

The Gibbs measure $\mu_{\Lambda_{L}}^{\eta}$ (defined in (1.3)) is naturally
associated with the following dynamics.

######  Definition 2.1 (Langevin dynamics).

Consider a box $\Lambda\subseteq\mathbb{Z}^{d}$, a realization of the random
field $\eta:\Lambda\to\mathbb{R}$ and let
$\left\\{B_{t}(x)\,:\,x\in\Lambda\right\\}$ be a collection of independent
standard Brownian motions. The Langevin dynamics
$\left\\{\phi_{t}(x)\,:\,x\in\Lambda\right\\}$ is then defined to be the
solution of the system of stochastic differential equations

(2.4) $\left\\{\begin{aligned} \mathrm{d}\phi_{t}(x)&=\sum_{e\ni
x}V^{\prime}\left(\nabla\phi_{t}(e)\right)\mathrm{d}t+\lambda\eta(x)\mathrm{d}t+\sqrt{2}\,\mathrm{d}B_{t}(x)&&(t,x)\in(0,\infty)\times\Lambda,\\\
\phi_{t}(x)&=0&&(t,x)\in(0,\infty)\times\partial\Lambda_{L},\\\
\phi_{0}(x)&=0&&x\in\Lambda.\end{aligned}\right.$

The Langevin dynamics (2.4) is stationary, reversible and ergodic with respect
to the Gibbs measure $\mu_{\Lambda}^{\eta}$.

####  2.2.2. De Giorgi-Nash-Moser regularity theory

The next statement we collect is the De Giorgi-Nash-Moser regularity estimate
(following the works of De Giorgi [22], Nash [43] and Moser [39, 40], see also
[31, Chapter 8] and [34, Chapter 3]) which states that the solutions of the
parabolic equation $\partial_{t}u-\nabla\cdot\mathbf{a}\nabla u=0$ are Hölder
continuous when the environment $\mathbf{a}$ is uniformly elliptic. In the
discrete setting, the result can be stated as follows.

######  Proposition 2.2 (De Giorgi-Nash-Moser regularity for solutions of
parabolic equation).

Fix an integer $L\geq 2$ and an environment $\mathbf{a}:[0,\infty)\times
E\left(\mathbb{Z}^{d}\right)\to[0,\infty)$ satisfying $c_{-}\leq\mathbf{a}\leq
c_{+}$. Then, there exist an exponent $\alpha:=\alpha(d,c_{+},c_{-})>0$ and a
positive constant $C:=C(d,c_{+},c_{-})<\infty$ such that, for any function
$u:(0,L^{2})\times\Lambda_{L}\to\mathbb{R}$ solution of the parabolic equation

$\partial_{t}u-\nabla\cdot\mathbf{a}\nabla u=0\hskip 5.69054pt\mbox{in}\hskip
5.69054pt(0,L^{2})\times\Lambda_{L},$

one has, for any time $t\in[L^{2}/2,L^{2}]$,

(2.5) $\sup_{e\in E\left(\Lambda_{L/2}\right)}\left|\nabla
u(t,e)\right|\leq\frac{C}{L^{\alpha}}\left(\frac{1}{L^{d+2}}\int_{0}^{L^{2}}\sum_{x\in\Lambda_{L}}\left|u(t,x)\right|^{2}\,dt\right)^{\frac{1}{2}}.$

Proposition 2.2 is a consequence of two standard results for solutions of
parabolic equations: the mean value inequality (see [24, Proposition 3.2]) and
the increase of oscillations (see [24, Proposition 3.2]). We include a short
proof of the result in Appendix B for completeness.

####  2.2.3. Fluctuations for random-field random surfaces and the Brascamp-
Lieb inequality

In this section, we record the main result of the article [21] which provides
an upper bound on the size of the fluctuations of the random surface in the
presence of a random-field in dimensions $d\geq 4$. The bound stated in (2.6)
is sharp, and a corresponding lower bound is included in the statement of [21,
Theorem 2] (as well as upper and lower bounds on the fluctuations of the
surface in the low dimensions $d\leq 3$).

######  Proposition 2.3 (Theorem 2 of [21]).

Fix $\lambda>0$, assume that $(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent
with $\mathbb{E}\left[\eta(x)^{2}\right]=1$ and a symmetric distribution. Then
there exists a constant $C>0$ such that, for any integer $L\geq 2$ and any
$x\in\Lambda_{L}$,

(2.6)
$\mathbb{E}\left[\left\langle\phi(x)\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}^{2}\right]\leq\left\\{\begin{aligned}
C\ln L~{}\mbox{if}~{}d=4,\\\ C~{}\mbox{if}~{}d\geq 5.\end{aligned}\right.$

######  Remark 2.4.

The statement of [21, Theorem 2] estimates the variance over the disorder of
the map $\left\langle\phi(x)\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}$. The
assumption that the distribution of the random field is symmetric ensures that
$\mathbb{E}[\left\langle\phi(x)\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}]=0$
for any $x\in\Lambda_{L}$. The inequality (2.6) can thus be deduced from the
upper bound on the variance stated in [21, Theorem 2].

Proposition 2.3 estimates the $L^{2}$-norm of _the thermal expectation_ of the
height of the random surface. It is also natural to consider the value of the
$L^{2}$-norm of the height of the surface under the annealed measure
$\mathbb{E}\mu_{\Lambda_{L}}^{\eta}$. This quantity can be decomposed
according to the formula

(2.7)
$\mathbb{E}\left[\left\langle\phi(x)^{2}\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}\right]=\mathbb{E}\bigg{[}\Big{\langle}\big{(}\phi(x)-\left\langle\phi(x)\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}\big{)}^{2}\Big{\rangle}_{\mu_{\Lambda_{L}}^{\eta}}\bigg{]}+\mathbb{E}\left[\left\langle\phi(x)\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}^{2}\right].$

The second term in the right-hand side is estimated by (2.6). An upper bound
on the first term is provided by the Brascamp-Lieb inequality [11, 12] which,
in our setting, reads as follows: for any dimension $d\geq 3$, any side length
$L\geq 1$, any vertex $x\in\Lambda_{L}$ and any realization of the random
field $\eta$,

(2.8)
$\left\langle\left(\phi(x)-\left\langle\phi(x)\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}\right)^{2}\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}\leq
C.$

A combination of Proposition 2.3 and the inequality (2.8) implies the
following estimate.

######  Lemma 2.5.

Fix $\lambda>0$, assume that $(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent
with $\mathbb{E}\left[\eta(x)^{2}\right]=1$ and a symmetric distribution. Then
there exists a constant $C>0$ such that, for any integer $L\geq 2$ and any
$x\in\Lambda_{L}$,

$\mathbb{E}\left[\left\langle\phi(x)^{2}\right\rangle_{\mu_{\Lambda_{L}}^{\eta}}\right]\leq\left\\{\begin{aligned}
C\ln L~{}\mbox{if}~{}d=4,\\\ C~{}\mbox{if}~{}d\geq 5.\end{aligned}\right.$

##  3\. Translation-covariant gradient Gibbs states in dimensions $d\geq 4$

This section is devoted to the proof of Theorem 1. The argument is decomposed
into two steps. In Section 3.1, we establish a general coupling statement
which allows to estimate the Wasserstein distance between the two measures
$\mu_{\Lambda_{0}}^{\eta}$, $\mu_{\Lambda_{1}}^{\eta}$ for two boxes
$\Lambda_{0},\Lambda_{1}$ of comparable size and large intersection
(Proposition 3.1). Section 3.2 then contains the proof of Theorem 1 based on
the result established in Proposition 3.1.

###  3.1. A general coupling statement

In this section, we construct a coupling between the two measures
$\mu_{\Lambda_{0}}^{\eta}$, $\mu_{\Lambda_{1}}^{\eta}$ so as to minimize the
$L^{\infty}$-norm of the difference of the gradients of the surfaces in the
bulk of the two boxes $\Lambda_{0}$ and $\Lambda_{1}$.

######  Proposition 3.1.

Fix a dimension $d\geq 4$, let $\lambda>0$, and assume that
$(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent with
$\mathbb{E}\left[\eta(x)^{2}\right]=1$ and a symmetric distribution. There
exists an exponent $\alpha:=\alpha(d,c_{+},c_{-})>0$ such that the following
holds. Let $L\geq 2$ be an integer and let
$\Lambda_{0},\Lambda_{1}\subseteq\mathbb{Z}^{d}$ be two boxes such that
$\Lambda_{L}\subseteq\Lambda_{0}\cap\Lambda_{1}$ and
$\Lambda_{0}\cup\Lambda_{1}\subseteq\Lambda_{4L}$. Then, for almost every
$\eta$, there exists a coupling $\nu^{\eta}$ between the measures
$\mu_{\Lambda_{0}}^{\eta}$ and $\mu_{\Lambda_{1}}^{\eta}$ such that

(3.1) $\mathbb{E}\left[\int_{\Omega\times\Omega}\sup_{e\in
E\left(\Lambda_{L/2}\right)}\left|\nabla\phi_{0}(e)-\nabla\phi_{1}(e)\right|^{2}\nu^{\eta}(d\phi_{0},d\phi_{1})\right]\leq
CL^{-\alpha}.$

###### Proof.

Fix an integer $L\geq 2$ and two boxes $\Lambda_{0},\Lambda_{1}$ satisfying
the assumptions of the proposition. We denote by $L_{0}$ and $L_{1}$ their
respective side lengths and note that $L_{0},L_{1}\leq 4L$. We then fix a
realization of the random field $\eta:\mathbb{Z}^{d}\to\mathbb{R}$, and let
$\phi_{0}$ and $\phi_{1}$ be two independent random surfaces distributed
according to the measures $\mu^{\eta}_{\Lambda_{0}}$ and
$\mu^{\eta}_{\Lambda_{1}}$ respectively. We consider the two Langevin dynamics
(driven by the same Brownian motions)

(3.2) $\left\\{\begin{aligned} \mathrm{d}\phi_{0,t}(x)&=\sum_{e\ni
x}V^{\prime}\left(\nabla\phi_{0,t}(e)\right)\mathrm{d}t+\lambda\eta(x)\mathrm{d}t+\sqrt{2}\mathrm{d}B_{t}(x)&&(t,x)\in(0,\infty)\times\Lambda_{0},\\\
\phi_{0,0}(x)&=\phi_{0}(x)&&x\in\Lambda_{0},\\\
\phi_{0,t}(x)&=0&&(t,x)\in(0,\infty)\times\partial\Lambda_{0},\end{aligned}\right.$

and

(3.3) $\left\\{\begin{aligned} \mathrm{d}\phi_{1,t}(x)&=\sum_{e\ni
x}V^{\prime}\left(\nabla\phi_{1,t}(e)\right)\mathrm{d}t+\lambda\eta(x)\mathrm{d}t+\sqrt{2}\mathrm{d}B_{t}(x)&&(t,x)\in(0,\infty)\times\Lambda_{1},\\\
\phi_{1,0}(x)&=\phi_{1}(x)&&x\in\Lambda_{1},\\\
\phi_{1,t}(x)&=0&&(t,x)\in(0,\infty)\times\partial\Lambda_{1}.\end{aligned}\right.$

Using that the Langevin dynamics (3.2) (resp. (3.3)) is invariant under Gibbs
measure $\mu^{\eta}_{\Lambda_{0}}$ (resp. $\mu^{\eta}_{\Lambda_{1}}$), we see
that the interface $\phi_{0,t}$ (resp. $\phi_{1,t}$) is distributed according
to the measure $\mu^{\eta}_{\Lambda_{0}}$ (resp. $\mu^{\eta}_{\Lambda_{1}}$)
for any time $t\geq 0$. We denote by $\mu$ the law of the joint process
$\left(\phi_{0,t},\phi_{1,t}\right)_{t\geq 0}$ and by
$\left\langle\cdot\right\rangle_{\mu}$ the expectation with respect to $\mu$.

Considering the difference between the equations (3.2) and (3.3) and using
that the driving Brownian motions are the same, we obtain that the map
$w(t,x):=\phi_{0,t}(x)-\phi_{1,t}(x)$ solves the equation

(3.4) $\mathrm{d}w(t,x)=\left(\sum_{e\ni
x}V^{\prime}\left(\nabla\phi_{0,t}(e)\right)-\sum_{e\ni
x}V^{\prime}\left(\nabla\phi_{1,t}(e)\right)\right)\mathrm{d}t\hskip
14.22636pt\mbox{in}\hskip
14.22636pt[0,\infty)\times\left(\Lambda_{0}\cap\Lambda_{1}\right).$

Using the definition (2.1) for the elliptic operator
$-\nabla\cdot\mathbf{a}\nabla$ and the assumption
$\Lambda_{L}\subseteq\Lambda_{0}\cap\Lambda_{1}$, we see that the identity
(3.4) implies that the function $w(t,x):=\phi_{0,t}(x)-\phi_{1,t}(x)$ solves
the parabolic equation

$\partial_{t}w-\nabla\cdot\mathbf{a}\nabla w=0\hskip 14.22636pt\mbox{in}\hskip
14.22636pt[0,\infty)\times\Lambda_{0}\cap\Lambda_{1},$

where the coefficient $\mathbf{a}$ is given by the formula

(3.5)
$\mathbf{a}(t,e):=\int_{0}^{1}V^{\prime\prime}\left(s\nabla\phi_{0,t}(e)+(1-s)\nabla{\phi_{1,t}}(e)\right)\,ds.$

Using the uniform convexity of the potential $V$ stated in (1.2) and applying
the De Giorgi-Nash-Moser regularity estimate stated in Proposition 2.2, we
obtain that there exists an exponent $\alpha:=\alpha(d,c_{+},c_{-})>0$ such
that

(3.6) $\displaystyle\sup_{e\in
E\left(\Lambda_{L/2}\right)}\left|\nabla\phi_{0,L^{2}}(e)-\nabla\phi_{1,L^{2}}(e)\right|^{2}$
$\displaystyle\leq\frac{C}{L^{\alpha}}\frac{1}{L^{d+2}}\int_{0}^{L^{2}}\sum_{x\in\Lambda_{L}}\left(\phi_{0,t}(x)-\phi_{1,t}(x)\right)^{2}\,dt$
$\displaystyle\leq\frac{C}{L^{\alpha}}\frac{1}{L^{d+2}}\int_{0}^{L^{2}}\sum_{x\in\Lambda_{L}}\phi_{0,t}(x)^{2}+\phi_{1,t}(x)^{2}\,dt.$

We then integrate over the measure $\mu$ and use the stationarity of the
Langevin dynamics (3.2) and (3.3). We obtain

(3.7) $\displaystyle\left\langle\sup_{e\in
E\left(\Lambda_{L/2}\right)}\left|\nabla\phi_{0,L^{2}}(e)-\nabla\phi_{1,L^{2}}(e)\right|^{2}\right\rangle_{\mu}$
$\displaystyle\leq\frac{C}{L^{\alpha}}\frac{1}{L^{d+2}}\int_{0}^{L^{2}}\sum_{x\in\Lambda_{L}}\left\langle\phi_{0,t}(x)^{2}\right\rangle_{\mu}+\left\langle\phi_{1,t}(x)^{2}\right\rangle_{\mu}\,dt$
$\displaystyle=\frac{C}{L^{\alpha}}\frac{1}{L^{d+2}}\int_{0}^{L^{2}}\sum_{x\in\Lambda_{L}}\left\langle\phi(x)^{2}\right\rangle_{\mu^{\eta}_{\Lambda_{0}}}+\left\langle\phi(x)^{2}\right\rangle_{\mu^{\eta}_{\Lambda_{1}}}\,dt$
$\displaystyle=\frac{C}{L^{\alpha}}\frac{1}{L^{d}}\sum_{x\in\Lambda_{L}}\left\langle\phi(x)^{2}\right\rangle_{\mu^{\eta}_{\Lambda_{0}}}+\left\langle\phi(x)^{2}\right\rangle_{\mu^{\eta}_{\Lambda_{1}}}.$

Averaging over the disorder and using Lemma 2.5, we obtain

$\mathbb{E}\left[\left\langle\sup_{e\in
E\left(\Lambda_{L/2}\right)}\left|\nabla\phi_{0,L^{2}}(e)-\nabla\phi_{1,L^{2}}(e)\right|^{2}\right\rangle_{\mu}\right]\leq\left\\{\begin{aligned}
C\left(\ln L_{0}+\ln L_{1}\right)L^{-\alpha}~{}\mbox{if}~{}d=4,\\\
CL^{-\alpha}~{}\mbox{if}~{}d\geq 5.\end{aligned}\right.$

Using the assumption $L_{0},L_{1}\leq 4L$ and reducing the value of the
exponent $\alpha$ to absorb the logarithm, we deduce that, in any dimension
$d\geq 4$,

(3.8) $\mathbb{E}\left[\left\langle\sup_{e\in
E\left(\Lambda_{L/2}\right)}\left|\nabla\phi_{0,L^{2}}(e)-\nabla\phi_{1,L^{2}}(e)\right|^{2}\right\rangle_{\mu}\right]\leq
CL^{-\alpha}.$

We then define $\nu^{\eta}$ to be the law of the pair
$\left(\phi_{0,L^{2}},\phi_{1,L^{2}}\right)$ to complete the proof of
Proposition 3.1.

∎

###  3.2. Proof of Theorem 1

This section is devoted to the proof of Theorem 1 building upon the result of
Proposition 3.1. We recall the definitions of the Wasserstein distance (2.2)
stated in Section 2.

###### Proof of Theorem 1.

We first observe that Lemma 2.5 and Proposition 3.1 imply the following
inequality: for each fixed integer $L\geq 1$ and each pair of boxes
$\Lambda_{0},\Lambda_{1}\subseteq\mathbb{Z}^{d}$ satisfying
$\Lambda_{L}\subseteq\Lambda_{0}\cap\Lambda_{1}$ and
$\Lambda_{0}\cup\Lambda_{1}\subseteq\Lambda_{4L}$

(3.9)
$\mathbb{E}\left[W_{\mathcal{X}}\left(\mu_{\Lambda_{0}}^{\eta},\mu_{\Lambda_{1}}^{\eta}\right)^{2}\right]\leq
CL^{-\alpha}.$

Let us then fix $n\in\mathbb{N}$. Applying (3.9) with the integer $L=2^{n}$,
the boxes $\Lambda_{0}=\Lambda_{2^{n}}$ and $\Lambda_{1}=\Lambda_{2^{n+1}}$,
we obtain

(3.10)
$\mathbb{E}\left[W_{\mathcal{X}}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu_{\Lambda_{2^{n+1}}}^{\eta}\right)^{2}\right]\leq
C2^{-\alpha n}.$

The inequality (3.10) and the comparison estimate (1.5) imply that the
sequence of measures $(\mu_{\Lambda_{2^{n}}}^{\eta})_{n\in\mathbb{N}}$ is
almost surely Cauchy in the metric space
$(\mathcal{P}_{2}\left(\mathcal{X}_{1}\right),W_{\mathcal{X}})$. Since this
metric space is complete, it implies that the sequence converges almost surely
and we denote the limiting measure by
$\mu^{\eta}_{\nabla}\in\mathcal{P}_{2}\left(\mathcal{X}_{1}\right)$. It
remains to verify that the measure $\mu^{\eta}_{\nabla}$ satisfies the
properties listed in the statement of Theorem 1.

We first verify the quantitative estimate (1.5). By the inequality (3.8) and
the gluing lemma (following the terminology of [50, Chapter 1]), there exists
a coupling $\pi^{\eta}_{\nabla}$ between the measures
$\left(\mu_{\Lambda_{2^{k}L}}^{\eta}\right)_{k\geq 0}$ such that, for any
$k\geq 0$,

(3.11)
$\mathbb{E}\left[\left\langle\sup_{e\in\Lambda_{2^{k-1}L}}\left|\chi_{2^{k}L}(e)-\chi_{2^{k+1}L}(e)\right|^{2}\right\rangle_{\pi^{\eta}_{\nabla}}\right]\leq
C(2^{k}L)^{-\alpha},$

where $\left\langle\cdot\right\rangle_{\pi^{\eta}_{\nabla}}$ denotes the
expectation with respect to $\pi^{\eta}_{\nabla}$ and $\chi_{2^{k}L}(e)$ is
the gradient at the edge $e$ for the measure $\mu_{\Lambda_{2^{k}L}}^{\eta}$.
The inequality (3.11) implies that, for almost every $\eta$, the gradients
$\chi_{2^{k}L}$ converge pointwise almost surely with respect to the measure
$\pi^{\eta}_{\nabla}$ to a random variable $\chi_{\infty}$ distributed
according to $\mu^{\eta}_{\nabla}$. The triangle inequality additionally shows

(3.12)
$\mathbb{E}\left[\left\langle\sup_{e\in\Lambda_{L/2}}\left|\chi_{L}(e)-\chi_{\infty}(e)\right|^{2}\right\rangle_{\pi^{\eta}_{\nabla}}\right]^{\frac{1}{2}}\leq\sum_{k=0}^{\infty}\mathbb{E}\left[\left\langle\sup_{e\in\Lambda_{2^{k-1}L}}\left|\chi_{2^{k}L}(e)-\chi_{2^{k+1}L}(e)\right|^{2}\right\rangle_{\pi^{\eta}_{\nabla}}\right]^{\frac{1}{2}}\leq
CL^{-\alpha},$

which implies (1.5).

We next verify the translation covariance (1.6) of the measure
$\mu^{\eta}_{\nabla}$. To this end, let us note that, by the definitions of
the measure $\mu_{\Lambda}^{\eta}$ and of the shift operator $\tau_{y}$, one
has the identity: for any point $y\in\mathbb{Z}^{d}$, any realization of the
random field $\eta$ and any side length $L\in\mathbb{N}$,

(3.13)
$\left(\tau_{y}\right)_{*}\mu_{\Lambda_{L}}^{\eta}=\mu_{y+\Lambda_{L}}^{\tau_{y}\eta}.$

We then select an integer $n$ large enough (depending on $y$) so that
$\Lambda_{2^{n-1}}\subseteq y+\Lambda_{2^{n}}$. Applying Proposition 3.1 with
the integer $L:=2^{n-1}$ and the boxes $\Lambda_{0}:=\Lambda_{2^{n}}$ and
$\Lambda_{1}:=y+\Lambda_{2^{n}}$, we obtain

(3.14)
$\mathbb{E}\left[W_{\mathcal{X}}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu_{y+\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]\leq
C2^{-\alpha n}.$

Using the inequality (3.14) and the translation invariance of the law of the
disorder $\eta$, we deduce that

(3.15)
$\displaystyle\mathbb{E}\left[W_{\mathcal{X}}\left(\mu^{\tau_{y}\eta}_{\nabla},\mu_{y+\Lambda_{2^{n}}}^{\tau_{y}\eta}\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[W_{\mathcal{X}}\left(\mu^{\eta}_{\nabla},\mu_{y+\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]$
$\displaystyle\leq
2\mathbb{E}\left[W_{\mathcal{X}}\left(\mu^{\eta}_{\nabla},\mu_{\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]+2\mathbb{E}\left[W_{\mathcal{X}}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu_{y+\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]$
$\displaystyle\leq C2^{-\alpha n}.$

Additionally, using the translation property (2.3) of the Wasserstein metric,
we may write

(3.16)
$\mathbb{E}\left[W_{\mathcal{X}}\left(\left(\tau_{y}\right)_{*}\mu^{\eta}_{\nabla},\left(\tau_{y}\right)_{*}\mu_{\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]\leq
Ce^{C|y|}\mathbb{E}\left[W_{\mathcal{X}}\left(\mu^{\eta}_{\nabla},\mu_{\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]\leq
Ce^{C|y|}2^{-\alpha n}.$

Combining the identity (3.13) with the inequalities (3.15), (3.16), we obtain,
for any integer $n\in\mathbb{N}$ (chosen large enough depending on $y$),

$\displaystyle\mathbb{E}\left[W_{\mathcal{X}}\left(\left(\tau_{y}\right)_{*}\mu^{\eta}_{\nabla},\mu^{\tau_{y}\eta}_{\nabla}\right)^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[W_{\mathcal{X}}\left(\left(\tau_{y}\right)_{*}\mu^{\eta}_{\nabla},\left(\tau_{y}\right)_{*}\mu_{\Lambda_{2^{n}}}^{\eta}\right)^{2}\right]+\mathbb{E}\left[W_{\mathcal{X}}\left(\mu_{y+\Lambda_{2^{n}}}^{\tau_{y}\eta},\mu^{\tau_{y}\eta}_{\nabla}\right)^{2}\right]$
$\displaystyle\leq C(1+e^{C|y|})2^{-\alpha n}.$

Take the limit $n\to\infty$ yields

(3.17)
$\mathbb{E}\left[W_{\mathcal{X}}\left(\left(\tau_{y}\right)_{*}\mu^{\eta}_{\nabla},\mu^{\tau_{y}\eta}_{\nabla}\right)^{2}\right]=0,$

which implies the translation-covariance property (1.6).

There remains to prove the Dobrushin–Lanford–Ruelle (DLR) equation (1.7). We
first note that, by the definition of the finite-volume measure
$\mu_{\Lambda_{2^{n}}}^{\eta}$, one has, for any function $f\in
C_{b}\left(\mathcal{X}\right)$ and any subset
$\Lambda\subseteq\Lambda_{2^{n}}$,

(3.18)
$\int_{\Omega}\mu^{\eta}_{\Lambda_{2^{n}}}(d\phi_{1})\int_{\Omega}f(\nabla\phi)\mu_{\Lambda}^{\eta,\phi_{1}}(d\phi)=\int_{\Omega}f(\nabla\phi)\mu^{\eta}_{\Lambda_{2^{n}}}(d\phi).$

One then observes that the mapping
$\phi_{1}\mapsto\int_{\Omega}f(\nabla\phi)\mu_{\Lambda}^{\eta,\phi_{1}}(d\phi)$
depends only on finitely many values of the (discrete) gradient of the map
$\phi_{1}$ (the values on which $f$ depends and the values of the gradient on
the edges which are in a finite neighborhood of the boundary
$\partial\Lambda$). It is additionally continuous and bounded. Consequently,
it belongs to the space $C_{b}(\mathcal{X})$. Taking the limit $n\to\infty$ in
both sides of the identity (3.18) and using the almost-sure weak convergence
of the sequence $\mu_{\Lambda_{2^{n}}}^{\eta}$ in the space
$\mathcal{P}\left(\mathcal{X}\right)$, we obtain

$\int_{\Omega}\mu^{\eta}_{\nabla}(d\chi)\int_{\Omega}f(\nabla\phi)\mu_{\Lambda}^{\eta,\psi(\chi)}(d\phi)=\int_{\Omega}f(\chi)\mu^{\eta}_{\nabla}(d\chi),$

where $\psi(\chi)$ is any function satisfying $\nabla\psi(\chi)=\chi$. The
proof of (1.7), and thus of Theorem 1, is complete. ∎

##  4\. Translation-covariant Gibbs states in dimensions $d\geq 5$

This section is devoted to the proof of Theorem 2. The argument is decomposed
in three steps, which correspond to Section 4.1, Section 4.2 and Section 4.3.
In Section 4.1, we prove a decorrelation estimate for the spatial average of
the height of the random surface in dimensions $d\geq 5$. In Section 4.2, we
combine the decorrelation inequality established in Section 4.1 with the
coupling estimate for the gradient of the surface proved in Proposition 3.1,
and establish an upper bound for the Wasserstein distance between the two
measures $\mu_{\Lambda_{0}}^{\eta}$ and $\mu_{\Lambda_{1}}^{\eta}$ where
$\Lambda_{0}$ and $\Lambda_{1}$ are two boxes of similar size and large
intersection. Section 4.3 contains the proof of Theorem 2, building upon the
results of Section 4.2.

###  4.1. A decorrelation estimate for the random surface

We first establish an upper bound on the $L^{2}$-norm of the spatially
averaged height of the random surface in a box $\Lambda_{0}$ under the
annealed measure $\mathbb{E}\mu_{\Lambda}^{\eta}$. We note that the arguments
used in the proof of Lemma 4.1 are similar to the ones of [20, Theorem 3.1].

######  Lemma 4.1.

Fix a dimension $d\geq 5$, let $\lambda>0$, and assume that
$(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent with
$\mathbb{E}\left[\eta(x)^{2}\right]=1$ and a symmetric distribution. There
exists a constant $C<\infty$ such that for any box $\Lambda_{0}$ of side
length $\ell\geq 1$ and any box $\Lambda$ satisfying
$\Lambda_{0}\subseteq\Lambda$,

(4.1)
$\mathbb{E}\left[\left\langle\left|\frac{1}{\left|\Lambda_{0}\right|}\sum_{x\in\Lambda_{0}}\phi(x)\right|^{2}\right\rangle_{\mu^{\eta}_{\Lambda}}\right]\leq\frac{C}{\ell^{d-4}}.$

The proof of Lemma 4.1 requires to use the following two results for which we
refer to the articles [33, 20]. The first one provides an estimate on the
covariance of the variables $\phi(x)$ and $\phi(y)$ when the surface $\phi$ is
distributed according to the measure $\mu_{\Lambda}^{\eta}$ and applies
uniformly over the realizations of the disorder $\eta$. Before stating the
result, we recall the notation $|\cdot|_{+}:=\max\left(|\cdot|,1\right).$

######  Lemma 4.2 (Proposition 2.3 (ii) of [20]).

Fix a dimension $d\geq 3$. There exists a constant
$C:=C(d,c_{+},c_{-})<\infty$ such that for any $\lambda>0$, any realization of
the disorder $\eta:\Lambda\to\mathbb{R}$, and any pair of points
$x,y\in\Lambda$,

(4.2)
$\left|\left\langle\phi(x)\phi(y)\right\rangle_{\mu^{\eta}_{\Lambda}}-\left\langle\phi(x)\right\rangle_{\mu^{\eta}_{\Lambda}}\left\langle\phi(y)\right\rangle_{\mu^{\eta}_{\Lambda}}\right|\leq\frac{C}{|x-y|_{+}^{d-2}}.$

The next statement we collect is a covariance inequality for functions of the
random field. Estimates of this nature were used in stochastic homogenization
by Naddaf and Spencer [42], and then by Gloria and Otto [32, 33] to obtain
optimal quantitative estimates on the first-order corrector. Lemma 4.3 below
can, in particular, be deduced from [33, Lemma 3]. In the setting of random-
field random surfaces, it was used by Cotar and Külske in [20, Proposition
2.4].

######  Lemma 4.3.

Assume that $(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent with
$\mathbb{E}[\eta(x)^{2}]=1$ and fix a box $\Lambda\subseteq\mathbb{Z}^{d}$.
For any pair $f,g$ of Lipschitz and differentiable functions of the collection
$\left(\eta_{v}\right)_{v\in\Lambda}\in\mathbb{R}^{\Lambda}$, one has the
estimate

$\left|\mathbb{E}\left[fg\right]-\mathbb{E}\left[f\right]\mathbb{E}\left[g\right]\right|\leq\sum_{z\in\Lambda}\sup_{\eta\in\mathbb{R}^{\Lambda}}\left|\frac{\partial
f}{\partial\eta(z)}\right|\times\sup_{\eta\in\mathbb{R}^{\Lambda}}\left|\frac{\partial
g}{\partial\eta(z)}\right|.$

Using the two previous lemmas, we are now ready to show Lemma 4.1.

###### Proof of Lemma 4.1.

We first prove the following inequality: for any box
$\Lambda\subseteq\mathbb{Z}^{d}$ and any pair of vertices $x,y\in\Lambda$,

(4.3)
$\left|\mathbb{E}\left[\left\langle\phi(x)\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}\right]\right|\leq\frac{C}{|x-y|_{+}^{d-4}}.$

Using that the inequality (4.2) of Lemma 4.2 is valid for any realization of
the random field, the estimate (4.3) is equivalent to the following upper
bound: for any box $\Lambda\subseteq\mathbb{Z}^{d}$ and any pair of points
$x,y\in\Lambda$,

(4.4)
$\left|\mathbb{E}\left[\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}\right]\right|\leq\frac{C}{|x-y|_{+}^{d-4}}.$

To prove the inequality (4.4), we note that the functions
$\eta\mapsto\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}$ and
$\eta\mapsto\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}$ depend
only on the values of the field $\eta$ in the box $\Lambda$ and are
differentiable. An explicit computation yields the formulae, for any point
$z\in\Lambda$,

$\frac{\partial\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}}{\partial\eta(z)}=\lambda\left\langle\phi(x)\phi(z)\right\rangle_{\mu_{\Lambda}^{\eta}}-\lambda\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}\left\langle\phi(z)\right\rangle_{\mu_{\Lambda}^{\eta}}$

and

$\frac{\partial\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}}{\partial\eta(z)}=\lambda\left\langle\phi(y)\phi(z)\right\rangle_{\mu_{\Lambda}^{\eta}}-\lambda\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}\left\langle\phi(z)\right\rangle_{\mu_{\Lambda}^{\eta}}.$

Applying Proposition 4.2, we obtain that the functions
$\eta\mapsto\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}$ and
$\eta\mapsto\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}$ are
Lipschitz and satisfy, for any vertex $z\in\Lambda$ and any realization of the
random field $\eta:\Lambda\to\mathbb{R}$,

(4.5)
$\left|\frac{\partial\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}}{\partial\eta(z)}\right|\leq\frac{C}{|x-z|_{+}^{d-2}}\hskip
8.53581pt\mbox{and}\hskip
8.53581pt\left|\frac{\partial\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}}{\partial\eta(z)}\right|\leq\frac{C}{|y-z|_{+}^{d-2}}.$

We can thus apply Lemma 4.2 and Lemma 4.3 with the functions
$f:\eta\mapsto\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}$ and
$g:\eta\mapsto\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}$, and
note that the symmetry assumption on the random field implies the identities
$\mathbb{E}[\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}]=\mathbb{E}[\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}]=0$.
We obtain

$\displaystyle\left|\mathbb{E}\left[\left\langle\phi(x)\right\rangle_{\mu_{\Lambda}^{\eta}}\left\langle\phi(y)\right\rangle_{\mu_{\Lambda}^{\eta}}\right]\right|$
$\displaystyle\leq
C\sum_{z\in\Lambda}\frac{1}{|x-z|_{+}^{d-2}}\frac{1}{|y-z|_{+}^{d-2}}$
$\displaystyle\leq
C\sum_{z\in\mathbb{Z}^{d}}\frac{1}{|x-z|_{+}^{d-2}}\frac{1}{|y-z|_{+}^{d-2}}$
$\displaystyle\leq\frac{C}{|x-y|_{+}^{d-4}}.$

The proof of the inequality (4.4), and thus of (4.3), is complete. We next
prove the inequality (4.1). To this end, we use the inequality (4.3) and write

(4.6)
$\displaystyle\mathbb{E}\left[\left\langle\left|\frac{1}{\left|\Lambda_{0}\right|}\sum_{y\in\Lambda_{0}}\phi(y)\right|^{2}\right\rangle_{\mu_{\Lambda}^{\eta}}\right]$
$\displaystyle=\mathbb{E}\left[\left\langle\left|\frac{1}{\left|\Lambda_{0}\right|}\sum_{y\in\Lambda_{0}}\phi(y)\right|^{2}\right\rangle_{\mu_{\Lambda}^{\eta}}\right]$
$\displaystyle=\frac{1}{\left|\Lambda_{0}\right|^{2}}\sum_{y,z\in\Lambda_{0}}\mathbb{E}\left[\left\langle\phi(y)\phi(z)\right\rangle_{\mu_{\Lambda}^{\eta}}\right]$
$\displaystyle\leq\frac{C}{\ell^{2d}}\sum_{y,z\in\Lambda_{0}}\frac{1}{|y-z|_{+}^{d-4}}$
$\displaystyle\leq\frac{C}{\ell^{d-4}}.$

The proof of Lemma 4.1 is complete. ∎

###  4.2. A general coupling statement

In this section, we establish a general coupling statement to estimate the
Wasserstein distance between the two measures $\mu_{\Lambda_{0}}^{\eta}$,
$\mu_{\Lambda_{1}}^{\eta}$where $\Lambda_{0},\Lambda_{1}$ are two boxes of
comparable size and large intersection.

######  Proposition 4.4.

Fix a dimension $d\geq 5$, let $\lambda>0$, and assume that
$(\eta(x))_{x\in\mathbb{Z}^{d}}$ are independent with
$\mathbb{E}[\eta(x)^{2}]=1$ and a symmetric distribution. There exists an
exponent $\alpha:=\alpha(d,c_{+},c_{-})>0$ such that the following statement
holds. Let $L\geq 2$ be an integer and let
$\Lambda_{0},\Lambda_{1}\subseteq\mathbb{Z}^{d}$ be two boxes such that
$\Lambda_{L}\subseteq\Lambda_{0}\cap\Lambda_{1}$ and
$\Lambda_{0}\cup\Lambda_{1}\subseteq\Lambda_{4L}$. Then, one has the estimate

$\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{0}}^{\eta},\mu_{\Lambda_{1}}^{\eta}\right)^{2}\right]\leq
CL^{-\alpha}.$

###### Proof.

We fix an integer $L\geq 2$, two boxes $\Lambda_{0},\Lambda_{1}$ satisfying
the assumptions of the proposition, and let $\nu^{\eta}$ be the coupling
between the measures $\mu_{\Lambda_{0}}^{\eta}$ and $\mu_{\Lambda_{1}}^{\eta}$
provided by Proposition 3.1. By Proposition A.1, we know that, for any
$x\in\mathbb{Z}^{d}$ and any $\ell\in\mathbb{N}$, there exists a function
$\mathbf{f}_{x+\Lambda_{\ell}}:\vec{E}\left(x+\Lambda_{\ell}\right)\to\mathbb{R}$
and a constant $C:=C(d)<\infty$ such that, for any
$\phi:(x+\Lambda_{\ell})\to\mathbb{R}$,

(4.7)
$\phi(x)=\sum_{e\in\vec{E}\left(x+\Lambda_{\ell}\right)}\mathbf{f}_{x+\Lambda_{\ell}}(e)\nabla\phi(e)+\frac{1}{\left|\Lambda_{\ell}\right|}\sum_{y\in(x+\Lambda_{\ell})}\phi(y),$

and, for any edge $e\in\vec{E}\left(x+\Lambda_{\ell}\right),$

(4.8) $\left|\mathbf{f}_{x+\Lambda_{\ell}}(e)\right|\leq C.$

Let $\ell$ be integer such that $\ell\leq L/4$ and fix a vertex
$x\in\Lambda_{L/4}$. Noting that $x+\Lambda_{\ell}\subseteq\Lambda_{L/2}$, we
can apply the identity (4.7) and the inequality (4.8) with the box
$\Lambda:=(x+\Lambda_{\ell})$. We obtain, for any pair of functions
$\phi_{0},\phi_{1}:\Lambda_{L}\to\mathbb{R}$,

(4.9) $\displaystyle\left|\phi_{0}(x)-\phi_{1}(x)\right|$ $\displaystyle\leq
C\left|\Lambda_{\ell}\right|\sup_{e\in
E\left(x+\Lambda_{\ell}\right)}\left|\nabla\phi_{0}(e)-\nabla\phi_{1}(e)\right|$
$\displaystyle\quad+\left|\frac{1}{\left|\Lambda_{\ell}\right|}\sum_{y\in(x+\Lambda_{\ell})}\phi_{0}(y)\right|+\left|\frac{1}{\left|\Lambda_{\ell}\right|}\sum_{y\in(x+\Lambda_{\ell})}\phi_{1}(y)\right|.$

Taking the expectation with respect to measure $\nu^{\eta}$ (using that the
marginals of this distribution are $\mu_{\Lambda_{0}}^{\eta}$ and
$\mu_{\Lambda_{1}}^{\eta}$) and then over the disorder $\eta$, we obtain

$\displaystyle\sup_{x\in\Lambda_{L/4}}\mathbb{E}\left[\int_{\Omega\times\Omega}\left|\phi_{0}(x)-\phi_{1}(x)\right|^{2}\nu^{\eta}(d\phi_{0},d\phi_{1})\right]$
$\displaystyle\leq
C\ell^{d}\mathbb{E}\left[\int_{\Omega\times\Omega}\sup_{e\in
E\left(\Lambda_{L/2}\right)}\left|\nabla\phi_{0}(e)-\nabla\phi_{1}(e)\right|^{2}\nu^{\eta}(d\phi_{0},d\phi_{1})\right]$
$\displaystyle\quad+3\sup_{x\in\Lambda_{L/4}}\mathbb{E}\left[\left\langle\left|\frac{1}{\left|\Lambda_{\ell}\right|}\sum_{y\in(x+\Lambda_{\ell})}\phi(y)\right|^{2}\right\rangle_{\mu_{\Lambda_{0}}^{\eta}}\right]+3\sup_{x\in\Lambda_{L/4}}\mathbb{E}\left[\left\langle\left|\frac{1}{\left|\Lambda_{\ell}\right|}\sum_{y\in(x+\Lambda_{\ell})}\phi(y)\right|^{2}\right\rangle_{\mu_{\Lambda_{1}^{\eta}}}\right].$

The first term in the right-hand side can be estimated thanks to Proposition
3.1, the second and third terms can be bounded from above by using Lemma 4.1.
We obtain

$\sup_{x\in\Lambda_{L/4}}\mathbb{E}\left[\int_{\Omega\times\Omega}\left|\phi_{0}(x)-\phi_{1}(x)\right|^{2}\nu^{\eta}(d\phi_{0},d\phi_{1})\right]\leq
C\ell^{d}L^{-\alpha}+C\ell^{4-d}.$

Optimising over the value of $\ell$ (for instance choosing
$\ell=L^{\alpha/(2d)}$) and reducing the value of the exponent $\alpha$, we
deduce that

(4.10)
$\sup_{x\in\Lambda_{L/4}}\mathbb{E}\left[\int_{\Omega\times\Omega}\left|\phi_{0}(x)-\phi_{1}(x)\right|^{2}\nu^{\eta}(d\phi_{0},d\phi_{1})\right]\leq
CL^{-\alpha}.$

A combination of Lemma 2.5 and (4.10) shows

$\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{0}}^{\eta},\mu_{\Lambda_{1}}^{\eta}\right)^{2}\right]\leq
CL^{-\alpha}.$

∎

###  4.3. Proof of Theorem 2

This section contains the proof of Theorem 2 based on the coupling statement
established in Section 4.2.

###### Proof of Theorem 2.

Fix an integer $n\in\mathbb{N}$. Applying Proposition 4.4 with $L=2^{n}$, the
boxes $\Lambda_{0}=\Lambda_{2^{n}}$ and $\Lambda_{1}=\Lambda_{2^{n+1}}$, we
obtain

(4.11)
$\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu_{\Lambda_{2^{n+1}}}^{\eta}\right)^{2}\right]\leq
C2^{-\alpha n}.$

The inequality (4.11) implies that the sequence of measures
$(\mu_{\Lambda_{2^{n}}}^{\eta})_{n\in\mathbb{N}}$ is almost surely Cauchy in
the metric space $(\mathcal{P}_{2}\left(\Omega_{1}\right),W_{\Omega})$. Since
this metric space is complete, it implies that it converges almost surely to a
probability distribution
$\mu^{\eta}\in\mathcal{P}_{2}\left(\Omega_{1}\right)$. It remains to verify
that the measure $\mu^{\eta}$ satisfies the properties listed in the statement
of Theorem 2. First, the inequality (4.11) implies the estimate: for any
$L\geq 2$,

(4.12)
$\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{L}}^{\eta},\mu^{\eta}\right)^{2}\right]\leq
CL^{-\alpha}.$

We then fix an integer $L\geq 2$ and let $n$ be the integer satisfying
$2^{n-1}<L\leq 2^{n}$. Applying Proposition 4.4 with the boxes
$\Lambda_{0}:=\Lambda_{L}$ and $\Lambda_{1}:=\Lambda_{2^{n}}$, we obtain

(4.13)
$\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu_{\Lambda_{L}}^{\eta}\right)^{2}\right]\leq
CL^{-\alpha}.$

Consequently, by the triangle inequality for the Wasserstein metric,

$\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{L}}^{\eta},\mu^{\eta}\right)^{2}\right]\leq
2\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu_{\Lambda_{L}}^{\eta}\right)^{2}\right]+2\mathbb{E}\left[W_{\Omega}\left(\mu_{\Lambda_{2^{n}}}^{\eta},\mu^{\eta}\right)^{2}\right]\\\
\leq CL^{-\alpha}.$

The proof of (4.12) is complete.

The proofs of the translation covariance (1.9) and the
Dobrushin–Lanford–Ruelle (DLR) equation (1.10) for the Gibbs states
$\mu^{\eta}$ are a notational modification of the corresponding proofs for the
gradient Gibbs states $\mu_{\nabla}^{\eta}$ established in Section 3.2. We
thus omit the details.

We finally prove the quantitative coupling estimate (1.8). Using (4.10) and
the gluing lemma (see [50, Chapter 1]), we obtain that there exists a coupling
$\pi^{\eta}$ between the measures
$\left(\mu_{\Lambda_{2^{k}L}}^{\eta}\right)_{k\geq 0}$ such that

(4.14) $\sup_{x\in
2^{k-1}L}\mathbb{E}\left[\left\langle\left|\phi_{2^{k}L}(x)-\phi_{2^{k+1}L}(x)\right|^{2}\right\rangle_{\pi^{\eta}}\right]\leq
C(2^{k}L)^{-\alpha},$

where $\left\langle\cdot\right\rangle_{\pi^{\eta}}$ denotes the expectation
with respect to $\pi^{\eta}$ and $\phi_{2^{k}L}(x)$ is the height of the
surface at the vertex $x$ for the measure $\mu_{\Lambda_{2^{k}L}}^{\eta}$. The
inequality (3.11) implies that, for almost every $\eta$, the random variables
$\phi_{2^{k}L}$ converge pointwise almost surely with respect to the measure
$\pi^{\eta}$ to a random variable $\phi_{\infty}$ distributed according to
$\mu^{\eta}$. Summing the inequality (4.14) over the integers $k\in\mathbb{N}$
and using the triangle inequality as in (3.12) completes the proof of (1.8). ∎

##  Appendix A Gradient of the Neumann Green’s function

In this section, we establish the pointwise estimate on the gradient of the
Neumann Green’s function which is used in the proof of Proposition 4.4. The
statement and proof are well-known in the literature and given here for the
reader’s convenience.

######  Proposition A.1.

Fix a dimension $d\geq 3$. There exists a constant $C:=C(d)<\infty$ such that
the following holds. Fix $x\in\mathbb{Z}^{d}$, $\ell\in\mathbb{N}$ and set
$\Lambda=x+\Lambda_{\ell}$. There exists a function
$\mathbf{f}_{\Lambda}:\vec{E}\left(\Lambda\right)\to\mathbb{R}$ satisfying the
following properties: for any function $\phi:\Lambda\to\mathbb{R}$,

(A.1)
$\sum_{e\in\vec{E}\left(\Lambda\right)}\mathbf{f}_{\Lambda}(e)\nabla\phi(e)=\phi(x)-\frac{1}{\left|\Lambda\right|}\sum_{x\in\Lambda}\phi(y)\hskip
14.22636pt\mbox{and}\hskip 14.22636pt\forall
e\in\vec{E}\left(\Lambda\right),~{}\left|\mathbf{f}_{\Lambda}(e)\right|\leq
C.$

###### Proof.

The usual version of the Sobolev inequality on $\mathbb{Z}^{d}$ (which follows
easily from the one on ${\mathbb{R}^{d}}$ by interpolation) states that there
exists a constant $C:=C(d)<\infty$ such that, for any function
$g:\Lambda\to\mathbb{R}$ satisfying $\sum_{x\in\Lambda}g(x)=0$,

$\left(\sum_{y\in\Lambda}\left|g(y)\right|^{\frac{2d}{d-2}}\right)^{\frac{d-2}{2d}}\leq
C\left(\sum_{e\in\vec{E}\left(\Lambda\right)}\left|\nabla
g(e)\right|^{2}\right)^{\frac{1}{2}}.$

In the discrete setting, we additionally have, for any function
$g:\Lambda\to\mathbb{R}$ and any vertex $z\in\Lambda$,

$|g(z)|\leq\left(\sum_{y\in\Lambda}\left|g(y)\right|^{\frac{2d}{d-2}}\right)^{\frac{d-2}{2d}}.$

A combination of the two previous displays yields, for any function
$g:\Lambda\to\mathbb{R}$ satisfying $\sum_{y\in\Lambda}g(y)=0$,

(A.2) $\sup_{z\in\Lambda}\left|g(z)\right|\leq
C\left(\sum_{e\in\vec{E}\left(\Lambda\right)}\left|\nabla
g(e)\right|^{2}\right)^{\frac{1}{2}}.$

Consider the variational problem

$\inf_{\begin{subarray}{c}g:\Lambda\to\mathbb{R}\\\
\sum_{y\in\Lambda}g(y)=0\end{subarray}}\frac{1}{2}\sum_{e\in\vec{E}\left(\Lambda\right)}\left|\nabla
g(e)\right|^{2}-g(x).$

From (A.2), the minimum exists and is uniquely attained. We denote it by
$G_{\Lambda}$ and set $\mathbf{f}_{\Lambda}=\nabla G_{\Lambda}$. By the first
variations, we have, for any function $g:\Lambda\to\mathbb{R}$ satisfying
$\sum_{y\in\Lambda}g(y)=0$,

(A.3) $\sum_{e\in\vec{E}\left(\Lambda\right)}\mathbf{f}_{\Lambda}(e)\nabla
g(e)=g(x).$

Consider a function $\phi:\Lambda\to\mathbb{R}$. Applying the identity (A.3)
with the function $g:=\phi-|\Lambda|^{-1}\sum_{y\in\Lambda}\phi(y)$ shows the
first identity of (A.1). We then use the inequality (A.2), combine it with the
identity (A.3) (applied with the function $g=G_{\Lambda}$) and obtain

(A.4) $\sup_{z\in\Lambda}|G_{\Lambda}(z)|\leq
C\left(\sum_{e\in\vec{E}\left(\Lambda\right)}\left|\nabla
G_{\Lambda}(e)\right|^{2}\right)^{\frac{1}{2}}=C|G_{\Lambda}(x)|^{\frac{1}{2}}\leq
C\left(\sup_{z\in\Lambda}|G_{\Lambda}(z)|\right)^{\frac{1}{2}}.$

A rearrangement of the inequality (A.4) yields
$\sup_{z\in\Lambda}|G_{\Lambda}(z)|\leq C$. Consequently, for any edge
$e=\\{y,z\\}\in E\left(\Lambda\right),$
$|\mathbf{f}_{\Lambda}(e)|=\left|G_{\Lambda}(y)-G_{\Lambda}(z)\right|\leq 2C$.
∎

##  Appendix B De Giorgi-Nash-Moser regularity

In this section, we prove the version of the De Giorgi-Nash-Moser regularity
estimate for solutions of parabolic equations stated in Proposition 2.2. The
result is again standard and well-known in the literature (see [31, Chapter
8], [34, Chapter 3]), and we present a short proof for completeness. We
mention that a closely related version of the result was established in the
discrete setting considered here and used to study the $\nabla\phi$-model
(without disorder) in [30, Appendix B].

We first collect two preliminary results: the mean value inequality and the
increase of oscillations for solutions of parabolic equations. In the
continuous setting, these results can be found in [40]. In the discrete
setting, we refer to Proposition 3.2 and Proposition 3.4 of [24]. In both
statements, we fix an environment
$\mathbf{a}:[0,\infty]\times\mathbb{Z}^{d}\to\mathbb{R}$ which is assumed to
satisfy the uniform ellipticity assumption $c_{-}\leq\mathbf{a}\leq c_{+}$.

######  Proposition B.1 (Mean value inequality, Proposition 3.2 of [24]).

There exists a constant $C:=C(d,c_{+},c_{-})<\infty$ such that, for any $L\geq
1$, $t\geq 0$, $x\in\mathbb{Z}^{d}$, and any function
$u:[t,t+4L^{2}]\times(x+\Lambda_{2L})\to\mathbb{R}$ solution of the parabolic
equation

$\partial_{t}u-\nabla\cdot\mathbf{a}\nabla
u=0~{}\mbox{in}~{}[t,t+4L^{2}]\times(x+\Lambda_{2L}),$

one has the estimate

$\sup_{[t+3L^{2},t+4L^{2}]\times(x+\Lambda_{L})}u^{2}\leq\frac{C}{L^{d+2}}\int_{t}^{t+4L^{2}}\sum_{y\in
x+\Lambda_{2L}}\left|u(t,y)\right|^{2}\,dt.$

In the next proposition, we use the notation $\operatorname*{osc}$ to refer to
the oscillation of a function $u$ defined by the formula

(B.1)
$\operatorname*{osc}_{[s,t]\times(x+\Lambda_{L})}u:=\sup_{[s,t]\times(x+\Lambda_{L})}u-\inf_{[s,t]\times(x+\Lambda_{L})}u.$

######  Proposition B.2 (Increase of oscillations, Proposition 3.4 of [24]).

There exists a constant $\theta:=\theta(d,c_{+},c_{-})>1$ such that, for any
$L\geq 1$, $t\geq 0$, $x\in\mathbb{Z}^{d}$ and any solution
$u:[t,t+4L^{2}]\times(x+\Lambda_{2L})\to\mathbb{R}$ of the parabolic equation

$\partial_{t}u-\nabla\cdot\mathbf{a}\nabla
u=0~{}\mbox{in}~{}[t,t+4L^{2}]\times(x+\Lambda_{2L}),$

one has the estimate

(B.2)
$\operatorname*{osc}_{[t,t+4L^{2}]\times(x+\Lambda_{2L})}u\geq\theta\operatorname*{osc}_{[t+3L^{2},t+4L^{2}]\times(x+\Lambda_{L})}u.$

###### Proof of Proposition 2.2.

The proof of Proposition 2.2 from Proposition B.1 and Proposition B.2 is
direct and is for instance given in [45, Theorem 1.31] or in [24, (2.4)]. As
Proposition 2.2 cannot be immediately obtained from the two aforementioned
results, we briefly rewrite the argument.

By iterating Proposition B.2, we obtain that there exists a constant
$C:=C(d,c_{+},c_{-})<\infty$ such that, for any time $t\in[L^{2}/2,L^{2}]$ and
any edge $e=\\{x,y\\}\in E\left(\Lambda_{L/2}\right)$,

$\left|\nabla
u(t,e)\right|=\left|u(t,x)-u(t,y)\right|\leq\operatorname*{osc}_{[t-4,t]\times(x+\Lambda_{2})}u\leq
C\theta^{-\frac{\ln L}{\ln
2}}\operatorname*{osc}_{[t-L^{2}/16,t]\times(x+\Lambda_{L/4})}u.$

Setting $\alpha:=\ln\theta/\ln 2>0$, using the definition of the oscillation
(B.1) and applying the mean value inequality stated in Proposition B.1, we
obtain

$\displaystyle\left|\nabla
u(t,e)\right|\leq\frac{C}{L^{\alpha}}\operatorname*{osc}_{[t-L^{2}/16,t]\times(x+\Lambda_{L/4})}u$
$\displaystyle\leq\frac{C}{L^{\alpha}}\sup_{[t-L^{2}/16,t]\times(x+\Lambda_{L/4})}|u|$
$\displaystyle\leq\frac{C}{L^{\alpha}}\left(\frac{1}{L^{d+2}}\int_{t-L^{2}/4}^{t}\sum_{y\in
x+\Lambda_{L/2}}\left|u(s,y)\right|^{2}\,ds\right)^{\frac{1}{2}}$
$\displaystyle\leq\frac{C}{L^{\alpha}}\left(\frac{1}{L^{d+2}}\int_{0}^{L^{2}}\sum_{y\in\Lambda_{L}}\left|u(s,y)\right|^{2}\,ds\right)^{\frac{1}{2}}.$

The proof of Proposition 2.2 is complete. ∎

## References

  * [1] S. Adams, R. Kotecký, and S. Müller.  Strict convexity of the surface tension for non-convex potentials.  arXiv preprint arXiv:1606.09541, 2016.
  * [2] M. Aizenman, M. Harel, and R. Peled.  Exponential decay of correlations in the $2D$ random field Ising model.  J. Stat. Phys., 180(1-6):304–331, 2020.
  * [3] M. Aizenman and R. Peled.  A power-law upper bound on the correlations in the $2D$ random field Ising model.  Comm. Math. Phys., 372(3):865–892, 2019.
  * [4] M. Aizenman and J. Wehr.  Rounding of first-order phase transitions in systems with quenched disorder.  Physical Review Letters, 62(21):2503, 1989.
  * [5] M. Aizenman and J. Wehr.  Rounding effects of quenched randomness on first-order phase transitions.  Comm. Math. Phys., 130(3):489–528, 1990.
  * [6] S. Armstrong and W. Wu.  ${C}^{2}$ regularity of the surface tension for the $\nabla\phi$ interface model.  arXiv preprint arXiv:1909.13325, 2019.
  * [7] M. Biskup and R. Kotecký.  Phase coexistence of gradient Gibbs states.  Probab. Theory Related Fields, 139(1-2):1–39, 2007.
  * [8] M. Biskup and H. Spohn.  Scaling limit for a class of gradient fields with nonconvex potentials.  Ann. Probab., 39(1):224–251, 2011.
  * [9] A. Bovier and C. Külske.  A rigorous renormalization group method for interfaces in random media.  Rev. Math. Phys., 6(3):413–496, 1994.
  * [10] A. Bovier and C. Külske.  There are no nice interfaces in (2+1)-dimensional SOS models in random media.  J. Stat. Phys., 83(3-4):751–759, 1996.
  * [11] H. J. Brascamp and E. H. Lieb.  On extensions of the Brunn-Minkowski and Prékopa-Leindler theorems, including inequalities for log concave functions, and with an application to the diffusion equation.  J. Functional Analysis, 22(4):366–389, 1976.
  * [12] H. J. Brascamp and E. H. Lieb.  Some inequalities for Gaussian measures and the long-range order of the one-dimensional plasma.  In Inequalities, pages 403–416. Springer, 2002.
  * [13] H. J. Brascamp, E. H. Lieb, and J. L. Lebowitz.  The statistical mechanics of anharmonic lattices.  Bull. Inst. Internat. Statist., 46(1):393–404 (1976), 1975.
  * [14] J. Bricmont and A. Kupiainen.  Phase transition in the $3d$ random field Ising model.  Comm. Math. Phys., 116(4):539–572, 1988.
  * [15] D. Brydges and H.-T. Yau.  Grad $\phi$ perturbations of massless Gaussian fields.  Comm. Math. Phys., 129(2):351–392, 1990.
  * [16] S. Chatterjee.  On the decay of correlations in the random field Ising model.  Comm. Math. Phys., 362(1):253–267, 2018.
  * [17] C. Cotar and J.-D. Deuschel.  Decay of covariances, uniqueness of ergodic component and scaling limit for a class of $\nabla\phi$ systems with non-convex potential.  Ann. Inst. Henri Poincaré Probab. Stat., 48(3):819–853, 2012\.
  * [18] C. Cotar, J.-D. Deuschel, and S. Müller.  Strict convexity of the free energy for a class of non-convex gradient models.  Comm. Math. Phys., 286(1):359–376, 2009.
  * [19] C. Cotar and C. Külske.  Existence of random gradient states.  Ann. Appl. Probab., 22(4):1650–1692, 2012.
  * [20] C. Cotar and C. Külske.  Uniqueness of gradient Gibbs measures with disorder.  Probab. Theory and Related Fields, 162(3-4):587–635, 2015.
  * [21] P. Dario, M. Harel, and R. Peled.  Random-field random surfaces.  arXiv preprint arXiv:2101.01711, 2021.
  * [22] E. De Giorgi.  Sulla differenziabilità e l’analiticità delle estremali degli integrali multipli regolari.  Mem. Accad. Sci. Torino. Cl. Sci. Fis. Mat. Nat. (3), 3:25–43, 1957\.
  * [23] T. Delmotte.  Inégalité de Harnack elliptique sur les graphes.  In Colloquium Mathematicum, volume 72, pages 19–37. Instytut Matematyczny Polskiej Akademii Nauk, 1997.
  * [24] T. Delmotte and J.-D. Deuschel.  On estimating the derivatives of symmetric diffusions in stationary random environment, with applications to $\nabla\phi$ interface model.  Probab. Theory Related Fields, 133(3):358–390, 2005.
  * [25] J.-D. Deuschel, G. Giacomin, and D. Ioffe.  Large deviations and concentration properties for $\nabla\phi$ interface models.  Probab. Theory Related Fields, 117(1):49–111, 2000.
  * [26] J. Ding and M. Wirth.  Correlation length of two-dimensional random field Ising model via greedy lattice animal.  arXiv preprint arXiv:2011.08768, 2020.
  * [27] J. Ding and J. Xia.  Exponential decay of correlations in the two-dimensional random field Ising model.  Invent. Math., pages 1–47, 2021.
  * [28] T. Funaki.  Stochastic interface models.  In Lectures on Probability Theory and Statistics, volume 1869 of Lecture Notes in Math., pages 103–274. Springer, Berlin, 2005.
  * [29] T. Funaki and H. Spohn.  Motion by mean curvature from the Ginzburg-Landau $\nabla\phi$ interface model.  Comm. Math. Phys., 185(1):1–36, 1997.
  * [30] G. Giacomin, S. Olla, and H. Spohn.  Equilibrium fluctuations for $\nabla\phi$ interface model.  Ann. Probab., 29(3):1138–1172, 2001.
  * [31] D. Gilbarg and N. S. Trudinger.  Elliptic partial differential equations of second order.  Classics in Mathematics. Springer-Verlag, Berlin, 2001.  Reprint of the 1998 edition.
  * [32] A. Gloria and F. Otto.  An optimal variance estimate in stochastic homogenization of discrete elliptic equations.  Ann. Probab., 39(3):779–856, 2011.
  * [33] A. Gloria and F. Otto.  An optimal error estimate in stochastic homogenization of discrete elliptic equations.  Ann. Appl. Probab., 22(1):1–28, 2012.
  * [34] Q. Han and F. Lin.  Elliptic partial differential equations, volume 1.  American Mathematical Soc., 2011.
  * [35] Y. Imry and S.-K. Ma.  Random-field instability of the ordered state of continuous symmetry.  Physical Review Letters, 35(21):1399, 1975.
  * [36] C. Külske and E. Orlandi.  A simple fluctuation lower bound for a disordered massless random continuous spin model in $D=2$.  Electron. Comm. Probab., 11:200–205, 2006.
  * [37] C. Külske and E. Orlandi.  Continuous interfaces with disorder: Even strong pinning is too weak in two dimensions.  Stochastic Process. Appl., 118(11):1973–1981, 2008.
  * [38] J. Miller.  Fluctuations for the Ginzburg-Landau $\nabla\phi$ interface model on a bounded domain.  Comm. Math. Phys., 308(3):591–639, 2011.
  * [39] J. Moser.  On Harnack’s theorem for elliptic differential equations.  Comm. Pure and Appl. Math., 14(3):577–591, 1961.
  * [40] J. Moser.  A Harnack inequality for parabolic differential equations.  Comm. Pure and Appl. Math., 17(1):101–134, 1964.
  * [41] A. Naddaf and T. Spencer.  On homogenization and scaling limit of some gradient perturbations of a massless free field.  Comm. Math. Phys., 183(1):55–84, 1997.
  * [42] A. Naddaf and T. Spencer.  Estimates on the variance of some homogenization problems, 1998, unpublished preprint.
  * [43] J. Nash.  Continuity of solutions of parabolic and elliptic equations.  American Journal of Mathematics, 80(4):931–954, 1958.
  * [44] S. Sheffield.  Random surfaces.  Astérisque, (304):vi+175, 2005.
  * [45] D. W. Stroock and W. Zheng.  Markov chain approximations to symmetric diffusions.  In Ann. Inst. Henri Poincaré Probab. Stat., volume 33, pages 619–649. Elsevier, 1997.
  * [46] A. C. van Enter and C. Külske.  Nonexistence of random gradient Gibbs measures in continuous interface models in $d=2$.  Ann. Appl. Probab., 18(1):109–119, 2008.
  * [47] A. C. van Enter and S. B. Shlosman.  First-order transitions for n-vector models in two and more dimensions: Rigorous proof.  Physical Review Letters, 89(28):285702, 2002.
  * [48] A. C. van Enter and S. B. Shlosman.  Provable first-order transitions for nonlinear vector and gauge models with continuous symmetries.  Comm. Math. Phys., 255(1):21–32, 2005.
  * [49] Y. Velenik.  Localization and delocalization of random interfaces.  Probab. Surv., 3:112–169, 2006.
  * [50] J. Villain.  Theory of one-and two-dimensional magnets with an easy magnetization plane. II. The planar, classical, two-dimensional magnet.  Journal de Physique, 36(6):581–590, 1975.
  * [51] C. Villani.  Optimal transport: old and new, volume 338.  Springer Science & Business Media, 2008.