Perturbed Divisible Sandpiles and Quadrature Surfaces

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 ℤd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {Z}^{d}$\end{document} (d ≥ 2) which continuously deforms occupied regions of the divisible sandpile model of Levine and Peres (J. Anal. Math. 111(1), 151–219 2010), by redistributing the total mass of the system onto 1m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\frac 1m$\end{document}-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 1m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\frac 1m$\end{document}. 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 1/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.


Background
In a recent work [1], the current authors introduced a new growth model on the lattice Z d (d ≥ 2) which redistributes a given initial mass on Z d onto a combinatorial free boundary. The growth rule asks vertices of 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 [13] 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 [15,16]). A QS, for a given source μ is the boundary of a domain D, that contains the support of μ, and has the property that the Poincaré balayage satisfies Bal : μ −→ H d−1 ∂D , (1.1) which is equivalent to for all h harmonic on a neighbourhood of 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 (wellbehaved) divisible sandpiles (DS), which redistribute 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 [9] for the single source, [10] for multi-source case, also [12] and [11]); 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 ( 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 ≈ 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 Fig. 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) n 1/d it is obvious that we are close to BS, and when F (n) = const, then we are having a DS. 1 In this paper we show, using compactness arguments that there exists m = F (n) ∞ 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 [13], for instance). It is however far from obvious, even for a two-point source (see Fig. 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 u m = mI {0<u m <1/m} − μ, where μ is the initial source, consisting of several point masses. A solution to this problem is not necessarily unique if μ is not a single point source. However, from free boundary technique we know that under certain good conditions (which are satisfied at least for a twosource point) the functions u m converge (for a subsequence) to a functionũ, as m tends to infinity, whereũ solves a free boundary problem of Bernoulli type 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 → ∞, as n → ∞, 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; see [7] Proposition 6.2. 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 Fig. 1 Numerical simulation of the sandpile with initial distribution concentrated at two points of Z 2 . On the left, each point (±47, 0) ∈ Z 2 carries mass 50 000, and the sandpile threshold m is set to 10. On the right we bring the two sources closer, by putting them at (±41, 0) ∈ Z 2 and keeping the rest of the parameters unchanged; the final configuration then evolves to a stadium shape. The blue region in both images depicts points of Z 2 having mass m = 10, and the gray area embraced by the blue region, carries no mass. The two marked red points in the gray area represent the support of the initial mass distribution 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.
It is worth exploring the idea of bridging free boundaries and particle dynamic processes that may define similar problems. Attempts to find and bridge these models should be valuable and interesting.

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 : hZ d → R be a given function.
where y ∼ h x means that x, y ∈ hZ d are lattice neighbours, i.e. ||x − y|| l 1 = h. When h = 1, we will simply write x ∼ y instead of x ∼ 1 y. Similar to the discrete Laplace operator, one defines discrete derivatives. Namely, for a unit vector e ∈ Z d , i.e. ||e|| l 1 = 1, the function 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 ϕ ∈ C ∞ 0 (R d ) one has 2d h ϕ → ϕ, and ∇ h e ϕ(x) → ∇ e ϕ(x) uniformly in R d , where and ∇ e are respectively the continuum Laplace operator and derivative in the direction e.
We next recall the definition of the fundamental solution to 1 in Z d . Namely, there exists a function g(x, y) : x is the Laplacian with respect to the x variable, and δ 0 is the characteristic function of 0 ∈ Z d . The existence and asymptotics of such g are well-known (see [4,8,17]) and we have where γ 0 is a constant, and |B 1 | is the volume of the unit ball in R d . The operator 1 also enjoys a maximum principle (see for instance [8,Exercise 1.4.7]).

Discrete Maximum Principle (DMP)
In this statement, and throughout the paper, the boundary of a non-empty set E ⊂ Z d , denoted by ∂E, is the set of all x ∈ Z d for which there exists x ∈ E such that x ∼ y.

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 μ 0 : Z d → R + (d ≥ 2), that is a nonnegative bounded function of finite support on Z d , and let n ≥ 0 be the total mass of the system, i.e. n = x∈Z d μ 0 (x). Fix also some threshold m > 0, and inductively define a triple (V k , u k , μ k ) for integer k ≥ 0 as follows. When k = 0 set V 0 = suppμ 0 and u 0 ≡ 0; we now define the passage from k to k + 1. At a time k call a particular site x ∈ Z d unstable if either of the following two (mutually exclusive) conditions hold: For an unstable site x we define an excess of mass denoted by E(x) and being equal to μ k (x) − m if x is unstable according to (a), and to μ k (x), if (b) holds for x. Then, a toppling of an unstable site x is the procedure of redistribution of the excess E(x) evenly among its 2d lattice neighbours.
Thus, at each discrete time k ≥ 0 we identify an unstable site, say x k and topple it. Then, otherwise. (1.4) The process terminates if there are no unstable sites. It follows easily from the definition of the process, that for any k ∈ Z + we have be any infinite toppling sequence, and u k and μ k be the odometer, and mass distribution at time k as defined above. Obviously, for each fixed x ∈ Z d the sequence {u k (x)} is non-decreasing, and hence there is a limit (odometer) u(x) := lim k→∞ u k (x), possibly infinite. Moreover, if u is finite at a given x and its lattice neighbours, then passing to the limit in Eq. 1.5 we see that μ k (x) is also convergent; let μ(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 Z d , and the final configuration (u, μ) 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 , μ 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 Z d of the singularly perturbed sandpile model defined in Section 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 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. Fig. 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δ 0 and fixed threshold m, converge to a ball under a scaling limit, as n → ∞, 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.

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 Section 1.3 will not stabilise after finitely many steps (modulo the trivial cases). Moreover, vertices in the set of visited sites On the boundary, darker colors indicate higher concentration of mass. One observes a change from a robust spherical shape to a polygonal shape of the BS as the threshold increases may become unstable infinitely many times during the lifetime of the process. This, as in the case of BS [1], brings us to the following class of the toppling sequences, which leaves a chance for the process to reach a stable state.
If not stated otherwise, in this section μ 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. 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 ∂V k (see the scanning procedure by hyperplanes as described in [1, Proposition 2.1]).
Let k ≥ 1 be given, and fix any x ∈ ∂V k \ V 0 , clearly u k (x) = 0. Assume y ∼ 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 ∼ y is not visited at time i, it follows that y ∈ ∂V i , in particular u i−1 (y) = 0, and hence y can topple only according to condition (a) for instability as defined in Section 1.3, and hence y must carry mass at least m. This means that at time k, when x ∈ ∂V k , the 1-neighbourhood of x has total mass at least m/(2d).
For each x ∈ Z d let B(x) be a 1-box centered at x, i.e. B(x) = {x} ∪ {y ∈ Z d : y ∼ x}. On one hand we proved that the total mass carried by vertices of B(x) where x ∈ ∂V k \V 0 is at least m/(2d). On the other hand, each x ∈ ∂V k \V 0 can be counted at most 2d times in each box B, and hence, we conclude that the total number of boundary points is bounded above by (2d) 2 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) ≥ u 1 (x) for any x ∈ Z d , which by symmetry will imply the desired result.
For i = 1, 2 by u i,k denote the odometer function and by μ i,k the distribution corresponding to the sequence T i after the k-th toppling occurs. Let 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 ≤ 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 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 ≤ 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 where the last equality follows from Eq. 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 = x k one has 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 ≤ k be the last time x has toppled prior to time k. Observe that i < k since x = x k . Then we have where the second inequality follows by inductive hypothesis. We thus have proved (2.3) for all x = x k . Consider the following inequality To complete the induction step suppose for a moment that (2.4) holds true. Then Rearranging the first and the last terms leads to where the last inequality is due to Eq. 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, If we are in case (a), then μ 1,k (x k ) = m in view of the toppling rule, and hence (2.4) follows in view of the stability of μ 2 . Next, if x k is set to topple because of (b), then u 1,k−1 (x k ) > 1 m n 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 From here, and the stability of μ 2 we obtain μ 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 where e i ∈ Z d is the i-th element of the standard basis. Let T be a hyperplane through the origin of Z d and with a normal collinear with some element of N . Then, any sandpile generated by initial distribution of the form nδ 0 , is symmetric with respect to T . Namely, if u is the odometer, and x ∈ Z d is arbitrary, for x * −the mirror reflection of x with respect to T (which is obviously in Z d in view of the choice of T ) one has 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 ± 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. Proof For each k ∈ Z + let u k be the odometer function under the additional restriction imposed by the barrier u b , and let μ 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) ≤ u b (x) for all k ∈ Z + and all x ∈ Z d , can be dropped. Then, the claim of the lemma will follow by the abelian property proved in Proposition 2.2.
Fix k ≥ 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, Assume we are in case (a). Then, according to the sandpile rule defined in Section 1.3 we need to topple the excess from m, i.e. μ k (x k ) − m. Following the restriction imposed by the barrier u b , the maximum amount of mass we can move out from i.e. we cannot move the entire excess. Set Then, invoking the k-th toppling, we will move 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 μ 0 we set u : Z d → R + , V ⊂ Z d , and μ : Z d → R + for its odometer, set of visited sites, and the final distribution respectively. In particular we have Let us also see what is the discrete PDE problem solved by odometer. If x ∈ Z d and 0 < u(x) ≤ 1 m n 2/d , then x has always toppled according to rule (a), and hence μ(x) = m. Next, if u(x) > 1 m n 2/d , then eventually x starts to topple according to (b), and hence we get μ(x) = 0. Finally, if x ∈ ∂{u > 0} then u(x) = 0, and hence x has never toppled, which implies μ ≤ m. Summarising the discussion, we conclude where I stands for the characteristic function of a set. Observe, that if the initial distribution μ 0 is already stable, then the odometer is identically zero.

The Smallest Super-Solution
We define a class of super-solutions to our sandpile model.

Definition 2.2 (Super-solutions)
Let μ 0 : Z d → R + be a given initial distribution with total mass n, and let m > 0 be fixed. Call u : Denote by W the set of super-solutions, which is non-empty since the odometer of the sandpile is in W. We have the following.
Proof Let us first show that for any w 1 , w 2 ∈ W one has u := min{w 1 , where the second inequality is in view of w 1 ∈ W. Similarly, if u(x) > 1 m n 2/d , then both w 1 (x), w 2 (x) > 1 m n 2/d , and arguing as above we obtain 1 u(x) + μ 0 (x) ≤ 0. We thus have condition (ii) for u, and (i) follows similarly, in a straightforward manner. In particular, we get that W is closed under taking minimum of finitely many elements. We next show that u * ∈ W. Let u 0 be the odometer of the sandpile. Since u 0 ∈ W, and hence 0 ≤ u * ≤ 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 ε > 0 and each Since ε > 0 is arbitrary, the claim of the lemma follows from Eq. 2.8 by a similar argument as we had above.
It follows from the definition of u and u * that u − u * is a sub-solution, hence from DMP we get that u ≤ u * in V * . Since both u and u * are 0 outside V * we obtain u ≤ u * on 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 ∈ V * then inequality of (i) follows from Eq. 2.9, otherwise, if x ∈ Z d \ V * we have u(x) = u * (x) = 0 which, combined with (i) for u * , and the inequality u ≤ u * on Z d , implies condition (i) of a super-solution for u. We next proceed to (ii). It is clear from Eq. 2.9 that it suffices to show that 1 u(x) + μ 0 (x) ≤ 0 if u(x) > 1 m n 2/d . But since u * ≥ u everywhere, for such x we have u * (x) > 1 m n 2/d and hence x / ∈ E. The latter and (2.9) imply 1 u(x) + μ 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 , and do the toppling following T , with an additional requirement that at any time k the cumulative emissions of mass from any x ∈ 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 ≤ u * which is due to definition of v 0 . Let x ∈ 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 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 μ corresponding to v 0 is not stable at x. This means that either μ(x) > m, or we have both μ(x) > 0 and v 0 (x) > 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 ∈ 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.

Analysis of Discrete Shapes for Single Sources
In this section we analyse the shape of the sandpile in the discrete space Z d generated by initial distribution nδ 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.

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 Z d , enjoys monotonicity in lattice directions. More precisely, denote by S the set of mirror symmetry hyperplanes of the unit cube 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. Then, for any X 1 , X 2 ∈ Z d , such that X 1 −X 2 is a non-zero vector orthogonal to T , we have 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 ∈ 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 (Z d ) * = Z d , in particular X 1 = (X 2 ) * . Set 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 ∈ V − . As in [1,Theorem 4.5], here as well we get 1 u T (x) ≤ 1 u(x), which, together with the equality u T 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 1 u T (x) ≤ 1 u(x * ). Observe that if u T (x) ≤ 1 m n 2/d then the last inequality gives the stability of the sandpile at x. Otherwise, if u T (x) > 1 m n 2/d we get u(x * ) > 1 m n 2/d and hence 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, completing the proof of the theorem.

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δ 0 and mass threshold equal to m. We will be using the following notation (cf. Fig. 2): • V n,m ⊂ Z d − the set of visited sites of the sandpile, • V n,0 ⊂ V n,m − the set where the odometer u > 1 m n 2/d ; in particular u is harmonic on V n,0 \ {0} and (thanks to the discrete monotonicity of u) the origin lies in V n,0 , • V n,1 := V n,m \ V n,0 − the set where the entire mass of the sandpile is concentrated. where we have used the fact that the ball B has no intersection with the region V n,0 . Since u = 0 on ∂V n,m we have w < 0 on 2 which, together with Eq. 3.1 implies where c d > 0 is some small constant depending on dimension d only. Rearranging the last inequality completes the proof of the lemma. Then, we have Proof We start with (i). For x ∈ B we have the representation 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 where T is the first exit time from B of the walk, g is the fundamental solution of 1 defined in Section 1.2, and E x stands for the expectation conditioned that the walk has started from x. Now fix any x, y ∈ B r/2 such that x ∼ y. Then, from Eq. 3.
where we have used the asymptotics (1.3) and that x ∼ y to estimate the difference of g-s, and the fact that x ∈ B r/2 in addition to those to bound the sum. With this estimate getting back to sum in Eq. 3.4 we obtain w∈B\{x,y} where the sum of the first term is estimated by a simple counting argument relying on the structure of B (see, e.g. [9, Lemma 5.2]), while bound on the sum involving the second summand follows from a trivial estimate |B| ≤ C d r d . Returning to Eq. 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 ∼ 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 Clearly, 1 v = −m everywhere and v ≤ 0 on ∂B. From here, we get and hence min We thus have 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 ∈ B r/2 one has
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 where x i ∈ B r/2 for 0 ≤ i ≤ k. Clearly, we can assume k |x − y|, by considering a path of the shortest length. Now, (i) follows from 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 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 . Since the origin does not lie in B r , we have 0 ≤ 1 u ≤ 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 )| ≤ Cmr 2 , consequently 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 [8,Theorem 1.7.2]), that Since u 1 ≤ 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 massregion. We are now in a position to also give a bound from below. 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 ∈ B r+1 ∩Z d such that u(z) ≥ n 2/d m . Now, applying Lemma 3.5 we obtain which completes the proof of the lemma.
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) ∈ Z d is inside V n,0 but is not an interior point of the set V n,0 . In particular, u(X R ) > 1 m n 2/d but for some y ∼ X R we have u(y) ≤ 1 m n 2/d . Let S R be the simplex with vertices at ±Re i , where i = 1, 2, ..., d, in other words S R is the ball of radius R in l 1 metric. Let us prove the following inclusions: Notice, that both the simplex and the cube in Eq. 3.7 are restricted to Z d . We take any x ∈ 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 ± e j , and obviously so is the simplex S R , and hence, it will be enough to prove (3.7) for X 0 = (x 1 , ..., x d ) ∈ Z d satisfying x d ≥ |x i |, i = 1, ..., d − 1. To accomplish that, we will use the directional monotonicity of the odometer.
in particular X 0 , X R ∈ H * . Next, consider the cone Since ν i · ν j ≥ 0 for all 1 ≤ i, j ≤ 2d − 1, we get C 0 ⊂ H * . We now translate C 0 to X 0 by setting C = X 0 + C 0 , and as X 0 ∈ H * we get C ⊂ H * . In view of the choice of the collection {ν i } 2d−1 i=1 and the points X 0 , X R , it is easy to see X R ∈ C. Finally, in the cone C we use the directional monotonicity given by Theorem 3.1 which implies that in C the odometer u attains its maximum at the vertex of the cone, i.e. at X 0 . Since X R ∈ C we obtain u(X 0 ) ≥ u(X R ), and hence X 0 ∈ V n,0 which completes the proof of the first inclusion of Eq. 3.7. The second inclusion, is much simpler, and follows easily by using monotonicity in the coordinate directions only.
With Eq. 3.7 at hand, the proof of the proposition will be complete, once we show that R 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 ∂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 In view of the monotonicity of u in the direction of e d , we have ⊂ V n,0 . Moreover, definitions of H * and S R imply that | | R d−1 with constants in the equivalence depending on dimension d only. Now, for a given X = (x, 0) ∈ where x ∈ Z d−1 , let t 1 ∈ Z + be the smallest integer such that X 1 = (x, t 1 ) / ∈ V n,0 and let t 2 ≥ t 1 be the smallest integer such that u vanishes at X 2 = (x, t 2 ). In view of the choice of X 1 and the discussion above, the cone X 1 + C 0 has no points of V n,0 . Moreover, X 2 ∈ C 0 and according to Lemma 3.6 we obtain that t 2 − t 1 r where r = 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 intersects the mass-region by r points, and since each point in V n,1 carries mass m and the total mass of the system is preserved, we get R d−1 rm n, from this, and the estimate of Lemma 3.6, we get R 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 the set ∂V n,m is a graph over Z d−1 ×{0}. Namely, any line (x, t) ∈ Z d−1 × Z + in this region intersects ∂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 ≤ i ≤ d .

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 centered at the origin, and let f : B r ∩ Z d → R be harmonic in B r . The following estimate for the discrete derivative of f ("difference estimate") is proved in [8, Theorem 1.7.1]; there exists a constant C > 0 independent of f and r such that |f (x) − f (0)| ≤ Cr −1 ||f || l ∞ for any x ∼ 0. Iterating this bound as in the proof of Corollary 3.4 we get for any x, y ∈ Z d \ 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 ∈ ∂V n,m and take r > 0 such that the ball B = B(x 0 , r)∩ Z d does not contain the origin. Next, we write u = u 1 + u 2 in B where 1 u 1 = 1 u in B and u 1 = 0 on ∂B r , and 1 u 2 = 0 in B and u 2 = u on ∂B r .
The latter, in view of discrete Harnack, implies that u 2 (x) ≤ Cr 2 m for all x ∈ B r/2 , with C > 0 depending on dimension only. This bound on u 2 and estimate (3.8) imply The last estimate coupled with Corollary 3.4 applied to u 1 , gives 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 ∈ ∂V n,m does not contain the origin, in particular we have 0 ≤ 1 u ≤ m in B ∩ Z d . With this in mind, we take r = 16C 0 n 1/d m and varying x 0 on the boundary of ∂V n,m , from (3.9) we get |u(x) − u(y)| ≤ Cn 1/d |x − y|, ∀x, y ∈ {z ∈ Z d : dist(z, V n,1 ) ≤ r/4} := V * . (3.10)

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 1 u 1 = 1 u − nδ 0 in V n,m and u 1 = 0 on ∂V n,m , and hence 1 u 2 = −nδ 0 in V n,m and u 2 = 0 on ∂V n,m .
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 for any x ∈ V n,0 with r 0 n 1/d ≤ |x| ≤ 2r 0 n 1/d and y ∼ x. For u 1 , we take any x ∈ V n,0 satisfying r 0 n 1/d ≤ |x| ≤ 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 where y ∼ x. Combining (3.11) and (3.12) leads to Lipschitz estimate where x ∼ y, x ∈ V n,0 with r 0 n 1/d ≤ |x| ≤ 2r 0 n 1/d .

Step 3. Interpolating between two regions
Consider the set This bound, combined with (3.10) completes the proof of the proposition.

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 Section 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δ 0 and threshold m, where m ≥ 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 ∈ Z d , and any x, y ∈ Z d satisfying |x|, |y| ≥ r 0 n 1/d , we have Proof For x ∈ Z d denote w(x) = ∇ 1 e u(x). We fix a radius r > 0, which is a large constant independent of n and m, a point x 0 ∈ Z d with |x 0 | ≥ r 0 n 1/d , and consider the problem in the ball B = B(x 0 , r) ∩ Z d . With this notation, we need to show that for any unit vector e ∈ 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 e ∈ Z d and suppress it from the subscript of ∇ 1 e . In view of the choice of r and Proposition 3.10, we have 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 Thus, it is left to handle the part with u 1 . To this end consider the set 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 1 w(x) = 0 on B \ E. Now, using the Green's representation of u 1 as in the proof of Proposition 3.10, for all x ∈ B r/2 we get (3.16) 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 [14] 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.

Scaling Limits
Here we prove that the sandpile shapes, generated from initial distribution concentrated at a single vertex of 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 ≥ 1 set h = n −1/d , and define the scaled odometer by u h (x) = h 2 u n (h −1 x) where x ∈ hZ d . Next, for 0 < h ≤ 1 and ξ = (ξ 1 , ..., ξ d ) ∈ hZ d define the half-open cube In order to study the scaling limit of the model, we need to extend each u h to a function defined on R d . We will use a standard extension of u h which preserves its discrete derivatives and hence h -Laplacian. Namely, for fixed 0 < h ≤ 1 define a function U h : R d → R + , where for each ξ ∈ hZ d and any Clearly for any ξ ∈ hZ d and any x ∈ C h (ξ ) we get Applying the scaling to the estimates of Proposition 3.10 and Proposition 3.11, we obtain that for any ρ > 0 there exists a constant C ρ independent of m ≥ 1, such that the following hold true: where e ∈ Z d is any unit vector. 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 ρ > 0 small, and for 0 < h ≤ 1 set Since the support of U h is uniformly bounded in 0 < h ≤ 1 thanks to Lemma which provides a Lipschitz extension for w h,e in the entire space R d . What we also obtain, in view of definition (4.7), is that and  (4.12) and comparing with Eq. 4.11, we see that what needs to be proved is that We will show that this set has measure 0, and hence can be ignored in Eq. 4.13. Indeed, observe that U ρ 0 trivially inherits, in the same form but now on 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 ∂{U ρ 0 = 1 m } is locally a graph of a Lipschitz function, and is hence rectifiable. In particular, the measure of ∂{U ρ 0 = 1 m } is 0, which completes the proof of Eq. 4.13, and hence (4.11) follows, completing the proof of Eq. 4.9.
We next proceed to the proof of Eq. 4.10. Fix any unit vector e ∈ Z d and any point Since U ρ 0 solves (4.9), then it is C 1 in the neighbourhood of x 0 , and hence it is enough to show that for any ϕ ∈ C ∞ 0 having a support in a small neighbourhood of x 0 , one has where the second equality is simply integration by parts, and what needs to be proved is the first one. For that one, we observe (4.14) and hence the equality in Eq. 4.10. Due to Eq. 4.4 the function W 0,ρ is C 0,1 in the neighbourhood of ∂{U ρ 0 > 0}, and is vanishing everywhere outside the support of U ρ 0 , and hence from Eq. 4.10 we obtain |∇U ρ 0 | = 0 on ∂{U ρ 0 > 0}. (4.15) Finally, applying diagonal argument as ρ → 0, we conclude the existence of a subsequence h k → 0 as k → ∞, and a function U 0 such that (a) U 0 ≥ 0 is compactly supported, and is C 1,1 in R d away from any open neighbourhood of the origin; (b) U h k → U 0 uniformly outside any open neighbourhood of the origin (c) Note, that in (c) we get the equation in R d instead of only {U 0 > 0}, as |∇U 0 | = 0 on ∂{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 h be the fundamental solution to h in hZ d . Then, h (u h − h ) = 0 in B r , where r > 0 is fixed small enough such that B r ⊂ {u h > 1/m}. The existence of such r follows from Proposition 3.7. Applying DMP we obtain, that |u h − h | ≤ C in B r , and hence the limit |U 0 − 0 | ≤ C in B r \ {0} and solves (U 0 − 0 ) = 0 in B r \ {0}. It follows that the origin is a removable singularity for the harmonic function U 0 − 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.
We shall now take the smallest super-solution u * and we only need to check that the infimum in the class of super-solutions does not degenerate to zero (i.e. u * ≡ 0) 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 V olume({0 < u * < k}) ≥ 1/m. 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 supersolution is actually a solution. Now in the set {u * > k}, we may always replace u * by h which is the solution of h = −δ 0 , with boundary values k on ∂{u * > k}. Also in the set {u * < k} we make a replacement with solution to h = kI {h>0} and the corresponding boundary values. Hence we may assume u * = mI {0<u<k} − δ 0 − μ * , where support (μ * ) is on the sphere ∂{u * > k}. Next let v * solve v * = −δ 0 + mI {0<u * <k} and v * = 0 on the boundary of the support of u * . Then by comparison principle v * ≤ u * , and hence v * ≤ −δ 0 + mI {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 V olume({0 < u * j < k}) ≤ 1 (this is the reverse of the previous inequality). Hence R * j ≈ R j + c d R 1−d j , 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 ≥ 3 and F j = a 2 max(log(R * j /|x|), 0), where a d is a normalization constant for d ≥ 2. Since F j = −δ 0 ≤ −δ 0 + mI {0<u * j <1/k} ≤ u * j , in B R * j and F j > 0 = u * j on ∂B R * j , we can apply comparison principle to deduce F j ≥ u * j on B R * j . But on the other hand for 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 Eq. A.1.
The conclusion is that if Eq. 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 Eq. A.1 with bounded support, it suffices to show that there is only one spherically symmetric solution to Eq. A.1 having bounded support. This we establish in the next lemma.