Perturbed divisible sandpiles and quadrature surfaces

TL;DR

This paper links quadrature surfaces with sandpile dynamics by introducing a new lattice Laplacian growth model, analyzing its scaling limits, and revealing how mass distribution shapes evolve from spherical to polygonal forms.

Contribution

It introduces a novel lattice Laplacian growth model that connects quadrature surfaces with divisible sandpiles and analyzes its scaling limits and shape transitions.

Findings

Scaling limit is a ball when the threshold m is fixed.

Mass redistributes onto an annular ring of thickness 1/m.

When m tends to infinity slowly, the limit is a ball with mass on its boundary.

Abstract

The main purpose of the present paper is to establish a link between quadrature surfaces (potential theoretic concept) and sandpile dynamics (Laplacian growth models). For this aim, we introduce a new model of Laplacian growth on the lattice $\mathbb{Z}^d$ $(d\geq 2)$ which continuously deforms occupied regions of the \emph{divisible sandpile} model of Levine and Peres, by redistributing the total mass of the system onto $\frac 1m$-sub-level sets of the odometer which is a function counting total emissions of mass from lattice vertices. In free boundary terminology this goes in parallel with singular perturbation, which is known to converge to a Bernoulli type free boundary. We prove that models, generated from a single source, have a scaling limit, if the threshold $m$ is fixed. Moreover, this limit is a ball, and the entire mass of the system is being redistributed onto an annular ring of thickness $\frac 1m$. By compactness argument we show that, when $m$ tends to infinity sufficiently slowly with respect to the scale of the model, then in this case also there is scaling limit which is a ball, with the mass of the system being uniformly distributed onto the boundary of that ball, and hence we recover a quadrature surface in this case. Depending on the speed of decay of $m$, the visited set of the sandpile interpolates between spherical and polygonal shapes. Finding a precise characterisation of this shape-transition phenomenon seems to be a considerable challenge, which we cannot address at this moment.

Equations

$$\mathrm{Bal}:\mu ⟶{\mathcal{H}}^{d-1}{\lfloor }_{\partial D},$$

$$\int h(x)d\mu ={\int }_{\partial D}h(x)d{\mathcal{H}}^{d-1},$$

$$\Delta u_m=m{\mathbb{I}}_{\{0<u_m<1/m\}}-\mu ,$$

$$\Delta \tilde{u}=c_0{\mathcal{H}}^{d-1}{\lfloor }_{\partial {\mathbb{I}}_{\{\tilde{u}>0\}}}-\mu ,(some{c}_0>0)$$

$${\Delta }^{h}f(x)=\frac{1}{2d{h}^2}\sum _{y{\sim }_{h}x}[f(y)-f(x)],$$

$${\nabla }_e^{h}f(x)=h^{-1}\left[f(x+he)-f(x)\right],$$

$$g(x,y)=\{\begin{array}{ll}-\frac{2}{\pi }\log |x-y|+{\gamma }_0+\mathrm{O}(|x-y{|}^{-2}),& \text{d=2},\\ \frac{2}{(d-2)|B_1|}|x-y{|}^{2-d}+\mathrm{O}(|x-y{|}^{-d}),& \text{d≥3},\end{array}$$

$${\mu }_{k+1}(y)=\{\begin{array}{ll}{\mu }_k(x)-\mathcal{E}(x),& \text{if y=x},\\ {\mu }_k(y)+\frac{1}{2d}\mathcal{E}(x),& \text{if y∼x},\\ {\mu }_k(y),& \text{otherwise}.\end{array}$$

$${\Delta }^1{u}_k(x)={\mu }_k(x)-{\mu }_0(x),x\in {\mathbb{Z}}^d.$$

$$u_2(x_k)\ge u_{1,k}(x_k),k=1,2,....$$

$$u_{1,k}(x_k)=u_{1,k-1}(x_k).$$

$$u_2(x_k)=u_2(x_i)\ge u_{1,i}(x_i)=u_{1,k-1}(x_i)=u_{1,k}(x_k),$$

$$u_2(x)\ge u_{1,k}(x).$$

$$u_{1,k}(x)=u_{1,i}(x)\le u_2(x),$$

$${\mu }_2(x_k)\le {\mu }_{1,k}(x_k).$$

$$\begin{gathered} \frac{1}{2d}\sum _{y\sim x_k}{u}_2(y)-u_2(x_k)={\Delta }^1{u}_2(x_k)={\mu }_2(x_k)-{\mu }_0(x_k)\stackrel{(\text{ref}mu2-dist-vs-mu1k)}{\le } \\ {\mu }_{1,k}(x_k)-{\mu }_0(x_k)={\Delta }^1{u}_{1,k}(x_k)=\frac{1}{2d}\sum _{y\sim x_k}{u}_{1,k}(y)-u_{1,k}(x_k). \end{gathered}$$

$$u_2(x_k)-u_{1,k}(x_k)\ge \frac{1}{2d}\sum _{y\sim x_k}\left(u_2(y)-u_{1,k}(y)\right)\ge 0,$$

$$u_2(x_k)=u_2(x_i)\ge u_{1,i}(x_i)=u_{1,i}(x_k)>\frac{1}{m}{n}^{\frac{2}{d}}.$$

$$\mathcal{N}:=\{e_i,e_i\pm e_j:1\le i\ne j\le d\},$$

$$u(x)=u(x_{\ast }).$$

$$u_b(x_k)<{\mu }_k(x_k)-m+u_k(x_k),$$

$${\epsilon }_0:={\mu }_k(x_k)-m+u_k(x_k)-u_b(x_k)>0.$$

$${\mu }_{k+1}(x_k)-m={\mu }_k(x_k)-m+u_k(x_k)-u_b(x_k)={\epsilon }_0>0,$$

$$V=\mathrm{supp}{\mu }_0\cup \{u>0\}\cup \partial \{u>0\}.$$

$${\Delta }^{1}u(x)+{\mu }_0(x)\le m,x\in {\mathbb{Z}}^d,$$

$${\Delta }^{1}u(x)+{\mu }_0(x)=m{\mathbb{I}}_{\left\{0<u\le \frac{1}{m}{n}^{2/d}\right\}}(x),x\in \{u>0\},$$

$${\Delta }^{1}u(x)+{\mu }_0(x)\le {\Delta }^1{w}_1(x)+{\mu }_0(x)\le m,$$

$$0\le w_{\epsilon }(x)-u_{\ast }(x)\le \epsilon \text{ for all }x\in {\mathbb{Z}}^d.$$

$${\Delta }^{1}u(x)+{\mu }_0(x)=m{\mathbb{I}}_E(x)\text{ in }{V}_{\ast }\text{and}u=0\text{ on }\partial V_{\ast }.$$

$${\Delta }^1{v}_0(x)\le {\Delta }^1{u}_{\ast }(x)$$

Full text

We prove that models, generated from a single source, have a scaling limit, if the threshold $m$ is fixed. Moreover, this limit is a ball, and the entire mass of the system is being redistributed onto an annular ring of thickness $\frac{1}{m}$. By compactness argument we show that, when $m$ tends to infinity sufficiently slowly with respect to the scale of the model, then in this case also there is scaling limit which is a ball, with the mass of the system being uniformly distributed onto the boundary of that ball, and hence we recover a quadrature surface in this case.

Depending on the speed of decay of $m$, the visited set of the sandpile interpolates between spherical and polygonal shapes. Finding a precise characterisation of this shape-transition phenomenon seems to be a considerable challenge, which we cannot address at this moment.

H. A. was supported by postdoctoral fellowship from Knut and Alice Wallenberg Foundation. H. Sh. was partially supported by Swedish Research Council.

1. 1 Introduction

2. 1.1 Background

3. 1.2 Discrete Laplacian

4. 1.3 Our model and the main results

5. 2 Basic properties of the model

6. 2.1 Termination and abelian property

7. 2.2 The smallest super-solution

8. 3 Analysis of discrete shapes for single sources

9. 3.1 Discrete monotonicity

10. 3.2 The size of the sandpile shapes

11. 3.3 Uniform Lipschitz estimate

12. 3.4 $C^{1,1}$ estimates

13. 4 Scaling limits

14. A Uniqueness of the scaling limit

1. Introduction

1.1. Background

In a recent work [1], the current authors introduced a new growth model on the lattice ${\mathbb{Z}}^{d}$ $(d\geq 2)$ which redistributes a given initial mass on ${\mathbb{Z}}^{d}$ onto a combinatorial free boundary. The growth rule asks vertices of ${\mathbb{Z}}^{d}$ lying in the interior of the visited sites of the model, and vertices on the boundary of the set of visited sites carrying mass larger than a prescribed threshold, to redistribute their entire mass evenly among their $2d$ lattice neighbours. This procedure creates a sequence of non-decreasing domains that was shown to converge to, what we termed, Boundary Sandpile (hereinafter BS for short). This local rule allows (and in fact forces) huge masses to accumulate on the free boundary of the moving front. In the case of a single source mass, the numerics indicates that shapes generated by this model do not converge to a sphere under a scaling limit, but to a shape somewhat reminiscent of the classical Abelian sandpile (see [2] for definition, and [12] for the scaling limit).

The initial motivation for BS was to find a sandpile dynamic that gives us the so-called Quadrature surfaces (QS) (see [14, 15]). A QS, for a given source $\mu$ is the boundary of a domain $D$, that contains the support of $\mu$, and has the property that the Poincaré balayage satisfies

[TABLE]

which is equivalent to

[TABLE]

for all $h$ harmonic on a neighbourhood of $\overline{D}$.

Our intention in this paper is to introduce a new sandpile dynamic that (contrary to BS) will correspond to a QS in its scaling limit. This would then parallel the theory of (well-behaved) divisible sandpiles (DS), which redistributes an initial mass by putting a prescribed amount of mass at each visited site, and only moving out the excess from this prescribed amount. In DS the limit shapes, even for multi-source masses (with some reasonable control over their distribution), have shown to be the so-called quadrature domains (QD) with prescribed density and given source (see [8] for the single source, [9] for multi-source case, also [11] and [10]); a well-established area in potential theory [5].

To establish a link between QS and sandpile dynamics, we have chosen an approach based on singular perturbation theory for free boundary value problems (see [3]). This approach suggests that if one considers a slight perturbation of DS starting with total mass $n$ at the origin, by putting larger mass $m$ on $(\frac{1}{m})$-sub-level sets of the odometer function, and letting $m=m_{n}$ tend to infinity as $n$ tends to infinity, then interesting shapes appear. By taking $m\approx n^{1/d}$ we obtain a shape close to that of BS, and by letting $m$ be fixed (but large) we obtain spherical shapes, see Figure 2. The question that arose was:

How fast/slow (relative to $n$) should $m$ grow, in order to reach a desirable sandpile?

More precisely, we are interested in finding various functions $F$, with $m=F(n)$, for which there is a new sandpile shape. When $F(n)\asymp n^{1/d}$ it is obvious that we are close to BS, and when $F(n)=\mathrm{const}$, then we are having a DS222The only difference between DS and this case is that mass is now redistributed to the sub-level sets, with constant amount..

In this paper we show, using compactness arguments that there exists $m=F(n)\nearrow\infty$ for which the sandpile shape converge to a sphere, and the entire mass is being uniformly distributed onto the boundary. This in general suggests that in analogy with constructing QD through DS, we can construct QS through sandpile dynamics too.

Remark 1.1**.**

(Technical remark for experts)* In this paper, when studying the scaling limit of the model, we only consider the case of a single source, and leave out the general case. The reason for this is several infeasible technical difficulties at the moment. One major problem for the case of general initial source is the uniqueness question for the QS. This may still not be an issue, if we can show that scaling limits of our problem are unique (as it was done in the case of Abelian Sandpile [12], for instance). It is however far from obvious, even for a two-point source (see Figure 1 for a numerical illustration), how such a uniqueness can be proven. In our case we have used the geometry of spheres to overcome this difficulty.*

On the other hand it will be apparent (from what follows in this paper) that we can prove that for each fixed $m$ subsequences of our model converge to solutions of a free boundary problem of the type

[TABLE]

where $\mu$ is the initial source, consisting of several point masses. A solution to this problem is not necessarily unique if $\mu$ is not a single point source. However, from free boundary technique we know that under certain good conditions (which is satisfied at least for a two-source point) the function $u_{m}$ converge (for a subsequence) to a function $\tilde{u}$, as $m$ tends to infinity, where $\tilde{u}$ solves a free boundary problem of Bernoulli type

[TABLE]

which amounts to a quadrature identity of type (1.2).

From here we want to deduce that by compactness we may consider a sequence of $m=m_{n}$ such that $m_{n}\to\infty$, as $n\to\infty$, in a way that we can achieve the quadrature identity (1.2) for the limit scenario. For single source case the analysis heavily depends on the uniqueness and stability of the limit shape, that is unavailable for multi-source case.

Remark 1.2**.**

(To non-expert readers)* Although it took a while to pull the right strings to get this result, we want to stress that the paper actually is using soft mathematics, coming from classical theory of free boundary problems, translated into the discrete language. Hence we expect the reading to be easy. We need also to stress that the novelty of the paper lies in the model itself, and connections between free boundary problems and sandpile dynamics. This we believe might have some potential to grow in the near future and to give rise to an independent research area. There are indeed many free boundary problems on one side and several particle dynamic processes on the other side that seem to be defining the same model from different perspectives. Attempts to bridge these models should be valuable and interesting.*

1.2. Discrete Laplacian

Before getting into the technical definitions and details, we fix some notation, and several basic facts which will be used throughout the text.

Let $h>0$ be fixed, and $f:h{\mathbb{Z}}^{d}\to{\mathbb{R}}$ be a given function. Its (normalised) $\Delta^{h}$-Laplacian is defined as

[TABLE]

where $y\sim_{h}x$ means that $x,y\in h{\mathbb{Z}}^{d}$ are lattice neighbours, i.e. $||x-y||_{l^{1}}=h$. When $h=1$, we will simply write $x\sim y$ instead of $x\sim_{1}y$. Similar to the discrete Laplace operator, one defines discrete derivatives. Namely, for a unit vector $e\in{\mathbb{Z}}^{d}$, i.e. $||e||_{l^{1}}=1$, the function

[TABLE]

is the discrete derivative of $f$ in the direction $e$ at the point $x\in h{\mathbb{Z}}^{d}$.

It is well-known that as the lattice spacing tends to 0, the discrete Laplacian, and the discrete gradient converge to their continuous analogues, in a sense that for any $\varphi\in C_{0}^{\infty}({\mathbb{R}}^{d})$ one has $2d\Delta^{h}\varphi\to\Delta\varphi$, and $\nabla_{e}^{h}\varphi(x)\to\nabla_{e}\varphi(x)$ uniformly in ${\mathbb{R}}^{d}$, where $\Delta$ and $\nabla_{e}$ are respectively the continuum Laplace operator and derivative in the direction $e$.

We next recall the definition of the fundamental solution to $\Delta^{1}$ in ${\mathbb{Z}}^{d}$. Namely, there exists a function $g(x,y):{\mathbb{Z}}^{d}\times{\mathbb{Z}}^{d}\to{\mathbb{R}}$ such that $\Delta^{1}_{x}g(x,y)=-\delta_{0}(y-x)$ for all $x,y\in{\mathbb{Z}}^{d}$, where $\Delta^{1}_{x}$ is the Laplacian with respect to the $x$ variable, and $\delta_{0}$ is the characteristic function of $0\in{\mathbb{Z}}^{d}$. The existence and asymptotics of such $g$ are well-known (see [4], [16], [7]) and we have

[TABLE]

where $\gamma_{0}$ is a constant, and $|B_{1}|$ is the volume of the unit ball in ${\mathbb{R}}^{d}$.

The operator $\Delta^{1}$ also enjoys a maximum principle (see for instance [7, Exercise 1.4.7]).

**Discrete maximum principle (DMP): ** Let $u,v:{\mathbb{Z}}^{d}\to{\mathbb{R}}$, $V\subset{\mathbb{Z}}^{d}$ be finite, and $\Delta^{1}u\geq\Delta^{1}v$ in the interior of $V$. Then $\max_{\partial V}(u-v)\geq\max_{V}(u-v)$.

In this statement, and throughout the paper, the boundary of a non-empty set $E\subset{\mathbb{Z}}^{d}$, denoted by $\partial E$, is the set of all $x\in{\mathbb{Z}}^{d}$ for which there exists $x\in E$ such that $x\sim y$.

1.3. Our model and the main results

We now give the formal definition of the sandpile process which is the focus of this paper. Throughout the paper we will always assume that the dimension $d$ of the space is at least two. Assume we have an initial distribution of mass $\mu_{0}:{\mathbb{Z}}^{d}\to{\mathbb{R}}_{+}$ $(d\geq 2)$, that is a non-negative bounded function of finite support on ${\mathbb{Z}}^{d}$, and let $n\geq 0$ be the total mass of the system, i.e. $n=\sum_{x\in{\mathbb{Z}}^{d}}\mu_{0}(x)$. Fix also some threshold $m>0$, and inductively define a triple $(V_{k},u_{k},\mu_{k})$ for integer $k\geq 0$ as follows. When $k=0$ set $V_{0}=\mathrm{supp}\mu_{0}$ and $u_{0}\equiv 0$; we now define the passage from $k$ to $k+1$. At a time $k$ call a particular site $x\in{\mathbb{Z}}^{d}$ unstable if either of the following two (mutually exclusive) conditions hold:

• (a)

$\mu_{k}(x)>m$,

• (b)

$0<\mu_{k}(x)\leq m$ and $u_{k}(x)>\frac{1}{m}{n^{\frac{2}{d}}}$.

For an unstable site $x$ we define an excess of mass denoted by $\mathcal{E}(x)$ and being equal to $\mu_{k}(x)-m$ if $x$ is unstable according to (a), and to $\mu_{k}(x)$, if (b) holds for $x$. Then, a toppling of an unstable site $x$ is the procedure of redistribution of the excess $\mathcal{E}(x)$ evenly among its $2d$ lattice neighbours.

Thus, at each discrete time $k\geq 0$ we identify an unstable site, say $x_{k}$ and topple it. Then, we set $V_{k+1}:=V_{k}\cup\{y\in{\mathbb{Z}}^{d}:\ y\sim x\}$, $u_{k+1}(x)=u_{k}(x)+\delta_{x}\mathcal{E}(x)$, $x\in{\mathbb{Z}}^{d}$, and

[TABLE]

The process terminates if there are no unstable sites. It follows easily from the definition of the process, that for any $k\in{\mathbb{Z}}_{+}$ we have

[TABLE]

Let $T=\{x_{k}\}_{k=0}^{\infty}\subset{\mathbb{Z}}^{d}$ be any infinite toppling sequence, and $u_{k}$ and $\mu_{k}$ be the odometer, and mass distribution at time $k$ as defined above. Obviously, for each fixed $x\in{\mathbb{Z}}^{d}$ the sequence $\{u_{k}(x)\}$ is non-decreasing, and hence there is a limit (odometer) $u(x):=\lim_{k\to\infty}u_{k}(x)$, possibly infinite. Moreover, if $u$ is finite at a given $x$ and its lattice neighbours, then passing to the limit in (1.5) we see that $\mu_{k}(x)$ is also convergent; let $\mu(x)$ be its limit in such case. With this notation, we fix the following:

Definition 1.1**.**

(Stabilizing sequence)* Call a sequence $T$ stabilizing, if the limit odometer $u$ is finite everywhere on ${\mathbb{Z}}^{d}$, and the final configuration $(u,\mu)$ is stable.*

For convenience, we will allow a toppling sequence $T$ to include stable sites as well. Thus, if a certain site $x$ included in $T$ is assigned to topple at time $k$, but is stable, then toppling $x$ does not change anything. In other words the toppling procedure for stable vertices is idle, and the triple $(V_{k},u_{k},\mu_{k})$ simply remains the same for time $k+1$.

Main results and organisation of the paper. Our main results concern general analysis in the discrete space ${\mathbb{Z}}^{d}$ of the singularly perturbed sandpile model defined in subsection 1.3, and the existence of the scaling limit of the model generated from a single source. More precisely, in Section 2 we show that the model is well-posed, in a sense that it reaches a stable configuration given that all unstable sites have been toppled infinitely many times (Proposition 2.1), moreover, this stable configuration does not depend on the order of topplings (Proposition 2.2). We then show that the odometer function is the smallest element in the class of all super-solutions to discrete PDE problem associated with the sandpile (Lemma 2.6). Specifying the analysis for initial distributions concentrated at a single point of ${\mathbb{Z}}^{d}$, we show that the odometer function is monotone in the lattice directions (Theorem 3.1). We then estimate the size of the annulus where the mass is concentrated (Lemma 3.6), and of the region which is free of mass (Proposition 3.7) (cf. Figure 2). Finally, we complete Section 3 by proving that the odometer is Lipschitz regular uniformly with respect to threshold $m$ (Proposition 3.10) and is $C^{1,1}$ with norm depending linearly on $m$ (Proposition 3.11).

In Section 4 we show that the sandpile shapes generated by an initial distribution of the form $n\delta_{0}$ and fixed threshold $m$, converge to a ball under a scaling limit, as $n\to\infty$, moreover, the limiting odometer solves a singular perturbation problem with Bernoulli type free boundary condition (Theorem 4.1). We also prove, that there is a scaling limit for a subsequence of the sandpiles, as the mass threshold $m$ tends to infinity slowly with $n$ (Theorem 4.3).

The appendix to the paper, proves a uniqueness result for solutions to singular perturbation problem with Bernoulli type free boundary condition. This result is being used for establishing the existence of the scaling limit of the model.

2. Basic properties of the model

2.1. Termination and abelian property

The aim of this section is to prove the well-posedness of the model (in a sense to be made precise below), and to present some basic properties of it.

It is not hard to realize that the process defined in subsection 1.3 will not stabilise after finitely many steps (modulo the trivial cases). Moreover, vertices in the set of visited sites may become unstable infinitely many times during the lifetime of the process. This, as in the case of BS [1], brings us to he following class of the toppling sequences, which leaves a chance for the process to reach a stable state.

Definition 2.1**.**

(Infinitive sequence)* Call a sequence $T=\{x_{k}\}_{k=1}^{\infty}\subset{\mathbb{Z}}^{d}$ infinitive (on a set $E\subset{\mathbb{Z}}^{d}$) if any $x\in{\mathbb{Z}}^{d}$ (correspondingly $x\in E$) appears in $T$ infinitely often.*

If not stated otherwise, in this section $\mu_{0}$ is any initial distribution of mass, with total mass $n>0$, and $m>0$ stands for the threshold of the sandpile process. Note, that $n$ and $m$ are not necessarily integers.

The next two propositions show that any initial distribution can be stabilized uniquely by an infinitive toppling sequence. Proofs follow closely to their corresponding results from our paper [1], without any additional difficulties, but we include some details here for the sake of completeness.

Proposition 2.1**.**

(Termination of the process)* If a toppling sequence $T$ is infinitive, then it is stabilizing according to Definition 1.1.*

**Proof. **As in the proof of [1, Proposition 2.1] it will be enough to show that the set of visited sites $V_{k}$ stays bounded uniformly in $k$, for which it suffices to obtain a uniform bound on the number of elements in $\partial V_{k}$ (see the scanning procedure by hyperplanes as described in [1, Proposition 2.1]).

Let $k\geq 1$ be given, and fix any $x\in\partial V_{k}\setminus V_{0}$, clearly $u_{k}(x)=0$. Assume $y\sim x$ is the generator of $x$, i.e. the first time $x$ was visited by the process, was due to the toppling of $y$, say at time $i<k$. Since $x\sim y$ is not visited at time $i$, it follows that $y\in\partial V_{i}$, in particular $u_{i-1}(y)=0$, and hence $y$ can topple only according to condition (a) for instability as defined in subsection 1.3, and hence $y$ must carry mass at least $m$. This means that at time $k$, when $x\in\partial V_{k}$, the 1-neighbourhood of $x$ has total mass at least $m/(2d)$.

For each $x\in{\mathbb{Z}}^{d}$ let $\mathrm{B}(x)$ be a 1-box centred at $x$, i.e. $\mathrm{B}(x)=\{x\}\cup\{y\in{\mathbb{Z}}^{d}:\ y\sim x\}$. On one hand we proved that the total mass carried by vertices of $\mathrm{B}(x)$ where $x\in\partial V_{k}\setminus V_{0}$ is at least $m/(2d)$. On the other hand, each $x\in\partial V_{k}\setminus V_{0}$ can be counted at most $2d$ times in each box $\mathrm{B}$, and hence, we conclude that the total number of boundary points is bounded above by $(2d)^{2}\frac{n}{m}$ . Since the bound is uniform in $k$, the proof of this proposition is complete. ∎

Proposition 2.2**.**

(Abelian property)* For any two infinitive toppling sequences $T_{1}$ and $T_{2}$ the corresponding final configurations are the same.*

**Proof. **The scheme of the proof follows that of [1, Proposition 2.2]. In view of the previous proposition, each $T_{i}$ is stabilizing and has well-defined odometer function, call it $u_{i}$. We will prove that $u_{2}(x)\geq u_{1}(x)$ for any $x\in{\mathbb{Z}}^{d}$, which by symmetry will imply the desired result.

For $i=1,2$ by $u_{i,k}$ denote the odometer function and by $\mu_{i,k}$ the distribution corresponding to the sequence $T_{i}$ after the $k$-th toppling occurs. Let $T_{1}=\{x_{k}\}_{k=1}^{\infty}$, we now show that

[TABLE]

Observe that (2.1) is enough for our purpose since the site $x_{k}$ appears infinitely often in $T_{1}$ and $u_{1,k}$ converges to $u_{1}$ pointwise as $k$ tends to infinity. Thus, in what follows we prove (2.1) which we will do by induction on $k$.

The base case of induction, i.e. when $k=1$, is trivial. We now assume that (2.1) holds for any time $1\leq i<k$, and prove it for time $k$. First we eliminate the trivial case, when $x_{k}$ in $T_{1}$ is stable at time $k$. Indeed, in such scenario, we get

[TABLE]

So, if $x_{k}$ has never toppled prior to time $k$, then $u_{1,k-1}(x_{k})=0$ and we are done. Otherwise, if $i\leq k-1$ is the last time $x_{k}$ has toppled in $T_{1}$, then clearly $x_{i}=x_{k}$, and hence, using the inductive hypothesis, we conclude

[TABLE]

where the last equality follows from (2.2). This completes the induction step for the case when $x_{k}$ was already stable.

We next proceed to the case when toppling $x_{k}$ in $T_{1}$ is not stable at time $k$. First we will prove that for any $x\neq x_{k}$ one has

[TABLE]

There are two possible sub-cases, either $x$ was never toppled in $T_{1}$ up to time $k$, which implies $u_{1,k}(x)=0$ and we get (2.3), or $x$ was toppled at some time before $k$. In the latter case let $i\leq k$ be the last time $x$ has toppled prior to time $k$. Observe that $i<k$ since $x\neq x_{k}$. Then we have

[TABLE]

where the second inequality follows by inductive hypothesis. We thus have proved (2.3) for all $x\neq x_{k}$. Consider the following inequality

[TABLE]

To complete the induction step suppose for a moment that (2.4) holds true. Then

[TABLE]

Rearranging the first and the last terms leads to

[TABLE]

where the last inequality is due to (2.3). This completes the induction. Thereby, to finish the proof of the proposition, we need to verify (2.4) which we do next.

Recall that $x_{k}$ is unstable, and hence there are two possible reasons for $x_{k}$ to topple in $T_{1}$ at time $k$. Namely,

• (a)

$\mu_{k-1,1}(x_{k})>m$

• (b)

$0<\mu_{k-1,1}(x_{k})\leq m$ and $u_{1,k-1}(x_{k})>\frac{1}{m}n^{\frac{2}{d}}$.

If we are in case (a), then $\mu_{1,k}(x_{k})=m$ in view of the toppling rule, and hence (2.4) follows in view of the stability of $\mu_{2}$. Next, if $x_{k}$ is set to topple because of (b), then $u_{1,k-1}(x_{k})>\frac{1}{m}n^{\frac{2}{d}}$. But this means that $x_{k}$ had already toppled prior to time $k$, and we let $i<k$ be the largest time before $k$, when it has toppled; in particular $x_{k}=x_{i}$. Applying the inductive hypothesis we arrive at

[TABLE]

From here, and the stability of $\mu_{2}$ we obtain $\mu_{2}(x_{k})=0$, completing the proof of 2.4), and hence the proposition is proved.

∎

A simple consequence of the abelian property, is the following symmetry for a point mass concentrated at the origin. Consider the set of vectors

[TABLE]

where $e_{i}\in{\mathbb{Z}}^{d}$ is the $i$-th element of the standard basis. Let $T$ be a hyperplane through the origin of ${\mathbb{Z}}^{d}$ and with a normal collinear with some element of $\mathcal{N}$. Then, any sandpile generated by initial distribution of the form $n\delta_{0}$, is symmetric with respect to $T$. Namely, if $u$ is the odometer, and $x\in{\mathbb{Z}}^{d}$ is arbitrary, for $x^{*}\ -$the mirror reflection of $x$ with respect to $T$ (which is obviously in ${\mathbb{Z}}^{d}$ in view of the choice of $T$) one has

[TABLE]

For coordinate directions, i.e. when $T$ has normal in the direction of some $e_{i}$, the proof follows by symmetrization of the toppling sequence (see [1, Corollary 2.4]). In case of the directions $e_{i}\pm e_{j}$, the claim follows by noticing that the symmetry with respect to $T$ is a composition of two reflections with respect to coordinate axes.

The next result, which is yet another easy corollary of the abelian property, will be used in the proof of Lemma 2.6.

Lemma 2.3**.**

Assume we are given a function $u_{b}:{\mathbb{Z}}^{d}\to{\mathbb{R}}_{+}$ and an infinitive toppling sequence $T=\{x_{k}\}_{k=0}^{\infty}$. Suppose in addition, that at any time $k\geq 0$ we require the total (cumulative) emissions of mass from $x_{k}$ to NOT exceed $u_{b}(x_{k})$. Under this restriction, if the limiting configuration is stable, then it coincides with the limiting (stable) configuration of the sandpile.

**Proof. **For each $k\in{\mathbb{Z}}_{+}$ let $u_{k}$ be the odometer function under the additional restriction imposed by the barrier $u_{b}$, and let $\mu_{k}$ be the mass distribution at time $k$. We will show, that for each $k$ toppling $x_{k}$ makes the site $x_{k}$ stable, meaning that the requirement $u_{k}(x)\leq u_{b}(x)$ for all $k\in{\mathbb{Z}}_{+}$ and all $x\in{\mathbb{Z}}^{d}$, can be dropped. Then, the claim of the lemma will follow by the abelian property proved in Proposition 2.2.

Fix $k\geq 0$, then if $x_{k}$ is stable, no toppling will be performed at time $k$, otherwise, there are two mutually exclusive reasons for instability of $x_{k}$ at time $k$. Namely,

• (a)

$\mu_{k}(x_{k})>m$

• (b)

$0<\mu_{k}(x_{k})\leq m$ and $u_{k}(x_{k})>\frac{1}{m}n^{\frac{2}{d}}$.

Assume we are in case (a). Then, according to the sandpile rule defined in subsection 1.3 we need to topple the excess from $m$, i.e. $\mu_{k}(x_{k})-m$. Following the restriction imposed by the barrier $u_{b}$, the maximum amount of mass we can move out from $x_{k}$ equals $u_{b}(x_{k})-u_{k}(x_{k})$. Now assume that

[TABLE]

i.e. we cannot move the entire excess. Set

[TABLE]

Then, invoking the $k$-th toppling, we will move $u_{b}(x_{k})-u_{k}(x_{k})$, meaning that $u_{k+1}(x_{k})=u_{b}(x_{k})$. Hence,

[TABLE]

and since no further toppling of $x_{k}$ will be allowed, the last expression shows that that $x_{k}$ remains unstable in the limiting configuration, which is a contradiction. We thus obtain, that in case (a) the usual toppling rule is not affected by the existence of a barrier $u_{b}$.

The same reasoning gives the claim in case (b) too, and hence the proof of the lemma is completed. ∎

Having a unique stable configuration for the sandpile, for fixed initial distribution $\mu_{0}$ we set $u:{\mathbb{Z}}^{d}\to{\mathbb{R}}_{+}$, $V\subset{\mathbb{Z}}^{d}$, and $\mu:{\mathbb{Z}}^{d}\to{\mathbb{R}}_{+}$ for its odometer, set of visited sites, and the final distribution respectively. In particular we have

[TABLE]

Let us also see what is the discrete PDE problem solved by odometer. If $x\in{\mathbb{Z}}^{d}$ and $0<u(x)\leq\frac{1}{m}n^{2/d}$, then $x$ has always toppled according to rule (a), and hence $\mu(x)=m$. Next, if $u(x)>\frac{1}{m}n^{2/d}$, then eventually $x$ starts to topple according to (b), and hence we get $\mu(x)=0$. Finally, if $x\in\partial\{u>0\}$ then $u(x)=0$, and hence $x$ has never toppled, which implies $\mu\leq m$. Summarising the discussion, we conclude

[TABLE]

[TABLE]

where $\mathbb{I}$ stands for the characteristic function of a set. Observe, that if the initial distribution $\mu_{0}$ is already stable, then the odometer is identically zero.

2.2. The smallest super-solution

We define a class of super-solutions to our sandpile model.

Definition 2.2**.**

(Super-solutions) Let $\mu_{0}:{\mathbb{Z}}^{d}\to{\mathbb{R}}_{+}$ be a given initial distribution with total mass $n$, and let $m>0$ be fixed. Call $u:{\mathbb{Z}}^{d}\to{\mathbb{R}}_{+}$ a super-solution to the sandpile if

• (i)

$\Delta^{1}u(x)+\mu_{0}(x)\leq m\qquad x\in{\mathbb{Z}}^{d}$,

• (ii)

$\Delta^{1}u(x)+\mu_{0}(x)\leq m\mathbb{I}_{\left\{0<u\leq\frac{1}{m}n^{2/d}\right\}}(x)\qquad x\in\{u>0\}$,

Denote by $\mathbb{W}$ the set of super-solutions, which is non-empty since the odometer of the sandpile is in $\mathbb{W}$. We have the following.

Lemma 2.4**.**

If $u_{*}(x)=\inf\limits_{w\in\mathbb{W}}w(x)$, $x\in{\mathbb{Z}}^{d}$, then $u_{*}\in\mathbb{W}$.

**Proof. **Let us first show that for any $w_{1},w_{2}\in\mathbb{W}$ one has $u:=\min\{w_{1},w_{2}\}\in\mathbb{W}$. Fix $x\in{\mathbb{Z}}^{d}$ and assume $0<u(x)\leq\frac{1}{m}n^{2/d}$. Then, $w_{i}(x)>0$ for $i=1,2$, and if $u(x)=w_{1}(x)$ then $w_{1}(x)\leq\frac{1}{m}n^{2/d}$. From this we get

[TABLE]

where the second inequality is in view of $w_{1}\in\mathbb{W}$. Similarly, if $u(x)>\frac{1}{m}n^{2/d}$, then both $w_{1}(x),w_{2}(x)>\frac{1}{m}n^{2/d}$, and arguing as above we obtain $\Delta^{1}u(x)+\mu_{0}(x)\leq 0$. We thus have condition (ii) for $u$, and (i) follows similarly, in a straightforward manner. In particular, we get that $\mathbb{W}$ is closed under taking minimum of finitely many elements.

We next show that $u_{*}\in\mathbb{W}$. Let $u_{0}$ be the odometer of the sandpile. Since $u_{0}\in\mathbb{W}$, and hence $0\leq u_{*}\leq u_{0}$, we get that $u_{*}$ has finite support, and we let $E=\{x_{1},...,x_{N}\}$ be the closure of the support of $u_{*}$, i.e. the support of $u_{*}$ along with its lattice boundary. Then, for any $\varepsilon>0$ and each $1\leq i\leq N$ there exists $w_{\varepsilon}^{i}\in\mathbb{W}$ such that $0\leq w_{\varepsilon}^{i}(x_{i})-u_{*}(x_{i})\leq\varepsilon$. Setting $w_{\varepsilon}:=\min\{w_{\varepsilon}^{1},...,w_{\varepsilon}^{N},u_{0}\}$, we obtain $w_{\varepsilon}\in\mathbb{W}$ and

[TABLE]

Since $\varepsilon>0$ is arbitrary, the claim of the lemma follows from (2.8) by a similar argument as we had above. ∎

Lemma 2.5**.**

Let $u_{*}$ be the smallest super-solution defined in Lemma 2.4. Then, $u_{*}$ satisfies Definition 2.2 with equality in (ii).

**Proof. **Let $V_{*}=\{u_{*}>0\}$ which is a finite set in ${\mathbb{Z}}^{d}$. If $V_{*}$ is empty we are done, otherwise, set $E=\{x:\ 0<u_{*}(x)\leq\frac{1}{m}n^{2/d}\}\subset V_{*}$, and let $u$ be the solution to

[TABLE]

It follows from the definition of $u$ and $u_{*}$ that $u-u_{*}$ is a sub-solution, hence from DMP we get that $u\leq u_{*}$ in $V_{*}$. Since both $u$ and $u_{*}$ are 0 outside $V_{*}$ we obtain $u\leq u_{*}$ on ${\mathbb{Z}}^{d}$. Let us now prove that $u$ is a super-solution; this will imply that $u=u_{*}$ and then (2.9) will complete the proof.

We start with part (i) of Definition 2.2. If $x\in V_{*}$ then inequality of (i) follows from (2.9), otherwise, if $x\in{\mathbb{Z}}^{d}\setminus V_{*}$ we have $u(x)=u_{*}(x)=0$ which, combined with (i) for $u_{*}$, and the inequality $u\leq u_{*}$ on ${\mathbb{Z}}^{d}$, implies condition (i) of a super-solution for $u$. We next proceed to (ii). It is clear from (2.9) that it suffices to show that $\Delta^{1}u(x)+\mu_{0}(x)\leq 0$ if $u(x)>\frac{1}{m}n^{2/d}$. But since $u_{*}\geq u$ everywhere, for such $x$ we have $u_{*}(x)>\frac{1}{m}n^{2/d}$ and hence $x\notin E$. The latter and (2.9) imply $\Delta^{1}u(x)+\mu_{0}(x)=0$, and the claim of the lemma follows. ∎

Lemma 2.6**.**

Let $u_{*}$ be as in Lemma 2.4, and let $u_{0}$ be the odometer of the sandpile. Then $u_{*}=u_{0}$.

**Proof. **If $u_{*}=0$, then the sandpile is stable, and we are done, since then $u_{0}=0$ as well. Now assume that $V_{*}:=\{u_{*}>0\}$ is not empty. Fix any infinitive toppling sequence on $V_{*}$, call it $T=\{x_{k}\}_{k=0}^{\infty}$, and do the toppling following $T$, with an additional requirement that at any time $k$ the cumulative emissions of mass from any $x\in{\mathbb{Z}}^{d}$ cannot exceed $u_{*}(x)$. Let the limiting odometer for this modified process be $v_{0}$. We claim that $v_{0}=u_{*}$. Observe, that this equality, combined with Lemma 2.3 and stability of $u_{*}$ completes the proof of the current lemma. In what follows we prove the desired equality.

It is enough to show that $v_{0}$ is a super-solution, since then the claim follows by minimality of $u_{*}$, and inequality $v_{0}\leq u_{*}$ which is due to definition of $v_{0}$. Let $x\in V_{*}$ be a point of touch, i.e. $u_{*}(x)=v_{0}(x)$; we claim that at these points $v_{0}$ is stable. Indeed, we have

[TABLE]

where we have used that $u_{*}(x)=v_{0}(x)$ and hence the characteristic functions coincide at $x$.

Now assume that $x$ is not a point of touch, meaning that $v_{0}(x)<u_{*}(x)$, and assume further that the sandpile configuration $\mu$ corresponding to $v_{0}$ is not stable at $x$. This means that either $\mu(x)>m$, or we have both $\mu(x)>0$ and $v_{0}(x)>\frac{1}{m}n^{2/d}$. Since $v_{0}(x)<u_{*}(x)$, we see that in both cases we could have toppled more mass from $x$, which means that in the toppling process we had not moved the maximal allowed excess mass, contradicting the toppling rule. We conclude that there can be no $x\in{\mathbb{Z}}^{d}$ which is not stable and at the same time $v_{0}(x)<u_{*}(x)$, hence the configuration corresponding to $v_{0}$ is stable, and the proof of the lemma is complete. ∎

3. Analysis of discrete shapes for single sources

In this section we analyse the shape of the sandpile in the discrete space ${\mathbb{Z}}^{d}$ generated by initial distribution $n\delta_{0}$ and threshold parameter $m$. Namely, we establish a certain discrete monotonicity property of the odometer function, estimate the size of regions in the set of visited sites with and without mass, and conclude the section by proving discrete Lipschitz (uniform with respect to the threshold $m$), and $C^{1,1}$ regularity estimates on the odometer function.

3.1. Discrete monotonicity

As in the case of BS [1], here as well the odometer function of the singularly perturbed sandpile started with initial mass concentrated at a single point of ${\mathbb{Z}}^{d}$, enjoys monotonicity in lattice directions. More precisely, denote by $\mathcal{S}$ the set of mirror symmetry hyperplanes of the unit cube $[0,1]^{d}$; clearly $\mathcal{S}$ consists of the hyperplanes $\{x_{i}=1/2\}$, $\{x_{i}=x_{j}\}$ and $\{x_{i}=-x_{j}\}$ where $i,j=1,...,d$ and $i\neq j$. In particular, there are $d^{2}$ elements in $\mathcal{S}$.

The next theorem is the analogue of what we had for BS (the result also holds true for the classical Abelian sandpile as is proved in [1]). The proof of the current version, is a simplification of the one we had in [1]. We will follow the same line of arguments as we had in [1, Theorem 4.5], and we will only outline the differences for this case, which are minor.

Theorem 3.1**.**

(Directional monotonicity)* For initial mass distribution $\mu_{0}=n\delta_{0}$ let $V\subset{\mathbb{Z}}^{d}$ be the set of visited sites, and let $u$ be the odometer function. Fix any hyperplane $T\in\mathcal{S}$. Then, for any $X_{1},X_{2}\in{\mathbb{Z}}^{d}$, such that $X_{1}-X_{2}$ is a non-zero vector orthogonal to $T$, we have*

[TABLE]

**Proof. **Translate $T$ to a position where it has equal distance from $X_{1}$ and $X_{2}$. Let $T_{0}$ be this unique translated copy of $T$. For a given $x\in{\mathbb{R}}^{d}$ denote by $x^{*}$ its mirror reflection with respect to the hyperplane $T_{0}$. Due to the choice of $T$ and $X_{1}$, $X_{2}$ it follows that this reflection preserves the lattice, and we we have $({\mathbb{Z}}^{d})^{*}={\mathbb{Z}}^{d}$, in particular $X_{1}=(X_{2})^{*}$.

Set $V^{*}:=\{x^{*}:\ x\in V\}$ for the reflection of $V$, and similarly define the reflected odometer function by $u^{*}(x):=u(x^{*})$, where $x\in{\mathbb{Z}}^{d}$. From the discussion above, we have that $u^{*}$ is defined on ${\mathbb{Z}}^{d}$ and $V^{*}\subset{\mathbb{Z}}^{d}$. Let $\mathcal{H}_{-}$ be the closed halfspace of ${\mathbb{R}}^{d}$ determined by $T_{0}$ and containing the origin, and let $\mathcal{H}_{+}={\mathbb{R}}^{d}\setminus\mathcal{H}_{-}$. Consider the set $V_{T}=V_{-}\cup V_{+}$ where $V_{-}:=\mathcal{H}_{-}\cap V$ and $V_{+}:=\mathcal{H}_{+}\cap V\cap V^{*}$. Our first goal is to show that $V_{T}=V$. By definition $V_{T}\subset V$, so we need to establish the reverse inclusion. Consider the function

[TABLE]

We claim that $u_{T}$ defines a super-solution in a sense of Definition 2.2. To see this, consider two cases.

Case 1. $x\in V_{-}$. As in [1, Theorem 4.5], here as well we get $\Delta^{1}u_{T}(x)\leq\Delta^{1}u(x)$, which, together with the equality $u_{T}(x)=u(x)$ gives

[TABLE]

Case 2. $x\in V_{+}$. If $u_{T}(x)=u(x)$ then we proceed as in Case 1, and get the stability at $x$. We will thus assume that $u_{T}(x)=u(x^{*})$. In this case, following again [1, Theorem 4.5], we obtain $\Delta^{1}u_{T}(x)\leq\Delta^{1}u(x^{*})$. Observe that if $u_{T}(x)\leq\frac{1}{m}n^{2/d}$ then the last inequality gives the stability of the sandpile at $x$. Otherwise, if $u_{T}(x)>\frac{1}{m}n^{2/d}$ we get $u(x_{*})>\frac{1}{m}n^{2/d}$ and hence

[TABLE]

completing the proof that $u_{T}$ satisfies condition (ii) of being a super-solution.

Part (i) follows from the fact that the support of $u_{T}$ is included in the support of $u$. We conclude that $u_{T}$ is a super-solution, and hence due to the minimality of Lemma 2.5, we get $u_{T}=u$ and in particular,

[TABLE]

completing the proof of the theorem. ∎

3.2. The size of the sandpile shapes

There are two distinct regions in the visited set of the sandpile process. Namely, the one which is free of mass, and the one where each vertex carries mass equal to the threshold $m$. The aim of this subsection, is to give estimates on the size of these two regions.

Recall, that we are dealing with sandpiles with initial distribution $n\delta_{0}$ and mass threshold equal to $m$. We will be using the following notation (cf. Figure 2):

• •

$V_{n,m}\subset{\mathbb{Z}}^{d}\ -$ the set of visited sites of the sandpile,

• •

$V_{n,0}\subset V_{n,m}\ -$ the set where the odometer $u>\frac{1}{m}n^{2/d}$; in particular $u$ is harmonic on $V_{n,0}\setminus\{0\}$ and (thanks to the discrete monotonicity of $u$) the origin lies in $V_{n,0}$,

• •

$V_{n,1}:=V_{n,m}\setminus V_{n,0}\ -$ the set where the entire mass of the sandpile is concentrated.

Lemma 3.2**.**

(Upper bound on the size of mass region)* Let $x_{0}\in V_{n,1}$ and $r>0$ be such that $B(x_{0},r)\cap V_{n,0}=\emptyset$. Then, there exists a dimension dependent constant $C_{d}>0$ such that for any $n>0$ large and any $1<m<n$ one has*

[TABLE]

**Proof. **Consider the function $w(x)=u(x)-u(x_{0})-m|x-x_{0}|^{2}$, for which we have $\Delta^{1}w(x)=0$ for all $x\in B(x_{0},r)\cap V_{n,1}=:D$, and hence by DMP we get

[TABLE]

Observe that

[TABLE]

where we have used the fact that the ball $B$ has no intersection with the region $V_{n,0}$. Since $u=0$ on $\partial V_{n,m}$ we have $w<0$ on $\Gamma_{2}$ which, together with (3.1) implies

[TABLE]

where $c_{d}>0$ is some small constant depending on dimension $d$ only. Rearranging the last inequality completes the proof of the lemma. ∎

Lemma 3.3**.**

Fix $r>0$ large, and a function $f:B\to{\mathbb{R}}_{+}$ bounded above by $m$, where $B=B(0,r)\cap{\mathbb{Z}}^{d}$. Let $u$ be the (unique) solution to

[TABLE]

Then, we have

• (i)

$|u(x)-u(y)|\leq Crm$, for any $x\in B_{r/2}=B(0,r/2)\cap{\mathbb{Z}}^{d}$ and $y\sim x$,

• (ii)

$-r^{2}m\leq u(0)\leq 0$,

with dimension dependent constant $C$.

**Proof. **We start with (i). For $x\in B$ we have the representation

[TABLE]

where $G$ is the Green’s kernel for $B$, i.e. the expected number of passages through $w$ of a random walk started at $x$ before escaping from $B$. More precisely, we have

[TABLE]

where $T$ is the first exit time from $B$ of the walk, $g$ is the fundamental solution of $\Delta^{1}$ defined in subsection 1.2, and $\mathbb{E}^{x}$ stands for the expectation conditioned that the walk has started from $x$.

Now fix any $x,y\in B_{r/2}$ such that $x\sim y$. Then, from (3.2) we have

[TABLE]

From (3.3) we get

[TABLE]

where we have used the asymptotics (1.3) and that $x\sim y$ to estimate the difference of $g$-s, and the fact that $x\in B_{r/2}$ in addition to those to bound the sum. With this estimate getting back to sum in (3.4) we obtain

[TABLE]

where the sum of the first term is estimated by a simple counting argument relying on the structure of $B$ (see, e.g. [8, Lemma 5.2]), while bound on the sum involving the second summand follows from a trivial estimate $|B|\leq C_{d}r^{d}$. Returning to (3.4) we are left to estimate only the first two sums on the r.h.s., but their contribution is bounded above by $C_{d}m$ in view of the representation (3.3), the fact that $x\sim y$, and then using that $g$ is symmetric with respect to coordinate axes and has bounded Laplacian.

We next proceed to the claim of (ii). The upper bound is a consequence of DMP; to establish the lower bound, consider the function

[TABLE]

Clearly, $\Delta^{1}v=-m$ everywhere and $v\leq 0$ on $\partial B$. From here, we get

[TABLE]

and hence

[TABLE]

We thus have

[TABLE]

The proof of the lemma is now complete. ∎

A trivial corollary of the 1-step Lipschitz estimate of Lemma 3.3 is the following.

Corollary 3.4**.**

Retain all notation and conditions of Lemma 3.3. Then, for any $x,y\in B_{r/2}$ one has

• (i)

$|u(x)-u(y)|\leq Crm|x-y|$* for any $x,y\in B_{r/2}$,*

• (ii)

$-Cr^{2}m\leq u(x)\leq 0$* for any $x\in B_{r/2}$,*

where $C>0$ is a constant depending on dimension only.

**Proof. **To see (i), take any path through $B_{r/2}$ connecting $x$ and $y$, namely

[TABLE]

where $x_{i}\in B_{r/2}$ for $0\leq i\leq k$. Clearly, we can assume $k\asymp|x-y|$, by considering a path of the shortest length. Now, (i) follows from

[TABLE]

where we have applied Lemma 3.3 to each summand.

For (ii) observe that the upper bound is again due to DMP, and the lower bound is in view of

[TABLE]

where we have used Lemma 3.3 (ii). The proof is now complete. ∎

The next result will be used in estimating the size of the mass region $V_{n,1}$.

Lemma 3.5**.**

(Quadratic upper bound)* Let $x_{0}\in\partial V_{n,m}$ be any. Take $r>0$ large so that the ball $B_{r}=B(x_{0},r)\cap{\mathbb{Z}}^{d}$ does not contain the origin of ${\mathbb{Z}}^{d}$. Then*

[TABLE]

**Proof. **We start with the splitting $u=u_{1}+u_{2}$ where

[TABLE]

and

[TABLE]

Since the origin does not lie in $B_{r}$, we have $0\leq\Delta^{1}u\leq m$ in $B_{r}$, and hence applying Lemma 3.3 (i) (to the function $u_{1}(x+x_{0})$) we get $|u_{1}(x_{0})|\leq Cmr^{2}$, consequently

[TABLE]

From here and using that fact that $u_{2}$ is harmonic in $B_{r}$ and non-negative on the boundary of $B_{r}$, we get, by discrete Harnack (see [7, Theorem 1.7.2]), that

[TABLE]

Since $u_{1}\leq 0$ in $B_{r}$ by DMP, the last inequality implies the claim of the current lemma. ∎

Recall that in Lemma 3.2 we established an upper bound on the thickness of the mass-region. We are now in a position to also give a bound from below.

Lemma 3.6**.**

(The size of mass-region)* Let $x_{0}\in\partial V_{n,m}$ and let $r>0$ be such that $B(x_{0},r)\subset V_{n,1}$, and $B(x_{0},r+1)\cap V_{n,0}\neq\emptyset$. Then*

[TABLE]

where equivalence holds with dimension dependent constants.

**Proof. **The upper bound on $r$ is due to Lemma 3.2, and we only need to prove the lower bound here. By assumption there exists $z\in B_{r+1}\cap{\mathbb{Z}}^{d}$ such that $u(z)\geq\frac{n^{2/d}}{m}$. Now, applying Lemma 3.5 we obtain

[TABLE]

which completes the proof of the lemma. ∎

Proposition 3.7**.**

(The size of mass-free region)* There exist dimension dependent positive constants $c_{1}<c_{2}$ such that for any fixed threshold $m>1$ and $n$ large enough one has*

[TABLE]

**Proof. **We first compare $V_{n,0}$ with sets having a very simple structure. Let $R>0$ be the largest integer such that the point $X_{R}=(0,...,0,R)\in{\mathbb{Z}}^{d}$ is inside $V_{n,0}$ but is not an interior point of the set $V_{n,0}$. In particular, $u(X_{R})>\frac{1}{m}n^{2/d}$ but for some $y\sim X_{R}$ we have $u(y)\leq\frac{1}{m}n^{2/d}$. Let $\mathcal{S}_{R}$ be the simplex with vertices at $\pm Re_{i}$, where $i=1,2,...,d$, in other words $\mathcal{S}_{R}$ is the ball of radius $R$ in $l_{1}$ metric. Let us prove the following inclusions:

[TABLE]

Notice, that both the simplex and the cube in (3.7) are restricted to ${\mathbb{Z}}^{d}$.

We take any $x\in\mathcal{S}_{R}$ and show that it is in $V_{n,0}$. Following (2.5), we know that the odometer function is symmetric with respect to coordinate axes, and hyperplanes through the origin with normals in the directions $e_{i}\pm e_{j}$, and obviously so is the simplex $\mathcal{S}_{R}$, and hence, it will be enough to prove (3.7) for $X_{0}=(x_{1},...,x_{d})\in{\mathbb{Z}}^{d}$ satisfying $x_{d}\geq|x_{i}|$, $i=1,...,d-1$. To accomplish that, we will use the directional monotonicity of the odometer.

For $1\leq i\leq d-1$, set $\nu_{i}=e_{d}-e_{i}$, $\nu_{i+d-1}=e_{d}+e_{i}$, and denote $\nu_{2d-1}=e_{d}$. For each $1\leq i\leq 2d-1$ consider the discrete halfspace

[TABLE]

and define $H_{*}:=H_{1}\cap...\cap H_{2d-1}$. In view of the choice of the vectors $\nu_{i}$ we have

[TABLE]

in particular $X_{0},X_{R}\in H_{*}$. Next, consider the cone

[TABLE]

Since $\nu_{i}\cdot\nu_{j}\geq 0$ for all $1\leq i,j\leq 2d-1$, we get $\mathcal{C}_{0}\subset H_{*}$. We now translate $\mathcal{C}_{0}$ to $X_{0}$ by setting $\mathcal{C}=X_{0}+\mathcal{C}_{0}$, and as $X_{0}\in H_{*}$ we get $\mathcal{C}\subset H_{*}$. In view of the choice of the collection $\{\nu_{i}\}_{i=1}^{2d-1}$ and the points $X_{0},X_{R}$, it is easy to see $X_{R}\in\mathcal{C}$. Finally, in the cone $\mathcal{C}$ we use the directional monotonicity given by Theorem 3.1 which implies that in $\mathcal{C}$ the odometer $u$ attains its maximum at the vertex of the cone, i.e. at $X_{0}$. Since $X_{R}\in\mathcal{C}$ we obtain $u(X_{0})\geq u(X_{R})$, and hence $X_{0}\in V_{n,0}$ which completes the proof of the first inclusion of (3.7). The second inclusion, is much simpler, and follows easily by using monotonicity in the coordinate directions only.

With (3.7) at hand, the proof of the proposition will be complete, once we show that $R\asymp n^{1/d}$ with constants depending only on $d$. In what follows we prove this, and hence the proposition. We will use the fact that $\partial V_{n,0}$ is locally a graph, along with bounds of the size of the mass-region. Here again, due to the symmetry, it will be enough to consider the region where $x_{d}\geq|x_{i}|$ for all $1\leq i\leq d-1$.

Let $\Pi\subset{\mathbb{Z}}^{d-1}\times\{0\}$ be the projection of the set $H_{*}\cap\mathcal{S}_{R}$ onto ${\mathbb{Z}}^{d-1}\times\{0\}$. In view of the monotonicity of $u$ in the direction of $e_{d}$, we have $\Pi\subset V_{n,0}$. Moreover, definitions of $H_{*}$ and $\mathcal{S}_{R}$ imply that $|\Pi|\asymp R^{d-1}$ with constants in the equivalence depending on dimension $d$ only. Now, for a given $X=(\overline{x},0)\in\Pi$ where $\overline{x}\in{\mathbb{Z}}^{d-1}$, let $t_{1}\in{\mathbb{Z}}_{+}$ be the smallest integer such that $X_{1}=(\overline{x},t_{1})\notin V_{n,0}$ and let $t_{2}\geq t_{1}$ be the smallest integer such that $u$ vanishes at $X_{2}=(\overline{x},t_{2})$. In view of the choice of $X_{1}$ and the discussion above, the cone $X_{1}+\mathcal{C}_{0}$ has no points of $V_{n,0}$. Moreover, $X_{2}\in\mathcal{C}_{0}$ and according to Lemma 3.6 we obtain that $t_{2}-t_{1}\asymp r$ where $r=\frac{1}{m}n^{1/d}$ and constants in the equivalence depend on dimension only. We thus see that each line in the direction of $e_{d}$ through points of $\Pi$ intersects the mass-region by $\asymp r$ points, and since each point in $V_{n,1}$ carries mass $m$ and the total mass of the system is preserved, we get

[TABLE]

from this, and the estimate of Lemma 3.6, we get $R\asymp n^{1/d}$ completing the proof of this proposition. ∎

Remark 3.8**.**

Observe, that putting together Lemma 3.6 and Proposition 3.7 we see that for large $n>1$ and $1<m<n^{1/d}$, the set of visited sites of the sandpile grows proportionally to $n^{1/d}$ uniformly in $m$.

Remark 3.9**.**

The proof of Proposition 3.7 shows, that in the cone

[TABLE]

the set $\partial V_{n,m}$ is a graph over ${\mathbb{Z}}^{d-1}\times\{0\}$. Namely, any line $(\overline{x},t)\in{\mathbb{Z}}^{d-1}\times{\mathbb{Z}}_{+}$ in this region intersects $\partial V_{n,m}$ by a single point. Thanks to the symmetry of the sandpile, this property also holds in the regions where the $i$-th coordinate of a point is the largest, for any $1\leq i\leq d$ .

3.3. Uniform Lipschitz estimate

In this section we prove that away from the origin the odometer function is Lipschitz uniformly in $m$, with the Lipschitz constant bounded above by the size of the model.

Assume $B_{r}$ is a ball of radius $r$ centred at the origin, and let $f:\overline{B_{r}}\cap{\mathbb{Z}}^{d}\to{\mathbb{R}}$ be harmonic in $B_{r}$. The following estimate for the discrete derivative of $f$ (“difference estimate“) is proved in [7, Theorem 1.7.1]; there exists a constant $C>0$ independent of $f$ and $r$ such that $|f(x)-f(0)|\leq Cr^{-1}||f||_{l^{\infty}}$ for any $x\sim 0$. Iterating this bound as in the proof of Corollary 3.4 we get

[TABLE]

Proposition 3.10**.**

Let $u$ be the odometer for the sandpile with initial distribution $n\delta_{0}$ and threshold $m$. Then, for any $r_{0}>0$ small there exists a constant $C=C(r_{0},d)$ such that

[TABLE]

for any $x,y\in{\mathbb{Z}}^{d}\setminus B(0,r_{0}n^{1/d})$.

**Proof. **The idea is to show that $u$ is Lipschitz in a neighbourhood of the mass-region $V_{n,1}$, and also in the neighbourhood of the origin. Then, we can conclude the Lipschitz estimate in between these two regions (where $u$ is harmonic) by DMP. For the clarity, we will split the proof into a few steps.

Step 1. Lipschitz bound near the mass-region

Fix any $x_{0}\in\partial V_{n,m}$ and take $r>0$ such that the ball $B=B(x_{0},r)\cap{\mathbb{Z}}^{d}$ does not contain the origin. Next, we write $u=u_{1}+u_{2}$ in $B$ where

[TABLE]

and

[TABLE]

Since $0=u(x_{0})=u_{1}(x_{0})+u_{2}(x_{0})$, by Lemma 3.3 we get

[TABLE]

The latter, in view of discrete Harnack, implies that $u_{2}(x)\leq Cr^{2}m$ for all $x\in B_{r/2}$, with $C>0$ depending on dimension only. This bound on $u_{2}$ and estimate (3.8) imply

[TABLE]

The last estimate coupled with Corollary 3.4 applied to $u_{1}$, gives

[TABLE]

Recall, that the thickness of the mass-region is bounded above by $C_{0}n^{1/d}m^{-1}$ according to Lemma 3.6. In view of Proposition 3.7 we see that the ball $B=B(x_{0},r)$ with $r=16C_{0}n^{1/d}m$ and $x_{0}\in\partial V_{n,m}$ does not contain the origin, in particular we have $0\leq\Delta^{1}u\leq m$ in $B\cap{\mathbb{Z}}^{d}$. With this in mind, we take $r=16C_{0}n^{1/d}m$ and varying $x_{0}$ on the boundary of $\partial V_{n,m}$, from we (3.9) we get

[TABLE]

Step 2. Lipschitz bound near the origin

Here again we will partition the solution into two parts, namely one with bounded Laplacian, and another one as the Green’s kernel. We write $u=u_{1}+u_{2}$ where

[TABLE]

and hence

[TABLE]

Notice that $u_{2}$ is the Green’s function of $V_{n,m}$ with the pole at the origin, multiplied by $n$. From Proposition 3.7 and [1, Lemma 5.1] we have

[TABLE]

for any $x\in V_{n,0}$ with $r_{0}n^{1/d}\leq|x|\leq 2r_{0}n^{1/d}$ and $y\sim x$. For $u_{1}$, we take any $x\in V_{n,0}$ satisfying $r_{0}n^{1/d}\leq|x|\leq 2r_{0}n^{1/d}$ and using the fact that $u_{1}$ has bounded Laplacian everywhere, from Lemma 3.3, choosing the radius of the ball to be $r_{0}n^{1/d}m^{-1}$, we obtain

[TABLE]

where $y\sim x$. Combining (3.11) and (3.12) leads to Lipschitz estimate

[TABLE]

where $x\sim y$, $x\in V_{n,0}$ with $r_{0}n^{1/d}\leq|x|\leq 2r_{0}n^{1/d}$.

Step 3. Interpolating between two regions

Consider the set

[TABLE]

It is left to establish a Lipschitz estimate for $u$ in $E$. Take any direction vector $\pm e_{i}$ and consider the function $w(x)=u(x+e_{i})-u(x)$ with $x\in E$. Clearly $\Delta^{1}w=0$ in $E$, and in view of estimates (3.10) and (3.13) we have $|w|\leq C_{r_{0}}n^{1/d}$ on $\partial E$. Invoking DMP we obtain $|w|\leq C_{r_{0}}n^{1/d}$ on $E$. This shows 1-step Lipschitz bound for $u$. Iterating it as in Corollary 3.4 leads to

[TABLE]

This bound, combined with (3.10) completes the proof of the proposition. ∎

3.4. $C^{1,1}$ estimates

Here we prove further regularity estimates for the odometer function. The goal is to show that discrete derivatives, as defined in subsection 1.2, are Lipschitz away from the origin, however here the Lipschitz constant would depend on $m$. We will see later, in Section 4, that these estimates, when properly scaled, force the gradient of odometer function to vanish on the boundary of (any) scaling limit of the set of visited sites.

Proposition 3.11**.**

Let $u$ be the odometer for the sandpile with initial distribution $n\delta_{0}$ and threshold $m$, where $m\geq 1$ is fixed and $n>1$ is large. Then, for any $r_{0}>0$ small there exists a constant $C=C(r_{0},d)$ such that for any unit vector $e\in{\mathbb{Z}}^{d}$, and any $x,y\in{\mathbb{Z}}^{d}$ satisfying $|x|,|y|\geq r_{0}n^{1/d}$, we have

[TABLE]

**Proof. **For $x\in{\mathbb{Z}}^{d}$ denote $w(x)=\nabla^{1}_{e}u(x)$. We fix a radius $r>0$, which is a large constant independent of $n$ and $m$, a point $x_{0}\in{\mathbb{Z}}^{d}$ with $|x_{0}|\geq r_{0}n^{1/d}$, and consider the problem in the ball $B=B(x_{0},r)\cap{\mathbb{Z}}^{d}$. With this notation, we need to show that

[TABLE]

for any unit vector $\widetilde{e}\in{\mathbb{Z}}^{d}$, since then the estimate of the proposition will follow by iteration, as in Corollary 3.4 for instance. We now fix a unit vector $\widetilde{e}\in{\mathbb{Z}}^{d}$ and suppress it from the subscript of $\nabla^{1}_{\widetilde{e}}$.

In view of the choice of $r$ and Proposition 3.10, we have

[TABLE]

with a constant $C=C(r_{0},d)$. As in the proof of Proposition 3.10, we do the splitting $w=u_{1}+u_{2}$ in $B$, where $u_{1}$ is the potential part of $w$ and $u_{2}$ is harmonic. By DMP, estimates (3.14), and (3.8) we have

[TABLE]

Thus, it is left to handle the part with $u_{1}$. To this end consider the set

[TABLE]

which is the $1$-discrete neighbourhood of the boundaries of $V_{n,m}$ and $V_{n,0}$ contained in the ball $B(x_{0},r)$. Observe, that $\Delta^{1}w(x)=0$ on $B\setminus E$. Now, using the Green’s representation of $u_{1}$ as in the proof of Proposition 3.10, for all $x\in B_{r/2}$ we get

[TABLE]

where the penultimate inequality is proved in the proof of Proposition 3.10, while the last one, with a constant depending on $r$ comes from the estimate on the number of points of $E$.

Putting together (3.15) and (3.16) we complete the proof of this proposition. ∎

Remark 3.12**.**

It should be remarked that some of the results and approaches used in this section, such as Lemmas 3.2, and 3.5, as well as Propositions 3.10, and 3.11 go in parallel with the theory of free boundary problems of the form (2.7) (with a more general r.h.s.) in continuous space (see the first chapter of [13] for instance). Although our approach in this section shares some similarities with the continuous analogues, the proofs however are much different due to the fact that we are working in a discrete space here.

4. Scaling limits

Here we prove that the sandpile shapes, generated from initial distribution concentrated at a single vertex of ${\mathbb{Z}}^{d}$, have a scaling limit (which is a ball), when the threshold $m$ is fixed, and the mass $n$ tends to infinity. Then, we show that there is also a scaling limit as $m$ tends to infinity along with $n$ but very slowly with respect to $n$. That shape is still a ball, but the entire mass is concentrated on the boundary.

We will need a few notation. For $n\geq 1$ set $h=n^{-1/d}$, and define the scaled odometer by $u_{h}(x)=h^{2}u_{n}(h^{-1}x)$ where $x\in h{\mathbb{Z}}^{d}$. Next, for $0<h\leq 1$ and $\xi=(\xi_{1},...,\xi_{d})\in h{\mathbb{Z}}^{d}$ define the half-open cube

[TABLE]

In order to study the scaling limit of the model, we need to extend each $u_{h}$ to a function defined on ${\mathbb{R}}^{d}$. We will use a standard extension of $u_{h}$ which preserves its discrete derivatives and hence $\Delta^{h}$-Laplacian. Namely, for fixed $0<h\leq 1$ define a function $U_{h}:{\mathbb{R}}^{d}\to{\mathbb{R}}_{+}$, where for each $\xi\in h{\mathbb{Z}}^{d}$ and any $x\in\mathrm{C}_{h}(\xi)$ we have set $U_{h}(x)=u_{h}(\xi)$. Clearly for any $\xi\in h{\mathbb{Z}}^{d}$ and any $x\in\mathrm{C}_{h}(\xi)$ we get

[TABLE]

Applying the scaling to the estimates of Proposition 3.10 and Proposition 3.11, we obtain that for any $\rho>0$ there exists a constant $C_{\rho}$ independent of $m\geq 1$, such that the following hold true:

• •

Uniform Lipschitz estimate

[TABLE]

• •

$C^{1,1}$-estimates

[TABLE]

where $e\in{\mathbb{Z}}^{d}$ is any unit vector.

Theorem 4.1**.**

(Scaling limit for fixed $m$)* There exists a compactly supported non-negative function $u_{0}\in C({\mathbb{R}}^{d}\setminus\{0\})$ which is spherically symmetric, and is $C^{1,1}$ outside any open neighbourhood of the origin, such that*

• (i)

$u_{h}\to u_{0}$* uniformly in ${\mathbb{R}}^{d}$ outside any open neighbourhood of the origin,*

• (ii)

$\Delta u_{0}=2dm\mathbb{I}_{\{0<u_{0}<\frac{1}{m}\}}-\delta_{0}$* in $\{u_{0}>0\}$ in the sense of distributions, where $\delta_{0}$ is the Dirac delta at the origin and $\Delta$ denotes the continuous Laplace operator,*

• (iii)

$|\nabla u_{0}|=0$* everywhere on $\partial\{u_{0}>0\}$.*

**Proof. **We first show that there is a convergent subsequence of $\{U_{h}\}$ whose limit satisfies the requirements of the theorem. Then, we conclude the proof by showing that any two convergent subsequences have the same limit. For the proof of the convergent subsequence we will follow our approach from [1, Theorem 5.3], with the only difference that here we also need to take care of the convergence of the gradient too.

Fix $\rho>0$ small, and for $0<h\leq 1$ set

[TABLE]

Since the support of $U_{h}$ is uniformly bounded in $0<h\leq 1$ thanks to Lemma 3.6 and Proposition 3.7, from (4.3) we get that $U_{h}$ is a bounded and Lipschitz function on $E_{\rho}(h)$, both uniformly in $0<h\leq 1$. We let $C_{\rho}>0$ be the Lipschitz constant of $U_{h}$ on $E_{\rho}(h)$, and define

[TABLE]

which, due to the mentioned properties of $U_{h}$, becomes a bounded and $C_{\rho}$-Lipschitz function on ${\mathbb{R}}^{d}$ which coincides with $U_{h}$ on the set $E_{\rho}(h)$. The extension defined by (4.5) is standard and well-known (see [6, Theorem 2.3] for instance). Moreover, by construction we have that the family $\{U_{h}^{\rho}\}_{0<h\leq 1}$ is uniformly bounded and is non-negative everywhere. Observe also that for $0<h\leq 1$, $\xi\in E_{2\rho}(h)$ and $x\in\mathrm{C}_{h}(\xi)$ (see (4.1)) by construction we have

[TABLE]

We also do a similar extension for the discrete derivatives. Namely, for a unit vector $e\in{\mathbb{Z}}^{d}$ set $w_{h,e}(x)=\nabla_{e}^{h}u_{h}(x)$ where $x\in h{\mathbb{Z}}^{d}$. Due to (4.4) we have that $w_{h,e}$ is $C_{\rho}$-Lipschitz on $E_{\rho}(h)$. Next, in analogy with (4.5) define

[TABLE]

which provides a Lipschitz extension for $w_{h,e}$ in the entire space ${\mathbb{R}}^{d}$. What we also obtain, in view of definition (4.7), is that

[TABLE]

where $\xi\in E_{2\rho}(h)$ and $x\in\mathrm{C}_{h}(\xi)$. Note that here we have used the fact that the extension $U_{h}$ preserves the discrete derivatives.

By construction, we have that $\{U_{h}^{\rho}\}$ and $\{W_{h,e}^{\rho}\}$ form a family of equicontinuous and uniformly bounded functions with respect to $h$, and hence, applying Arzelà-Ascoli we extract a subsequence $h_{k}\to 0$ as $k\to\infty$, such that $U_{h_{k}}^{\rho}\to U_{0}^{\rho}$ and $W_{h_{k},e}^{\rho}\to W_{0,e}^{\rho}$ locally uniformly in ${\mathbb{R}}^{d}$ as $k\to\infty$, where $U_{0}^{\rho}$, and $W_{0,e}^{\rho}$ for any $e$, are Lipschitz functions defined on ${\mathbb{R}}^{d}$. We now show that

[TABLE]

and

[TABLE]

where $\Delta$ and $\nabla_{e}$ are the Laplacian and $e$-directional derivative respectively.

We start with (4.9). Following the definition of a weak solution to (4.9), we need to check that

[TABLE]

for any $\varphi\in C_{0}^{\infty}(B)$ where $B\subset\{U_{0}^{\rho}>0\}\setminus\overline{B_{2\rho}}$ is a ball. Using discrete integration by parts, we have

[TABLE]

and comparing with (4.11), we see that what needs to be proved is that

[TABLE]

Observe, that for any $x\in\{U_{0}^{\rho}>0\}\setminus\overline{B_{2\rho}}$ where $U_{0}^{\rho}(x)\neq\frac{1}{m}$, the convergence in (4.13) is pointwise in view of the convergence of $U_{h_{k}}^{\rho}$ and estimate (4.6). Now, if $U_{0}^{\rho}(x)=\frac{1}{m}$ and $x$ is an interior point of the set $\{U_{0}^{\rho}=\frac{1}{m}\}$, then we have the convergence in (4.13) in the neighbourhood of $x$, thanks to (4.12) and the fact that $\Delta U_{0}^{\rho}=0$ in a neighbourhood of $x$. Thus, in (4.13), it is left to cover the case when $x\in\partial\{U_{0}^{\rho}=\frac{1}{m}\}$. We will show that this set has measure 0, and hence can be ignored in (4.13). Indeed, observe that $U_{0}^{\rho}$ trivially inherits, in the same form but now on ${\mathbb{R}}^{d}$, the directional monotonicity of odometers established in Theorem 3.1. But then, following Remark 3.9 (see also the proof of Proposition 3.7, and [1, Theorem 5.3 (v)]), we get that the set $\partial\{U_{0}^{\rho}=\frac{1}{m}\}$ is locally a graph of a Lipschitz function, and is hence rectifiable. In particular, the measure of $\partial\{U_{0}^{\rho}=\frac{1}{m}\}$ is 0, which completes the proof of (4.13), and hence (4.13) follows, which completes the proof of (4.9).

We next proceed to the proof of (4.10). Fix any unit vector $e\in{\mathbb{Z}}^{d}$ and any point $x_{0}\in\{U_{0}^{\rho}>0\}\setminus\overline{B_{2\rho}}$. Since $U_{0}^{\rho}$ solves (4.9), then it is $C^{1}$ in the neighbourhood of $x_{0}$, and hence it is enough to show that for any $\varphi\in C_{0}^{\infty}$ having a support in a small neighbourhood of $x_{0}$, one has

[TABLE]

where the second equality is simply integration by parts, and what needs to be proved is the first one. For that one, we observe

[TABLE]

and hence the equality in (4.10).

Due to (4.4) the function $W_{0,\rho}$ is $C^{0,1}$ in the neighbourhood of $\partial\{U_{0}^{\rho}>0\}$, and is vanishing everywhere outside the support of $U_{0}^{\rho}$, and hence from (4.10) we obtain

[TABLE]

Finally, applying diagonal argument as $\rho\to 0$, we conclude the existence of a subsequence $h_{k}\to 0$ as $k\to\infty$, and a function $U_{0}$ such that

• (a)

$U_{0}\geq 0$ is compactly supported, and is $C^{1,1}$ in ${\mathbb{R}}^{d}$ away from any open neighbourhood of the origin;

• (b)

$U_{h_{k}}\to U_{0}$ uniformly outside any open neighbourhood of the origin

• (c)

$\Delta U_{0}=2dm\mathbb{I}_{\{0<U_{0}<1/m\}}-\delta_{0}$ in ${\mathbb{R}}^{d}$.

Note, that in (c) we get the equation in ${\mathbb{R}}^{d}$ instead of only $\{U_{0}>0\}$, as $|\nabla U_{0}|=0$ on $\partial\{U_{0}>0\}$. The only thing which needs a proof, is the equation in (c) near the origin. That, again, can be handled as in [1, Theorem 5.3]. Namely, for each $h>0$ let $\Phi_{h}$ be the fundamental solution to $\Delta^{h}$ in $h{\mathbb{Z}}^{d}$. Then, $\Delta^{h}(u_{h}-\Phi_{h})=0$ in $B_{r}$, where $r>0$ is fixed small enough such that $B_{r}\subset\{u_{h}>1/m\}$. The existence of such $r$ follows from Proposition 3.7. Applying DMP we obtain, that $|u_{h}-\Phi_{h}|\leq C$ in $B_{r}$, and hence the limit $|U_{0}-\Phi_{0}|\leq C$ in $B_{r}\setminus\{0\}$ and solves $\Delta(U_{0}-\Phi_{0})=0$ in $B_{r}\setminus\{0\}$. It follows that the origin is a removable singularity for the function harmonic $U_{0}-\Phi_{0}$, implying the equality in (c) near the origin.

We now apply Lemma A.1 from Appendix, which states that there is a unique solution to (c), and that solution is spherically symmetric. With this at hand, we conclude the existence of the scaling limit of the sequence $U_{h}$, without passing to a subsequence, since we get that any sequence of scales, contains a subsequence, along which $U_{h}$ converges to the same limit.

The proof of the theorem is now complete. ∎

Remark 4.2**.**

Observe, that convergence in (4.13) implies that the set of a visited sites of the sandpile, after scaling, converges to the support of the limiting odometer. In particular, the family of odometers $\{u_{h}\}$ do not degenerate as $h\to 0$.

Now, a simple compactness argument allows us to take limits when $m\to\infty$. This in particular, gives a sandpile dynamics, which redistributes all mass on the free boundary, which has a circular shape.

Theorem 4.3**.**

(Scaling limit as $m\to\infty$)* Let $u_{n,m}$ be the odometer of the sandpile with initial distribution $n\delta_{0}$, and threshold $1<m<n$. Then there is a slowly increasing function $F:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$, with $F(+\infty)=+\infty$, such that the odometers $u_{n,F(n)}$, after rescaling by $n^{-2/d}$, converge to a spherically symmetric function, as $n\to\infty$.*

**Proof. **Thanks to Theorem 4.1, for each $m$ fixed, there exists a scaling limit for $\{u_{n,m}\}$ as $n\to\infty$. Let $u_{0,m}$ be this limit, which we know is spherically symmetric. Next, according to Lemma A.2, we get that $u_{0,m}$ converges uniformly, away from the origin, to some $u_{0}$, a spherically symmetric function. Finally, for each given $m$ we choose $n=n(m)$ large enough, such that

[TABLE]

We may obviously choose $n=n(m)$ as increasing function of $m$. Taking $m$ to infinity, and using the convergence of $u_{0,m}$ to $u_{0}$, we have a final spherical configuration. Now reverting the relation $m\mapsto n$, we get a function $n\mapsto F(n)$ which satisfies the conditions in the theorem. This completes the proof. ∎

Appendix A Uniqueness of the scaling limit

For any $m,k>0$ fixed, let $u$ (being compactly supported)333This assumption is not any restriction to our theory, as the support of odometer functions are compact. It should be remarked that one can show existence of unbounded solutions to this particular free boundary problem. solve

[TABLE]

We want to prove that such a $u$ is unique and spherically symmetric. To see this, we show that if $u$ is not spherically symmetric, then we may generate two distinct spherically symmetric solutions by starting from infimum/supremum over all rotations of $u$. Then, we prove that given the solution is spherically symmetric, then it is uniquely determined by coefficients $m$ and $k$. These two components put together settle the uniqueness of solutions to (A.1).

We shall now take the smallest super-solution solution and we only need to check that the infimum in this class does not degenerate to zero and that it solves the problem. Obviously we may only consider spherically symmetric super-solutions with shrinking radius of support. Hence integration by parts (using that $u_{*}$ vanishes outside a large ball) implies $Volume(\{0<u_{*}<k\})\geq 1$. This bound, coupled with the spherical symmetry of $u_{*}$ implies that the support of $u_{*}$ must contain a ball of a fixed radius, in particular the support cannot shrink too much. So it remains to show that the smallest super-solution is actually a solution. Now in the set $\{u_{*}>k\}$, we may always replace $u_{*}$ by $h$ which is the solution of $\Delta h=-\delta_{0}$, with boundary values $k$ on $\partial\{u_{*}>k\}$. Also in the set $\{u_{*}<k\}$ we make a replacement with solution to $\Delta h=k\mathbb{I}_{\{h>0\}}$ and the corresponding boundary values. Hence we may assume $\Delta u_{*}=m\mathbb{I}_{\{0<u<k\}}-\delta_{0}-\mu_{*}$, where $support(\mu_{*})$ is on the sphere $\partial\{u_{*}>k\}$. Next let $v_{*}$ solve $\Delta v_{*}=-\delta_{0}+m\mathbb{I}_{\{0<u_{*}<k\}}$ and $v_{*}=0$ on the boundary of the support of $u_{*}$. Then by comparison principle $v_{*}\leq u_{*}$, and hence $\Delta v_{*}\leq-\delta_{0}+m\mathbb{I}_{\{0<v_{*}<k\}}$, implying that $v_{*}$ is also a super-solution to our problem. Now $u_{*}$ being infimum in the class implies $v_{*}=u_{*}$, and hence $u_{*}$ solves (A.1).

A similar argument works by taking the largest sub-solution, among all solutions with bounded support. Here again we may consider spherically symmetric sub-solutions. We need to show that the largest sub-solution in the class stays bounded. As in previous case we may replace the maximizing sequence of sub-solutions $u_{j}^{*}$ with exact solutions in the sets $\{u_{j}^{*}>k\}$ and $\{u_{j}^{*}<k\}$. To see that the support of such sub-solutions should stay uniformly bounded we argue as follows. Let $R_{j}^{*}$ be the radius of the support of $u^{*}_{j}$, and $R_{j}$ be such that $u_{j}^{*}=k$ on $|x|=R_{j}$. Using sub-solutions properties, we can conclude that $Volume(\{0<u^{*}_{j}<k\})\leq 1$ (this is the reverse of the previous inequality). Hence $R_{j}^{*}\approx R_{j}+c_{d}R_{j}^{1-d}$, for some dimensional constant $c_{d}>0$.

Let now $F$ denote the fundamental solution, and set $F_{j}=\max(F-a_{d}(R_{j}^{*})^{2-d},0)$, for $d\geq 3$ and $F_{j}=a_{2}\max(\log(R_{j}^{*}/|x|),0)$, where $a_{d}$ is a normalization constant for $d\geq 2$. Since

[TABLE]

and

[TABLE]

we can apply comparison principle to deduce $F_{j}\geq u^{*}_{j}$ on $B_{R_{j}^{*}}$. But on the other hand for $|x|=R_{j}\approx R_{j}^{*}-c_{d}R_{j}^{1-d}$ we have

[TABLE]

once $R_{j}^{*}\approx R_{j}+c_{d}R_{j}^{1-d}$ is large enough. Hence a contradiction. This implies that the support of sub-solutions should stay bounded. Now taking the largest sub-solution, and using the fact that their supports are uniformly bounded, we get yet another solution to (A.1).

The conclusion is that if (A.1) has a solution with bounded support, which is not spherically symmetric, then we can produce two distinct spherically symmetric solutions. Hence to show uniqueness of solutions to (A.1) with bounded support, it suffices to show that there is only one spherically symmetric solution to (A.1) having bounded support. This we establish in the next lemma.

Lemma A.1**.**

For any $m,k>0$ there is a unique spherically symmetric solution to (A.1) having bounded support.

**Proof. **We shall give a computational proof, which is elementary, but tedious. First we consider the case when $d>2$. Let $\omega_{d}>0$ be the normalising constant of the Green’s kernel444It is well-known that $\omega_{d}=(d(d-2)|B_{1}|)^{-1}$ where $|B_{1}|$ is the volume of the unit ball in ${\mathbb{R}}^{d}$, but the actual value of this constant will be of no relevance to our proof.. The spherical symmetry of $u$ implies that it should be of the form

[TABLE]

where $a_{i}$ and $r_{j}$ are constants. The aim is to show that these constant are uniquely determined from (A.1), which we do next.

By a straightforward computation we have

[TABLE]

where (A.2) is due to the condition $u=0$ on $|x|=r_{2}$, (A.3) comes from $|\nabla u|=0$ on $|x|=r_{2}$, (A.4) is in view of $u=k$ on $|x|=r_{1}$ and (A.5) is the continuity of $u$ on $|x|=r_{1}$, finally (A.6) is the continuity of the gradient on $|x|=r_{1}$.

From (A.6) and (A.3) we have

[TABLE]

and hence

[TABLE]

Subtracting (A.3) from (A.5), to eliminate $a_{2}$, we obtain

[TABLE]

Plugging the value of $r_{2}$ from (A.8), we get

[TABLE]

To complete the proof we need to see that (A.10) has a unique solution when $0<r_{1}<\infty$. Rearranging (A.10), consider the function

[TABLE]

where $0<x<\infty$. We need to show that $f$ has a unique zero in $(0,\infty)$.

Observe, that $f(0+)=-\infty$, and $f(+\infty)=k>0$, and hence the continuity of $f$ gives the existence of a zero for $f$. We are left with establishing the uniqueness of this zero. To this end, computing the derivative of $f$ we get

[TABLE]

We want to show that $f^{\prime}(x)>0$ everywhere, which will complete the proof. Ignoring the $x$ in front of $f^{\prime}$ in the last expression and denoting $y:=(d-1)\omega_{d}x^{-d}$, we obtain that $f^{\prime}(x)=xF(y)$, where

[TABLE]

and hence it is enough to show that $F>0$ everywhere. The latter follows easily from the fact that $F(0)=0$ and $F^{\prime}>0$ in $(0,\infty)$. We conclude the proof of uniqueness in the case when $d\geq 3$.

We now treat the case of $d=2$. Set $\omega_{2}=-1/(2\pi)$ which is the normalising constant for the Green’s kernel. Again, the exact value of $\omega_{2}$ is not important for the proof, and we will only use that $\omega<0$. Relying on the spherical symmetry of $u$, as above, we see that $u$ has to be of the form

[TABLE]

where $a_{i}$ and $r_{j}$ are unknown constants. The goal is to show that these constants are being determined uniquely given the properties of $u$. As a above, a direct computation leads to

[TABLE]

where (A.12) is due to the condition $u=0$ on $|x|=r_{2}$, (A.13) comes from $|\nabla u|=0$ on $|x|=r_{2}$, (A.14) is in view of $u=k$ on $|x|=r_{1}$ and (A.15) is the continuity of $u$ on $|x|=r_{1}$, finally (A.16) is the continuity of the gradient on $|x|=r_{1}$.

Using (A.13) we get $a_{3}=-\frac{m}{2}r_{2}^{2}$, which together with (A.16) implies

[TABLE]

Subtracting (A.15) from (A.12) we obtain

[TABLE]

From the above relation between $r_{1}$ and $r_{2}$ we have

[TABLE]

Plugging the values of $a_{3}$ and $r_{2}$ we get

[TABLE]

and thus need to show that the last equation has a unique solution in $(0,\infty)$. To this end, set $y:=-\frac{2}{m}\frac{\omega_{2}}{r_{1}^{2}}$. With this notation, we need to prove that the function

[TABLE]

has a unique zero in the interval $(0,\infty)$. Observe, that $F(0+)=-k+\frac{\omega_{2}}{2}<0$ as $\omega_{2}<0$, and $F(+\infty)=+\infty$, hence the existence. For uniqueness, computing the derivative of $F$, we see

[TABLE]

which is always positive in the range $y\in(0,\infty)$, and the uniqueness follows.

The proof of the lemma is now complete. ∎

As we saw in Section 4, the scaling limit of the sandpile generated by a single source at the origin, is determined by the following PDE

[TABLE]

where the support of $u$ bounded and contains a ball of some fixed radius, in both cases uniformly with respect to $m$. From Lemma A.1 we have that $u$ is spherically symmetric, moreover, if we let $x_{m}$ be the radius of the ball $\{u\leq\frac{1}{m}\}$, then following (A.11) and (A.17) we have that $x_{m}$ is the unique zero in the range $(0,\infty)$ of the function

[TABLE]

for dimension $d>2$, and for dimension 2, rearranging (A.17), and plugging the value of $\omega_{2}=-\frac{1}{2\pi}$, we see that $x_{m}$ is the unique solution in the range $(0,\infty)$ of the equation

[TABLE]

Our final result shows that the limiting shapes of the sandpiles, corresponding to each $m$ converge, as $m\to\infty$. For this it will be enough to prove that the radii $x_{m}$ converge. This we do next; the proof is elementary and straightforward.

Lemma A.2**.**

The radii $x_{m}$ converge as $m\to\infty$.

**Proof. **Due to their relation with the sandpile, the set $\{x_{m}\}_{m\geq 1}$ is bounded away from zero and infinity, as we outlined above. We first consider the case when $d>2$. Rearranging the last two terms in (A.19) we get, for $x=x_{m}$, that

[TABLE]

where we have used the fact that $\{x_{m}\}$ is bounded away from 0 and infinity. Since $f(x_{m})=0$, from the last expression, taking $m\to\infty$, we see that

[TABLE]

We now settle the case of $d=2$. Defining $y=\pi mx^{2}$ in (A.20), we get

[TABLE]

Now, using that $x_{m}$, the solution to (A.20) is bounded from 0 and infinity, we do asymptotic expansion in the last expression, when $m\to\infty$ (and hence $y\to\infty$), which gives

[TABLE]

and hence, as $y=\pi mx_{m}^{2}$, we obtain that

[TABLE]

The proof of the lemma is now complete. ∎