We prove a shape theorem for internal diffusion limited aggregation on mated-CRT maps, a family
of random planar maps which approximate Liouville quantum gravity (LQG) surfaces.
The limit is an LQG harmonic ball, which we constructed in a companion paper.
We also prove an analogous result for the divisible sandpile.

§ INTRODUCTION

Left: IDLA on the Tutte embedding of a finite $\gamma = \sqrt{2}$ mated-CRT map with boundary, with $10^4$ vertices.
Walkers are added until the first time one of them hits a boundary vertex.
Cells are colored according to the first time one of their vertices is hit by a random walk.
Right: Illustration of the graph upon which IDLA is run. The finite mated-CRT map with boundary can be constructed from a pair of conditioned Brownian motions via a similar procedure as in Definition <ref>. Away from the boundary, it locally looks like the infinite mated-CRT map considered in this paper, but its Tutte embedding is easier to define and simulate. See Remark <ref> for further discussion.

In this article we study the large-scale behavior of a random growth model (IDLA) on a random graph (the mated-CRT map)
which approximates a random fractal surface embedded in the plane (Liouville quantum gravity, LQG). We establish, for each $\gamma \in (0,2)$, the convergence in probability of the IDLA cluster on the $\gamma$-mated-CRT map to a so-called $\gamma$-LQG harmonic ball (Theorem <ref>). Harmonic balls are random subsets of the plane characterized by a mean-value property for harmonic functions with respect to the LQG measure (Definition <ref>). In a companion work <cit.>, harmonic balls are constructed via Hele-Shaw flow and it is shown that they are neither Lipschitz domains nor LQG metric balls.

§.§ Background

We now briefly and somewhat informally describe the aforementioned objects, delaying the more technical definitions to Section <ref> below. Minimal prior knowledge of LQG is needed to read this article.

§.§.§ Internal DLA

Internal diffusion limited aggregation (Internal DLA or IDLA) was first introduced to model the physical process of `material removal'
by Meakin-Deutch in <cit.>. It was also independently introduced by Diaconis-Fulton in  <cit.>  as a special case of a commutative algebra.

In its simplest form, IDLA produces an increasing sequence of subsets of the vertex set of a graph $G$ (often taken to be a lattice such as $\Z^2$) as follows. Fix a root vertex $o$ for $G$. Send out an explorer starting from $o$ and allow it to travel according to a simple random walk on $G$ until it reaches an unoccupied vertex, at which point it stops and occupies that vertex. This process is repeated, each explorer starts at the root and walks randomly until it reaches an unoccupied vertex. For $n \geq 1$, the time $n$ IDLA cluster is the collection of occupied vertices after $n $ iterations of this.

Following remarkable work by several groups of mathematicians, the large $n$ behavior of the IDLA clusters is now well understood on several families of graphs: $\Z^d$
Cayley graphs of groups with polynomial <cit.> and  exponential <cit.> growth, supercritical percolation clusters <cit.>, Sierpinski gasket graphs <cit.>, comb lattices <cit.>, and cylinders <cit.>. In most of these cases, the scaling limit of the IDLA clusters is described by a metric ball in the corresponding ambient space (the main result of paper will give an example where the scaling limit of IDLA is not a metric ball).

We would like to examine what happens when $G$ is taken to be a random graph embedded in the plane. Specifically, instead of taking some deterministic lattice
which approximates Euclidean space, we take $G$ to be a graph which approximates a Liouville quantum gravity (LQG) surface.

§.§.§ Liouville quantum gravity

LQG was introduced by Polyakov in the 1980s in the context of bosonic string theory <cit.> as a model of a `random two-dimensional Riemannian manifold'. LQG is too rough to be defined as a manifold in the rigorous sense, although the following loose description can be made precise via various regularization procedures.
For $\gamma \in (0,2)$ and a domain $D\subset\C$, a $\gamma$-LQG surface parameterized by $D$ is the two-dimensional Riemannian manifold with Riemannian metric tensor $e^{\gamma h} \, (dx^2+dy^2)$, where $dx^2+dy^2$ is the Euclidean metric tensor and $h$ is a variant of the Gaussian free field (GFF) on $D$.

Of particular relevance to this work is the associated volume form, the $\gamma$-Liouville measure or $\gamma$-LQG measure. This is a Radon measure
on $D$ which is (informally) given by

\begin{equation} \label{eq:lqg-measure}
\mu_h(dz) = e^{\gamma h(z)} \, d z,
\end{equation}

where $dz$ denotes two-dimensional Lebesgue measure <cit.>. We will give a precise definition of $\mu_h$ and the relevant variants of the GFF in Section <ref>. For now,	the reader only needs to know that $\mu_h$ is a random non-atomic, locally finite Radon measure on $D$ which assigns positive mass to every open subset of $D$ and is mutually singular with respect to Lebesgue measure.
The works <cit.> showed that an LQG surface can also be endowed with a canonical random metric, but we will not need this metric in the present paper.

§.§.§ Random planar maps and the mated-CRT map

Left: A sample construction of the mated-CRT map. We draw the graphs of $C-L$ and $ R$ in the interval $[0, 12 {\eps}]$
where $C$ is a large constant chosen so that the graphs do not intersect. Each vertical strip containing the interval $[a-{\eps}, a]$ for $a \in {\eps} \Z$  corresponds to a vertex
in the mated-CRT map. The adjacency condition (<ref>) for $L$ (resp. $R$) corresponds to two vertices $a,b \in {\eps} \Z$ being adjacent if there is a horizontal line segment above the graph of $C-L$ (resp. below the graph of $ R$) which intersects the graph only in the vertical strips which contain the intervals $[a-{\eps}, a]$ and $[b-{\eps},b]$.
In the figure, we have illustrated, for each pair $(a,b)$ for which this adjacency condition holds, the lowest (resp. highest) such horizontal line segment in green.
Right: A planar embedding of the mated-CRT map on the left, under which it is a triangulation. Edges arising from the adjacency condition for $L$ (resp. $R$) which do not join consecutive vertices are shown in red (resp. blue). Edges joining consecutive vertices are shown in black. A similar illustration appeared as <cit.>.

LQG is widely expected, and in some cases rigorously known, to describe the scaling limit of
random planar maps. A planar map is a graph embedded in $\C$ in such a way that no two edges cross, viewed modulo orientation preserving homeomorphisms $\C \rightarrow \C$.
Uniform random planar maps (including uniform triangulations, quadrangulations, etc.) converge to $\sqrt{8/3}$-LQG surfaces in the Gromov-Hausdorff sense <cit.> and, in the case of triangulations, when embedded via the so-called Cardy embedding <cit.>. Similar convergence results are expected to hold for various types of non-uniform random planar maps toward $\gamma$-LQG with $\gamma\not=\sqrt{8/3}$. For example, random planar maps sampled with probability proportional to the Ising model partition function are expected to converge to $\sqrt 3$-LQG (see, e.g., <cit.> for more on such conjectures).

Mated-CRT maps are a one-parameter family of random planar maps, indexed by $\gamma\in(0,2)$, whose connection to LQG is better understood than for most other random planar maps. In particular, it was proven by Gwynne-Miller-Sheffield <cit.> that the re-scaled counting measure on vertices of the $\gamma$-mated-CRT  map converges to the $\gamma$-LQG measure when the map is embedded into $\C$ via the Tutte embedding (i.e., the embedding where the position of each vertex is the average of the positions of its neighbors). Building on this, it was also shown by Berestycki-Gwynne in <cit.> that random walk on the mated-CRT map converges to so-called Liouville Brownian motion, see Section <ref>
for the definition.

The mated-CRT map is built via a pair of correlated linear Brownian motions as in Figure <ref>.

Fix $\gamma \in (0,2)$ and consider a two dimensional Brownian motion $(L_t,R_t)_{t\in\R}$ started at the origin with correlation coefficient given by $ - \cos ( \pi \gamma^2/ 4)$.  For ${\eps} > 0$, the $\gamma$-mated-CRT map with cell size ${\eps}$ associated to $(L_t, R_t)_{t\in\R}$ is the graph ${\Geps}$ whose vertex set is ${\eps} \Z = {\VGeps}$, with two vertices $a < b$ connected by an edge if and only if

\begin{equation} \label{eq:edge-map}
\left( \inf_{t \in [a - {\eps},a]} X_t \right) \vee \left( \inf_{t \in [b - {\eps}, b]} X_t \right)\le \inf_{t \in [a, b - {\eps}]} X_t
\end{equation}

where $X$ can be either $L$ or $R$. There are two edges connecting $a$ and $b$ if $|b-a| > {\eps}$ and (<ref>)
is satisfied for both $L$ and $R$.  See Figure <ref> for an illustration and a description of the planar map structure of $\mathcal {\Geps}$.

By Brownian scaling, the law of ${\Geps}$ (as a graph) does not depend on ${\eps}$. However, ${\eps}$ is a convenient scaling parameter when talking about limits, as we will see just below.
The definition of the mated-CRT map is a semi-discrete generalization of so-called mating-of-trees bijections between random planar maps and random walk excursions.
For instance, the Mullin bijection <cit.> is a discrete analog of Definition <ref> which produces a planar map decorated by a spanning tree from a two-dimensional random walk excursion.  See <cit.> for a survey of mating-of-trees bijections and their applications.

The main reason why the mated-CRT map is more tractable than other random planar maps is that, due to the results of <cit.>, it has an a priori embedding into $\C$ defined in terms of LQG and Schramm-Loewner evolution (SLE). We will now discuss this embedding.

§.§.§ SLE/LQG embedding

In this paper, we work with an embedding of the mated-CRT map constructed via a space-filling Schramm-Loewner evolution (SLE) curve. We remark, however, that by <cit.>, our results extend
to the so-called Tutte (or harmonic) embedding of the mated-CRT map — see Figure <ref> and Remark <ref>.

Let $h$ be the random distribution on $\C$ corresponding to a $\gamma$-quantum cone (a slight modification of the whole-plane GFF).
Let $\mu_h$ be the $\gamma$-LQG measure associated to $h$.
Let $\eta$ be a whole-plane space-filling SLE, independent of $h$, with parameter $\kappa = 16/\gamma^2 \in (4,\infty)$.
We parameterize $\eta$ so that $\eta(0) = 0$ and $\mu_h(\eta([a,b])) = b-a$ for every $a,b\in\R$ with $a < b$.
We will review the definitions of $h$ and $\eta$ in Section <ref> below.
For now, the unfamiliar reader can just think $\eta$ as a random continuous, space-filling curve in $\C$ from $\infty$ to $\infty$ parameterized by a locally finite Radon measure.

For $\gamma \in (0,2)$ and ${\eps} > 0$, let ${\Geps}$ be the ${\eps}$-mated CRT map with vertex set ${\VGeps}$ as defined in Definition <ref>.
The SLE/LQG embedding is defined as follows. Each vertex $a \in {\VGeps} = {\eps}\Z$ corresponds to a cell,

\begin{equation} \label{eq:sle-cell}
\Heps{a} := \eta([a-{\eps}, a]) .
\end{equation}

Two distinct vertices $a,b \in {\VGeps}$ are connected by one (resp. two) edges if and only if $\Heps{a} \cap \Heps{b}$ has one
(resp. two) connected components which are not singletons.  This provides an embedding
of ${\Geps}$ into $\C$ by sending each vertex $a \in {\VGeps}$ to the point $\eta(a) \in \Heps{a}$ then drawing edges between adjacent vertices in such a way that no two edges cross (the precise locations of the edges is unimportant for our purposes).

The equivalence of the above SLE/LQG description of the mated-CRT map and the Brownian motion description, Definition <ref>,
is a consequence of the main result of the seminal work <cit.>. See, in particular, <cit.>. For an illustration of the SLE/LQG embedding, see Figure <ref>.

Top left: A space-filling SLE curve $\eta$ for $\kappa \geq 8$ divided into cells $\eta([a - {\eps}, a])$ for a collection of $a \in {\eps} \Z$. Top right: The
same curve with a blue path showing the order in which cells are traversed by $\eta$.	Bottom left: A point in each cell, corresponding to a vertex in ${\Geps}$, is displayed in red
and red edges are drawn to adjacent neighbors. Bottom right: The same as bottom left but with the black edges removed — this illustrates the embedding of ${\Geps}$ into $\C$.
A similar figure has appeared previously as <cit.>.

§.§.§ Harmonic balls

In this article we show that IDLA on the SLE/LQG embedding of the mated-CRT map converges to a harmonic ball.

Fix $\gamma \in (0,2)$ and let $\mu_h$ be the LQG measure (<ref>) where $h$ is a variant of the Gaussian free field.
A $\gamma$-LQG harmonic ball is a domain $\Lambda \subset \C$ which satisfies the following mean-value property for harmonic functions:

\begin{equation}
\int_{\Lambda} f(y) d \mu_h(y) = \mu_h(\Lambda) f(0), \quad \forall \mbox{$f:\C \to \R$ harmonic in a neighborhood of $\cl(\Lambda)$}.
\end{equation}

In a companion paper <cit.>, it is shown that if $\gamma\in(0,2)$ and $h$ is a suitable variant of the GFF on $\C$, then there is a unique growing family of LQG harmonic balls $\{\Lambda_t\}_{t>0}$, parameterized so that $\mu_h(\Lambda_t)=t$. These harmonic balls are constructed directly in the continuum via a certain optimization problem involving the measure $\mu_h$. This optimization problem, as well as several important properties of harmonic balls proved in <cit.>, are used in our proof. A statement of the facts about LQG harmonic balls which will be needed later appears in Theorem <ref> below. It is also shown in <cit.> that LQG harmonic balls are neither Lipschitz domains nor LQG metric balls.

§.§ Main results

Fix $\gamma \in (0,2)$ and for ${\eps} > 0$ let ${\Geps}$ be the ${\eps}$-mated CRT map as in Definition <ref>,
with vertices ${\VGeps}$ embedded into $\C$ via the SLE/LQG embedding as in Definition <ref>.
For a set $X \subset \C$ and $\delta > 0$ we denote the $\delta$-outer and $\delta$-inner neighborhoods of $X$ by

\begin{equation} \label{eq:outer-inner-neighborhoods}
\begin{aligned}
\Bdeltap(X) &:= \{ x \in \C : \dist(x, X) \leq \delta\} \\
\Bdeltam(X) &:= \{ y \in X : \dist(y, X^c) \geq \delta\}
\end{aligned}
\end{equation}

where $\dist$ denotes Euclidean distance.

§.§.§ Internal DLA

Our notation for the IDLA cluster on ${\Geps}$ is as follows. Begin with $A_{{\eps}}(0) = \emptyset$. For $m\geq 1$, inductively let $A_{{\eps}}(m)$ be the union of $A_{{\eps}}(m-1)$ and the first point at which a simple random walk on ${\Geps}$ started at the origin $ 0 \in {\VGeps} $
exits $A_{{\eps}}(m-1)$. By filling in the mated-CRT map cells corresponding to vertices in $A_{\eps}$ (Definition <ref>) we may map a subset $A_{{\eps}} \subset {\VGeps}$ to a subset of $\C$ via

\begin{equation} \label{eq:vertex-scaling}
\overline A_{{\eps}} := \bigcup_{a \in A_{{\eps}}} \Heps{a}.
\end{equation}

Let $\{\Lambda_t\}_{t\geq 0}$ be the growing family of $\gamma$-LQG harmonic balls with $\mu_h(\Lambda_t) = t$, as discussed just after Definition <ref>, with $h$ the field corresponding to the $\gamma$-quantum cone, as in Definition <ref>.
For each $\delta > 0$ and $t > 0$, it holds with probability tending to 1 as ${\eps}\to 0$ that the union of the cells corresponding to the time $\lfloor t{\eps}^{-1}\rfloor$-IDLA cluster approximates $\Lambda_t$ in the following sense:

\[
\Bdeltam(\Lambda_t)	\subset \overline A_{{\eps}}(\lfloor t {\eps}^{-1} \rfloor) \subset \Bdeltap(\Lambda_t),
\]

where $B_{\delta}^{\pm}$ are as in (<ref>).

One may  also prove a version of Theorem <ref> for the disk version of the mated-CRT map under the Tutte embedding, as depicted in Figure <ref>.
The quantum disk is the canonical LQG surface with the topology of the disk <cit.>.
In <cit.>, the authors introduce a disk version of the mated-CRT map with a distinguished boundary and define its Tutte embedding (the definition of the Tutte embedding is easier for planar maps with boundary, which is why one wants to consider the disk version of the mated-CRT map).
They also show that the Tutte embedding of the disk version of the mated-CRT map is close, in the uniform distance, to its SLE/LQG embedding.
Using the fact that the generalized function representing the quantum disk behaves like a GFF away from its boundary, one can extend the results of <cit.> to define a growing family of harmonic balls on the quantum disk started at a uniformly sampled interior point and stopped at the first time when they hit the boundary.
Via another absolute continuity argument, one can extend the result of Theorem <ref> to get that the IDLA clusters on the mated-CRT map with the disk topology, stopped at the first time they hit the boundary, converge to harmonic balls on the quantum disk under the SLE/LQG embedding.
We will not provide the details of these absolute continuity arguments in this paper, but see <cit.> and Theorem <ref> below for similar arguments.
Using the results of <cit.>, one gets that IDLA clusters on the disk version of the mated-CRT map also converge to harmonic balls under the Tutte embedding.

§.§.§ Divisible sandpile

Our proof of convergence of IDLA proceeds by first establishing an analogous result for an auxiliary particle system which may be thought of as the `expected value' of IDLA.

The divisible sandpile is a deterministic diffusion process on a graph $G$,
which, to the best of our knowledge, was first introduced by Zidarov <cit.>
and (independently) studied by Levine in his thesis <cit.>. Start with masses (real numbers) on vertices of a graph and redistribute the mass according to the following rule. Whenever a vertex has mass $t$ larger than one, it is unstable and topples by distributing the excess mass, $(t-1)$, to each of its neighbors equally. This process, called stabilizing, continues until every vertex has mass less than or equal to one. Although stabilizing may take infinite time, if we have a finite amount of mass on an infinite graph, the end configuration exists and does not depend on the order in which unstable sites are toppled <cit.>. From another perspective, the divisible sandpile is simply the Jacobi iterative algorithm, , <cit.>, for solving the Poisson equation.

In one variant of the model, called the single-source divisible sandpile, we start with a pile of mass $t > 0$ at a root vertex $o$ in $G$ and stabilize. The final collection of vertices which contain some mass is the cluster. In each IDLA convergence theorem recalled in Section <ref>, it has been shown that the limit shape of the single-source divisible sandpile cluster coincides with that of IDLA. This relationship was first observed and established on $\Z^d$ by Levine-Peres in <cit.>. See  <cit.> for a survey relating IDLA and the divisible sandpile.

In our setting, the IDLA cluster also shares the same limit shape as the divisible sandpile. Denote by $D_{{\eps}}(t)$ the divisible sandpile cluster starting with mass $t > 0$ at the origin in ${\Geps}$.

Theorem <ref> holds with the divisible sandpile cluster
started with mass $t {\eps}^{-1}$ at the origin,  $\overline D_{{\eps}}(t {\eps}^{-1})$, in place of the IDLA cluster
started with $\lfloor t {\eps}^{-1} \rfloor$ explorers at the origin,
$\overline A_{{\eps}}(\lfloor t {\eps}^{-1} \rfloor)$.

The divisible sandpile has an explicit description in terms of a discrete optimization problem, which we recall in Section <ref> below.
In particular, from this optimization problem it is easy to deduce that the divisible sandpile cluster enjoys a certain discrete mean-value property; the clusters are
discrete harmonic balls <cit.>. From this, it is no surprise that divisible sandpile clusters converge to continuum harmonic balls.

We would also like to point out that the divisible sandpile appears to be the only way to find discrete harmonic balls. In fact, this property of the divisible sandpile was used to establish logarithmic fluctuations of IDLA around its limit shape in <cit.>.

§.§ Open problems

We suggest several directions for further study. The most pressing problem is to extend our convergence result to other
random planar maps. This problem is mentioned in the companion work <cit.>.

Show that the scaling limit of IDLA and the divisible sandpile on random planar maps, other than mated-CRT maps, is also described by $\gamma$-LQG harmonic balls.
For example, on a uniform random planar map establish convergence of IDLA to a $\sqrt{8/3}$-LQG harmonic ball.

Our proofs of Theorems <ref> and <ref> rely heavily on the convergence of random walk on the mated-CRT map to Liouville Brownian motion. This fact has not been established for random planar maps besides the mated-CRT map. We speculate that a solution to Problem <ref> which does not use the convergence of random walk might rely on special symmetries of random planar maps, for example, the peeling process for uniform planar maps <cit.>.

The next two questions concern IDLA on the mated-CRT map. While we have established a qualitative convergence theorem much more quantitative results are known on $\Z^d$. Lawler <cit.> used fine estimates on the Green's function in $\Z^d$ to establish subdiffusive fluctuations of IDLA around its limit shape. Later, Asselah-Gaudillère <cit.> and independently Jerison-Levine-Sheffield <cit.> established logarithmic fluctuations. We do not expect logarithmic fluctuations in our setting, mainly because the cell sizes in the mated-CRT map have  a `multi-fractal' behavior. See Figure <ref>.

Provide bounds on the rate of convergence of IDLA on the mated-CRT map.

A close up of the interface in Figure <ref>.

After proving logarithmic fluctuations in $\Z^d$, in a remarkable work, Jerison-Levine-Sheffield showed that the scaling limit of the fluctuations themselves exists.
Namely, fluctuations of IDLA on $\Z^d$ converge weakly to a variant of the Gaussian free field <cit.>. Darrow extended
the analysis of Jerison-Levine-Sheffield on $\Z^2$ to the case of multiple sources <cit.>.

Do the fluctuations of IDLA on the mated-CRT map around its limit shape have a scaling limit?

§.§ Paper and proof outline

Our proof uses three main inputs; the pairing of IDLA and the divisible sandpile,
regularity properties of the mated-CRT map, and the convergence of random walk on the mated-CRT map to Liouville Brownian motion.

The connection between IDLA and the divisible sandpile was implicit in the first work on IDLA by Lawler-Bramson-Griffeath <cit.>
and made explicit in the thesis of Levine <cit.>. In particular, following the arguments of <cit.>,
a lower bound on the divisible sandpile cluster implies the same lower bound for the IDLA cluster.
The divisible sandpile, being deterministic and enjoying a variational description, , the least action principle for its odometer (see Section <ref> below for definitions), is more tractable than IDLA.

In the cases considered so far, , <cit.>,
convergence of the divisible sandpile is proved by explicitly bounding its odometer via sharp estimates on stopped Green's functions.  These sharp estimates are unavailable in our setting and hence we proceed with a somewhat more robust method. Roughly, we use the fact that the mated-CRT map has an elliptic Harnack inequality and random walk converges to Brownian motion to show that the divisible sandpile odometer converges to its continuum counterpart, the odometer for Hele-Shaw flow (the solution to the optimization problem used to construct harmonic balls).

We start by recalling in Section <ref> several mated-CRT map estimates and the invariance principle for random walk
proved in <cit.>. In this section we also record some basic properties of harmonic balls and Hele-Shaw flow on an LQG surface which we proved in a companion paper <cit.>.

In Section <ref>, we use the estimates proved in <cit.> to establish a
general tightness lemma (Lemma <ref>) for bounded functions with bounded graph Laplacian
on the mated-CRT map. We then use this lemma to establish the convergence of solutions to the discrete Dirichlet problem on the mated-CRT map to solutions of the continuum Dirichlet problem (Lemma <ref>) and convergence of the discrete Green's kernel to the continuum Green's function (Lemma <ref>).

In Section <ref>, we define the odometer for the divisible sandpile model on the mated-CRT map and state some basic properties of it which are proved in the Appendix, Section <ref>.
We then prove convergence of the divisible sandpile odometer to its continuum counterpart. Our proof method is similar to, but easier than, that used to show convergence of the Abelian sandpile odometer <cit.>.
In Section <ref> we also observe that the convergence of the divisible sandpile odometer implies a lower bound on the divisible sandpile cluster: , for fixed $t , \delta > 0$ it holds with probability tending to 1 as ${\eps}\to 0$ that

\begin{equation} \label{eq:div-lower-outline}
\Bdeltam(\Lambda_t) \subset \overline D_{{\eps}}( t {\eps}^{-1} ).
\end{equation}

In Section <ref> we apply the technique of Lawler-Bramson-Griffeath <cit.> together with the divisible sandpile lower bound (<ref>) and random walk estimates proved in <cit.> to prove an analogous lower bound for the IDLA cluster.
In fact, the lower bound holds for the stopped IDLA cluster where random walkers which exit (a mated-CRT map approximation of) $\Lambda_t$ are stopped.
By the Abelian property of IDLA, Lemma <ref>, the IDLA cluster can be built by first generating the stopped IDLA cluster and then restarting
the stopped random walkers.
Thus, an upper bound of the form

\begin{equation} \label{eq:idla-upper-outline}
\overline A_{{\eps}}(\lfloor t {\eps}^{-1} \rfloor) \subset \Bdeltap(\Lambda_t)
\end{equation}

will follow if one can show that the restarted random walks cannot travel further than a Euclidean distance of $\delta$ away from $\Lambda_t$.

In fact, if the limit shape has the correct size, , has volume equal to the total number of walkers,
and has a measure zero boundary, then by the lower bound and conservation of mass, the number of restarted walkers can be made arbitrarily small.  These two limit shape properties are true for an LQG harmonic ball, as proved in the companion paper <cit.>.  Therefore, the proof of the upper bound then reduces to controlling the behavior of order $\delta {\eps}^{-1}$ restarted random walks near the boundary of $\Lambda_t$.  We do this by combining the arguments
of <cit.> together with those from LQG theory.

Specifically, it was shown by Duminil-Copin-Lucas-Yadin-Yehudayoff <cit.> that if a graph is regular enough, namely is volume doubling, and has a metric which can be controlled by the Euclidean metric, then a small number of walkers cannot spread very far. The mated-CRT map does not enjoy these properties, even approximately, and thus we cannot apply the Duminil-Copin-Lucas-Yadin-Yehudayoff argument directly. Instead, we use methods from LQG theory to provide weak harmonic measure estimates for random walk on the mated-CRT map. A modification of the iterative argument of Duminil-Copin-Lucas-Yadin-Yehudayoff then yields the upper bound.  See the beginning of Section <ref> for a detailed outline of the argument.

In Section <ref> we use the harmonic measure estimates proved in Section <ref> to prove an upper bound on the divisible sandpile cluster. The argument here is a simplified, discrete version of the upper bound for LQG harmonic balls from <cit.>.

All of the aforementioned results were proved conditional on $t$ being small. In Section <ref> we combine the upper and lower bounds together with a scaling argument to remove this constraint on $t$ and complete the proof.

§.§ Notation and conventions

* Inequalities/equalities between functions/scalars are interpreted pointwise.

* For a set $D \subset \C$, $\partial D$ denotes its topological boundary, $\cl(D) = D \cup \partial D$ its closure,
and $\inte(D)$ its interior.

* For a set $A \subset {\VGeps}$, $\partial A$ is the set of vertices in ${\VGeps} \setminus A$ which are joined by an edge to
a vertex in $A$ and $\cl(A) = A \cup \partial A$ is its closure.
The set $\inte(A)$ are the vertices in $A$ which do not share an edge with a vertex in $A^c$.

* For two sets $X, Y \subset \C$, we write $X \Subset Y$ if $\cl(X) \subset Y$.

* $B_r(x)$ denotes the open ball of Euclidean radius $r > 0$ centered at $x \in \C$, when $x$ is omitted, the ball is centered at 0.

* For $0 < r_1 < r_2$, denote an open annulus centered at $z \in \C$ by

\begin{equation} \label{eq:annulus}
\A_{r_1, r_2}(z) = B_{r_2}(z) \setminus \cl(B_{r_1}(z))
\end{equation}

and $\A_{r_1, r_2} := \A_{r_1,r_2}(0)$.

* Let $\{E^{r}\}_{z > 0}$ be a one-parameter family of events.
We say that $E^r$ occurs with polynomially high probability as $r \to 0$
if there exists $p>0$ such that $\P[E^r] \geq 1 - O(r^p)$.

* For two sets $X, Y \subset \C$, we define $\dist(X,Y) = \inf_{x \in X, y \in Y} \dist(x,y)$ where $\dist$ denotes the Euclidean distance between two points.

* For a set $D \subset \C$, $C(D)$ denotes the set of continuous functions on  $D$.

§.§ Acknowledgments

A.B. thanks Charlie Smart and Bill Cooperman for useful discussions.
A.B. was partially supported by NSF grant DMS-2202940 and a Stevanovich fellowship.
E.G. was partially supported by a Clay research fellowship.

§ PRELIMINARIES

In this section we recall several preliminary results on the Gaussian free field, Liouville quantum gravity, space-filling SLE, harmonic balls, random walk on the mated-CRT map, and Liouville Brownian motion. Of note, we record the convergence of random walk on the mated-CRT map to Liouville Brownian motion <cit.>
and the regularity properties of $\gamma$-LQG harmonic balls established in <cit.>.

§.§ Gaussian free field

In this paper, we work with a particular kind of Gaussian free field $h$ which arises naturally from the scaling limit of random planar maps.
This field corresponds to an LQG surface called a $\gamma$-quantum cone. For our purposes we just define the field and give its scaling properties
— a complete definition of the surface is given in <cit.>.

We start by defining the whole-plane Gaussian free field (GFF) $h^{\C}$ which is the centered Gaussian random generalized
function on $\C$ with covariances

\begin{equation} \label{eq:cov-kernel}
\Cov(h^\C(x), h^\C(y)) := \log \frac{\max(|x|,1) \max(|y|,1)}{|x-y|} ,\quad\forall x,y\in\C .
\end{equation}

The GFF $h^{\C}$ is not defined pointwise as the covariance kernel in (<ref>) diverges to $\infty$ as $x\to y$.
However, for $z\in\C$ and $r>0$, one can define the average of $h^{\C}$ over the circle of radius $r$ centered at $z$, which we denote by $h^{\C}_r(z)$ <cit.>. 	The whole plane GFF is sometimes defined modulo additive constant. Our choice of covariance in (<ref>) corresponds to fixing this additive constant so that $h^\C_1(0) = 0$ (see, e.g., <cit.>).

The circle average of $h^{\C}$, $t \to h_{e^{-t}}(0)$ is a Brownian motion, <cit.>.
The field corresponding to the $\gamma$-quantum cone is constructed by enforcing its circle-average to correspond to a particular conditioned Brownian motion.
Fix $\gamma \in (0,2)$ and let

\begin{equation} \label{eq:Q}
Q := \frac{2}{\gamma} + \frac{\gamma}{2}.
\end{equation}

The circle average embedding of the $\gamma$-quantum cone is the distribution $h$ defined as follows.
Let $B$ be a standard linear Brownian motion and let $\hat{B}$ be an independent standard linear Brownian motion conditioned so that $\hat{B}_t + (Q-\gamma) t > 0$ for all $t > 0$.

\begin{equation}
A_t := \begin{cases}
B_t + \gamma t , \quad &t \geq 0 \\
\hat{B}_{-t} + \gamma t ,\quad &t < 0 .
\end{cases}
\end{equation}

The distribution $h$ is defined so that $t \to h_{e^{-t}(0)}$ has the same law as the process $A$ and $h - h_{|\cdot|}$ is independent
from $h_{|\cdot|}$ and has the same law as the analogous process for the whole-plane GFF, $h^{\C}$.

We note that the field $h$ in Definition <ref> has the property that $\sup \{ r > 0 : h_r(0) + Q \log r = 0 \} = 1$. Furthermore, it is immediate from Definition <ref> that $h$ restricted to the unit disk agrees in law with the corresponding restriction of a whole-plane GFF plus $ - \gamma \log |\cdot|$, normalized so that its circle average over the disk is zero:

\begin{equation} \label{eq:gff}
h|_{B_1} \overset{d}{=} ( h^{\C} - \gamma \log |\cdot| )|_{B_1} .
\end{equation}

The law of the whole-plane GFF, viewed modulo additive constant, is invariant under complex affine transformations of $\C$. This implies the following scaling property for $h^{\C}$,

\begin{equation} \label{eq:gff-scaling}
h^{\C} \overset{d}{=} h^{\C}(x\cdot +y) - h^{\C}_{|x|}(y), \quad \forall x \in \C \setminus \{0\}, \quad \forall  y \in \C.
\end{equation}

The distribution $h$ has a similar scale invariance. For $b > 0$, let

\begin{equation} \label{eq:cone-scale-radius}
R_b := \sup \{ r > 0 : h_r(0) + Q \log r = \frac{1}{\gamma} \log b\}.
\end{equation}

Then, by <cit.>,

\begin{equation} \label{eq:quantum-cone-scaling}
h(R_b \cdot) + Q \log R_b - \frac{1}{\gamma} \log b \overset{d}{=} h.
\end{equation}

§.§ Liouville quantum gravity

Liouville quantum gravity (LQG) is a one-parameter family of random fractal surfaces which were introduced by Polyakov in the 1980s in the context of bosonic string theory <cit.>.
We give some basic properties of LQG and refer the interested reader to the introductory texts <cit.>.

Let $\mu_h$ denote the $\gamma$-LQG area (Liouville) measure associated to $h$, the $\gamma$-quantum cone defined in (<ref>). One of the several ways of defining $\mu_h$ is as the a.s. weak limit

\begin{equation} \label{eq:measure-lim}
\mu_h = \lim_{\epsilon \rightarrow 0} \epsilon^{\gamma^2/2} e^{\gamma h_\epsilon(z)} \,dz
\end{equation}

where $dz$ denotes Lebesgue measure and $h_\epsilon(z)$ is the circle average <cit.>.
In fact, the measure $\mu_{\tilde h}$ can be constructed for any random generalized function $\tilde h$ of the form  $h + f$ where $f$ is a possibly random continuous function
and $h$ is the whole-plane GFF, defined in (<ref>). In particular, this includes the $\gamma$-quantum cone from Definition <ref>.
For later use, we record some basic properties of the LQG measure.
[LQG measure]
The LQG area measure $\mu_h$ satisfies the following properties.

* Radon measure. A.s., $\mu_h$ is a non-atomic Radon measure.

* Locality. For every deterministic open set $U \subset \C$,
$\mu_h(U)$ is given by a measurable function of
$h \vert_{U}$.

* Weyl scaling. A.s., $e^{\gamma f} \cdot \mu_h = \mu_{h + f}$
for every continuous function $f: \C \to \R$.

* Conformal covariance. A.s., the following is true. Let $U , \tilde U \subset \C$ be open and let $\phi$ be a conformal
map from $\tilde U$ to $U$. Then, with $Q$ as in (<ref>),

\begin{equation} \label{eq:measure-coord}
\mu_{h \circ \phi + Q \log |\phi'|}(A) = \mu_h(\phi(A)) \quad \mbox{for all Borel measurable $A \subset \tilde U$}.
\end{equation}

The first three properties in Fact <ref> are immediate from the definition (<ref>). The conformal covariance property was proven to hold a.s. for a fixed conformal map in <cit.> and extended to all conformal maps simultaneously in <cit.>.

§.§ Harmonic balls

A harmonic ball is a domain which satisfies the mean-value property for harmonic functions with respect to the LQG measure — Definition <ref>.
These balls are constructed in a companion work via the following optimization problem.
For each $t > 0$ and ball $B_r = B_r(0)$, let

\begin{equation} \label{eq:least-super-solution}
\overline w_t^{B_r} = \inf\{ w \in C(\cl(B_r)) : \Delta w \leq \mu_h \mbox{ in $B_r$}  \mbox{ and }  w \geq  - t G_{B_r}(0, \cdot) \mbox{ in $\cl(B_r)$}\},
\end{equation}

where $C(\cl(B_r))$ denotes the set of continuous functions on the closed ball, $\Delta w$ is interpreted in the distributional sense, $G_{B_r}$ is the Green's function for standard Brownian motion on $B_r$ with zero boundary conditions, and the infimum is pointwise.
The definition in (<ref>) is the obstacle problem formulation of Hele-Shaw flow <cit.> and may be thought of as the continuum
analogue of the odometer for the divisible sandpile, defined below in Section <ref>.

On an event of probability one, there exists a unique family of harmonic balls $\{\Lambda_t\}_{t > 0}$
satisfying the following properties:

* For each $t > 0$, $\mu_h(\Lambda_t) = t$, $\mu_h(\partial \Lambda_t) = 0$, and $\Lambda_t$ is equal to the interior of its closure.

* The domains $\Lambda_t$ are bounded, connected, open, contain the origin, increase continuously in $t$ (in the Hausdorff topology), and satisfy $\cap_{t > 0} \Lambda_t = \{0\}$.

* If, for some $r > 0$, $\Lambda_t \subset B_r$, then, for all $\rho > r$,

\[
\Lambda_t = \inte \left( \cl \left( \{ x \in B_\rho : \overline w_t^{B_\rho} > - t G_{B_\rho}(0, \cdot)\} \right) \right).
\]

Moreover, $\overline w_t^{B_\rho}$ is continuous and

\[
\begin{cases}
\Delta \overline w_t^{B_\rho} = \mu_h |_{\Lambda_t}  \quad \mbox{on $B_\rho$} \\
\overline w_t^{B_\rho} = 0 \quad \mbox{on $\partial B_\rho$}.
\end{cases}
\]

Strictly speaking, <cit.> considers the case where $h$ is not the quantum cone but rather a whole-plane GFF plus a log-singularity.
Transferring the proof to the quantum cone is a standard absolute continuity argument which we now provide.

On an event of probability one, there exists a unique family of harmonic balls $\{\Lambda_t\}_{0 < t < T}$
which satisfies the properties of Theorem <ref> (with $t$ restricted to be in $(0,T)$), where

\[
T = \sup \{ t > 0 : \Lambda_t \subset B_{1/2}\}.
\]

By (<ref>), the restrictions to the unit disk of the quantum cone field $h$ and the whole-plane GFF plus $-\gamma \log |\cdot|$ agree in law.
Since harmonic balls depend locally on the field <cit.>, this equality in law together with  <cit.> directly implies the existence of a family of harmonic balls $\{ \Lambda_t \}_{0 < t < T}$ satisfying the properties in Theorem <ref>.
Uniqueness of such a family is guaranteed by <cit.>
(that proposition is stated for a family of harmonic balls defined for all times, but the proof still works if the harmonic balls are only defined up to some fixed time).

We will now extend the above construction to all time using the scale invariance property of the $\gamma$-quantum cone.
The following is the quantum cone analog of <cit.> (which gives a similar statement for the LQG harmonic balls associated with the whole-plane GFF).

For $b >0$, define the random radius $R_b > 0$ as in (<ref>).
A.s., there exists a unique family of harmonic balls
$\{\Lambda_{t}\}_{0 < t < b T^b}$
such that $\{R_b^{-1} \Lambda_{b t}\}_{0 < t < T^b}$
has the same law as $\{\Lambda_t\}_{0 < t < T}$,
$T$ is as in Lemma <ref>,

\[
T^b = \sup\{ t > 0 : R_b^{-1} \Lambda_{b t} \subset B_{1/2} \}.
\]

In fact, we have the equality of joint laws

\[
\left( h, \{\Lambda_t\}_{0 < t < T} \right) \overset{d}{=}\left( h(R_b \cdot) + Q \log R_b - \frac{1}{\gamma} \log b, \{R_b^{-1} \Lambda_{b t}\}_{0 < t < T^b} \right).
\]

Fix $b > 0$ and let

\[
h^b := h(R_b \cdot) + Q \log R_b - \frac{1}{\gamma} \log b,
\]

so that, as in (<ref>), $h^b \overset{d}{=} h$. By Lemma <ref>, a.s. there exists a unique family of harmonic balls $\{\Lambda^b_t\}_{0 < t < T^b}$ which satisfies the properties of Theorem <ref> for the measure $\mu_{h^b}$ (with $t$ restricted to be in $(0,T^b)$), where

\[
T^b = \sup\{ t > 0 : \Lambda^b_{t} \subset B_{1/2} \}.
\]

\[
\Lambda_{b t}  :=   R_b  \Lambda^b_t, \quad \forall 0 < t < T^b
\]

and note that by the properties of $\Lambda^b_t$, the LQG coordinate change formula, and Weyl scaling (Fact <ref>), a.s.,

\[
t = \mu_{h^b}(\Lambda^b_t) = b^{-1} \mu_h(\Lambda_{ b t}), \quad \forall 0 < t < T^b.
\]

For similar reasons, the domains $\{\Lambda_{t}\}_{0 < t < b T^b}$ are a.s. harmonic balls and satisfy parts (a) and (b) of Theorem <ref> for the measure $\mu_h$
(with $t$ constrained to be less than $b T^b$).
These properties are enough to apply <cit.> (which works if the harmonic balls are only defined up to some fixed time) and conclude that $\{\Lambda_{t}\}_{0 < t < b T^b}$ is the unique such family. As $h^{b}$ has the same law as $h$, this shows the equality of joint laws.

By Lemma <ref>, a.s. for every integer $b \geq 0$,
there exists a unique family of harmonic balls $\{\Lambda_{t}\}_{0 < t < b T^b}$
which satisfies the conditions of the theorem statement up to time $b T^b$,
where $T^b$ is as in Lemma <ref>. Since $T^b \overset{d}{=} T$ for each $b>0$, where $T$ is as in Lemma <ref>, and $T$ is a strictly positive random variable, we deduce that $\lim_{b \to \infty} b T^b = \infty$ in probability. This complete the proof.

By applying the proof of Lemma <ref> to the infinite family of harmonic balls given by Theorem <ref>, we have the following.

For each $b > 0$, we have the equality of the joint laws

\[
\left( h, \{\Lambda_t\}_{t > 0} \right) \overset{d}{=}\left( h(R_b \cdot) + Q \log R_b - \frac{1}{\gamma} \log b, \{R_b^{-1} \Lambda_{b t}\}_{t > 0} \right).
\]

where the random radius $R_b > 0$ is as in (<ref>).

§.§ Whole plane space-filling SLE

The Schramm-Loewner evolution ($\text{SLE}_{\kappa}$) for $\kappa > 0$ was introduced by Schramm in <cit.>
as the only conformally invariant family of random curves which satisfy the domain Markov property. There are three phases of  $\text{SLE}_{\kappa}$: the curves are simple when $\kappa \in (0,4]$,  have self-intersections but are not space-filling for $\kappa \in (4,8)$, and are space-filling when $\kappa \geq 8$.

In the embedding of the mated-CRT map, we work with the whole plane space-filling SLE $\eta$, a variant of SLE introduced in
<cit.>.  This is a random space-filling curve $\eta$ in $\C$
from $\infty$ to $\infty$ which fills all of $\C$ and does not intersect the interior of its past.

For context, we informally describe how one constructs space-filling SLE. When $\kappa \geq 8$, whole plane space-filling SLE is simply ordinary SLE. In this case the
cells $\eta([x- {\eps}, x])$ are topological rectangles and their boundaries are the union of four $\text{SLE}_{\underline{\kappa}}$ curves for $\underline{\kappa} = 16/\kappa$.
When $\kappa \in (4,8)$,  $\eta$ is built by taking ordinary $\text{SLE}_{\kappa}$ and `filling in' bubbles which are disconnected from infinity.
This results in a somewhat complicated geometry. In particular, for $\kappa \in (4,8)$ it is possible for the intersection of two cells to be an uncountable, totally disconnected Cantor-like set. Pairs of cells which intersect in this manner do not correspond to edges of the mated-CRT map in Definition <ref>.

In this paper, we mainly just use the fact that the law of $\eta$ is scale-invariant:

\begin{equation} \label{eq:eta-scale-invariant}
r \eta \overset{d}{=} \eta, \quad \forall r > 0, \quad \mbox{when viewed as curves modulo time parameterization}.
\end{equation}

For a detailed treatment of SLE see <cit.> and for a more in-depth exposition of whole plane space-filling SLE see <cit.>.

§.§ Mated-CRT map

We defined most of what we need concerning the mated-CRT map in Section <ref>. We provide some additional notation which we will use in the sequel.

We fix the convention that when referring to points in $\C$, we use the letters $x,y,z$; while vertices in the mated-CRT map correspond to letters $a,b$.
Sets in $\C$ are $X, Y, D$, while sets of vertices in the mated-CRT map are $A,B$.

For ${\eps} > 0$ and $\gamma \in (0,2)$ let ${\Geps}$ be the mated-CRT map with vertex set ${\VGeps}$ embedded in the plane via a whole plane space-filling SLE curve $\eta$ with $\kappa = 16/\gamma^2$
as in Definitions <ref> and <ref>.

It is important to note that while a.s. each vertex in the mated-CRT map has a finite number of neighbors, the number of neighbors may be arbitrarily large.
Write $\deg^{{\eps}}(a)$ for the degree of $a \in {\VGeps}$. 	We will need a way to pass between vertices in the map and points in $\C$:  for $x \in \C$, let

\begin{equation} \label{eq:plane-to-embedding}
a_x^{{\eps}} := \mbox{(smallest $a \in {\eps} \Z$ such that $x \in \Heps{a}$)}.
\end{equation}

For a set $D \subset \C$, we define

\begin{equation} \label{eq:cells-domain-restriction}
{\Geps}(D) := \mbox{(subgraph of ${\Geps}$ induced by $\{a \in {\eps} \Z: \Heps{a} \cap D \neq \emptyset\}$)}
\end{equation}

and write ${\VGeps}(D)$ for the vertices in ${\Geps}(D)$.
Several times throughout the paper we will work with the sets ${\VGeps}(D) \setminus {\VGeps}(\partial D)$ rather than ${\VGeps}(D)$ to account
for the coarseness of the cells in the mated-CRT map. Indeed, even if $D$ is a smooth, bounded open set, typically ${\VGeps}(D)$ overlaps with ${\VGeps}(\partial D)$ — see Figure <ref>.

A set $D\subset \C$ and the set of mated-CRT map cells corresponding to the vertices in ${\Geps}(D)$, as in (<ref>).
The boundary of the domain $\partial D$ is drawn in red and the cells of ${\Geps}(D)$ are filled in pink. A point $x \in \C$ and the cell containing it, $\Heps{a^{{\eps}}_x}$, is drawn
in light blue. A similar figure appeared in  <cit.>.

§.§ Cell size and volume estimates

In this section, we recall some estimates proved in the appendix of <cit.>.
The following estimate gives worst-case upper and lower bounds for the LQG mass of Euclidean balls.

For each $\beta^+ \in (0, (2 -\gamma)^2/2)$ and $\beta^- > (2+\gamma)^2/2$
it holds with polynomially high probability as $\delta \to 0$ that

\begin{equation} \label{eq:ball-volume}
\delta^{\beta^-} \leq \mu_h(B_\delta(z)) \leq \delta^{\beta^+} ,\quad\forall z \in B_{1-\delta} .
\end{equation}

For convenience, we fix $\beta^{\pm}$ in a fashion depending only on $\gamma$ so that when we refer to Lemma
<ref> we may just refer to a single choice of $\beta^{\pm}$.

The next lemma gives similar worst-case upper and lower bounds for the size of the cells in the mated-CRT map.

Fix a small parameter $\zeta \in (0,1)$ and $\rho \in (0,1)$. With polynomially high probability as ${\eps} \to 0$, the Euclidean diameter of each mated-CRT map cell $\Heps{a}$ which intersects $B_\rho$ satisfies

\begin{equation}
{\eps}^{\frac{2}{(2 - \gamma)^2} + \zeta} \leq \mathrm{diam}(\Heps{a}) \leq {\eps}^{\frac{2}{(2 + \gamma)^2} - \zeta}.
\end{equation}

The following lemma is a quantitative version of the important fact that the counting measure on vertices in the mated-CRT map converges weakly to the LQG area measure.

Fix $\rho \in (0,1)$. There exists $\alpha = \alpha(\gamma) > 0$ and $\beta = \beta(\gamma) > 0$
such that with polynomially high probability as ${\eps} \to 0$, the following is true.
Let $D \subset B_{\rho}$ and let $f: B_{{\eps}^{\alpha}}(D) \to [0,\infty)$
be a non-negative function which is ${\eps}^{-\beta}$-Lipschitz continuous and satisfies $\|f\|_{\infty} \leq {\eps}^{-\beta}$.
For $a \in {\VGeps}(D)$, let $w^{{\eps}}_a$ be an arbitrary point of $\Heps{a} \cap \cl(D)$.
If we let $\mu_h$ be the $\gamma$-LQG area measure induced by $h$, then

\begin{equation}
{\eps}^{-1} \int_D f(z)d \mu_h(z) - {\eps}^{-1 + \alpha} \leq
\sum_{a \in {\VGeps}(D)} f(w^{\eps}_a) \leq
{\eps}^{-1} \int_{B_{{\eps}^\alpha}(D)} f(z) d \mu_h(z)  + {\eps}^{-1 + \alpha}
\end{equation}

simultaneously for every choice of $D$ and $f$ as above.

An analogous result to the above holds when cells are weighted by their degrees, in which case, since the mated-CRT map is a planar triangulation,
there is an extra factor of 6 in the limit.

For ${\eps}> 0$, let $\mu_h^{{\eps}}$ be the measure whose restriction to each cell $\Heps{a}$ for $a \in {\eps} \Z$
is equal to $\deg^{{\eps}}(a) \mu_h |_{\Heps{a}}$. Then as ${\eps} \to 0$, we have $\mu^{{\eps}}_h \to 6 \mu_h$
in probability w.r.t. the vague topology.

§.§ Liouville Brownian motion and random walk

Liouville Brownian motion (LBM) is the natural diffusion associated with $\gamma$-LQG and was constructed in <cit.>. Roughly, LBM is obtained from ordinary Brownian motion by changing time
so that $\mu_h$ is invariant. That is, LBM started at $z\in\C$ is defined as $B^z_{\phi^{-1}(t)}$, where $B^z_t$ is Brownian motion started at $z $ and

\[
\phi(t) = \lim_{{\eps} \to 0} \int_{0}^{t} {\eps}^{\gamma^2/2} e^{\gamma h_{{\eps}}(B_s^z)} ds,
\]

where $h_{{\eps}}$ is the average of $h$ over the circle of radius ${\eps}$.

In this paper we consider LBM rescaled by its (annealed) median exit time from a ball,

\begin{equation} \label{eq:rescaled-lbm}
X^z_t = B^z_{\phi^{-1}(t m_0)}, \quad \forall z \in \C,
\end{equation}

\begin{equation} \label{eq:mo}
m_0 = \mbox{median}(\inf\{t > 0 : |B_{\phi^{-1}(t)}| \geq 1/2\}).
\end{equation}

For $z \in \C$ and ${\eps} > 0$, let $X^{z, {\eps}} : \N \to {\VGeps}$ be a simple random walk on ${\Geps}$ started from $a_z^{{\eps}}$ (defined in (<ref>)).
On occasion we will also write, for $a \in {\VGeps}$, the walk started from $a$ as $X^{a, {\eps}}$.

LBM was shown in <cit.> to describe the scaling limit of random walk on the mated-CRT map. To state this precisely, we introduce the scaling factor

\begin{equation} \label{eq:meps}
{\meps} := (\mbox{median exit time of $\eta(X^{0, {\eps}})$ from $B_{1/2}$}),
\end{equation}

where $B_{1/2}$ is the Euclidean ball of radius $1/2$ centered at 0. By (1.8) in <cit.>,

\begin{equation} \label{eq:meps-scaling}
C^{-1} {\eps}^{-1} \leq {\meps} \leq C {\eps}^{-1}, \quad \mbox{for all sufficiently small ${\eps} > 0$},
\end{equation}

where $C > 1$ is a deterministic constant depending only on $\gamma$. We also extend the domain of definition of the embedded
walk $\N \ni j \to \eta(X^{z, {\eps}}_j)$ from $\N$ to $[0, \infty)$ by piecewise linear interpolation.

For each compact subset $K \subset \C$ and each $z \in K$,  the conditional law of the embedded linearly interpolated walk $(\eta(X^{z, {\eps}}_{m_{{\eps}} t}))_{t \geq 0}$
given $(h, \eta)$ converges in probability to the law of rescaled $\gamma$-LBM, defined by (<ref>), started from $z$ associated with $h$ with respect to the Prokhorov topology induced by the local uniform metric on curves $[0,\infty) \to \C$.

§.§ Discrete potential theory

Recall the notation for simple random walk from Section <ref>.
We denote the Green's function for simple random walk killed upon exiting a set $D \subset \C$ by

\begin{equation} \label{eq:greens-function}
\begin{aligned}
&\Gr{{\eps}}{D}(a,b) \\
&:= \E[\mbox{number of times that $X^{a, {\eps}}$ hits $b$ before exiting ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} | h, \eta]
\end{aligned}
\end{equation}

and its normalized version, the Green's kernel by

\begin{equation} \label{eq:greens-kernel}
\gr{{\eps}}{D}(a,b) := \frac{\Gr{{\eps}}{D}(a,b)}{\deg^{\eps}(b)}.
\end{equation}

By defining

\[
\Gr{{\eps}}{D}(a, \cdot) = \Gr{{\eps}}{D}(\cdot, a) = 0, \quad \forall a \not \in {\VGeps}(D) \setminus {\VGeps}(\partial D)
\]

we extend the domain of definition of both the Green's function and kernel to all of ${\VGeps}$.

For $a \in {\VGeps}$, the exit time of simple random walk on the mated-CRT map from a domain $D \subset \C$ is defined by

\begin{equation}  \label{eq:exit-time}
\taueps{D}{a} := \mbox{(first exit time of $X^{a, {\eps}}$ from ${\VGeps}(D) \setminus {\VGeps}(\partial D)$)}
\end{equation}

and its mean is

\begin{equation} \label{eq:expected-exit-time}
\Qeps{D}(a):= \E[\taueps{D}{a} | h, \eta]  = \sum_{b \in {\VGeps}(D) \setminus {\VGeps}(\partial D)} \Gr{{\eps}}{D}(a,b) .
\end{equation}

We will also need to consider the normalized expected exit time, for  $a \in {\VGeps}$

\begin{equation} \label{eq:normalized-exit-time}
\qeps{D}(a)= \sum_{b \in {\VGeps}(D) \setminus {\VGeps}(\partial D)} \gr{{\eps}}{D}(a,b).
\end{equation}

Let $\Delta^{{\eps}}$ denote the degree normalized graph Laplacian on $\mathcal{G}^{{\eps}}$,

\begin{equation} \label{eq:degree-normalized-laplacian}
\begin{aligned}
\Delta^{{\eps}} u(a) &= \frac{1}{\deg^{{\eps}}(a)} \sum_{b \sim a} (u(b) - u(a)) \quad \forall a \in {\VGeps}, \\
&\mbox{for functions $u: {\VGeps} \to \R$},
\end{aligned}
\end{equation}

where the sum $b \sim a$ is over the sites $b \in {\VGeps}$ which are joined by an edge to $b$.
As $\Delta^{{\eps}}$ is the generator of simple random walk on ${\Geps}$, for $D \subset \C$, by, , <cit.>,

\begin{equation} \label{eq:discrete-exit-laplacian}
\begin{cases}
\Delta^{{\eps}} \Qeps{D}(\cdot) = -1 \quad &\mbox{in ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} \\
\Qeps{D}(\cdot) = 0 \quad &\mbox{otherwise}
\end{cases}
\end{equation}

\begin{equation} \label{eq:discrete-normalized-exit-laplacian}
\begin{cases}
\Delta^{{\eps}} \qeps{D}(\cdot) = -1/\deg^{{\eps}}(\cdot) \quad &\mbox{in ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} \\
\qeps{D}(\cdot) = 0 \quad &\mbox{otherwise}
\end{cases}
\end{equation}

and for each $a \in {\VGeps}(D) \setminus {\VGeps}(\partial D)$,

\begin{equation} \label{eq:discrete-green-laplacian}
\begin{cases}
\Delta^{{\eps}} \Gr{{\eps}}{D}(\cdot,a) = -1\{\cdot = a\} \quad &\mbox{in ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} \\
\Gr{{\eps}}{D}(\cdot,a) = 0 \quad &\mbox{otherwise}
\end{cases}
\end{equation}

\begin{equation} \label{eq:discrete-normalized-green-laplacian}
\begin{cases}
\Delta^{{\eps}} \gr{{\eps}}{D}(a,\cdot) = -\frac{1\{\cdot = a\}}{\deg^{{\eps}}(a)} \quad &\mbox{in ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} \\
\gr{{\eps}}{D}(a,\cdot) = 0 \quad &\mbox{otherwise}.
\end{cases}
\end{equation}

For a subset $A$ of  ${\VGeps}$, we overload notation and
define $\Gr{{\eps}}{A}, 	\gr{{\eps}}{A}, \Qeps{D},$ and $\qeps{A}$
as above with the instances of ${\VGeps}(D) \setminus {\VGeps}(\partial D)$
replaced by $A$.  Each of (<ref>), (<ref>), (<ref>),
and (<ref>) hold with this same substitution.

We also write $\Delta$ for the continuum Laplacian on $\C$. The operator $\Delta^{{\eps}}$ shares several basic properties with its continuum counterpart $\Delta$, , linearity and monotonicity. Of importance to us is the following maximum principle for $\Delta^{{\eps}}$. Recall that a function $u: {\VGeps} \to \R$ is subharmonic (resp. superharmonic) on $A \subset {\VGeps}$ if $\Delta^{{\eps}} u \geq 0$
(resp. $\Delta^{{\eps}} u \leq 0$) on $A$ and is harmonic if it is both subharmonic and superharmonic.  When we need to emphasize the distinction between $\C$ and ${\VGeps}$, we qualify functions as $\Delta$-harmonic or $\Delta^{{\eps}}$-harmonic.

For every connected subset $A \subset {\VGeps}$, if $u: \cl(A) \to \R$ is subharmonic on $A$, then

\begin{equation} \label{eq:maximum-principle}
\max_{a \in A} u(a) \leq \max_{b \in \partial A} u(b).
\end{equation}

This is standard. Consider a point $a \in A$ and apply the definition of subharmonicity to find a
neighbor $\cl(A) \ni b \sim a$ for which $u(b) \geq u(a)$.

In this article, we will use discrete Green's function and exit times estimates established in <cit.>.
Fix a $\rho \in (0,1)$ — several of our estimates will be stated only on $B_\rho$. The reason for working on $B_\rho$ is that $h|_{B_1}$ agrees in law with the corresponding restriction of a whole-plane GFF plus $\gamma\log(1/|\cdot|)$, see (<ref>). Consequently the constants in several of these estimates will depend on $\rho$.
First, we recall an upper bound on the Green's kernel.

There exists $\beta = \beta(\gamma) > 0$
and $C = C(\rho, \gamma) > 0$ such that
with polynomially high probability as ${\eps} \to 0$,

\[
\gr{{\eps}}{B_{\rho}}(a,b)
\leq C \log \left(  \frac{1}{|\eta(a) - \eta(b)|} \right) + C
\]

simultaneously for all $a,b \in {\VGeps}(B_{\rho})$
with $|\eta(a) - \eta(b)| \geq {\eps}^{\beta}$.

We will also need an upper bound on the Green's function. Recall that the Green's function is the Green's kernel multiplied by the degree.
Since the vertex degrees on the mated-CRT map are unbounded, we cannot directly apply the previous lemma.
Nevertheless, the following crude upper bound will suffice.

For each $\zeta \in (0,1)$, it holds
with polynomially high probability as ${\eps} \to 0$ that

\[
\sup_{a,b \in {\VGeps}(B_{\rho})} \Gr{{\eps}}{B_{\rho}}(a,b)
\leq {\eps}^{-\zeta}.
\]

We also need a lower bound on both the Green's kernel and function.

There exists $\omega = \omega(\gamma) > 0$ and $C = C(\rho, \gamma) > 0$
such that with polynomially high probability as ${\eps} \to 0$ the following holds.
For each $a \in {\VGeps}(B_{\rho - {\eps}^{\omega}})$ and
each $r \in [{\eps}^{\omega}, \dist(\eta(a), \partial B_{\rho})]$,

\[
\begin{aligned}
\Gr{{\eps}}{B_r(\eta(a))}(a,b)  \geq \gr{{\eps}}{B_r(\eta(a))}(a,b) &\geq C^{-1} \log \left( \frac{r}{|\eta(a) - \eta(b)|} \wedge {\eps}^{-1} \right) - 1, \\
\quad \forall b \in {\VGeps}(B_{r/3}(\eta(a))).
\end{aligned}
\]

We conclude this subsection by formulating an upper bound on the expected exit time function in terms of the coarseness
of the ${\eps}$-mated CRT map and a certain continuum quantity. To that end, we consider the following random variable (from (5.1) in <cit.>)

\begin{equation} \label{eq:largest-exit-time-rv}
\Mrho{D} := \sup_{z \in D} \int_{D} \left( \log \left( \frac{1}{|z-w|} \right) + 1 \right) d\mu_h(w), \quad \mbox{for $D \subset \C$}
\end{equation}

and recall the following bounds.

Almost surely,

\[
\sup_{D \subseteq B_{\rho}} \Mrho{D} < \infty
\]

and with polynomially high probability as $\delta \to 0$

\[
\sup_{y \in B_{\rho}} \Mrho{B_{\delta}(y)} \leq \delta^q,
\]

for some deterministic $q = q(\gamma) > 0$.

Step 1: $	\mathcal{M}(D) < \infty$.

As $\mu_h$ is a locally finite measure, it suffices to bound the log integrand in (<ref>).
By <cit.>, there are constants $c_0, c_1$ depending only on $\rho$ and $\gamma$ such that for each $A > 1$,
it holds with probability at least $1 - c_0 e^{-c_1 A}$ that

\begin{equation} \label{eq:intermediate-bound}
\sup_{z \in D} \int_{D}  \log \left( \frac{1}{|z-w|} \right) d\mu_h(w) \leq A ( \mu_h(D) + e^{-c_1 A}).
\end{equation}

By the Borel-Cantelli lemma, this shows $\sup_{D \subseteq B_{\rho}} \Mrho{D} < \infty$ almost surely.

Step 2: $\sup_{y \in B_{\rho - \delta}} \Mrho{B_{\delta}(y)} \leq \delta^q$.

This follows by a partition of the integral into dyadic annuli as in the proof of <cit.>. An explicit proof is given in
the proof of Corollary 5.2 in <cit.> – this bound is exactly (5.3) there.

By <cit.>, we can bound the expected exit time in terms of this random variable.

There exists deterministic $\alpha = \alpha(\gamma) > 0$ and $C = C(\rho, \gamma) > 0$ such that with polynomially high probability as ${\eps} \to 0$
it holds simultaneously for every Borel set $D \subset B_{\rho}$ with $\mu_h(D) \geq {\eps}^{\alpha}$
that for all $a \in {\VGeps}(D)$

\[
\E[(\taueps{D}{a})^N| h, \eta]  \leq N! C^N {\eps}^{-N} \left( \sup_{z \in D} \int_{B_{{\eps}^{\alpha}}(D)}  \left ( \log \left( \frac{1}{|z-w|} \right) + 1  \right) d \mu_h(w) \right)^N
\quad \forall N \in \N,
\]

in particular, for all $D \subset B_{\rho}$

\[
\qeps{D}(a) \leq \Qeps{D}(a)\leq C {\eps}^{-1} \Mrho{B_{{\eps}^{\alpha}}(D)}, \quad \forall a \in {\VGeps}(D).
\]

By combining the previous two lemmas, we have the following bounded on the expected exit time.

There exists deterministic $\alpha = \alpha(\gamma) > 0$ and $C = C(\rho, \gamma) > 0$ such that with polynomially high probability as ${\eps} \to 0$
it holds simultaneously for every Borel set $D \subset B_{\rho}$ with $\mu_h(D) \geq {\eps}^{\alpha}$

\[
\qeps{D} \leq \Qeps{D} < \infty.
\]

Moreover,	for some deterministic $q = q(\gamma) > 0$, it holds with probability tending to 1 as ${\eps} \to 0$ and then $\delta \to 0$ that

\[
\qeps{B_{\delta}(y)} \leq \Qeps{B_{\delta}(y)} \leq {\eps}^{-1} \times \delta^q,
\]

for every $y \in \C$ such that $B_{\delta}(y) \subset B_{\rho}$.

Apply Lemma <ref> with some fixed $\rho' \in (\rho, 1)$ to Lemma <ref> and absorb the constants into the factor of $\delta^q$.

§ CONVERGENCE OF THE DIRCHLET PROBLEM

In this section we use the estimates recalled in the previous section to prove uniform convergence of both the Green's kernel and the solution to the Dirichlet problem to their continuum counterparts.  Throughout this section, we fix $\rho \in (0,1)$ such that $B_{\rho} \Subset B_1$.

§.§ Scaling

Recall the definition of ${\meps}$  and recall from (<ref>) that ${\meps} \asymp {\eps}^{-1}$.
For a set $D \subset B_{\rho}$, we rescale a function, $u^{{\eps}}: {\VGeps}(D) \to \R $
by defining $\overline{u}^{{\eps}}: D \to \R$ in a piecewise constant fashion

\begin{equation} \label{eq:discrete-scaling}
\overline{u}^{{\eps}}(x) := {\meps}^{-1}  u^{{\eps}}(a^{{\eps}}_x), \quad \forall x \in D,
\end{equation}

where $a^{{\eps}}_x$ is as in (<ref>).  That is, $\overline{u}^{{\eps}}$ is constant within the cells of the mated-CRT map.

Throughout this section, we work with functions which are of order ${\eps}^{-1}$. We do this because this is the order of the expected exit time function,
Lemma <ref>, and the divisible sandpile odometer in Section <ref>.

§.§ Tightness

In this subsection we prove a general tightness result concerning bounded functions with bounded Laplacian.

For every connected set $D \subset B_{\rho}$
and $\sigma(h,\eta)$-measurable sequence  $\mathcal{E} \ni {\eps} \to 0$,
if a sequence of functions $u^{{\eps}}: {\VGeps}(D) \to \R $ indexed by ${\eps} \in \mathcal E$ satisfies

\[
\begin{aligned}
|\Delta^{{\eps}} u^{{\eps}}| \leq C &\quad \mbox{on ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} \\
u^{{\eps}} \leq C {\eps}^{-1} &\quad \mbox{on ${\VGeps}(D)$}
\end{aligned}
\]

uniformly over ${\eps} \in \mathcal{E}$ for some $\sigma(h,\eta)$-measurable $C > 0$, then the following is true.
Almost surely, there exists a $\sigma(h,\eta)$-measurable subsequence $\mathcal{E}_0 \subset \mathcal{E}$ and a function $\overline{u} \in C(D)$ such that for every compact subset $O \Subset D$,
we have that  $\overline{u}^{{\eps}} \to \overline{u}$ uniformly in $O$ as $\mathcal{E}_0 \ni {\eps} \to 0$.

Moreover, if $D$ is simply connected, $u^{{\eps}}$ is discrete harmonic in ${\VGeps}(D) \setminus {\VGeps}(\partial D)$, and if $u^{{\eps}}$ has Hölder continuous boundary data, in the sense that there exists $\chi \in (0,1]$ and
$C' > 0$ so that

\begin{equation} \label{eq:holder-boundary-data}
|u^{{\eps}}(a) - u^{{\eps}}(b)| \leq C'( {\eps} \vee |\eta(a) - \eta(b)|)^{\chi}, \quad \forall a,b \in {\VGeps}(\partial D), \quad \forall {\eps}  \in (0,1)
\end{equation}

then the convergence occurs uniformly in $\cl(D)$.

Lemma <ref> will follow from the Arzéla-Ascoli theorem and the following
Hölder continuity of $\Delta^{{\eps}}$-harmonic functions on the mated-CRT map.

There exists $\xi = \xi(\rho, \gamma) \in (0,1)$ and $A = A(\rho, \gamma) > 0$
such that the following holds with polynomially high probability as ${\eps} \to 0$.
Let $D \subset B_{\rho}$ be a connected domain and let $f^{{\eps}}: {\VGeps}(D) \to \R$ be
$\Delta^{{\eps}}$-harmonic on ${\VGeps}(D) \setminus {\VGeps}(\partial D)$. Then,

\begin{equation}
|f^{{\eps}}(a) - f^{{\eps}}(b)| \leq A \left( \sup_{a' \in {\VGeps}(D)} |f^{{\eps}}(a')| \right) \left( \frac{{\eps} \vee |\eta(a) - \eta(b)|}{\dist(\eta(a), \partial D)} \right)^\xi, \quad \forall a,b \in {\VGeps}(D).
\end{equation}

Moreover, if $D$ is simply connected and if $f^{{\eps}}$ has Hölder continuous boundary data in the sense of (<ref>), then

\begin{equation} \label{eq:boundary-estimate}
|f^{{\eps}}(a) - f^{{\eps}}(b)| \leq \max \left\{ C, A \left( \sup_{a' \in {\VGeps}(D)} |f^{{\eps}}(a')| \right) \right \}( {\eps} \vee |\eta(a) - \eta(b)|)^{\xi}, \quad \forall a,b \in {\VGeps}(D).
\end{equation}

The conclusions of Lemmas <ref>, <ref>, and <ref>
hold with probability approaching one as $\mathcal{E} \ni {\eps} \to 0$. In particular, by the Borel-Cantelli lemma, almost surely, there exists a deterministic subsequence $\mathcal{E}_1 \subset \mathcal{E}$ such that the conclusions of all of these lemmas hold for all sufficiently small ${\eps} \in \mathcal{E}_1$. We henceforth restrict to such an event.

We first prove interior convergence: let $O \Subset D \subset B_{\rho}$ be given. By the Arzéla-Ascoli theorem, and the boundedness assumption on $\overline{u}^{{\eps}}$, it suffices to prove the following estimate:
almost surely, for each small $k > 0$, there exists possibly random choices of $\delta = \delta(k) > 0$ and ${\eps}_1 = {\eps}_1(k, \delta) \in \mathcal{E}_1$ so that for all ${\eps} \in (0, {\eps}_1) \cap \mathcal{E}_1$,

\begin{equation} \label{eq:equicontinuity}
|x-y| \leq \delta \implies |\overline{u}^{{\eps}}(x) - \overline{u}^{{\eps}}(y)| \leq k, \quad \forall x,y \in O.
\end{equation}

We now prove (<ref>). Let $k > 0$ be given, fix $\delta > 0$ to be determined below, in a manner depending on $k$. By Lemma <ref>,

\begin{equation} \label{eq:cell-size-estimate}
\sup_{a \in {\VGeps}(B_{\rho})} \diam(\Heps{a}) \leq \delta^2, \quad \forall {\eps} \in \mathcal{E}_1 \cap (0, {\eps}_1)
\end{equation}

where ${\eps}_1$ depends on $\delta$ through Lemma <ref>.

As $O \Subset D$, there exists $\delta_1 > 0$ so that $O + B_{\delta_1} \Subset D$.
Suppose $\delta \in (0,  \min(\delta_1^{2}, \delta_1)/100 )$ and let $x \in O$ be given.  Let $\underline{\delta} = \delta^{1/2} + \delta$ and
$f^{{\eps}}:{\VGeps}(\cl(B_{\underline{\delta}}(x))) \to \R$ be the harmonic extension of $u^{{\eps}}$,

\[
\begin{cases}
\Delta^{{\eps}} f^{{\eps}} = 0 \quad \mbox{on ${\VGeps}(B_{\underline{\delta}}(x)) \setminus {\VGeps}(\partial B_{\underline{\delta}}(x))$} \\
f^{{\eps}} = u^{{\eps}} \quad \mbox{on ${\VGeps}(\partial B_{\underline{\delta}}(x))$} .
\end{cases}
\]

By the maximum principle (Lemma <ref>) and our uniform boundedness assumption on $u^{{\eps}}$,

\begin{equation} \label{eq:harmonic-extension-bounded}
\sup_{a' \in {\VGeps}(\cl(B_{\delta}(x)))} |f^{{\eps}}(a')| \leq \sup_{a' \in {\VGeps}(\cl(B_{\delta}(x)))} |u^{{\eps}}(a')| \leq  C {\eps}^{-1}  , \quad \forall {\eps} \in (0, {\eps}_1) \cap \mathcal{E}.
\end{equation}

By the assumption $|\Delta^{{\eps}} u^{{\eps}}| \leq C$ and (<ref>), the functions $(u^{{\eps}} - f^{{\eps}}) \pm C \times \Qeps{B_{\underline{\delta}}(x)}$ are super (resp. sub) harmonic on ${\VGeps}(B_{\underline{\delta}}(x)) \setminus {\VGeps}(\partial B_{\underline{\delta}}(x))$ and identically equal to 0 on ${\VGeps}(\partial B_{\underline{\delta}}(x))$ for each ${\eps}$. Therefore, by the maximum principle,

\[
|u^{{\eps}}-f^{{\eps}}| \leq  C \times \Qeps{B_{\underline{\delta}}(x)}  \quad \mbox{on ${\VGeps}(\cl(B_{\underline{\delta}}(x)))$}
\]

for all ${\eps} \in (0, {\eps}_1) \cap \mathcal{E}_1$.
Thus,  we have, for such ${\eps}$, by the triangle inequality,

\begin{equation} \label{eq:harmonic-decomposition}
\sup_{a,b \in {\VGeps}(B_{\delta}(x))} |u^{{\eps}}(a) - u^{{\eps}}(b)|  \leq \sup_{a,b \in {\VGeps}(B_{\delta}(x))} |f^{{\eps}}(a) - f^{{\eps}}(b)| + 2 C \times \sup_{a\in  {\VGeps}(B_{\delta}(x))} \Qeps{B_{\underline{\delta}}(x)}(a) .
\end{equation}

By Lemma  <ref>, we may take ${\eps}_1$ smaller, depending only on $\delta$, so that

\begin{equation} \label{eq:second-term-bound}
2 C \times \sup_{a\in  {\VGeps}(B_{\delta}(x))} \Qeps{B_{\underline{\delta}}(x)}(a)   \leq {\eps}^{-1} \delta^q, \quad \forall {\eps} \in (0, {\eps}_1) \cap \mathcal{E}_1
\end{equation}

for a deterministic $q  = q(\gamma) > 0$.

It remains to bound the first term in (<ref>): with $A = A(\rho, \gamma) > 0$ and $\xi \in (0,1)$ as in Lemma <ref>,

\begin{align*}
&\sup_{a,b \in {\VGeps}(B_{\delta}(x))} |f^{{\eps}}(a) - f^{{\eps}}(b)| \\
&\leq  A \left( \sup_{a' \in {\VGeps}(B_{\delta}(x))} |f^{{\eps}}(a')| \right) \sup_{a,b \in {\VGeps}(B_{\delta}(x))} \left( \frac{{\eps} \vee |\eta(a) - \eta(b)|}{\dist(\eta(a), \partial B_{\underline{\delta}}(x))} \right)^\xi \quad \mbox{(Lemma \ref{lemma:discrete-holder-continuity})} \\
&\leq A C {\eps}^{-1} \sup_{a,b \in {\VGeps}(B_{\delta}(x))} \left( \frac{{\eps} \vee |\eta(a) - \eta(b)|}{\dist(\eta(a), \partial B_{\underline{\delta}}(x))} \right)^\xi \quad \mbox{(by \eqref{eq:harmonic-extension-bounded})}  \\
&\leq A C {\eps}^{-1} \left( \frac{\delta + \delta^2}{\delta^{1/2} - \delta^2} \right)^\xi \quad \mbox{(by \eqref{eq:cell-size-estimate} and since $\underline\delta = \delta^{1/2} + \delta$)}  \\
&\leq {\eps}^{-1} k/2 \quad \mbox{(for small $\delta = \delta(k, \xi)$)},  \stepcounter{equation}\tag{\theequation}\label{eq:first-term-bound}
\end{align*}

for all ${\eps} \in (0, {\eps}_1) \cap \mathcal{E}_1$.

By plugging (<ref>) and (<ref>) into (<ref>), we get that for sufficiently small $\delta = \delta(k, \xi)$,

\[
\sup_{y \in {\VGeps}(B_{\delta}(x))} |u^{{\eps}}(y) - u^{{\eps}}(x)|  \leq {\eps}^{-1} k, \quad \forall {\eps} \in (0, {\eps}_1) \cap \mathcal{E}_1.
\]

By combining this inequality with (<ref>) and the definition (<ref>) of $\overline{u}^{{\eps}}$, we have (<ref>).

We now assume the conditions of the moreover clause. Since we have assumed $u^{{\eps}}$ is harmonic, we can apply (<ref>) of Lemma <ref> to get the analogue of (<ref>) up to the boundary: for each $k > 0$, there exists $\delta = \delta(k,C',\chi) > 0$ such that for each small enough ${\eps} \in \mathcal E_1$,

\begin{equation} \label{eq:up-to-boundaryequicontinuity}
|x-y| \leq \delta \implies |\overline{u}^{{\eps}}(x) - \overline{u}^{{\eps}}(y)| \leq k, \quad \forall x,y \in {\VGeps}(D).
\end{equation}

And so, we may conclude using the Arzéla-Ascoli theorem.

§.§ Uniform convergence of the Green's kernel

For $z \in \C$, let

\[
X^{z, {\eps}} : \N \to {\VGeps}
\]

be a simple random walk on ${\Geps}$ started from the vertex $a^{{\eps}}_z \in {\VGeps}$, defined in (<ref>),
stopped upon exiting ${\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})$.
Let $X^{z}$ be LBM started at $z$ stopped upon exiting $B_{\rho}$
and  let

\[
\hat{X}^{z, {\eps}}: [0,\infty) \to \C
\]

be the process
$t \to \eta(X^{z, {\eps}}_{{\meps} t})$ extended from ${\meps}^{-1} \N$ to $[0,\infty)$ by piece-wise linear interpolation.
The goal of this section is to prove the following.

Fix a deterministic choice of $r>0$ and $z\in\C$ such that $B_r(z) \subset B_\rho$ and
$y \in B_r(z)$. Let $g^{{\eps}}(\cdot) := {\eps}^{-1} \gr{{\eps}}{B_r(z)}(a^{{\eps}}_y , \cdot)$ and define $\overline{g}^{\eps}$ as in (<ref>).
For each small $\delta > 0$ and $O \subset \cl(B_r(z)) \setminus \{y\}$ which lies at positive distance from $y$,
with probability tending to 1 as ${\eps} \to 0$,

\begin{equation} \label{eq:green-uniform}
\sup_{x \in O} \left|\overline{g}^{{\eps}}(x)  - \frac{m_0}{6} \times G_{B_r(z)}(y , x)\right| \leq \delta,
\end{equation}

where $G_{B_r(z)}(y,x)$ is the continuum Green's function for standard Brownian motion killed upon exiting $B_r(z)$.

Lemma <ref> will be shown as a consequence of the following result and Lemma <ref>.

For each deterministic choice of $z\in\C$ and $r > 0$ such that $B_r(z) \subset B_{\rho}$, each $y \in B_r(z)$, and each bounded, Hölder-continuous, $\sigma(h, \eta)$-measurable function $f: B_{r}(z) \to \R$,

\[
\begin{aligned}
{\meps}^{-1}  \sum_{x \in {\VGeps}(B_r(z))} \Gr{{\eps}}{B_r(z)}(a^{{\eps}}_{y}, x) f(\eta(x)) \to m_0 \times \int_{B_r(z)} G_{B_r(z)} (y, x) f(x) d \mu_h(x),
\end{aligned}
\]

in probability as ${\eps} \to 0$, where $G_{B_r(z)}$ is the continuum Green's function, and $a^{{\eps}}_{y}$ is as in (<ref>).

This follows from convergence of random walk on ${\Geps}$ to LBM and the fact that the discrete and continuum Green's functions can be represented in terms of these objects.

The conclusions of Lemmas <ref> and <ref>
hold with probability approaching one as ${\eps} \to \infty$. Restrict to the event that these occur for all small enough ${\eps}$.

Let the deterministic parameters $y \in B_r(z) \subset B_{\rho}$ and $f: B_r(z) \to \C$ a bounded and Hölder continuous function which is $\sigma(h, \eta)$-measurable be given. Consider the bounded, continuous
functional on the space of continuous curves $w:   [0,\infty)  \to B_{\rho}$,

\[
\phi^{f}_T(w) := \int_0^T  f(w_t)  dt, \quad \forall T > 0 .
\]

By boundedness and
Theorem <ref>,

\begin{equation} \label{eq:rw-to-lbm}
\lim_{{\eps} \to 0} \E[ \phi^{f}_T(\hat{X}^{y, {\eps}}) | h , \eta] = \E\left[ \int_0^T f(X^y_t) dt | h, \eta\right], \quad \forall T > 0 ,
\end{equation}

where the convergence is in probability.

The term on the left in (<ref>) can be used to control its discrete counterpart.
Indeed, as $f$ is Hölder continuous and bounded and the cell size of the mated-CRT map is controlled by Lemma <ref>, for all ${\eps}$ small depending only on $h, \eta$,

\[
{\meps}^{-1} \sum_{t=0}^{\lceil T {\meps} \rceil} f(\eta(X^{y, {\eps}}_t))  - C {\eps}^{q} \leq \phi^{f}_T(\hat{X}^{y, {\eps}}) \leq {\meps}^{-1} \sum_{t=0}^{\lceil T {\meps} \rceil} f(\eta(X^{y, {\eps}}_t))  + C {\eps}^{q}
\]

for all $T > 0$, where $C = C(f) > 0$, $q = q(\gamma, f) > 0$, and $a^{{\eps}}_y$ is as in (<ref>).  Thus, by (<ref>),

\begin{equation} \label{eq:green-trunc-weak-convergence}
\lim_{{\eps} \to 0} {\meps}^{-1}\E\left[ \sum_{t=0}^{\lceil T {\meps} \rceil } f(\eta(X^{y, {\eps}}_t))  | h , \eta\right] =  \E\left[ \int_0^T f(X^y_t) dt | h, \eta\right], \quad \forall T > 0,
\end{equation}

in probability as ${\eps} \to 0$. By Lemmas <ref>, <ref>  and the boundedness of $f$, with $\taueps{B_r(z)}{a^{{\eps}}_y}$ the first exit time as in (<ref>),

\[
|{\meps}^{-1}\ \sum_{t=0}^{\lceil T {\meps} \rceil } f(\eta(X^{y, {\eps}}_t)) |
< {\meps}^{-1} \sup_{x \in B_r(z)} |f(x)| \taueps{B_r(z)}{a^{{\eps}}_y}
< \infty, \quad \forall T  >0, \quad \forall {\eps} \in (0,1),
\]

and ${\meps}^{-1} \E[\taueps{B_r(z)}{a^{{\eps}}_y}  | h, \eta]$
is a.s. bounded uniformly in ${\eps}$.

Therefore, as (<ref>) holds for each $T > 0$, we may take $T \to \infty$ and interchange the limits to see that

\begin{align*}
{\meps}^{-1} \sum_{b \in {\VGeps}(B_r(z))} \Gr{{\eps}}{B_r(z)}(a^{{\eps}}_y, b) f(\eta(b))
&= {\meps}^{-1}\E\left[ \sum_{t=0}^{\infty } f(\eta(X^{y, {\eps}}_t))  | h , \eta\right]  \\
&\overset{{\eps} \to 0}{\rightarrow}  \E\left[ \int_0^\infty f(X^y_t) dt | h, \eta\right] \\
&= m_0 \times \int_{B_r(z)} G_{B_r(z)} (y, x) f(x) d \mu_h(x),
\end{align*}

completing the proof.

We now use the integral convergence of the Green's function (Lemma <ref>) together with Lemma <ref> and tightness to show uniform convergence of the Green's kernel away from its pole.

Let a sequence of ${\eps}$s tending to zero be given. Fix $y \in B_r(z) \subset B_{\rho}$ and let $O \Subset B_{r}(z) \setminus \{y\}$ be open.
Denote by  $\mathcal{P}(B_r(z))$ the set of polynomials with rational coefficients in $B_r(z)$.

There exists a deterministic subsequence $\mathcal{E}$ along which a.s. the convergence in Lemma  <ref>
occurs and the convergence in Lemma <ref>
occurs for each $f \in \mathcal{P}(B_r(z))$.
We may further assume, by the Borel-Cantelli lemma, that along this subsequence the conclusions of Lemma <ref>
and Lemma <ref> hold.  Restrict to the subsequence $\mathcal{E}$ and the event of the preceding two sentences.

By Lemma <ref>,  $\gr{{\eps}}{B_{\rho}}(a^{{\eps}}_y , \cdot)$ is uniformly bounded in ${\VGeps}(O)$ by a random constant along $\mathcal{E}$, and by (<ref>) is $\Delta^{{\eps}}$-harmonic in ${\VGeps}(O)$. Thus, by Lemma <ref>, there is a random subsequence along which $\overline g^{\eps}$ converges uniformly in $O$ to some continuous function $\overline g^* : O \to \R$ (possibly depending on the subsequence).

Fix such a subsequence $\mathcal{E}_1 \subset \mathcal{E}$, and subsequential limit $\overline{g}^*$. Let $\mu_h^{{\eps}}$ be the measure on $\C$ whose restriction to each cell $\eta([x-{\eps}, x])$ for $x \in {\eps} \Z$
is equal to $\deg^{{\eps}}(x)$ times $\mu_h |_{\eta([x-{\eps},x])}$. As we have assumed a.s. convergence of the
measure $\mu_h^{{\eps}}$ to $6 \mu_h$ along $\mathcal{E}_1$ with respect to the local Prokhorov topology, we have that as $\mathcal{E}_1 \ni {\eps}  \to 0$,

\[
\int_{O} f(x) \overline{g}^{{\eps}}(x) d \mu^{{\eps}}_h(x)
\to 6 \times \int_{O} f(x) \overline{g}^*(x) d \mu_h(x),
\quad \forall \mbox{$f:\mathcal{P}(B_r(z)) \to \R$},
\]

almost surely.
By unpacking the definitions of $\overline{g}^{{\eps}}$ and $\mu_h^{{\eps}}$,
and using Lemma <ref> to control the approximation error, we have that a.s. as $\mathcal E \ni {\eps} \to 0$,

\[
\begin{aligned}
&\left | \int_{O} f(x) \overline{g}^{{\eps}}(x) d \mu^{{\eps}}_h(x) - 	{\meps}^{-1} \sum_{b \in {\VGeps}(O)} f(\eta(x)) \Gr{{\eps}}{B_{r}(z)}(a^{{\eps}}_y, b) \right|  \leq C {\eps}^{q}, \\
&\quad \mbox{for constants $C = C(f)$ and $q = q(\gamma, f) > 0$,} \\
& \quad \forall \mbox{$f:\mathcal{P}(B_r(z)) \to \R$}.
\end{aligned}
\]

Combining the previous two intended equations together with Lemma <ref>
shows that $\overline{g}^*(x) = \frac{m_0}{6} \times G_{B_{r}(z)}(y, x)$. Hence, (<ref>) holds when $O\Subset B_r(z) \setminus \{y\}$.

After identifying the limit, we now argue that the convergence occurs up to $\partial B_r(z)$. By a diagonal argument, we may arrange so that
$\overline g^{\eps}$ converges uniformly to $ \frac{m_0}{6} \times G_{B_{r}(z)}(y, x)$ along $\mathcal{E}$ on each compact subset of $B_r(z) \setminus \{y\}$.
Let $s > 0$ be given and let $r_1 \in (0, r)$ be such that $y$ lies at positive distance from $\cl(\A_{r_1, r}(z))$ and

\[
\sup_{x \in \cl(\A_{r_1, r}(z))} G_{B_{r}(z)}(y, x) < s/4
\]

Indeed, this is possible as $G_{B_{r}(z)}(y, \cdot)$ is smooth up to the boundary away from its pole and is identically zero on $\partial B_r(z)$.
By the maximum principle, Lemma <ref>, and local uniform convergence we have that

\[
\sup_{b \in {\VGeps}(\A_{r_1, r}(z))} \gr{{\eps}}{B_{\rho}}(a^{{\eps}}_y , b)
\sup_{b \in {\VGeps}(\partial B_{r_1}(z))} \gr{{\eps}}{B_{\rho}}(a^{{\eps}}_y , b)
< s/2,
\]

for all ${\eps}$ sufficiently small. This implies convergence up to the boundary, $\partial B_r(z)$.

§.§ Convergence of the Dirichlet problem

We follow a similar strategy as in the previous subsection to prove the following.

Fix a deterministic choice of $z\in\C$ and $r>0$ such that $B_r(z) \subset B_{\rho}$, a bounded, deterministic Hölder-continuous function $\psi: B_{r}(z) \to \R$, and a deterministic
Hölder continuous function $\phi: \cl(B_r(z)) \to \R$.
The solution to the discrete Dirichlet problem

\[
\begin{cases}
\Delta^{{\eps}} f^{{\eps}}(a) = \psi(\eta(a)) \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
f^{{\eps}}(b) = {\meps} \times  \phi(\eta(b)) \quad \mbox{for $b \in {\VGeps}(\partial B_r(z))$}
\end{cases}
\]

converges in probability in the uniform topology on $\cl(B_r(z))$ to the solution of the continuum Dirichlet problem,

\[
\begin{cases}
\Delta \overline{f} = m_0 \psi \times \mu_h  \quad\mbox{on $B_r(z)$} \\
\overline{f}  = \phi \quad \mbox{on $\partial B_r(z)$}.
\end{cases}
\]

Moreover, a.s. the limiting function $\overline{f}$ is Hölder continuous in $\cl(B_{r}(z))$.

We first prove a pointwise convergence result concerning harmonic functions.

For each deterministic choice of $z\in\C$, $r>0$ such that $B_r(z) \subset B_{\rho}$ and each deterministic Hölder continuous function $\psi: \cl(B_{r}(z)) \to \R$
the following occurs. For $y \in \C$, recall that $\taueps{B_r(z)}{a_y^{{\eps}}}$
is the first time simple random walk started at $a_y^{{\eps}} \in {\VGeps}$ exits ${\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$. For $y \in B_{r}(z)$, let

\[
\tau_{B_r(z)}^{y} = \min \{t \geq 0 : X^{y}_t \not \in B_r(z) \}
\]

denote the first time LBM started at $y$ exits $B_r(z)$.

\[
\lim_{{\eps} \to 0}  \E \left[ \psi \left( \eta \left( X^{a^{{\eps}}_y,{\eps}}_{\taueps{B_r(z)}{a_y^{{\eps}}}} \right) \right)  | h, \eta \right]  = \E \left[ \psi \left( X^y_{\tau_{B_r(z)}^{y}} \right) | h, \eta \right], \quad \forall y \in B_{r}(z),
\]

where the convergence is in probability as ${\eps} \to 0$.

The proof is nearly identical to that of Lemma <ref>, but for completeness we sketch the proof.

Let the parameters $y \in B_{r}(z)$ and $\psi: \cl(B_{r}) \to \R$, a Hölder continuous function, be given.
Consider the bounded
functional on the space of continuous curves $w: C([0,\infty)) \to \cl(B_{r}(z))$ defined by

\[
\phi^{\psi}(w) := \psi(w(\tau)),
\]

where $\tau$ is the first time $w(t) \not \in B_r(z)$.
Recall from the beginning of Section <ref> that $\hat X^{y,{\eps}}$ is the piece-wise linear interpolation
of random walk started at $a^{{\eps}}_y$ on ${\VGeps}$ stopped upon exiting ${\VGeps}(B_\rho) \setminus {\VGeps}(\partial B_{\rho})$.
By Theorem <ref>,

\begin{equation} \label{eq:rw-to-lbm-harmonic}
\lim_{{\eps} \to 0} \E \left[ \phi^{\psi}(\hat{X}^{y, {\eps}}) | h , \eta \right] = \E \left[ \psi(X^y_{\tau^{y}}) | h, \eta \right],
\end{equation}

in probability as ${\eps} \to 0$.
We conclude by observing that the term on the left-hand-side  of (<ref>) is asymptotically close to its discrete counterpart.
Indeed, as $\psi$ is Hölder continuous on $\cl(B_{r}(z))$ and the cell-size of the mated-CRT map can be controlled by Lemma <ref>,

\[
|\phi^{\psi}(\hat{X}^{y, {\eps}}) -  \psi(\eta(X^{y}_{\tau^{{\eps}, y}}))| \leq C \times {\eps}^{q},
\]

in probability as ${\eps} \to 0$, for some $C = C(\psi) > 0$ and $q = q(\gamma, \psi) > 0$.

This pointwise convergence is used together with tightness and integral convergence of the Green's function to prove uniform convergence of solutions
to the discrete Dirichlet problem.

Decompose $f^{{\eps}} = f_1^{{\eps}} + f_2^{{\eps}}$ where

\[
\begin{cases}
\Delta^{{\eps}} f_1^{{\eps}}(a) = 0 \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
f_1^{{\eps}}(b) = {\meps} \phi(\eta(b)) \quad \mbox{for $b \in {\VGeps}(\partial B_r(z))$}
\end{cases}
\]

\[
\begin{cases}
\Delta^{{\eps}} f_2^{{\eps}}(a) = \psi(\eta(a)) \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
f_2^{{\eps}}(b) = 0 \quad \mbox{for $b \in {\VGeps}(\partial B_r(z))$}.
\end{cases}
\]

Step 1: Pointwise convergence.

Let a deterministic sequence $\mathcal{E}$ be given. Recall from <cit.> that both $f_1^{{\eps}}$ and $f_2^{{\eps}}$
can be expressed in terms of the Green's function of random walk.  Hence, by Lemma <ref> and Lemma <ref> respectively,
there is a deterministic subsequence $\mathcal{E}_1 \subset \mathcal{E}$ along which a.s.		$\overline{f}_1^{{\eps}}$ and $\overline{f}_2^{{\eps}}$
converge pointwise at every rational point in $B_r(z)$ to
$\overline{f}_1$ and $\overline{f}_2$ where

\[
\begin{cases}
\Delta \overline{f}_1 = 0 \quad\mbox{on $B_r(z)$} \\
\overline{f}_1  = \phi \quad \mbox{on $\partial B_r(z)$}
\end{cases}
\]

\[
\begin{cases}
\Delta \overline{f}_2 = m_0 \mu_h \times \psi \quad\mbox{on $B_r(z)$} \\
\overline{f}_2  = 0 \quad \mbox{on $\partial B_r(z)$}
\end{cases}
\]

Step 2: Local uniform convergence.

We now show local uniform convergence in $B_r(z)$. By the maximum principle (Lemma <ref>), $f_1^{{\eps}} \leq {\meps} \sup_{x \in B_r(z)} |\phi(x)|$.
Moreover, as

\[
f_2^{{\eps}} \pm \left( \sup_{x \in B_r(z)} |\psi(x)| \right) \times \Qeps{B_r(z)}
\]

is super (resp. sub) harmonic in ${\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$ and zero on ${\VGeps}(\partial B_r(z))$, the maximum principle also implies that $|f_2^{{\eps}}| \leq \left(\sup_{x \in B_r(z)} |\psi(x)| \right) \Qeps{B_r(z)}$
By the previous two sentences and Lemma <ref> we have that $\overline{f}_1^{{\eps}}$ and $\overline{f}_2^{{\eps}}$
are bounded by a random constant along a deterministic subsequence of $\mathcal{E}_1$. Thus, by Lemma <ref> there is a random subsubsequence $\mathcal{E}_2 \subset \mathcal{E}_1$ along which $\overline{f}_1^{{\eps}}$ converges uniformly in $\overline{B}_r(z)$ and  $\overline{f}_2^{{\eps}}$ locally uniformly in $B_r(z)$ to continuous functions $\overline{f}_1^* : \overline{B}_r(z) \to \R$ and $\overline{f}_2^* : \overline{B}_r(z) \to \R$.
By Step 1 we have that  $\overline{f}^*_1 = \overline{f}_1$
and $\overline{f}^*_2 = \overline{f}_2$.

Step 3: Uniform convergence.

We now argue that $\overline{f}_2^{{\eps}}$ converges uniformly in $\cl(B_r(z))$ to $\overline{f_2}$ by a standard barrier argument.

\begin{equation}
M_{\psi} = \sup_{a' \in B_r(z)} |\psi(a')|.
\end{equation}

and recall from (<ref>) that the function $\Qeps{B_r(z)}: {\VGeps}(\cl(B_r(z))) \to \R$

\[
\begin{cases}
\Delta^{{\eps}} \Qeps{B_r(z)}(a) = -1 \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
\Qeps{B_r(z)}(b) = 0 \quad \mbox{for $b \in {\VGeps}(\partial B_r(z))$}.
\end{cases}
\]

By the same argument as Steps 1 and 2, $M_{\psi} \times \Qeps{B_r(z)}$ converges locally uniformly in $B_r(z)$
along $\mathcal{E}_2$
to $M_{\psi} \times \overline{\mathfrak{Q}}_{B_r(z)}$ where

\[
\begin{cases}
\Delta \overline{\mathfrak{Q}}_{B_r(z)} = -1 \quad\mbox{on $B_r(z)$} \\
\overline{\mathfrak{Q}}_{B_r(z)}  = 0 \quad \mbox{on $\partial B_r(z)$}.
\end{cases}
\]

By, for example, <cit.>, $\overline{\mathfrak{Q}}_{B_r(z)}$ is Hölder continuous in $\cl(B_r(z))$.
Since $f^{{\eps}}_2 \pm M_{\psi} \times \Qeps{B_r(z)}$ is super (resp. sub) harmonic in ${\VGeps}(B_{r}(z)) \setminus {\VGeps}(\partial B_{r}(z))$ and is $0$ on ${\VGeps}(\partial B_{r}(z))$, we deduce convergence up to the boundary using the maximum principle exactly as in the proof of  Lemma <ref>.

§.§ Convergence of the normalized expected exit time

Recall that for a domain $D \subset B_{\rho}$, the normalized exit time is defined as $\qeps{D}(a)= \sum_{b \in {\VGeps}(D) \setminus {\VGeps}(\partial D)} \gr{{\eps}}{D}(a,b)$
and by (<ref>) solves the following discrete Dirichlet problem,

\[
\begin{cases}
\Delta^{{\eps}} \qeps{D}(\cdot) = -(\deg^{{\eps}})^{-1}(\cdot) \quad &\mbox{in ${\VGeps}(D) \setminus {\VGeps}(\partial D)$} \\
\qeps{D}(\cdot) = 0 \quad &\mbox{otherwise}.
\end{cases}
\]

As the function $(\deg^{{\eps}})^{-1}$ is not Hölder continuous, Lemma <ref> does not directly imply
that the normalized expected exit time converges. Nevertheless, we can use tightness and convergence of the Green's kernel to show the following.

Fix a deterministic choice of $z \in \C$ and $r > 0$ such that $B_r(z) \subset B_{\rho}$. As ${\eps} \to 0$, $\overline{\qeps{B_r(z)}}$ converges
uniformly in $\cl(B_r(z))$
to $\overline{\mathfrak q}_{B_r(z)}$
in probability

\[
\begin{cases}
\Delta \overline{\mathfrak q}_{B_r(z)} = -m_0/6 \times \mu_h \quad \mbox{on $B_r(z)$} \\
\overline{\mathfrak q}_{B_r(z)} = 0 \quad \mbox{ on $\partial B_{r}(z)$}.
\end{cases}
\]

Fix a deterministic choice of $z \in \C$ and $r > 0$ such that $B_r(z) \subset B_{\rho}$. We split the proof into steps.
Step 1: Identify a subsequence.

Let a sequence of ${\eps}$s tending to $0$ be given and fix a deterministic subsequence $\mathcal{E}$ along which the convergence result in Lemma <ref>
occurs in $\cl(B_r(z))$ and Lemmas <ref> and <ref> occur almost surely.
Also suppose that along this subsequence the convergence in Lemma <ref> occurs for
$\phi \equiv 0$ and for $\psi$ within a countable family of bump functions $\mathcal F(B_r(z))$ with the following property.
For each rational $\delta' >0$ and rational $x' \in B_r(z)$ such that $B_{\delta'}(x') \subset B_r(z)$, there is a function in $\mathcal F(B_r(z))$ which takes values in $[0,1]$, is 1 on $B_{\delta'}(x')$, and is 0 on $B_r(z) \setminus B_{2 \delta'}(x')$.

By Lemma <ref>
and our choice of subsequence, we have that that $\overline{\qeps{B_r(z)}}$
is uniformly bounded by a random constant along $\mathcal{E}$. Since $\deg^{{\eps}} \geq 1$, by (<ref>),
this function also has bounded $\Delta^{{\eps}}$-Laplacian. Thus, by Lemma <ref>, there
is a random subsequence $\mathcal{E}_0$ along which $\overline{\qeps{B_r(z)}}$ converges. In the remainder of the proof we use convergence
of the Green's kernel away from its pole to identify the subsequential limit.
Step 2: Convergence away from the singularity.

Let  $x \in B_r(z) \cap \Q^2$ and a rational $\delta >0$ such that $B_{\delta}(x) \Subset B_r(z)$ be given.
We would like to apply the uniform convergence of the Green's kernel away from its pole (Lemma <ref>).
So, we `carve out' the singularity at the pole of $\gr{{\eps}}{B_r(z)}(\cdot, \cdot)$ by considering,

\[
\begin{aligned}
\mathfrak q^{{\eps}, \delta^-}_{B_r(z)}(\cdot) &= \sum_{b \in {\VGeps}(B_r(z)) \setminus ({\VGeps}(\partial B_r(z)) \cup {\VGeps}(B_\delta(x)))} \gr{{\eps}}{B_r(z)}(\cdot, b) \\
\mathfrak q^{{\eps}, \delta^+}_{B_r(z)}(\cdot) &= \sum_{b \in {\VGeps}(B_\delta(x))} \gr{{\eps}}{B_r(z)}(\cdot, b).
\end{aligned}
\]

Observe that the same argument leading to (<ref>)
shows that

\begin{equation} \label{eq:trunc-laplacian-delta-}
\begin{cases}
\Delta^{{\eps}} \mathfrak q^{{\eps}, \delta^-}_{B_r(z)}(a) &= -(\deg^{{\eps}}(a))^{-1}  \\
&\times 1\{a \in {\VGeps}(B_r(z)) \setminus ({\VGeps}(\partial B_r(z)) \cup {\VGeps}(B_\delta(x))) \}\\
& \quad\mbox{on ${\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
\mathfrak q^{{\eps}, \delta^-}_{B_r(z)}  = 0 &\quad \mbox{otherwise.}
\end{cases}
\end{equation}

\begin{equation}\label{eq:trunc-laplacian-delta+}
\begin{cases}
\Delta^{{\eps}} \mathfrak q^{{\eps}, \delta^+}_{B_r(z)}(a) = -(\deg^{{\eps}}(a))^{-1} 1\{a \in {\VGeps}(B_\delta(x)) \} \quad\mbox{on ${\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
\mathfrak q^{{\eps}, \delta^+}_{B_r(z)}  = 0 \quad \mbox{otherwise}.
\end{cases}
\end{equation}

Apply the same argument as in the last paragraph of Step 1 to a countable set of rationals $\delta' > 0$ and $x' \in \Q^2$
to see that there is a random subsequence  of $\mathcal{E}_0$ (independent of $\delta$ and $x$)  along which $\overline{\qeps{B_r(z)}}$, $\overline{\mathfrak q}^{{\eps}, \delta^-}_{B_r(z)}$, $\overline{\mathfrak q}^{{\eps}, \delta^+}_{B_r(z)}$ converge uniformly in $\overline{B}_r(z)$ to continuous functions $\overline{\mathfrak q}^{*}_{B_r(z)}$, $\overline{\mathfrak q}^{*, \delta^-}_{B_r(z)}$, $\overline{\mathfrak q}^{*, \delta^+}_{B_r(z)}$ (with the limit possibly depending on the subsequence).

Fix such a subsequence $\mathcal{E}_1 \subset \mathcal{E}_0$ and limits $\overline{\mathfrak q}^{*}_{B_r(z)}$, $\overline{\mathfrak q}^{*, \delta^-}_{B_r(z)}$, $\overline{\mathfrak q}^{*, \delta^+}_{B_r(z)}$ and observe that

\begin{equation} \label{eq:asymptotic-equal}
\overline{\mathfrak q}^{*}_{B_r(z)} = \overline{\mathfrak q}^{*, \delta^-}_{B_r(z)} + \overline{\mathfrak q}^{*, \delta^+}_{B_r(z)}.
\end{equation}

Let $g^{{\eps}}(\cdot) := {\eps}^{-1} \gr{{\eps}}{B_r(z)}(a^{{\eps}}_x , \cdot)$ so that by Lemma <ref>,

\[
g^{{\eps}}(\cdot) \to (m_0/6) G_{B_r(z)}(x,\cdot)
\]

uniformly in $\cl(B_r(z)) \setminus B_{\delta}(x)$
as $\mathcal{E}_1 \ni {\eps} \to 0$.  As $g^{{\eps}}$ is zero on $\partial B_r(z)$,
it satisfies the Hölder continuity estimate (<ref>) from Lemma <ref> in $B_r(z) \setminus B_{\delta}(x)$.
This allows us to use $\overline{g}^{{\eps}}$ as a test function in Lemma <ref> so that

\[
\begin{aligned}
\int_{B_{r}(z) \setminus B_{\delta}(x)} \overline{g}^{{\eps}}(y) d \mu_h(y)
&\to \frac{m_0}{6} \int_{B_{r}(z) \setminus B_{\delta}(x)} G_{B_r(z)}(x, y) d \mu_h(y).
\end{aligned}
\]

as	$\mathcal{E}_1 \ni {\eps}  \to 0$.

By unpacking the definition of $\overline{g}^{{\eps}}(y)$ and using Lemma <ref>,
we see that the left-hand-side of the above converges to the same limit as $\overline{\mathfrak q}^{{\eps}, \delta^-}_{B_r(z)}(x)$.
This implies

\begin{equation} \label{eq:asymptotic-away}
\overline{\mathfrak q}^{*, \delta^-}_{B_r(z)}(x) = \frac{m_0}{6} \int_{B_{r}(z) \setminus B_{\delta}(x)} G_{B_r(z)}(x, y) d \mu_h(y).
\end{equation}

Step 3: Remove the singularity.

We claim that we may conclude once we show that

\begin{equation} \label{eq:asymptotic-small}
\lim_{\delta \to 0} \sup_{y \in B_r(z)} \overline{\mathfrak{q}}^{*, \delta^+}_{B_r(z)}(y) = 0.
\end{equation}

Indeed, (<ref>), (<ref>), and (<ref>) together imply that

\[
\overline{\mathfrak q}^{*}_{B_r(z)}(x) = \frac{m_0}{6} \int_{B_{r}(z)} G_{B_r(z)}(x, y) d \mu_h(y).
\]

Since $\Delta G_{B_r(z)}(x, \cdot) = -\delta_x(\cdot)$
in $B_r(z)$, the integral above
is equal to $\overline{\mathfrak{q}}_{B_r(z)}(x)$. This completes the proof as $x$ was an arbitrary rational and $\overline{\mathfrak{q}}^{*}_{B_r(z)}$ is continuous.

It remains to show (<ref>). Recall the definition of smooth bump functions $\mathcal{F}(B_r(z))$ from Step 1.
Let $\phi_{\delta} \in \mathcal{F}(B_r(z))$ be a smooth, positive bump function which is 1 on $B_{\delta}(x)$ and 0 on $B_{r}(z) \setminus B_{2 \delta}(x)$
and consider
the solution to the Dirichlet problem, $f^{{\eps}}_{\delta}: {\VGeps}(\cl(B_r(z))) \to \R$ defined by

\[
\begin{cases}
\Delta^{{\eps}} f^{{\eps}}_{\delta}(a) = \phi_{\delta}(\eta(a)) \quad \mbox{for $a \in {\VGeps}(B_{r}(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
f^{{\eps}}_{\delta}(b) = 0  \quad \mbox{for $b \in {\VGeps}(\partial B_r(z))$}.
\end{cases}
\]

Since $q^{{\eps}, \delta^+}_{B_r(z)} \pm f^{{\eps}}_{\delta}$ is sub (resp. super) harmonic in ${\VGeps}(B_{r}(z)) \setminus {\VGeps}(\partial B_{r}(z))$ and is $0$ on ${\VGeps}(\partial B_{r}(z))$, by the maximum principle,

\begin{equation} \label{eq:bound-q-by-f}
|\mathfrak q^{{\eps}, \delta^+}_{B_r(z)}| \leq |f^{{\eps}}_{\delta}|, \quad \mbox{in ${\VGeps}(B_{r}(z))$}.
\end{equation}

By Lemma <ref>, $\overline{f}^{{\eps}}_{\delta}$ converges uniformly in $B_r(z)$ to $\overline{f}_{\delta}$ where

\[
\begin{cases}
\Delta \overline{f}_{\delta} = m_0 \phi_{\delta} \times \mu_h \quad \mbox{on $B_r(z)$} \\
\overline{f}_{\delta} = 0  \quad \mbox{on $\partial B_r(z)$}.
\end{cases}
\]

We can represent

\begin{align*}
\overline{f}_{\delta}(\cdot) &= m_0 \int_{B_r(z)}  \phi_{\delta}(y) G_{B_r(z)}(\cdot, y) d \mu_h(y) \\
&\leq m_0 \int_{B_{2 \delta}(z)}  G_{B_1}(\cdot, y) d \mu_h(y) \\
&\leq m_0 \times \Mrho{B_{2 \delta}(z)}.
\end{align*}

Hence, by the second half of Lemma <ref>
there exists $q = q(\gamma) > 0$ such that

\[
\overline{f}_{\delta} \leq \delta^q  \quad \mbox{for all $\delta$ sufficiently small},
\]

which by (<ref>), implies (<ref>), completing the proof.

§ CONVERGENCE OF THE DIVISIBLE SANDPILE ODOMETER

Fix $\rho \in (0,1)$  and $t > 0$. We briefly recall the single-source divisible sandpile model, introduced in Section <ref>, and define its odometer.
In the model, we start with a continuous amount of mass $(t {\eps}^{-1})$ at the origin $0 \in {\VGeps}$ and zero elsewhere. At each time step, vertices $a \in {\VGeps}$ which have mass $\sigma(a) > 1$, where $\sigma$ is the current configuration, are unstable and topple, distributing the excess mass $(\sigma(a) - 1)$
equally among the neighbors of $a$.
Since ${\Geps}$ is an infinite graph eventually each site $a$ will have mass at most $1$.
Denote the final configuration of mass by $s^{{\eps}}_{t}: {\VGeps} \to [0, 1]$.

The odometer, which we denote by $v^{{\eps}}_{t}: {\VGeps} \to [0, \infty)$, tracks the total amount of mass each site has emitted during this process.
In particular,  by definition of the graph Laplacian, we have that

\begin{equation} \label{eq:final-sandpile-definition}
s^{{\eps}}_{t} = t {\eps}^{-1} \delta_0 +  \deg^{{\eps}} \times \Delta^{{\eps}} \left( \frac{v^{{\eps}}_{t}}{\deg^{{\eps}}} \right).
\end{equation}

(The two $\deg^{{\eps}}$ factors in the above expression are needed to ensure the above expression matches the prior description of toppling.)

By definition, we see that the divisible sandpile cluster satisfies

\begin{equation} \label{eq:cluster-definition}
D_{{\eps}}(t {\eps}^{-1})  = \{ x \in {\VGeps} : s^{{\eps}}_{t} > 0 \}
= \cl( \{ x \in {\VGeps} : v_t^{{\eps}} > 0 \}).
\end{equation}

The reason for the closure in (<ref>) is that there are some sites with mass
in the cluster which have not toppled. These sites, by definition, must have a neighbor which has toppled.

A fundamental property in the study of the sandpile is the least action principle,

\begin{equation} \label{eq:least-action-principle}
\frac{v^{{\eps}}_{t}}{\deg^{{\eps}}} = \min \{ w:{\VGeps} \to [0,\infty) : \Delta^{{\eps}} w + {\deg^{{\eps}}(0)}^{-1} t {\eps}^{-1} \delta_0 \leq (\deg^{{\eps}})^{-1}\},
\end{equation}

where the minimum is pointwise. See, , <cit.> for a proof. This property is closely related to (<ref>).
In fact, one sees from the definition that if $D_{{\eps}}(t {\eps}^{-1})  \subset {\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})$,
then we can express, using (<ref>), the solution to
(<ref>) as

\begin{equation} \label{eq:odometer-to-normalized-odometer}
\frac{v^{{\eps}}_{t}(a)}{\deg^{{\eps}}(a)} =
\begin{cases}
\wteps(a) + (t {\eps}^{-1}) \times \gr{{\eps}}{B_\rho}(0, a) \quad \mbox{for $a \in {\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})$} \\
0 \quad \mbox{otherwise}
\end{cases}
\end{equation}

\begin{equation} \label{eq:discrete-least-supersolution}
\begin{aligned}
\wteps := &\min \{ w: {\VGeps}(\cl(B_{\rho})) \to \R : \Delta^{{\eps}} w \leq (\deg^{{\eps}})^{-1} \mbox{ on ${\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})$} \\
&\quad \mbox{ and $w \geq -(t {\eps}^{-1}) \times  \gr{{\eps}}{B_\rho}(0, \cdot)$}\},
\end{aligned}
\end{equation}

where the minimum is pointwise. And, in general,

\begin{equation} \label{eq:lower-bound-odometer}
\frac{v^{{\eps}}_{t}}{\deg^{{\eps}}} \geq \wteps + (t {\eps}^{-1}) \times \gr{{\eps}}{B_\rho}(0, \cdot) \quad \mbox{on ${\VGeps}(B_{\rho})$}.
\end{equation}

Let $\overline \wteps$ be the re-scaled version of $\wteps$ as in (<ref>).
We will show the convergence of $\overline \wteps$ to its continuum counterpart defined by (<ref>). In order to ensure the divisible sandpile cluster is contained in $B_{\rho}$, we consider
$t \in (0, \TT)$ where

\begin{equation} \label{eq:cluster-stopping-time}
\TT = \sup \{t > 0 : \Lambda_t \subset B_{\rho/3} \},
\end{equation}

where $\Lambda_t$ is the LQG harmonic ball as in Theorem <ref>.
We take $B_{\rho/3}$ rather than $B_{\rho}$ in the above definition
to allow some room for error in the discrete estimates of Section <ref>.

Recall the median exit time $m_0$ defined in (<ref>).
For each $t > 0$, on the event $\{t < \TT\}$ we have that $\overline{\wteps}$ converges uniformly in
$B_{\rho}$ as ${\eps} \to 0$ to $\frac{m_0}{6} \overline{w}^{B_\rho}_t$
in probability, where $\overline w_t^{B_\rho}$ is as in Theorem <ref>.

The extra factor of $\frac{m_0}{6}$ in Theorem <ref> is a result of the factor ${\meps}$
appearing in the scaling definition (<ref>) and the fact that since the mated-CRT map is a planar triangulation,
$\E[\deg^{{\eps}}(0)] = 6$.

A result similar to that of Theorem <ref> was proved on $\Z^d$ by Levine and Peres in <cit.>
using different methods. Levine and Peres used precise estimates on the lattice Green's function in $\Z^d$ which are unavailable in our setting.
Our approach is similar to (but easier than) the method used to prove convergence of the Abelian sandpile — see <cit.>. (Also see <cit.> for a systematic approach to proving convergence of finite-difference schemes.)

Before proceeding we point out that convergence of the odometer implies a lower bound
on the divisible sandpile cluster.

For each $\delta \in (0, \rho)$ and $t > 0$, on the event $\{t < \TT\}$,
it holds except on an event of probability tending to 0 as ${\eps} \to 0$  that

\[
\Bdeltam(\Lambda_t) \subset \overline{D}_{{\eps}}(t {\eps}^{-1}).
\]

Fix $\delta \in (0, \rho)$ and $t > 0$ and restrict to the event that $\{t < \TT\}$.
Let a deterministic sequence of ${\eps}$s approaching zero be given. Recall from Theorem <ref> that

\begin{equation} \label{eq:recall-cluster-def}
\Lambda_t = \inte \left( \overline{\{ x \in B_\rho : \overline w_t^{B_\rho} > - t G_{B_\rho}(0, \cdot)\}} \right),
\end{equation}

where $G_{B_\rho}$ is the Green's function for $\Delta$ on $B_{\rho}$.
Since $\overline w_t^{B_\rho}$ is continuous and $G_{B_\rho}(0,\cdot)$ blows up at the origin, we can pick $r \in (0,1)$ so that

\begin{equation} \label{eq:green-singularity-blow-up}
\inf_{x \in B_{2 r}} t G_{B_\rho}(0, x) > \sup_{y \in B_{r}} -\overline w_t^{B_\rho}(y) .
\end{equation}

Since $\gr{{\eps}}{B_\rho}$ is superharmonic, by the maximum principle, Lemma <ref>,

\begin{equation} \label{eq:singularity-blow-up}
\inf_{a \in {\VGeps}(B_{r'})}\gr{{\eps}}{B_\rho}(0, a)
\inf_{b \in {\VGeps}(\partial B_{r'})}\gr{{\eps}}{B_\rho}(0, b), \quad \forall r' \in (0, \rho).
\end{equation}

Let $g^{{\eps}}(\cdot) := {\eps}^{-1} \gr{{\eps}}{B_\rho}(0 , \cdot)$ and define $\overline{g}^{\eps}$ as in (<ref>).
By Theorem <ref>, there exists a deterministic subsequence, $\mathcal{E}$, along which a.s.		$\overline{\wteps}$  converges uniformly in $B_{\rho}$ to $\frac{m_0}{6} \overline w_t^{B_\rho}$.
Further, by Lemma <ref> we may arrange so that, a.s. along $\mathcal{E}$,
$\overline{g}^{{\eps}}$ converges uniformly in $\A_{r, \rho}$ to $\frac{m_0}{6} \times G_{B_\rho}(0 , \cdot)$.  By (<ref>) and the continuity of the two functions involved, there exists $s > 0$ so that

\begin{equation} \label{eq:strictly-above-in-interior}
\begin{aligned}
\overline w_t^{B_\rho}(x) &> - t G_{B_\rho}(0, x) + s  \\
&\mbox{for all $x \in \Lambda_t$ such that $\dist(x, \partial \Lambda_t) \geq \delta/4$}.
\end{aligned}
\end{equation}

We may further assume that a.s., the cell-size estimates of Lemma <ref> hold for all sufficiently small ${\eps}\in \mathcal{E}$.

Therefore, by the uniform convergence of $\overline{\wteps}$ in $B_{\rho}$ and the
uniform convergence of $\overline{g}^{{\eps}}$ in $\A_{r, \rho}$ 	we have by (<ref>)

\[
\wteps(a)> \sup_{b \in {\VGeps}(\partial B_{2 r})}  -(t {\eps}^{-1}) \times \gr{{\eps}}{B_\rho}(0, b) , \quad \forall a \in {\VGeps}(B_{r}), \quad \forall   {\eps} \in\mathcal E \mbox{ sufficiently small}
\]

and thus by (<ref>)

\[
\wteps(a)> - (t {\eps}^{-1}) \times \gr{{\eps}}{B_\rho}(0, a) , \quad \forall a \in {\VGeps}(B_{r}), \quad \forall   {\eps} \in\mathcal E \mbox{ sufficiently small.}
\]

Also, by uniform convergence of $\overline{\wteps}$ in $B_{\rho}$ and
uniform convergence of $\overline{g}^{{\eps}}$ in $\A_{r, \rho}$ 	and (<ref>) we have

\[
\begin{aligned}
\wteps(a) &> - (t {\eps}^{-1}) \times \gr{{\eps}}{B_\rho}(0, a)  \\
&\mbox{$\forall a \in {\VGeps}(\Lambda_t \setminus B_{r})$ such that $\dist(\eta(a), \partial \Lambda_t) \geq \delta/2$} \\
&\mbox{$\forall {\eps} \in \mathcal{E} $ sufficiently small.}
\end{aligned}
\]

Combining the previous two indented equations with the definition of $\overline{D}_{{\eps}}$, (<ref>), and (<ref>)  completes the proof.

§.§ Basic properties of the divisible sandpile

For completeness, we recall some basic properties of the discrete least supersolution.
Since these results are standard, the proofs are given in Appendix <ref>.
We also note that continuum analogs of these properties are proven in <cit.>.

Recall the definition of $\qeps{B_\rho}$ from (<ref>).
For all $t > 0$ and ${\eps} > 0$,

\[
\max(-\qeps{B_{\rho}}, -(t {\eps}^{-1}) \times \gr{{\eps}}{B_\rho}(0, \cdot)) \leq \wteps \leq 0
\]

and  $\Delta^{{\eps}} \wteps \leq (\deg^{{\eps}})^{-1}$ on ${\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})$.

For the statement of the next lemma, we denote the discrete cluster by

\begin{equation} \label{eq:discrete-cluster}
\Lambdateps = \{ a \in {\VGeps}(B_{\rho}) :  \wteps(a) > - (t {\eps}^{-1}) \gr{{\eps}}{B_{\rho}}(0, a)\}.
\end{equation}

Observe that by (<ref>) for $t \in (0,\TT)$,

\begin{equation} \label{eq:cluster-and-lambdat}
D_{{\eps}}(t {\eps}^{-1})  = \cl(\Lambdateps).
\end{equation}

For all $t > 0$ and ${\eps} > 0$, we have that $\Delta^{{\eps}} \wteps = (\deg^{{\eps}})^{-1}$ on $\Lambdateps$
and $0 \leq \Delta^{{\eps}} \wteps \leq (\deg^{{\eps}})^{-1}$ on ${\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})$.
Moreover,  $\Lambdateps$ is connected and if $ t {\eps}^{-1}  > 1$,  then $0 \in \Lambdateps$.

We conclude with a discrete conservation of mass lemma.

For all $t > 0$ and ${\eps} > 0$, we have that $|\Lambdateps| \leq   t {\eps}^{-1} $. Moreover,
if $\Lambdateps \subset \inte({\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho}))$,

\[
\sum_{a \in {\VGeps}(B_{\rho}) \setminus {\VGeps}(\partial B_{\rho})} \Delta^{{\eps}} \wteps(a) \deg^{{\eps}}(a) =  t {\eps}^{-1} .
\]

§.§ Convergence of the odometer along subsequences

In this subsection we prove tightness of the re-scaled odometer function $\overline \wteps$ defined in (<ref>) and (<ref>). In subsequent subsections we show that, for $t \in (0, \TT)$, each subsequential limit is given uniquely by (<ref>).

Almost surely, for each $t > 0$ and each deterministic sequence of ${\eps}$s converging to 0,
there exists a $\sigma(h,\eta)$-measurable subsequence
and a Hölder continuous function $\overline{w}^*_t \in C(\cl(B_\rho))$ such that $\overline{\wteps} \to \overline{w}^*_t$
uniformly in $\cl(B_{\rho})$.

We will prove Lemma <ref> by combining Lemma <ref>, the results of the previous subsection,
and the following lemma.

A.s., for each  deterministic sequence of ${\eps}$s converging to 0, there exists a deterministic subsequence $\mathcal{E}$ and a random constant $C$ so that for each $t > 0$,

\[
\sup_{x \in B_{\rho}} |\overline{\wteps}(x)| \leq C, \quad \forall {\eps} \in \mathcal{E}.
\]

By Lemma <ref>,

\begin{equation} \label{eq:odometer-bounds}
-\qeps{B_{\rho}} \leq \wteps \leq 0, \quad \forall t > 0, \forall {\eps} > 0
\end{equation}

and by Lemma <ref> and (<ref>)

\begin{equation}
-\qeps{B_{\rho}} \geq -C_0 {\meps}  \Mrho{B_{\rho}}
\end{equation}

for some deterministic constant $C_0 = C_0(\rho, \gamma)$
with polynomially high probability as ${\eps} \to 0$, where $\Mrho{B_{\rho}}$ is as in (<ref>).
By Lemma  <ref>,  $\Mrho{B_{\rho}}$ is finite almost surely.
Hence, the previous two indented inequalities complete the proof by the definition of $\overline{\wteps}$ and the Borel-Cantelli lemma.

We now use the tightness lemma,  Lemma <ref>, to prove subsequential convergence.

By Lemma <ref> and Lemma <ref>, for each $t>0$, the assumptions of Lemma <ref>
are satisfied along a deterministic subsequence for the functions $\wteps$. Hence, $\overline{\wteps}$ converges a.s. locally uniformly in $B_{\rho}$ along a random subsequence. To see that the convergence occurs uniformly in $\overline{B_{\rho}}$, we use the standard barrier argument given in Step 3 of the proof of Lemma <ref>. To implement the barrier argument we use that $0 \leq \Delta^{{\eps}} w^{{\eps}}_t \leq 1$.

§.§ Properties of the subsequential limit

Convergence of the discrete Dirichlet problem, Lemma <ref>, is used
to establish some properties of each subsequential limit of the odometer.
We will use this in the next subsection to uniquely identify the limit.

A.s., for each deterministic sequence of ${\eps} \to 0$,
there exists a deterministic subsequence so that a subsequential limit $\overline{w}^*_t$ of $\overline \wteps$ along this subsequence (as given by Lemma <ref>) satisfies the following:

* $0 \geq \overline{w}^*_t \geq -t \frac{m_0}{6} G_{B_{\rho}}(0, \cdot)$ on $\cl(B_{\rho})$;

* $\Delta \overline{w}^*_t \geq 0$ on $B_{\rho}$;

* $\Delta \overline{w}^*_t \leq \frac{m_0}{6} \mu_h$ on $B_{\rho}$ and $\Delta \overline{w}^*_t = \frac{m_0}{6} \mu_h$ on $\{x \in B_{\rho} : \overline{w}^*_t(x) > -t \frac{m_0}{6} G_{B_{\rho}}(0, x)\}$.

By Lemma <ref>, we may take a deterministic subsequence along which,
the convergence of 	$\overline g^{\eps}$ from that lemma, with $y = 0$ and $B_r(z) = B_{\rho}$, occurs almost surely.
Take another subsequence
for which the convergence of Lemma <ref>
occurs a.s. for all $z \in \Q^2$ and rational $r > 0$,  $\psi \equiv 0$, and
all $\phi \in \mathcal{F}(B_r(z))$, the set of polynomials with rational coefficients
restricted to $\cl(B_r(z))$. Further suppose that along this subsequence the convergence in Lemma <ref> holds for all
$z \in \Q^2$ and rational $r > 0$.  By Lemma <ref>, there is a random subsequence $\mathcal E$ of the preceding subsequence such that $\lim_{\mathcal{E} \ni {\eps} \to 0} \overline{\wteps} \to \overline{w}^*_t$.
Step 1: Proof of (1).

By Lemma <ref>,

\begin{equation}
0 \geq \wteps \geq -(t {\eps}^{-1}) \gr{{\eps}}{B_{\rho}}(0, \cdot) \quad \mbox{on ${\VGeps}(\cl(B_{\rho}))$} , \quad \forall {\eps}>0.
\end{equation}

This inequality persists in the limit by the uniform convergence of $\wteps$ and the convergence of the rescaled version of $m_{{\eps}} {\eps}^{-1} \gr{B_{\rho}}{{\eps}}$ to $\frac{m_0}{6} G_{B_{\rho}}$, for $\mathcal{E} \ni {\eps} \to 0$ as given by Lemma <ref>.
Step 2: Proof of (2).

Let $\delta > 0$, a rational $r > 0$, and a $z \in \Q^2$ such that $B_r(z) \Subset B_{\rho}$ be given.
For each ${\eps} \in \mathcal{E}$ consider the discrete harmonic function $f_\delta^{\eps}$ defined by

\[
\begin{cases}
\Delta^{{\eps}} f^{{\eps}}_{\delta}(a) = 0 \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
f^{{\eps}}_{\delta}(b)  = m_{{\eps}} \phi^{\delta}(\eta(b)) + \delta \quad \mbox{for $b  \in {\VGeps}(\partial B_r(z))$}
\end{cases}
\]

where $\phi^{\delta} \in \mathcal{F}(B_r(z))$ is a rational polynomial satisfying

\begin{equation} \label{eq:approximation-by-delta}
\sup_{x \in \cl(B_r(z))} |\phi^{\delta}(x) - \overline{w}^*_t(x)| < \delta/4.
\end{equation}

By Lemma <ref>,  $\wteps$ is $\Delta^{{\eps}}$-subharmonic.
Thus, by the maximum principle, Lemma <ref>, we have

\[
\sup_{a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))}(\wteps - f^{{\eps}}_{\delta})(a) \leq \sup_{b  \in {\VGeps}(\partial B_r(z))}
(\wteps - f^{{\eps}}_{\delta})(b).
\]

As $\overline{\wteps}$ converges uniformly to $\overline{w}^*_t$ in $\cl(B_{\rho})$, by definition of $f_\delta^{\eps}$
and (<ref>), we have

\[
\sup_{b  \in {\VGeps}(\partial B_r(z))}
(\wteps - f^{{\eps}}_{\delta})(b) < 0, \quad \forall \mbox{ small ${\eps} \in \mathcal{E}$}.
\]

By combining the previous two indented equations, we have

\begin{equation} \label{eq:lower-bound}
f^{{\eps}}_{\delta}(a) \geq \wteps(a), \quad \forall a \in {\VGeps}(B_r(z)), \quad \forall \mbox{ small ${\eps} \in \mathcal{E}$}.
\end{equation}

By Lemma <ref>, we have that $\overline f^{{\eps}}_{\delta}$ converges uniformly in
$\cl(B_r(z))$ to a harmonic function $\overline f^{*}_{\delta}: \cl(B_r(z)) \to \R$ which satisfies

\begin{equation} \label{eq:harmonic-function-same-bdry-values}
\begin{cases}
\Delta \overline f^{*}_{\delta} = 0 \quad\mbox{on $B_r(z)$} \\
\overline f^{*}_{\delta}  =  \phi^{\delta} + \delta \quad \mbox{on $\partial B_r(z)$}.
\end{cases}
\end{equation}

\begin{align*}
\overline{w}^*_t(z) &\leq \overline f^{*}_{\delta}(z) \qquad \mbox{(by \eqref{eq:lower-bound} and uniform convergence)} \\
&= \int_{\partial B_{r}(z)} \overline{f}^*_{\delta} \qquad \mbox{(by harmonicity)} \\
&= \int_{\partial B_{r}(z)} (\phi^{\delta} + \delta) \qquad \mbox{(by \eqref{eq:harmonic-function-same-bdry-values})}.
\end{align*}

As the choice of rational $r > 0$ and $z \in \Q^2$ was arbitrary and the above inequality holds for any $\delta > 0$, by (<ref>),
and the continuity of $\overline{w}^*_t(z)$, we have

\[
\overline{w}^*_t(z) \leq \int_{\partial B_{r}(z)} \overline{w}^*_t, \quad \mbox{for all $r > 0$ and $z \in \C$ such that $B_r(z) \subset B_{\rho}$},
\]

completing the proof.
Step 3: Proof of (3).

We will show that

\begin{equation} \label{eq:limit-admissible}
\Delta \overline{w}^*_t \leq \frac{m_0}{6} \mu_h  \quad \mbox{on $B_{\rho}$},
\end{equation}

and then

\begin{equation} \label{eq:limit-equal}
\Delta \overline{w}^*_t = \frac{m_0}{6} \mu_h  \quad \mbox{on $\{\overline{w}^*_t > -t \frac{m_0}{6} G_{B_{\rho}}\}$},
\end{equation}

appealing to the argument of Step 2.

Let rational $r > 0$ and $z \in \Q^2$ such that $B_r(z) \subset B_{\rho}$ and $\delta > 0$ be given. Recall that for each ${\eps} > 0$,
$\qeps{B_r(z)}$ satisfies

\[
\begin{cases}
\Delta^{{\eps}} \qeps{B_r(z)}(a) = -(\deg^{{\eps}}(a))^{-1} \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$} \\
\qeps{B_r(z)}(b)  = 0 \quad \mbox{for $b  \in {\VGeps}(\partial B_r(z))$}.
\end{cases}
\]

By Lemma <ref>,

\begin{equation} \label{eq:difference-super-harmonic}
\Delta^{{\eps}} (\qeps{B_r(z)} + \wteps)(a) \leq 0 \quad\mbox{for $a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))$}.
\end{equation}

By Lemma <ref>, $\overline{\qeps{B_r(z)}}$ converges uniformly
in $\cl(B_r(z))$ to $\overline{\mathfrak q}_{B_r(z)}$

\begin{equation} \label{eq:limit-equation-laplacian}
\begin{cases}
\Delta \overline{\mathfrak q}_{B_r(z)} = -m_0/6 \times \mu_h \quad \mbox{on $B_r(z)$} \\
\overline{\mathfrak q}_{B_r(z)} = 0 \quad \mbox{ on $\partial B_{r}(z)$}.
\end{cases}
\end{equation}

As $(\qeps{B_r(z)} + \wteps)$ is discrete superharmonic (<ref>), the same argument as in Step 2 shows that

\begin{equation} \label{eq:limit-equation-sub-harmonic}
\Delta (\overline{\mathfrak q}_{B_r(z)} + \overline{w}^{*}_t) \leq 0 \quad\mbox{on $B_r(z)$}.
\end{equation}

Combining this with (<ref>) shows (<ref>).

We now show (<ref>). Fix rational $r > 0$ and $z \in \Q^2$ so that

\[
B_r(z) \subset \{x \in B_{\rho}: \overline{w}^*_t(x) > -t \frac{m_0}{6} G_{B_{\rho}}(x)\}.
\]

By the continuity and uniform convergence of the two functions involved, for all ${\eps} \in \mathcal{E}$ sufficiently small, we have that

\begin{equation} \label{discrete-cluster-contained}
\wteps(a) > -(t {\eps}^{-1}) \gr{B_{\rho}}{{\eps}}(0,a), \quad \forall a \in {\VGeps}(B_{r}(z)).
\end{equation}

The inequality (<ref>)
and Lemma <ref> together show that

\begin{equation} \label{eq:difference-harmonic}
\Delta^{{\eps}} (\qeps{B_r(z)} + \wteps)(a) = 0 \quad \forall a \in {\VGeps}(B_r(z)) \setminus {\VGeps}(\partial B_r(z))
\end{equation}

for all ${\eps} \in \mathcal{E}$ sufficiently small.
Again applying the argument of Step 2 twice, we see that

\begin{equation} \label{eq:limit-equation-harmonic}
\Delta (\overline{\mathfrak q}_{B_r(z)} + \overline{w}^{*}_t) = 0 \quad\mbox{on $B_r(z)$},
\end{equation}

completing the proof by (<ref>).

§.§ Proof of Theorem <ref>

Note that by the definition of $\TT$ and Theorem <ref>,

\[
\Lambda_t = \inte\left( \cl \left( \{x \in B_\rho : \overline{w}^{B_\rho}_t > - G_{B_{\rho}}(0, \cdot) \} \right) \right), \quad \forall t \in (0, \TT).
\]

By (1) and (3) of Lemma <ref> we have that each possible subsequential limit $\frac{6}{m_0} \overline{w}^*_t$ as in Lemma <ref> is admissible in
and hence

\begin{equation} \label{eq:admissible-supersolution}
\frac{m_0}{6} \overline{w}^{B_\rho}_t \leq \overline{w}^*_t, \quad \forall t \in (0, \TT).
\end{equation}

For the other direction, we first note that by (<ref>),

\[
\Lambda_t \subset   \inte \left( \cl \left( \{x \in B_\rho : \overline{w}^*_t > - \frac{m_0}{6} G_{B_{\rho}}(0, \cdot) \} \right) \right), \quad \forall t \in (0, \TT).
\]

Therefore, for each $t \in (0,\TT)$, by Lemma <ref> and (3) of Theorem <ref>,

\[
\begin{cases}
\Delta(\overline{w}^*_t - \frac{m_0}{6} \overline{w}^{B_\rho}_t) \geq 0 \quad \mbox{on $B_\rho$} \\
\overline{w}^*_t - \frac{m_0}{6} \overline{w}^{B_\rho}_t = 0 \quad \mbox{on $\partial B_\rho$}
\end{cases}
\]

and hence by the maximum principle in the continuum,

\begin{equation} \label{eq:admissible-subsolution}
\frac{m_0}{6} \overline{w}^{B_\rho}_t \geq \overline{w}^*_t
\end{equation}

which, together with (<ref>) shows that $\overline{w}^*_t = \overline{w}^{B_\rho}_t$.

§ IDLA LOWER BOUND

We use the convergence of the divisible sandpile and the argument of Lawler-Bramson-Griffeath <cit.>
to prove the following lower bound on the IDLA cluster. As in the previous section we fix a $\rho \in (0,1)$ so that $B_{\rho} \Subset B_1$.
Recall the notation for the subset of $\C$ corresponding to a subset of ${\VGeps}$, (<ref>).

Recall the time $\TT$ from (<ref>).
For each $\delta \in (0, \rho)$ and $t > 0$, on the event $\{t < \TT\}$,
it holds except on an event of probability tending to 0 as ${\eps} \to 0$  that

\[
\Bdeltam(\Lambda_t) \subset \overline{A}_{{\eps}}(\lfloor t {\eps}^{-1} \rfloor)
\]

where $\Lambda_t$ is as in Theorem <ref>.
Moreover, this lower bound holds even when walkers are stopped upon
exiting ${\VGeps}(\Lambda_t)$.

Our proof of Proposition <ref> relies on the following result of Lawler-Bramson-Griffeath <cit.>.
Recall Remark <ref>.

Suppose  $O_1^{{\eps}} \subset O_2^{{\eps}} \subset {\VGeps}(B_{\rho})$ are $\sigma(h,\eta)$-measurable domains containing the origin
such that for some deterministic $q > 0$ and random $\sigma(h,\eta)$-measurable finite random variable $\mathfrak{s} > 0$,

\begin{equation} \label{eq:mean-value-ineq}
|O_2^{{\eps}}| \Gr{{\eps}}{O_2^{{\eps}}}(0, a) \geq (1 + \mathfrak{s}) \sum_{b \in O_2^{{\eps}}} \Gr{{\eps}}{O_2^{{\eps}}}(b, a), \quad \forall a \in O_1^{{\eps}},
\end{equation}

\begin{equation} \label{eq:L-lower-bound}
\frac{ \sum_{b \in O^{{\eps}}_2} \Gr{{\eps}}{O_2^{{\eps}}}(b, a) }{ {\Gr{{\eps}}{O_2^{{\eps}}}(a, a)}} \geq {\eps}^{-q}, \quad \forall a \in O_1^{{\eps}},
\end{equation}

with probability approaching one as ${\eps}$ goes to zero,  then

\begin{equation} \label{eq:lower-bound-arbitrary}
O_1^{{\eps}} \subset A_{{\eps}}(|O_2^{{\eps}}|)
\end{equation}

with probability approaching one as ${\eps}$ goes to zero.
Moreover, the bound (<ref>) holds even when walkers are stopped upon exiting $O_2^{{\eps}}$.

We will give the proof of Lemma <ref> at the end of this section, following the argument of <cit.>. We will first prove Proposition <ref> by verifying the two conditions of Lemma <ref>.
Specifically, we will show that the mean-value inequality (<ref>) is satisfied for the discrete divisible sandpile cluster $\Lambdateps$.

For each $\delta \in (0, \rho)$ and $t > 0$ on the event $\{t < \TT\}$,
there exists a $\sigma(h,\eta)$-measurable finite random variable $\mathfrak{s} = \mathfrak{s}(\delta) > 0$ such that, except on an event of probability tending to 0 as ${\eps} \to 0$ ,

\begin{equation} \label{eq:mean-value-ineq-verify}
|\Lambdateps| \Gr{{\eps}}{\Lambdateps}(0, a) \geq (1 + \mathfrak{s}) \sum_{b \in \Lambdateps} \Gr{{\eps}}{\Lambdateps}(b, a), \quad  \mbox{for all $a \in \Lambdateps$ s.t. $\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta$}.
\end{equation}

We will also prove the lower bound in  (<ref>).

There exists a deterministic $q > 0$ such that, for each $\delta \in (0, \rho)$ and $t > 0$ on the event $\{t < \TT\}$,
it holds except on an event of probability tending to 0 as ${\eps} \to 0$  that

\begin{equation} \label{eq:L-lower-bound-verify}
\frac{\sum_{b \in \Lambdateps} \Gr{{\eps}}{\Lambdateps}(b, a)  }{ {\Gr{{\eps}}{\Lambdateps}(a, a)}} \geq {\eps}^{-q},
\quad  \mbox{for all $a \in \Lambdateps$ s.t. $\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta$}.
\end{equation}

These lemmas imply Proposition <ref>.

We restrict to the event that $\{t < \TT\}$ for $t > 0$.
By Lemmas <ref> and  <ref> (with $\delta/4$ instead of $\delta$) we may apply Lemma <ref>
with the sets $O^{{\eps}}_1 = \{a \in \Lambdateps:\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta/4 \}$ and $O_2^{{\eps}} = \Lambdateps$ which shows, together
with the cell-size estimates given by
Lemma <ref>

\[
\Bdeltamm{\delta/2}(\overline{\Lambdateps}) \subset \overline{A}_{{\eps}}(|\Lambdateps|),
\]

except on an event of probability going to zero as ${\eps} \to 0$,  with the bound holding even when walkers are stopped upon exiting $\Lambdateps$.
By Proposition <ref> and (<ref>), (and the cell-size estimates again) we have

\[
\Bdeltamm{\delta/2}(\Lambda_t) \subset \overline{\Lambdateps},
\]

except on an event of probability going to zero as ${\eps} \to 0$. This completes the proof by the upper bound on $|\Lambdateps|$ given by Lemma <ref>.

We now prove Lemmas <ref>, <ref>, and <ref>.
We start with Lemma <ref>, which relies on the estimates established in <cit.> as recalled in Section <ref>.

We restrict to the event that $\{t < \TT\}$ for $t > 0$.
By the fact $\Lambdateps \subset {\VGeps}(B_{\rho})$ and Lemma <ref>, we have

\begin{equation}
\Gr{{\eps}}{\Lambdateps}(a,a) \leq \Gr{{\eps}}{B_{\rho}}(a,a) \leq {\eps}^{-1/2}, \quad \forall a \in {\VGeps}(B_{\rho})
\end{equation}

except on an event of probability going to zero as ${\eps} \to 0$.
Therefore, it suffices to show that there exists a deterministic constant $C > 0$ independent of ${\eps}$ such that

\begin{equation} \label{eq:lower-bound-exit}
\sum_{b \in \Lambdateps} \Gr{{\eps}}{\Lambdateps}(b,a) \geq C \delta {\eps}^{-1}, \quad \forall
a \in \Lambdateps:\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta
\end{equation}

except on an event of probability going to zero as ${\eps} \to 0$. Indeed, we just need the desired estimate for ${\eps}$ small
and we may absorb the factor of $C \delta$ into ${\eps}^{1/2}$ for ${\eps}$ small enough, depending on $C$ and $\delta$.

Let $a \in \Lambdateps$ such that $\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta$ be given
and observe that

\[
\sum_{b \in \Lambdateps}  \Gr{{\eps}}{\Lambdateps}(b, a) \geq
\sum_{b \in {\VGeps}(B_{\delta/2}(a))} \Gr{{\eps}}{B_{\delta/2}(a)}(b,a).
\]

By Lemma <ref> and Lemma <ref>

\[
\sum_{b \in {\VGeps}(B_{\delta/2}(a))}  \Gr{{\eps}}{B_{\delta/2}(a)}(b,a)   \geq C \delta {\eps}^{-1},
\]

for some deterministic constant $C = C(\rho, \gamma)$ except on an event of probability going to zero as ${\eps} \to 0$.

We next prove the mean-value inequality using the basic properties of the discrete obstacle problem
and the convergence afforded by Theorem <ref>.

We restrict to the event that $\{t < \TT\}$ for $t > 0$.
As the function $\deg^{{\eps}}$ is reversible for simple random walk on ${\VGeps}$, by, for example, <cit.>,

\[
\Gr{{\eps}}{\Lambdateps}(a,b)/\deg^{{\eps}}(b) = 	\Gr{{\eps}}{\Lambdateps}(b,a)/\deg^{{\eps}}(a), \quad \forall a,b \in \Lambdateps.
\]

Therefore, we can rewrite the desired inequality as

\[
|\Lambdateps| \Gr{{\eps}}{\Lambdateps}(a, 0) \frac{\deg^{{\eps}}(a)}{\deg^{{\eps}}(0)} \geq (1 + \mathfrak{s}) \sum_{b \in \Lambdateps} \Gr{{\eps}}{\Lambdateps}(a, b)  \frac{\deg^{{\eps}}(a)}{\deg^{{\eps}}(b)} .
\]

By the definitions in Section <ref> and Remark <ref>, this shows that it suffices to establish

\begin{equation} \label{eq:mean-value-ineq-verify-changed}
|\Lambdateps| \gr{{\eps}}{\Lambdateps}(a, 0) \geq (1 + \mathfrak{s}) \qeps{\Lambdateps}(a) , \quad \forall
a \in \Lambdateps:\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta,
\end{equation}

for some finite $\sigma(h,\eta)$-measurable random variable $\mathfrak{s} = \mathfrak{s}(\rho, \gamma, \delta)> 0$ not depending on ${\eps}$,
except on an event of probability going to zero as ${\eps} \to 0$.

\begin{equation} \label{eq:odometer-tilde}
\begin{aligned}
v_t^{{\eps}}(a) &= \wteps(a) + (t {\eps}^{-1}) \gr{{\eps}}{B_{\rho}}(a, 0) \\
\tilde v_t^{{\eps}}(a) &= |\Lambdateps| \gr{{\eps}}{\Lambdateps}(a, 0) - \qeps{\Lambdateps}(a)
\end{aligned}
\end{equation}

and note that on $\Lambdateps$

\begin{align*}
&\Delta^{{\eps}} (v_t^{{\eps}} - \tilde v_t^{{\eps}}) \\
&= (\deg^{{\eps}})^{-1} + \deg^{{\eps}}(0)^{-1} t {\eps}^{-1} 1\{\cdot = 0\} - \Delta^{{\eps}} \tilde v_t^{{\eps}}
\quad \mbox{(by \eqref{eq:discrete-normalized-green-laplacian} and Lemma \ref{lemma:discrete-laplacian-bound})} \\
&= (\deg^{{\eps}})^{-1} + \deg^{{\eps}}(0)^{-1} t {\eps}^{-1} 1\{\cdot = 0\} -  (\deg^{{\eps}}(0)^{-1}  |\Lambdateps| 1\{\cdot = 0 \}
+ (\deg^{{\eps}})^{-1}) \\
&\quad \mbox{(by \eqref{eq:discrete-normalized-green-laplacian} and \eqref{eq:discrete-normalized-exit-laplacian})} \\
&\geq 0  \quad\mbox{($|\Lambdateps| \leq t {\eps}^{-1}$ by Lemma \ref{lemma:discrete-conservation-of-mass})}.
\end{align*}

Since $\wteps = \gr{{\eps}}{B_{\rho}}(\cdot, 0) = \qeps{\Lambdateps} = 0$ on $\partial \Lambdateps$,
this implies that

\[
\begin{cases}
\Delta^{{\eps}} (v_t^{{\eps}} - \tilde v_t^{{\eps}}) \geq 0 \quad \mbox{on $\Lambdateps$} \\
(v_t^{{\eps}} - \tilde v_t^{{\eps}}) = 0 \quad \mbox{on $\partial \Lambdateps$}.
\end{cases}
\]

By the maximum principle, Lemma <ref>, this implies that

\[
\tilde v_t^{{\eps}} \geq v_t^{{\eps}} \quad \mbox{on $\Lambdateps$}.
\]

By uniform convergence of $\overline v^{{\eps}}_t$,  (scaling defined in (<ref>))
which comes from Theorem <ref> and Lemma <ref>, for some finite
$\sigma(h,\eta)$-measurable random variable $C = C(\delta) > 0$, we have that

\[
\tilde v_t^{{\eps}}(a) \geq	v_t^{{\eps}}(a) > C {\eps}^{-1} , \quad \forall  a \in \Lambdateps:\dist(\eta(a), \partial \eta(\Lambdateps)) > \delta,
\]

except on an event of probability going to zero as ${\eps} \to 0$.
Also, by Lemma <ref>

\[
\sup_{a \in \Lambdateps}  \qeps{\Lambdateps}(a) \leq C_1 {\eps}^{-1},
\]

with probability approaching one as ${\eps} \to 0$ for some finite $\sigma(h,\eta)$-measurable random variable $C_1>0$ not depending on ${\eps}$. By combining the previous two intended equations with the definition (<ref>) of $\tilde v_t^{{\eps}}$,
we have (<ref>)
with $\mathfrak{s} := (C/C_1)$.

We conclude by including, for completeness, the arguments of <cit.>.

Recall that for each ${\eps} > 0$, the IDLA cluster is formed by running independent random walks
$X^{0, {\eps}}(i) := X^{{\eps}}(i)$ started at the origin on ${\Geps}$ and stopped upon exiting the cluster.
We let each independent random walk evolve forever even after it has left the occupied cluster
and introduce the stopping times:

\begin{align*}
\sigma^i &= \min \{ t: X^{{\eps}}_t(i) \not \in A_{{\eps}}(i -1) \} \\
&= \mbox{the time it takes the $i$th walker to leave the occupied cluster} \\
\tau_b^i &= \min \{ t: X^{{\eps}}_t(i) = b \}  \\
&= \mbox{the time it takes the $i$th walker to hit site $b$} \\
T^i &= \min\{ t : X^{{\eps}}_t(i) \not \in O_2^{{\eps}} \} \\
&= \mbox{the time it takes the $i$th walker to leave $O_2^{{\eps}}$}.
\end{align*}

For $a\in {\VGeps}$, let $E_a({\eps})$ denote the event that $a$ does not belong to the cluster $A_{{\eps}}(|O_2^{{\eps}}|)$.
We claim that to prove (<ref>), it suffices to show that
with probability approaching one as ${\eps} \to 0$,

\begin{equation} \label{eq:exponential-bound}
\P[E_a({\eps}) | h, \eta] \leq \exp(-c {\eps}^{-c}), \quad \forall a \in O_1^{{\eps}},
\end{equation}

for some deterministic $c > 0$ which is independent of ${\eps}$ and $a$.
Indeed, since the cells of ${\VGeps}$ have $\mu_h$-mass ${\eps}$ and intersect only along their boundaries there exists a random $C  = C(\rho, \gamma) > 0$ so that

\[
|{\VGeps}(B_{\rho})| \leq C {\eps}^{-1},
\]

with polynomially high probability as ${\eps} \to 0$.  Therefore, (<ref>) can be obtained from (<ref>) and a union bound over all $a \in O_1^{{\eps}}$.

Our aim now is to show (<ref>). Fix an ${\eps}$ for which (<ref>) holds and a $b \in O_2^{{\eps}}$.
Consider the random variables

\begin{align*}
M &= \sum_{i \leq |O_2^{{\eps}}|}  1\{ \tau_b^i \leq T^i \}  \\
&= \mbox{number of walks that visit $b$ before leaving $O_2^{{\eps}}$} \\
L &= \sum_{i \leq |O_2^{{\eps}}|} 1\{ \sigma^i \leq \tau_b^i < T^i  \} \\
&= \mbox{number of walks that visit $b$ before leaving $O_2^{{\eps}}$ but after exiting the current cluster}  .
\end{align*}

Observe that for all $k \geq 0$,

\begin{equation} \label{eq:upper-bound}
\begin{aligned}
\P[E_b({\eps}) | h , \eta]
&\leq \P[M = L | h, \eta] \\
&\leq \P[\{M \leq k \}\cup \{L \geq k\} | h , \eta]  \\
&\leq \P[M \leq k| h, \eta] + \P[L \geq k | h , \eta].
\end{aligned}
\end{equation}

We will choose $k$ below to minimize this bound.  First note that by, for example, <cit.>

\begin{equation} \label{eq:markov-chain-renewal}
\P^a[\tau_b < T | h ,\eta] =  \frac{\Gr{{\eps}}{O_2^{{\eps}}}(a, b)}{\Gr{{\eps}}{O_2^{{\eps}}}(b, b)}
\end{equation}

for all $a \in {\VGeps}$ where $\P^a[\tau_b < T | h ,\eta]$ denotes the probability
that an independent walker started at $a \in {\VGeps}$ hits site $b$ before it leaves  $O^{{\eps}}_2$.

As $M$ is a sum of i.i.d. random variables, by (<ref>)

\begin{equation} \label{eq:expected-value-of-m}
\E[M | h, \eta] = |O_2^{{\eps}}| \P[ \tau_b^{1} \leq T^{1} | h , \eta] = |O_2^{{\eps}}| \times \frac{\Gr{{\eps}}{O_2^{{\eps}}}(0, b)}{\Gr{{\eps}}{O_2^{{\eps}}}(b, b)}.
\end{equation}