# Limit shape of the leaky Abelian sandpile model with multiple layers

Théo Ballu, Cédric Boutillier, Sevak Mkrtchyan, Kilian Raschel

► To cite this version:

Théo Ballu, Cédric Boutillier, Sevak Mkrtchyan, Kilian Raschel. Limit shape of the leaky Abelian sandpile model with multiple layers. 2025. hal-04905111

HAL Id: hal-04905111

https://hal.science/hal-04905111v1

Preprint submitted on 22 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Limit shape of the leaky Abelian sandpile model with multiple layers

Théo Ballu¹, Cédric Boutillier²,³, Sevak Mkrtchyan⁴, and Kilian Raschel¹,⁵

¹Univ Angers, CNRS, LAREMA, SFR MATHSTIC, F-49000 Angers, France

²Sorbonne Université, CNRS, Laboratoire de Probabilités, Statistique et Modélisation, LPSM, UMR 8001, F-75005 Paris, France

³Institut Universitaire de France

⁴University of Rochester, Department of Mathematics, Rochester, NY, United States

⁵CNRS, International Research Laboratory France-Vietnam in Mathematics and its Applications, Vietnam Institute for Advanced Study in Mathematics, Hanoi, Vietnam

January 22, 2025

## Abstract

In this paper we study a triple generalization of the Leaky Abelian Sandpile Model (LASM) of Alevy and Mkrtchyan, originally analyzed in the case of the square lattice in dimension two. First, we work in any dimension. Second, each site can hold several different stacks of sand, one for each of a certain given number of different layers or colors. Third, when a stack of one color at a site topples, it can send sand not only to its nearest neighbors in equal amounts, but to all possible locations and colors, according to a fixed but arbitrary mass distribution. Stacks of different colors can topple according to different distributions and different leakiness parameters, however the toppling rule should be site-independent. We obtain three main results. First, in this generality, when the LASM is started with $N$ grains of sand in one color at the origin, the final stable configuration, after scaling down by $\log N$, converges to a limit shape as $N$ goes to infinity. Second, when the leakiness parameter converges to infinity and the toppling distribution has finite range, the limit shape converges to a polytope. Third, when the leakiness parameters converge to one, which means the leakiness disappears, the limit shape of the sandpile converges to an ellipsoid. From a technical point of view, we rely on a strong relation between the Green function for random walk and the shape of the sandpile. Finally, the limit shape exhibits interesting duality properties, which we also investigate.

## 1 Introduction and main results

The Abelian Sandpile Model (ASM) is a model introduced by Bak, Tang and Wiesenfeld in 1987 [5] as a “toy model” for studying self-organized criticality. In its basic form it can be thought of as a cellular automaton evolving on an arbitrary graph $G = (V, E)$, where $V$ is the set of vertices and $E$ the set of edges. A sandpile configuration on $G$ is a function $s: V \to \mathbb{Z}_{\ge 0}$. Given a configuration $s$, it evolves in discrete time according to the following principle: if the size $s(x)$ of the sandpile at a site $x$ is at least as large as the vertex degree of $x$, then $x$ “topples”, sending one grain of sand to each neighboring vertex. It was observed by Bak, Tang and Wiesenfeld that this model organizes itself into a critical state in the following sense. If the model is started with a random

<footnote>Keywords: Abelian sandpile model, killed random walk, Green function, limit shape
2020 Mathematics Subject Classification: 60G50; 60K35; 82B41</footnote>
configuration on a bounded portion of the square lattice $\mathbb{Z}^2$ and evolved until it stabilizes, the final stable configuration is critical, which means: if a small random addition of sand is made at a location, this change propagates throughout the lattice at all scales.

It is thus of interest to understand what can be said of this final, critical configuration. More precisely, what can be said about the final configuration as $N$ grows if the sandpile is started with the delta initial condition of $N$ grains of sand in one location? Simulations suggest that the region visited by the sandpile converges to a limit shape, and the final configuration has intricate local structure; see Figure 1. While there is an understanding of the local structure in terms of Apollonian circle packings [23, 30], few results about the limit shape have been obtained. It was shown in [24] that the final configuration is bounded between two circles. The existence of a limit shape was established in [29] while in [1] it was shown that the shape is Lipschitz.

A modified, dissipative or leaky version of the model was considered by Richard Kenyon and Lionel Levine. On the square lattice $\mathbb{Z}^2$ the leaky model was studied in [2]. A new parameter $1 \le m \in \mathbb{R}$ was introduced, which controls how leaky the system is. The function $s$ now takes nonnegative real values and if a site $x$ has at least $4m$ grains of sand, then it topples, sending 1 grain of sand to each neighbor, with the remaining $4(m-1)$ grains disappearing from the system. The quantity $1 - \frac{1}{m}$ can naturally be interpreted as a mass (for an underlying massive Laplacian operator), so the case $m=1$ is sometimes called massless. In the massless case $m=1$ this specializes to the original ASM. In [2] the precise limit shape of the uniform nearest neighbor Leaky Abelian Sandpile Model (LASM) on the 2-dimensional square lattice was computed for any $m > 1$. It was also shown, that the limit shape converges to an $L_1$ ball as $m \to \infty$ and to a circle as $m \to 1$. It is worth noting, that even though the limit shape as $m \to 1$ converges to a circle, the limit shape when $m=1$, i.e. in the case of the original ASM, is not expected to be a circle, so the limit shape, as a function of $m$ should be discontinuous at 1.

In the present paper we study the LASM in a vastly more general setup. First, instead of working in dimension 2, we work in arbitrary dimension. Second, each site can hold several different stacks of sand, one for each of $p$ different layers or colors. Third, when a stack of one color at a site topples, it can send sand not only to the nearest neighbors in equal amounts, but to all possible locations and colors, and according to a fixed but arbitrary mass distribution. Stacks of different colors can topple according to different distributions and different leakiness parameters, however the toppling rule should be site-independent. The only restrictions we place are some irreducibility and aperiodicity conditions, to ensure that all sites can be reached in an aperiodic way, and a moment condition to ensure mass doesn't quickly escape to infinity (see Assumption 1).

Our main result is a proof that, in this generality, when the LASM is started with $N$ grains of sand in one color at the origin, the final stable configuration, after scaling down by $\log N$, converges to a limit shape (Theorem 31). To state the result, we need to introduce the function $\Gamma^{-1}: \mathbb{S}^{d-1} \to \mathbb{R}^d$ as in Proposition 14; this function actually parametrizes the exponential decay of a Green function at infinity, see Proposition 5.

**Theorem 1** (Theorem 31). In direction $u \in \mathbb{S}^{d-1}$, for fixed $N$, the shape of the final configuration of the sandpile lies between two radii called $r_{N,u}$ and $R_{N,u}$. We have

$$ \lim_{N \to +\infty} \frac{r_{N,u}}{\log N} = \lim_{N \to +\infty} \frac{R_{N,u}}{\log N} = \frac{1}{\Gamma^{-1}(u) \cdot u} $$

and these limits are uniform in $u$. This means that the limit shape of the LASM is delimited by the curve

$$ C = \left\{ \frac{u}{\Gamma^{-1}(u) \cdot u} \mid u \in \mathbb{S}^{d-1} \right\}, \qquad (1) $$
and this limit holds for the Hausdorff distance (and in a stronger uniform sense, to be introduced in Definition 22).

See Figure 1 for various illustrations of Theorem 1. Moreover, in Theorem 31 we prove that the deviation of the unscaled final configuration from the limit shape scaled up by $\log N$ is only of constant order. However, the exact formula (1) for the limit shape is not easy to compute explicitly, since it involves inverting the gradient of a given spectral radius, see Proposition 14. As a further result (Proposition 32), we identify the limit curve with the dual of a curve which is the level set of the spectral radius mentioned above. For example, if there is only one color ($p = 1$), then the spectral radius reduces to the Laplace transform of the increments of a random walk, whose level sets are well understood. Back to the multicolor setting, as the boundary of a convex shape, it follows that the limit shape itself bounds a convex region.

Note, that the fact that the sandpile is started with $N$ grains of sand in one location and one color is not essential. If the sandpile is started from finitely many locations, with different number of grains in each of those location and in any colors, and $N$ is the largest number of grains of sand in any one color and location, then the limit shape will be exactly the same as that in Theorem 1.

We study two limiting regimes. First, we allow the leakiness parameter to converge to infinity at the same rate in all colors and study the case when the toppling distribution has finite range. In this regime we show that the limit shape converges to a polytope (Theorem 39), which depends only on the support of the toppling distribution, and not the distribution itself. Let $\mathcal{X}$ be the polytope defined in (38).

**Theorem 2 (Theorem 39).** The rescaled shape of the sandpile, $(\log m) \cdot C$, where $C$ is defined in (1), converges to the convex hull of $\mathcal{X}$ as $m$ goes to $+\infty$. The convergence is uniform in the sense of Definition 22. In particular, in the uncolored case, the limit shape of the sandpile when $m \to \infty$ is the convex hull of the support of the toppling distribution.

See the second row of Figure 1 for an illustration of Theorem 2. The limit polytope can be described as the set of points which have the same first-passage-time for the underlying walk of the toppling distribution (Theorem 40). This limiting polytope is the equivalent of the $L_1$ ball that was observed in [2].

The second limit we study is when the leakiness parameters converge to 1, which means the leakiness disappears. We show that in this limit the limit shape of the sandpile converges to an ellipsoid (Theorem 43).

**Theorem 3 (Theorem 43).** If the drift of the model is zero (in a sense to be defined in Section 6.3), then the rescaled shape of the sandpile, $\sqrt{m-1} \cdot C$, where $C$ is defined in (1), converges to the ellipsoid

$$ \{s \in \mathbb{R}^d \mid 2s^T \sigma^{-1} s \le 1\} $$

as $m$ tends to 1, where $\sigma^{-1}$ can be interpreted as the inverse of a covariance matrix. The convergence is uniform in the sense of Definition 22.

Note, however, that this does not identify the limit shape of the model when the leakiness parameters are in fact equal to 1, which corresponds to the ordinary, non-leaky ASM. It follows from our work that the limit shape as a function of the leakiness parameters is continuous away from 1, but at 1 it is expected to be discontinuous.

All our limit shape results hold for uniform convergence in spherical coordinates (Definition 22), which by our Proposition 24 implies convergence in Hausdorff distance (recalled in Definition 23).

To obtain the limit shape results, we establish a link between the limit shape of the LASM and the Green function of a related killed random walk, showing that the limit shape can be bounded
between two level curves of the Green function, thus reducing the study of the limit shape to the study of the corresponding Green function. See our Proposition 21, which proves that if the Green function is above (resp. below) a threshold, then a given site is (resp. is not) in the shape of the final configuration. This connection between the sandpile and the Green function is crucial for at least two reasons. First, from a technical point of view, we can use existing results on the asymptotics of the Green function. Second, the Green function is ubiquitous in random walk theory, and more generally in probability theory, so it has been studied using various approaches. In Section 3 we review some of the literature and mention five interpretations of the exponential decay, all of which are relevant to the sandpile model.

To give some more details, by applying a Doob transformation, the killed random walk is transformed into a non-killed random walk, the asymptotics of whose Green function is classical [28] when $p = 1$ and was recently obtained in [13, 6] in the multicolor setting. The exponential decay of the Green function established in [6] implies that the distance between the two level curves of the Green function that bound the boundary of the region visited by the sandpile is of constant order, whereas their scale is of logarithmic order in the number of the initial grains of sand, establishing the limit shape.

Note, that in the regime when the leakiness converges to 1 fast enough, the distance between the two level curves grows, becoming comparable to their sizes, making it impossible to establish a limit shape result. This is partly the justification why the link between these Green functions and the limit shape breaks in that regime, making the discontinuity mentioned above possible.

Instead of multiple colors, the model we study can be equivalently described in terms of a single-color model on a larger graph with a site-dependent toppling distribution. More precisely, consider a graph $\mathcal{G}$ on which $\mathbb{Z}^d$ acts freely, but not transitively, and for which the quotient $F = \mathcal{G}/\mathbb{Z}^d$ has a finite number of vertices. On $\mathcal{G}$ consider a vertex-dependent toppling distribution which is invariant under the action of $\mathbb{Z}^d$. Such a model is equivalent to a colored model on $\mathbb{Z}^d$ with the number of colors equal to the size of the fundamental domain $F$. Note, that adding colors to the model on $\mathcal{G}$ would not change anything, as such a model could still be equivalently formulated as a colored model on $\mathbb{Z}^d$ by just adding extra colors.

The model with $p$ colors on $\mathbb{Z}^d$ can also equivalently be described as a single-color model on the vertices of the graph $\mathbb{Z}^d \times \{1, \dots, p\}$. Also note, that if a lattice can be embedded in the square lattice, then our results extend to such a lattice as well. For example, our results apply to the 2-dimensional hexagonal lattice as well.

In Appendix B we connect some of our results to the geometry of an amoeba, which is naturally associated with the model. In particular, we show that the limit shape (1) is the dual of the connected component of the complement of the amoeba.

## 2 Leaky sandpile model and the associated random walk

In the following, $d, p, N \in \mathbb{N} := \{1, 2, \dots\}$. A sandpile configuration is a function $s: \mathbb{Z}^d \times \{1, \dots, p\} \to \mathbb{R}_{\ge 0}$. To each $(x, i, j) \in \mathbb{Z}^d \times \{1, \dots, p\}^2$, we associate $c(x, i, j) \ge 0$ and to each $i \in \{1, \dots, p\}$, we associate a leakiness parameter $m_i \ge 1$.

**Model 1 (Leaky Abelian Sandpile Model (LASM) on $\mathbb{Z}^d \times \{1, \dots, p\})$.** We assume that for all $i \in \{1, \dots, p\}$, $\sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(y, i, j) < +\infty$.

* The initial configuration is $s_0 := N\delta_{(0,i_0)}$. It consists of $N$ grains of sand at the origin of $\mathbb{Z}^d \times \{1, \dots, p\}$.
<figure><img src="image_2.png" /><figcaption>Figure 1: We consider the following sandpile model with four colors. Color 1 sends 1 chip to location (0,0,0) to each of the colors 2,3,4, $m = 4/3$ so 1 chip is lost. Color 2 sends 2 chips to color 1, to the locations (0,0,-1) and (0,0,1), $m = 2$, so no chip is lost. Color 3 sends 2 chips to color 1, to the locations (0,-1,0) and (0,1,0), $m = 2$, so no chip is lost. Color 4 sends 2 chips to color 1, to the locations (-1,0,0) and (1,0,0), $m = 2$, so no chip is lost. For each of the three rows, the leftmost picture shows the convex hull of the final configuration; the second picture shows a unit sphere at each point of the boundary of the visited region, the darker the sphere, the larger the height; the four right pictures show the middle horizontal slice of the final configuration in each of the 4 colors. The first row represents the above sandpile model with initial configuration of $10^{34}$ chips. The second row represents the same sandpile with the exception that $m = 10^8$ for color 1, and with initial configuration of $10^{300}$ chips. The last row represents the same sandpile with the exception that $m = 3$ for color 1, so it is not a leaky model, and with initial configuration of $10^6$ chips.</figcaption></figure>

* A sandpile configuration $s$ is said to be stable at site $(x, i)$ if:

$$s(x, i) \le m_i \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(y, i, j).$$

We say a configuration is stable if it is stable at every site. Otherwise, we say it is unstable.

* As long as a configuration is unstable at site $(x, i)$, it can topple, which means:

$$s(x, i) \leftarrow s(x, i) - m_i \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(y - x, i, j), \\ \forall (y, j) \in \mathbb{Z}^d \times \{1, \dots, p\}, \quad s(y, j) \leftarrow s(y, j) + c(y - x, i, j).$$

* Sites topple until the configuration eventually becomes stable.

Each time a site topples, $(m_i - 1) \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(y, i, j)$ grains of sand disappear, hence the word *leaky*. Starting from $s_0$, the configuration will eventually stabilize; in the case where all the $m_i$ are strictly greater than 1, there will be at most

$$\frac{N}{\min_{i \in \{1, \dots, p\}} (m_i - 1) \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(y, i, j)}$$
topplings. Besides, when several sites are unstable, the order in which they topple does not modify the stable configuration, hence the word Abelian.

We define a killed random walk associated with the sandpile model. This random walk is a Markov-additive process on $\mathbb{Z}^d \times \{1, \dots, p\}$. Markov-additive processes are Markov chains on the state-space $\mathbb{Z}^d \times E$ (here $E = \{1, \dots, p\}$), whose jumps are invariant along the $\mathbb{Z}^d$ component, that is, a Markov chain $(X_n)_{n \ge 0}$ such that

$$ \forall n \ge 0, \quad \forall i, j \in \{1, \dots, p\}, \quad \forall x, y \in \mathbb{Z}^d, \qquad \mathbb{P}(X_{n+1} = (y, j) | X_n = (x, i)) = \mathbb{P}(X_1 = (y-x, j) | X_0 = (0, i)). \quad (2) $$

Therefore, transitions of a Markov-additive process are fully described by the probabilities $\mu_{i,j}(x)$ to jump from $(0, i)$ to $(x, j)$ for $i, j \in \{1, \dots, p\}$ and $x \in \mathbb{Z}^d$. In the case where $p = 1$, the Markov-additive process is simply a sum of independent random variables with distribution $\mu = \mu_{1,1}$. Generalities on Markov-additive processes and their applications to queuing theory can be found in Chapter XI of [4].

**Model 2** (Killed Random Walk (KRW) associated with the sandpile model). For $x \in \mathbb{Z}^d$ and $i, j \in \{1, \dots, p\}$, let

$$ \mu_{i,j}(x) = \frac{c(x, i, j)}{m_i \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(y, i, j)}. \quad (3) $$

We define a killed Markov-additive process $(X_n)_{n \ge 0}$ on $(\mathbb{Z}^d \times \{1, \dots, p\}) \cup \{\mathbf{k}\}$ by the transitions

$$ \mathbb{P}((x, i) \to (y, j)) = \mu_{i,j}(y-x) \\ \mathbb{P}((x, i) \to \mathbf{k}) = 1 - \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mu_{i,j}(y) = 1 - \frac{1}{m_i}, $$

where the state $\mathbf{k}$ is the cemetery state of the walk after it is killed.

The measure $\sum_{j=1}^p \mu_{i,j}$ is a probability measure if and only if $m_i = 1$, that is, if there is no killing on color $i$. Otherwise, it is a sub-probability measure.

In order to apply previously known results on the Green function of Markov-additive processes, we will need the following assumptions on the Killed Random Walk (KRW).

**Assumption 1.**

1. The KRW is killed at least on one color, that is, there exists $i \in \{1, \dots, p\}$ such that $m_i > 1$.

2. The KRW is irreducible, that is:

$$ \forall (x, i, j) \in \mathbb{Z}^d \times \{1, \dots, p\}^2, \quad \exists n \in \mathbb{N}, \quad \exists (x_1, i_1), \dots, (x_{n-1}, i_{n-1}) \in \mathbb{Z}^d \times \{1, \dots, p\}, \\ \mu_{i,i_1}(x_1)\mu_{i_1,i_2}(x_2 - x_1) \cdots \mu_{i_{n-1},j}(x - x_{n-1}) > 0. $$

3. The KRW is aperiodic, that is:

$$ \forall (x, i) \in \mathbb{Z}^d \times \{1, \dots, p\}, \quad \text{gcd}\{n \ge 1 \mid \mathbb{P}_{(x,i)}(X_n = (x,i)) > 0\} = 1. $$

4. The KRW has finite exponential moments, that is:

$$ \forall i, j \in \{1, \dots, p\}, \quad \forall \gamma > 0, \quad \sum_{x \in \mathbb{Z}^d} e^{\gamma \|x\|} \mu_{i,j}(x) < +\infty. $$
# 3 Asymptotics of the Green function

## 3.1 The main result

**Definition 4.** The Green function of the KRW from Model 2 is

$$G((x, i), (y, j)) = \sum_{n=0}^{+\infty} \mathbb{P}_{(x,i)}(X_n = (y, j)).$$

As a consequence of (2), we have $G((x, i), (y, j)) = G((0, i), (y - x, j))$, so in the following, we will only consider values of the Green function of the form $G((0, i), (z, j))$, where the space coordinate of the first variable is 0. The main objective of this section is to prove the following asymptotic result. For $x \in \mathbb{R}^d$ we denote by $\|x\|$ its Euclidean norm and by $\hat{x} = \frac{x}{\|x\|} \in \mathbb{S}^{d-1}$ its direction.

**Proposition 5.** There exist continuous functions $\chi_{i,j}: \mathbb{S}^{d-1} \to \mathbb{R}_{>0}$ such that, when $\|x\| \to +\infty$,

$$G((0, i), (x, j)) \sim \chi_{i,j}(\hat{x}) e^{-\|x\|\Gamma^{-1}(\hat{x})} \|x\|^{-\frac{d-1}{2}},$$

where $\Gamma$ is the function introduced in Proposition 14.

An important feature of the previous proposition is that, although these asymptotics depend on the direction $\hat{x}$, they are uniform in $\hat{x}$.

First of all we recall some relevant and inspiring literature.

## 3.2 Literature on the Green function

To provide some context for the asymptotics given in Proposition 5, and to introduce natural ideas and tools, we briefly review the literature on the asymptotics of the Green function. Let us start with the following general asymptotic statement of the Green function $G(x, y) = G(0, y - x)$ for random walks with drift on $\mathbb{Z}^d$: as $\|x\| \to \infty$,

$$G(0, x) \sim \chi(\hat{x}) \exp(-\|x\|\nu(\hat{x})) \|x\|^{-\frac{d-1}{2}}. \quad (4)$$

The quantity $\chi(\hat{x})$ is the prefactor and $\nu(\hat{x}) \ge 0$ can be interpreted as the exponential decay to zero of the Green function. To illustrate the numerous results in the probabilistic and theoretical physics literature, let us mention five different interpretations of the exponential decay $\nu(\hat{x})$ and the underlying questions raised by these interpretations. In this paper we take inspiration from all of these perspectives. When the Green function has such asymptotics, it is said to have an Ornstein-Zernike decay (see e.g. [3]).

**A homeomorphism proved by Hennequin.** In [28, Thm 2.2], Ney and Spitzer derive the asymptotics (4). See [12, Thm 2.4] for a similar result when $\|x\| \to \infty$ in the direction of the drift. In particular, Ney and Spitzer compute in [28] the exponential decay of the Green function when $\|x\| \to \infty$ in terms of a homeomorphism between the sphere and level sets of the Laplace transform of the increment distribution found by Hennequin [18]. Since then, several alternative proofs of this asymptotic result (4) have been given, for example in [34, Thm 25.15].

In our more general context of walks on $\mathbb{Z}^d$ with $p$ colors or layers, i.e. $\mathbb{Z}^d \times \{1, \dots, p\}$, the asymptotics of the Green function were obtained independently in [13, Prop. 3.27] and [6, Thm 1.5]. Again, the exponential rate of decay to zero of the Green function is expressed in terms of a
homeomorphism between the unit sphere and the level set of a given spectral radius. In all the above references, the asymptotics are derived in the case of zero mass (i.e., non-killed random walks). As we shall see, the classical Doob transformation can transfer the results in the presence of a non-zero mass (random walks with killing).

**Second rate function in large deviation theory.** In [7, Cor. 5.7] it was shown that the exponential decay $\nu(\hat{x})$ in (4) can be reformulated in terms of the second rate function from large deviation theory. The latter is defined as follows: if $\Lambda$ denotes the classical rate deviation function associated to the random walk, then the second rate function is, for $x \in \mathbb{Z}^d$,

$$D(x) = \inf_{s>0} \frac{\Lambda(sx)}{s}.$$

This gives a different formulation and an interesting connection with large deviation theory. In the particular case of the simple random walk in dimension 2, the expression of the exponential decay in terms of the second large deviation function was independently obtained in [25, Prop. 8].

**Dual curves and Wulff crystals.** Another equivalent formulation of the exponential decay $\nu(\hat{x})$ in (4) is provided in the paper [7], see Theorem 1 there. It is proved that the function $D(x)$ above describes actually the dual set of the convex level set of the increment Laplace transform.

In the physics literature, the above dual set is often called a Wulff crystal or a Wulff shape. For instance, in the paper [25], the author studies the surface tension of the Ising model at criticality. In the setting of the killed two-dimensional simple random walk, he derives the logarithmic asymptotics of the two-point correlation function (the Green function), see [25, Thm 1]. In particular, he analyses the variation of the exponential decay of the Green function as a function of the killing rate, and studies the associated Wulff shapes. Wulff shapes actually arise in a variety of contexts in mathematical physics, such as Ising model and percolation theory, see [9].

**Inverse correlation length.** In the paper [3] the authors consider the Green function of a large class of “self-repulsive in average” models, including the killed random walk and other self-repulsive polymer models, including in particular the self-avoiding walk. Under certain assumptions they prove that the Green function has Ornstein-Zernike behaviour, which means that the asymptotics of the Green function are given by (4). In particular, the quantities $\chi$ and $\nu$ in (4) can be studied as a function of mass. The exponential decay $\nu(\hat{x})$ is interpreted as the inverse correlation length. Note that Wulff shapes also appear in [3], not to describe the rate of convergence to zero, but as a tool in the pathwise analysis of the Green function.

**Anisotropic norms.** In [26] the authors consider the simple random walk in any dimension and obtain four different asymptotic regimes for the associated Green functions, depending on how (simultaneously) the point $x$ goes to infinity and the mass tends to 0 or to $\infty$. One of the regimes corresponds to Ornstein-Zernike behaviour; the authors show that other relevant regimes appear. In [26] the quantity $x \mapsto \nu(\hat{x})\|x\|$ is interpreted as an anisotropic norm. A particularly interesting feature of the paper is the analysis of this norm as the mass goes to zero and infinity. In the zero-mass asymptotics the norm recovers an isotropy; on the contrary, in the infinite mass regime it converges to a “non-smooth” norm.

In dimension 2, similar regimes are obtained in [2, Thm 1.4] for the Green function of the killed simple random walk.
**Lattice Green functions and differential equations.** Although less related to asymptotic estimates, let us mention here a line of research started in the seventies, see e.g. [19], the main objective of which is to compute certain lattice Green functions. Indeed, as shown in [19], it is possible, for example, to compute the simple cubic lattice Green function in terms of certain special functions called Heun functions. More generally, there are several works which compute Green functions in closed form, see [15, 16]. These expressions can be hypergeometric functions, elliptic functions, Heun functions, classical integrals, diagonals of rational functions, etc. Sometimes they are expressed using their differential operators, and the need for computer algebra techniques can be particularly useful [17]. This line of research is still active, see e.g. [10]. Although this has not yet been done, the use of analytic combinatorics in several variables techniques [31] is promising to deduce from these expressions their asymptotic behaviour.

### 3.3 Proof of the main asymptotics result on the Green function

In this section, we derive the asymptotics of the Green function $G$ stated in Proposition 5 from the asymptotics of a non-killed walk studied in [6, 13]. To do so, we will perform a Doob transform to turn the killed random walk into a non-killed one with a drift.

The (real) Laplace transform of the walk is the function $L\mu: \mathbb{R}^d \to \mathcal{M}_p(\mathbb{R})$ defined by:

$$ \forall t \in \mathbb{R}^d, \quad (L\mu)_{i,j}(t) := \sum_{x \in \mathbb{Z}^d} e^{t \cdot x} \mu_{i,j}(x). $$

Assumption 1 on the finite exponential moments ensures that it is well defined on the whole of $\mathbb{R}^d$.

**Lemma 6** (Lem. 2.10 in [6]). Let $t \in \mathbb{R}^d$. We assume that 1 is an eigenvalue of the Laplace transform $L\mu(t)$ with an eigenvector $\varphi_t$ whose all entries are positive. We define the measures $(\mu_t)_{i,j}$ by

$$ \forall i, j \in \{1, \dots, p\}, \quad \forall x \in \mathbb{Z}^d, \quad (\mu_t)_{i,j}(x) = \frac{(\varphi_t)_j}{(\varphi_t)_i} e^{t \cdot x} \mu_{i,j}(x). $$

Then for all $i \in \{1, \dots, p\}$, $\sum_{j=1}^p (\mu_t)_{i,j}$ is a probability measure on $\mathbb{Z}^d$. This means that the Markov-additive process with transitions $\mu_t$ is a non-killed process.

The proof of Lemma 6 is the same as that of [6, Lem. 2.10]. We briefly recall the main idea below to keep our paper self-contained. In the work [32], the Doob transform technique is applied to statistical mechanics by relating two models, namely random rooted spanning forests and random spanning trees, and providing various applications.

*Proof.* The positivity of $\varphi_t$ ensures that $\sum_{j=1}^p (\mu_t)_{i,j}$ is non-negative, and

$$ \sum_{x \in \mathbb{Z}^d} \sum_{j=1}^p (\mu_t)_{i,j}(x) = \sum_{j=1}^p \sum_{x \in \mathbb{Z}^d} \frac{(\varphi_t)_j}{(\varphi_t)_i} e^{t \cdot x} \mu_{i,j}(x) = \frac{1}{(\varphi_t)_i} (L\mu(t)\varphi_t)_i = \frac{1}{(\varphi_t)_i} (\varphi_t)_i = 1. \quad \square $$

Assumption 1 on irreducibility and aperiodicity ensures that $L\mu(t)$ is a primitive matrix for every $t \in \mathbb{R}^d$. Therefore, if the spectral radius of $L\mu(t)$ is 1, then the Perron-Frobenius theorem ensures that there exists an eigenvector $\varphi_t$ like in Lemma 6, which is unique up to multiplication. This leads to the following definitions.

**Definition 7.** Let $t \in \mathbb{R}^d$. The spectral radius of $L\mu(t)$ is denoted $\rho(t)$. We define, in $\mathbb{R}^d$, $\mathcal{T} := \rho^{-1}([0, 1])$ and $\partial \mathcal{T} := \rho^{-1}(\{1\})$.
**Definition 8.** Let $t \in \partial T$ and $\varphi_t$ as in Lemma 6. The Doob transform of parameter $t$ of the process $(X_n)_{n \ge 0}$ is the Markov-additive process $(X_n^t)_{n \ge 0}$ whose jump matrix is $\mu_t$ defined by

$$ \forall i, j \in \{1, \dots, p\}, \quad \forall x \in \mathbb{Z}^d, \quad (\mu_t)_{i,j}(x) = \frac{(\varphi_t)_j}{(\varphi_t)_i} e^{t \cdot x} \mu_{i,j}(x). $$

It is a non-killed process because Lemma 6 ensures that every $\sum_{j=1}^p (\mu_t)_{i,j}$ is a probability measure, and no longer a sub-probability measure.

**Lemma 9.** We have $\rho(0) < 1$, that is, $0 \in \mathcal{T}$.

In the proof below, inequalities between vectors are meant coordinate by coordinate.

*Proof.* Let $M = L\mu(0) = (\sum_{x \in \mathbb{Z}^d} \mu_{i,j}(x))_{1 \le i,j \le p}$. It is a sub-stochastic matrix, so $\rho(M) \le 1$, which leaves us to prove that $\rho(M) \ne 1$. At least one of the rows of $M$ has a sum strictly smaller than 1, because our running assumption implies that there exists $i \in \{1, \dots, p\}$ such that $m_i > 1$. If we denote by $\mathbb{1}$ the vector of $\mathbb{R}^p$ with all coordinates equal to 1, it means that $M\mathbb{1} \le \mathbb{1}$ and $M\mathbb{1} \ne \mathbb{1}$. Let $x \in \mathbb{R}^d$ be a left eigenvector of $M$ associated with the eigenvalue $\rho(M)$, that is $x^\top M = \rho(M)x^\top$. The Perron-Frobenius theorem applied to $M^\top$ ensures that $x > 0$. If $\rho(M) = 1$, we would have

$$ x^\top (M \mathbb{1} - \mathbb{1}) = x^\top \mathbb{1} - x^\top \mathbb{1} = 0. $$

But $M\mathbb{1} - \mathbb{1} \le 0$ with at least one non-zero coefficient and $x > 0$, so $x^\top (M\mathbb{1} - \mathbb{1}) < 0$, a contradiction. $\square$

The Green function $G$ of the KRW and the Green function $G_t$ of the Doob transform of parameter $t$ are related as follows.

**Proposition 10** (Prop. 2.13 in [6]). Let $t \in \partial T$, $x \in \mathbb{Z}^d$ and $i, j \in \{1, \dots, p\}$. The Green function of the Doob transform is given by

$$ G_t((0,i), (x,j)) = \frac{(\varphi_t)_j}{(\varphi_t)_i} e^{t \cdot x} G((0,i), (x,j)), $$

where $\varphi_t$ is the Perron-Frobenius eigenvector introduced in Lemma 6.

Proposition 10 is stated (and proved) in [6] as Prop. 2.13. The main argument is briefly recalled below.

*Proof.* Let $x \in \mathbb{Z}^d$, $i, j \in \{1, \dots, p\}$ and $n \ge 0$. We have

$$ \mathbb{P}_{(0,i)}(X_n = (x,j)) = \sum_{(0,i)=(x_0,i_0) \to (x_1,i_1) \to \dots \to (x_n,i_n) = (x,j)} \mu_{i,i_1}(x_1) \mu_{i_1,i_2}(x_2 - x_1) \cdots \mu_{i_{n-1},j}(x - x_{n-1}), $$

where the sum is taken on all paths from $(0, i)$ to $(x, j)$ of length $n$ in the graph associated with the Markov chain. Similarly,

$$ \begin{aligned} \mathbb{P}_{(0,i)}(X_n^t = (x,j)) &= \sum_{(0,i) \to (x_1,i_1) \to \dots \to (x_j)} \frac{(\varphi_t)_{i_1}}{(\varphi_t)_i} e^{t \cdot x_1} \mu_{i,i_1}(x_1) \cdots \frac{(\varphi_t)_j}{(\varphi_t)_{i_{n-1}}} e^{t \cdot (x - x_{n-1})} \mu_{i_{n-1},j}(x - x_{n-1}) \\ &= \sum_{(0,i) \to (x_1,i_1) \to \dots \to (x_j)} \frac{(\varphi_t)_j}{(\varphi_t)_i} e^{t \cdot x} \mu_{i,i_1}(x_1) \mu_{i_1,i_2}(x_2 - x_1) \cdots \mu_{i_{n-1},j}(x - x_{n-1}) \\ &= \frac{(\varphi_t)_j}{(\varphi_t)_i} e^{t \cdot x} \mathbb{P}_{(0,i)}(X_n = (x,j)). \end{aligned} $$

Summing over $n \ge 0$ gives the result announced in Proposition 10. $\square$
In the following results we study the spectral radius $\rho$ in order to obtain useful properties of the set $\mathcal{T}$ introduced in Definition 7.

**Proposition 11** (Prop. 2.14 in [6]). The function $\rho$ is:

* convex (and even logarithmically convex);

* norm-coercive, i.e., $\rho(t) \xrightarrow{\|t\| \to +\infty} +\infty$.

We recall the proof given in [6], as we think it gives an interesting insight into the model.

*Proof.* The proof of the convexity relies on the result of [22]. If $\mu_{i,j}(x) > 0$, then $t \mapsto e^{t \cdot x} \mu_{i,j}(x)$ is logarithmically convex. Therefore, the entries of $L\mu$ are 0 if $\mu_{i,j} = 0$, or logarithmically convex functions, as sums of logarithmically convex functions. Besides, by irreducibility and aperiodicity, for $n$ large enough, all entries of $(L\mu)^n$ are positive. If all entries are logarithmically convex, then the result of [22] ensures that $\rho$ is logarithmically convex. To deal with the 0 entries, we take a closer look at the proof of [22] and notice that if $\mathrm{Tr}((L\mu)^n) > 0$ for $n$ large enough, which is the case, then the result holds.

We move to the proof of norm-coercivity. When a function is convex, its norm-coercivity is equivalent to coercivity in every direction. Hence we only need to prove that for every $u \in \mathbb{S}^{d-1}$, $\rho(ru) \to +\infty$ as $r \to +\infty$. According to the Perron-Frobenius theorem, see for example [33, Cor. 1],

$$ \rho(ru)^n \geq \min_{i \in \{1, \dots, p\}} \sum_{j=1}^{p} (L\mu(ru)^n)_{i,j}, $$

so it is enough to prove that for every $u \in \mathbb{S}^{d-1}$, there exists $n \in \mathbb{N}$ such that for every $i, j \in \{1, \dots, p\}$, $(L\mu(ru)^n)_{i,j} \to +\infty$ as $r \to +\infty$. Let $u \in \mathbb{S}^{d-1}$ and $i, j \in \{1, \dots, p\}$. Let $x \in \mathbb{Z}^d$ such that $x \cdot u > 0$. Then as a consequence of aperiodicity and irreducibility, there exists $n_{i,j} \in \mathbb{N}$ such that for every $n \geq n_{i,j}$, $\mathbb{P}_{(0,i)}(X_n = (x, j)) > 0$. Let $n = \max_{i,j \in \{1, \dots, p\}} n_{i,j}$. Then

$$ (L\mu(ru)^n)_{i,j} \geq e^{rx \cdot u} \mathbb{P}_{(0,i)}(X_n = (x, j)) \xrightarrow{r \to +\infty} +\infty, $$

which proves the norm-coercivity.

<figure><img src="image_3.png" /></figure>

**Proposition 12.** For every $t \in \partial \mathcal{T}$, $\nabla \rho(t) \neq 0$.

*Proof.* The function $\rho$ is convex, so if $\nabla \rho(t) = 0$ for some $t \in \partial \mathcal{T}$, then $\rho$ has a global minimum in $t$. But $\rho(t) = 1$ and $\rho(0) < 1$, hence a contradiction.

<figure><img src="image_4.png" /></figure>

The gradient of $\rho(t)$ is the right notion for the drift of the Doob transform of parameter $t$, see [6, Def. 1.3, Prop. 2.16]. Therefore, Proposition 12 means that the Doob transform of parameter $t \in \partial \mathcal{T}$ is non-centered, allowing to use previously known results on non-killed walks with a drift.

**Proposition 13** (Prop. 2.16 in [6]). For every $t \in \partial \mathcal{T}$, the Hessian $\mathbf{H}_{\rho}(t)$ of $\rho$ is positive-definite.

The next result is crucial to get the exponential decay of the Green function in Proposition 5. With the interpretation of the gradient of $\rho$ as the drift of the process, it means that the Doob transform allows us to change the drift in any direction of $\mathbb{R}^d$.
**Proposition 14** (Thm 2.17 in [6]). The function

$$ \Gamma: \partial \mathcal{T} \longrightarrow \mathbb{S}^{d-1} \\ t \longmapsto \frac{1}{\|\nabla \rho(t)\|} \nabla \rho(t) $$

is a homeomorphism between $\partial \mathcal{T}$ and the sphere $\mathbb{S}^{d-1}$ of $\mathbb{R}^d$.

**Proposition 15.** Let $u \in \mathbb{S}^{d-1}$. Then $\Gamma^{-1}(u) \cdot u > 0$.

*Proof.* It comes down to proving that for every $t \in \partial \mathcal{T}$, $t \cdot \nabla \rho(t) > 0$. By convexity of $\rho$, for every $t \in \partial \mathcal{T}$, $\mathcal{T}$ is included in the affine half-space

$$ H_t = \{x \in \mathbb{R}^d \mid (x - t) \cdot \nabla \rho(t) \le 0\}. $$

Besides, since the Hessian is positive-definite by Proposition 13, $\partial H_t \cap \mathcal{T} = \{t\}$. But $t \neq 0$ because $0 \in \mathcal{T} \setminus \partial \mathcal{T}$, so $0 \notin \partial H_t \cap \mathcal{T}$, hence $0 \in H_t \setminus \partial H_t$, i.e., $t \cdot \nabla \rho(t) > 0$. $\square$

The previous proposition ensures that the decay of Proposition 5 is actually exponential in every direction. This is not true for non-killed walks studied in [6], where the decay in the direction of the drift is only polynomial.

Now we can complete the proof of Proposition 5.

*Proof of Proposition 5.* Let us fix $t_0 \in \partial \mathcal{T}$. We remind that $G_{t_0}$ is the Green function of a non-killed walk, with a non-zero drift. Therefore, under Assumption 1, [6, Thm 1.5] shows that there exist continuous functions $\tilde{\chi}_{i,j}: \mathbb{S}^{d-1} \to \mathbb{R}_{>0}$ and $c: \mathbb{S}^{d-1} \to \mathbb{R}^d$ such that

$$ G_{t_0}((0, i), (x, j)) \sim \tilde{\chi}_{i,j} \left( \frac{x}{\|x\|} \right) e^{-\|x\|c\left(\frac{x}{\|x\|}\right) \cdot \frac{x}{\|x\|}} \|x\|^{-\frac{d-1}{2}}. $$

Hence using Proposition 10

$$ G((0, i), (x, j)) = \frac{(\varphi_{t_0})_i}{(\varphi_{t_0})_j} e^{-t_0 \cdot x} G_{t_0}((0, i), (x, j)) \sim \frac{(\varphi_{t_0})_i \tilde{\chi}_{i,j}}{(\varphi_{t_0})_j} \left( \frac{x}{\|x\|} \right) e^{-\|x\|(t_0+c(\frac{x}{\|x\|}) \cdot \frac{x}{\|x\|}) \|x\|^{-\frac{d-1}{2}}}. \quad (5) $$

It remains to identify the exponential decay constant $t_0 + c(\frac{x}{\|x\|})$ in (5). To do so we fix a direction $u \in \mathbb{S}^{d-1}$ and work with positive multiples of $u$, that is $x = \|x\|u$. We choose $t_u = \Gamma^{-1}(u)$, so that

$$ G((0, i), (\|x\|u, j)) = \frac{(\varphi_t)_i}{(\varphi_t)_j} e^{-\|x\|\Gamma^{-1}(u) \cdot u} G_{t_u}((0, i), (\|x\|u, j)). \quad (6) $$

We apply once again [6, Thm 1.5] to get the asymptotics of $G_{t_u}((0, i), (\|x\|u, j))$ as $\|x\| \to +\infty$. This time, our choice of $t_u$ makes the drift of the Doob transform in the direction of $u$. Therefore, there is no exponential decay in $G_{t_u}((0, i), (\|x\|u, j))$, which only decreases polynomially. Identifying the exponential decay in (5) and (6), we get $t_0 + c(\frac{x}{\|x\|}) = \Gamma^{-1}(\frac{x}{\|x\|})$, and thus (5) gives the result. $\square$

Let us make two remarks on Proposition 5. First, in the proof, we could not directly choose $t = t_u$ instead of a fixed $t_0$, because in that case, the Green function $G_{t_u}$ changes with the direction of $x$. Therefore, [6, Thm 1.5] would only yield the asymptotics in a fixed direction $u$, while we want the asymptotics to be uniform in $u$.

Moreover, using Proposition 15 and a continuity-compactness argument, we have

$$ \min_{u \in \mathbb{S}^{d-1}} \Gamma^{-1}(u) \cdot u > 0. \quad (7) $$

Denoting this minimum by $\gamma$, we see that the Green function decreases faster than $e^{-\gamma \|x\|}$.
## 4 Odometer functionand thresholds of the Green function

The main goal of this section is to provide a link between the Green function of the killed random walk and the shape of the sandpile. In Proposition 21, generalizing the idea of [2], we obtain two thresholds for the Green function $G((0, i_0), (x, i))$: one above which the point $(x, i)$ is inside the sandpile and one under which it is outside the sandpile. To do so, we introduce the odometer function, a classical function in sandpiles, and related operators.

**Definition 16.** We define the following operator $T: \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\}) \to \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$, which can be interpreted as a massive Laplacian:

$$ (Tv)(x, i) = \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} \frac{c(x - y, j, i)}{m_j \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(z, j, k)} v(y, j) - v(x, i) \quad (8) $$

$$ = \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}((y,j) \to (x,i))v(y,j) - v(x,i). \quad (9) $$

If we see $T$ as a square matrix with rows and columns indexed by $\mathbb{Z}^d \times \{1, \dots, p\}$, we have

$$ \begin{aligned} T((x, i), (y, j)) &= T1_{(y,j)}(x, i) = \frac{c(x - y, j, i)}{m_j \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(z, j, k)} - 1_{(x,i)=(y,j)} \\ &= \mathbb{P}((y,j) \to (x,i)) - 1_{(x,i)=(y,j)}. \end{aligned} $$

Let us check that the formula (8) for $Tv$ does indeed define an operator $\ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\}) \to \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$. One has

$$ \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} |\mathbb{P}((y,j) \to (x,i))v(y,j)| \le \|v\|_{\infty} \sum_{j=1}^{p} \sum_{y \in \mathbb{Z}^d} \mathbb{P}((0,j) \to (x-y,i)) \le p \|v\|_{\infty}, $$

therefore $Tv$ is well defined in $\ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$ and $\|Tv\|_{\infty} \le (p+1)\|v\|_{\infty}$.

**Definition 17.** The transpose of the Green function defines an operator $G^\top: \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\}) \to \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$ by

$$ G^\top v(x, i) = \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} G((y,j), (x,i))v(y,j). $$

Let us check that $G^\top$ does indeed define an operator $\ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\}) \to \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$:

$$ \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} |G((y,j), (x,i))v(y,j)| \le \|v\|_{\infty} \sum_{j=1}^{p} \sum_{z \in \mathbb{Z}^d} G((0,j), (z,i)). $$

According to Proposition 5 and to (7), there exist constants $\gamma, C > 0$ such that for every $(z, i, j)$, $G((0, j), (z, i)) \le Ce^{-\gamma\|z\|}$. Therefore,

$$ \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} |G((y,j), (x,i))v(y,j)| \le Cp\|v\|_{\infty} \sum_{z \in \mathbb{Z}^d} e^{-\gamma\|z\|} < +\infty. $$

This shows that $G^\top v \in \ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$ and $\|G^\top v\|_{\infty} \le Cp\|v\|_{\infty} \sum_{z \in \mathbb{Z}^d} e^{-\gamma\|z\|}$.

The following proposition expresses $T$ as an inverse of the Green function $G$.
**Proposition 18.** The products $G^\top T$ and $TG^\top$ are well defined and equal to the negative identity
$$-I_{\mathbb{Z}^d \times \{1, \dots, p\}}.$$

In particular, $T^{-1} = -G^\top$ shows that $T^{-1}$ has negative coefficients. Therefore, if $v \ge 0$, then
$T^{-1}v \le 0$.

*Proof.* The products are well defined because the columns of these matrices are in $\ell^\infty(\mathbb{Z}^d \times \{1, \dots, p\})$.
Indeed, for $T$ we only need to bound probabilities by 1, and for $G^\top$, we use the bound $G((0, j), (z, i)) \le Ce^{-\gamma \|z\|} \le C$ for constants $\gamma, C > 0$, which is a consequence of Proposition 5.

By definition of $T$ and $G$,

$$ (TG^\top)_{(x,i),(y,j)} = \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}((z,k) \to (x,i)) G((y,j), (z,k)) - G((y,j), (x,i)). \quad (10) $$

Besides, conditioning according to $X_{n-1}$, we get:

$$
\begin{align*}
G((y,j), (x,i)) &= \mathbb{1}_{(y,j)=(x,i)} + \sum_{n=1}^{+\infty} \mathbb{P}_{(y,j)}(X_n = (x,i)) \\
&= \mathbb{1}_{(y,j)=(x,i)} + \sum_{n=1}^{+\infty} \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}_{(y,j)}(X_n = (x,i) \mid X_{n-1} = (z,k)) \mathbb{P}_{(y,j)}(X_{n-1} = (z,k)) \\
&= \mathbb{1}_{(y,j)=(x,i)} + \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}((z,k) \to (x,i)) \sum_{m=0}^{+\infty} \mathbb{P}_{(y,j)}(X_m = (z,k)) \\
&= \mathbb{1}_{(y,j)=(x,i)} + \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}((z,k) \to (x,i)) G((y,j), (z,k)).
\end{align*}
$$

Putting it into (10), we get the expected $(TG^\top)_{(x,i),(y,j)} = -\mathbb{1}_{(y,j)=(x,i)}$.

Similarly,

$$ (G^\top T)_{(x,i),(y,j)} = \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} G((z,k), (x,i)) \mathbb{P}((y,j) \to (z,k)) - G((y,j), (x,i)). \quad (11) $$

Besides, conditioning according to $X_1$, we get:

$$
\begin{align*}
G((y,j), (x,i)) &= \mathbb{1}_{(y,j)=(x,i)} + \sum_{n=1}^{+\infty} \mathbb{P}_{(y,j)}(X_n = (x,i)) \\
&= \mathbb{1}_{(y,j)=(x,i)} + \sum_{n=1}^{+\infty} \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}_{(y,j)}(X_n = (x,i) \mid X_1 = (z,k)) \mathbb{P}_{(y,j)}(X_1 = (z,k)) \\
&= \mathbb{1}_{(y,j)=(x,i)} + \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}((y,j) \to (z,k)) \sum_{m=0}^{+\infty} \mathbb{P}_{(z,k)}(X_m = (x,i)) \\
&= \mathbb{1}_{(y,j)=(x,i)} + \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} \mathbb{P}((y,j) \to (z,k)) G((z,k), (x,i)).
\end{align*}
$$

Putting it into (11), we get the expected $(G^\top T)_{(x,i),(y,j)} = -\mathbb{1}_{(y,j)=(x,i)}$. $\square$
**Definition 19 ([2]).** The odometer function is the function $u: \mathbb{Z}^d \times \{1, \dots, p\} \to \mathbb{R}$ defined by
$$u(x, i) = \text{total mass of sand emitted from } (x, i) \text{ before stabilization.}$$

We say that a site $(x, i)$ is in the shape of the final configuration if $u(x, i) > 0$.

The odometer function counts the grains sent to the neighbours and those that disappeared because of the leakiness. Note that it implicitly depends on the initial configuration.

One might notice that in the shape of the final configuration, we do not count the sites that received grains of sand but did not emit any. Note, that those sites are at distance at most 1 (in the graph associated with the Markov chain) of sites that were counted, so in the case when the support of $\mu$ is finite, adding them will not modify the limit shape of the sandpile as $N$ goes to infinity. Furthermore, it might happen that sites received sand during the stabilization process but sent it all and have an amount of sand equal to zero in the final configuration. With our definition, those sites have a positive odometer and thus belong to the final shape.

The following proposition links the operator $T$ and the odometer function to the final configuration of the sandpile. It is a generalization of [2, Eq. (3.2)] to our context.

**Proposition 20.** Let $s_0$ be any initial configuration with a finite number of grains of sand, and denote by $f$ the final configuration obtained after stabilization. Then

$$Tu = f - s_0.$$

In particular, if we start the sandpile with $N$ grains at $(0, i_0)$, then

$$Tu = f - N\delta_{(0,i_0)}.$$

*Proof.* According to (8),

$$(Tu)(x, i) = \sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} \frac{c(x - y, j, i)}{m_j \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(z, j, k)} u(y, j) - u(x, i).$$

The definition 1 of the LASM shows that when the site $(y, j)$ topples, the proportion of sand leaving $(y, j)$ that is sent to $(x, i)$ is

$$\frac{c(x - y, j, i)}{m_j \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(z, j, k)}.$$

Therefore, the total mass received by $(x, i)$ is

$$\sum_{(y,j) \in \mathbb{Z}^d \times \{1, \dots, p\}} \frac{c(x - y, j, i)}{m_j \sum_{(z,k) \in \mathbb{Z}^d \times \{1, \dots, p\}} c(z, j, k)} u(y, j).$$

Hence

$$(Tu)(x, i) = \text{total mass received by } (x, i) - \text{total mass emitted from } (x, i) \\ = f(x, i) - s_0(x, i). \quad \square$$

The result below gives the connection between the level sets of the Green function and the final configuration. It is an adaptation of [2, Lem. 3.3] to our general setting.

**Proposition 21.** There exist constants $\alpha, \beta > 0$ such that:
• if $G((0, i_0), (x, i)) > \frac{\alpha}{N}$, then $(x, i)$ is in the shape of the final configuration;

• if $G((0, i_0), (x, i)) < \frac{\beta}{N}$, then $(x, i)$ is not in the shape of the final configuration.

*Proof.* By Propositions 18 and 20,

$$T\left(-G((0, i_0), \cdot) + \frac{1}{N}u\right) = \frac{1}{N}f.$$

The finite family $\left(m_i \sum_{(y,j)} c(y, i, j)\right)_{i \in \{1, \dots, p\}}$ is bounded by a constant $a$, and in the final configuration, no site has more than $a$ grains of sand. Therefore,

$$0 \le T\left(-G((0, i_0), \cdot) + \frac{1}{N}u\right) \le \frac{a}{N}1$$

where $1$ is the constant function equal to 1. Multiplying by $T^{-1} = -G^T$, which has negative coefficients, we get

$$0 \ge \frac{1}{N}u(x, i) - G((0, i_0), (x, i)) \ge -\frac{a}{N}G^T1(x, i).$$

Moreover,

$$G^T1(x, i) = \sum_{(y,j)} G((y,j), (x,i)) = \sum_{j=1}^{p} \sum_{z \in \mathbb{Z}^d} G((0,j), (z,i)),$$

and this sum is finite because of the exponential decay of the Green function. Moreover, it does not depend on $x$ and there is a finite number of $i \in \{1, \dots, p\}$, so it is bounded uniformly in $(x, i)$, hence the existence of a constant $\alpha$ such that:

$$0 \ge \frac{1}{N}u(x, i) - G((0, i_0), (x, i)) \ge -\frac{\alpha}{N} \qquad (12)$$

If $(x, i)$ is not in the shape of the final configuration, then $u(x, i) = 0$ and from (12), we get

$$G((0, i_0), (x, i)) \le \frac{\alpha}{N},$$

whereas if it is in the shape, then $(x, i)$ has emitted sand at least once, thus $u(x, i) \ge m_i \sum_{(y,j)} c(y, i, j)$ and (12) leads to

$$G((0, i_0), (x, i)) \ge \frac{1}{N}u(x, i) \ge \frac{1}{N} \min_{i \in \{1, \dots, p\}} m_i \sum_{(y,j)} c(y, i, j),$$

$$\text{so we can choose the constant } \beta := \min_{i \in \{1, \dots, p\}} m_i \sum_{(y,j)} c(y, i, j). \quad \square$$

# 5 Limit shape when the number of grains tends to infinity

In this section, we prove Theorem 1, which is our main result. The existence of a limit shape for the stable configuration as the number of grains of sand $N$ goes to infinity, after normalization by $\log N$, is proved. We give an explicit formula for the limit shape that uses the homeomorphism $\Gamma$ of Proposition 14 and also describe it as the dual curve of the set $\partial T$ introduced in Definition 7.

The general idea is to use the exponential decay of the Green function from Proposition 5 to show that the two level sets of Proposition 21 are close enough, so that after rescaling by $\log N$, they converge to the same curve.
## 5.1 Notions of convergence for sets

We first define the notion of convergence of sets that we use. Given that the exponential decay of the Green function depends on the direction, we use convergence in spherical coordinates.

**Definition 22.** Let $(E_N)_{N \ge 0}$ be a sequence of subsets of $\mathbb{R}^d$ and $C$ a curve which can be described in spherical coordinates as

$$C = \{r_u u \mid u \in \mathbb{S}^{d-1}\},$$

with $r: u \in \mathbb{S}^{d-1} \mapsto r_u \in \mathbb{R}_+$. We say that the limit shape of the sets $E_N$ is delimited by $C$ if, for every $u \in \mathbb{S}^{d-1}$, there exist two sequences $(r_{u,N})_{N \ge 0}$ and $(R_{u,N})_{N \ge 0}$ such that, for every $N \in \mathbb{N}$,

$$\{ru \mid u \in \mathbb{S}^{d-1}, r \le r_{u,N}\} \subset E_N \subset \{Ru \mid u \in \mathbb{S}^{d-1}, R \le R_{u,N}\}$$

and for every $u \in \mathbb{S}^{d-1}$,

$$r_{u,N} \xrightarrow[N \to +\infty]{} r_u \quad \text{and} \quad R_{u,N} \xrightarrow[N \to +\infty]{} r_u.$$

If these two limits are uniform in $u \in \mathbb{S}^{d-1}$, that is

$$\sup_{u \in \mathbb{S}^{d-1}} |r_{u,N} - r_u| \xrightarrow[N \to +\infty]{} 0 \quad \text{and} \quad \sup_{u \in \mathbb{S}^{d-1}} |R_{u,N} - R_u| \xrightarrow[N \to +\infty]{} 0,$$

we say that the convergence is uniform.

We now relate the convergence of sets from Definition 22 to the classical Hausdorff distance convergence, which we first recall.

**Definition 23.** Let $X, Y \subset \mathbb{R}^d$ be non-empty. The Hausdorff distance between $X$ and $Y$ is

$$d_H(X, Y) = \max \left( \sup_{x \in X} d(x, Y), \sup_{y \in Y} d(y, X) \right),$$

where $d(a, B) = \inf_{b \in B} \|a - b\|$.

The Hausdorff distance restricted to closed, bounded, non-empty subsets of $\mathbb{R}^d$ defines a distance in the sense of metric spaces.

**Proposition 24.**

1. If a sequence of sets $(E_N)_{N \ge 0}$ converges uniformly to

$$C = \{r_u u \mid u \in \mathbb{S}^{d-1}\}$$

in the sense of Definition 22, then we have $d_H(E_N, E) \xrightarrow[N \to +\infty]{} 0$ where

$$E = \{ru \mid u \in \mathbb{S}^{d-1}, 0 \le r \le r_u\}$$

is the set delimited by the curve $C$.

2. The conclusion is not necessarily true if the convergence to the limit shape is not uniform.

3. The converse of the assertion 1 is false; Hausdorff convergence does not even imply convergence in the non-uniform sense of Definition 22.
<figure><img src="image_5.png" /><figcaption>Figure 2: Sets that converge to the unit ball for the Hausdorff distance, but not in the sense of Definition 22.</figcaption></figure>

*Proof.* Let us assume the uniform convergence with the notations of Definition 22. Let $x = ru \in E_N$ with $r \ge 0$ and $u \in \mathbb{S}^{d-1}$. If $r \le r_u$, then $d(x, E) = 0$; whereas if $r \ge r_u$, then $r \le R_{u,N}$, so

$$d(x, E) \le d(ru, r_u u) = |r - r_u| \le R_{u,N} - r_u.$$

Therefore,

$$\sup_{x \in E_N} d(x, E) \le \max \left( 0, \sup_{u \in \mathbb{S}^{d-1}} (R_{u,N} - r_u) \right) \le \sup_{u \in \mathbb{S}^{d-1}} |R_{u,N} - r_u|.$$

Similarly, let $y = ru \in E$, which means $u \in \mathbb{S}^{d-1}$ and $r \le r_u$. If $r \le r_{u,N}$, then $y \in E_N$ so $d(y, E_N) = 0$; whereas if $r > r_{u,N}$, then $d(y, E_N) \le d(y, r_{u,N}u) = r - r_{u,N} \le r_u - r_{u,N}$. Therefore,

$$\sup_{y \in E_N} d(y, E_N) \le \max \left( 0, \sup_{u \in \mathbb{S}^{d-1}} (r_u - r_{u,N}) \right) \le \sup_{u \in \mathbb{S}^{d-1}} |r_u - r_{u,N}|.$$

In conclusion,

$$d_H(E_N, E) \le \max \left( \sup_{u \in \mathbb{S}^{d-1}} |R_{u,N} - r_u|, \sup_{u \in \mathbb{S}^{d-1}} |r_u - r_{u,N}| \right) \xrightarrow{N \to +\infty} 0.$$

We now give a counterexample in the non-uniform case. Let $(u_N)_{N \ge 0}$ be an injective sequence of $\mathbb{S}^{d-1}$ and

$$E_N = \bar{B}(0, 1) \cup [0, 2u_N]$$

the closed unit ball pricked with a toothpick at $u_N$. Then, with notations from Definition 22, we can choose $R_{u,N} = r_{u,N} = 1$ if $u \ne u_N$ and $R_{u_N,N} = r_{u_N,N} = 2$. The sequence $(u_N)_{N \ge 0}$ being injective, for fixed $u$, we have $R_{u,N} = r_{u,N} = 1$ for $N$ large enough. Therefore, $(E_N)_{N \ge 0}$ converges (non-uniformly) to the unit ball $\bar{B}(0, 1)$ in the sense of Definition 22. However, $d_H(E_N, \bar{B}(0, 1)) = 1$ for every $N \in \mathbb{N}$.

We finally give an example of convergence with respect to the Hausdorff distance where there is no convergence in the sense of Definition 22. Let $E_N$ be a circular sector as in Figure 2 with $\theta_N = \pi/N$. Then it is easy to see that $d_H(E_N, \bar{B}(0, 1)) \le \pi/2N$, so $E_N$ converges to the unit ball for the Hausdorff distance, but for every $N \in \mathbb{N}$, the radius in direction $(1, 0)$ is 0 which does not converge to 1. $\square$

## 5.2 Convergence of the sandpile to a limit shape

Let us fix $i_0, j \in \{1, \dots, p\}$. In order to use some tools from classical analysis, we first extend continuously the function $x \mapsto G((0, i_0), (x, j))$ from $\mathbb{Z}^d$ to $\mathbb{R}^d$, by saying that the value at $x \in \mathbb{R}^d$ is
a convex combination of the values at neighboring lattice points. To do so, we define $\theta : [0, +\infty) \to [0, 1]$ by

$$ \theta(s) = \begin{cases} 1-s & \text{if } s \in [0, 1], \\ 0 & \text{if } s > 1, \end{cases} $$

and

$$ G((0, i_0), (x, j)) = \frac{\sum_{y \in \mathbb{Z}^d} \theta(\|x - y\|_{\infty}) G((0, i_0), (y, j))}{\sum_{y \in \mathbb{Z}^d} \theta(\|x - y\|_{\infty})} $$

We will denote $u \in \mathbb{S}^{d-1}$ and

$$ g_u : r \mapsto G((0, i_0), (ru, j)) \quad (13) $$

the Green function in direction $u$. Note that the function $(u, r) \in \mathbb{S}^{d-1} \times [0, +\infty) \mapsto g_u(r)$ is continuous.

**Definition 25.** We define radii $r_{N,u}$ and $R_{N,u}$ by:

$$ r_{N,u} = \inf \left\{ r > 0 \mid g_u(r) \le \frac{\alpha}{N} \right\} \quad \text{and} \quad R_{N,u} = \sup \left\{ r > 0 \mid g_u(r) \ge \frac{\beta}{N} \right\}. $$

In other words, using the following generalized inverses

$$ g_{1,u}^{-1} : y \mapsto \inf \{r > 0 \mid g_u(r) \le y\} \quad \text{and} \quad g_{2,u}^{-1} : y \mapsto \sup \{r > 0 \mid g_u(r) \ge y\}, $$

we have

$$ r_{N,u} = g_{1,u}^{-1} \left( \frac{\alpha}{N} \right) \quad \text{and} \quad R_{N,u} = g_{2,u}^{-1} \left( \frac{\beta}{N} \right). $$

The generalized inverses of $g_u$ are well defined, at least for $y > 0$ small enough, because $g_u(r) \to 0$ as $r \to +\infty$. Besides, by continuity of $g_u$, for $i \in \{1, 2\}$,

$$ g_u \left( g_{i,u}^{-1}(y) \right) = y. \quad (14) $$

The need for these definitions of the radii comes from the fact that $r \mapsto G((0, i_0), (ru, j))$ is not necessarily decreasing (think of a model with long jumps and no short jumps), hence the use of generalized inverses. Besides, in order to use Proposition 21, we cannot use the same generalized inverse for both radii.

If $r < r_{N,u}$, then $G((0, i_0), (ru, j)) > \frac{\alpha}{N}$ and if $r > R_{N,u}$, then $G((0, i_0), (ru, j)) < \frac{\beta}{N}$. Therefore, Proposition 21 shows the following.

**Proposition 26.** In direction $u$, for fixed $N$, the shape of the final configuration of the sandpile lies between radii $r_{N,u}$ and $R_{N,u}$.

To get a limit shape after normalization by $\log N$, we need to prove that $\frac{r_{N,u}}{\log N}$ and $\frac{R_{N,u}}{\log N}$ have a common limit, uniformly in $u$. According to Proposition 5, there exist constants $C_u, \gamma_u$ (where $C_u$ depends on $i_0, j$ and $u$; $\gamma_u = \Gamma^{-1}(u) \cdot u$ only depends on $u$, and both are continuous functions of $u$ and positive by Proposition 15) such that, when $r$ goes to infinity,

$$ g_u(r) = C_u e^{-\gamma_u r} r^{-\frac{d-1}{2}} (1 + o(1)). $$

Besides, these asymptotics are uniform in $u$, which means that there is a function $\varepsilon$ that does not depend on $u$ such that $\varepsilon(r) \xrightarrow[r \to +\infty]{} 0$ and $|o(1)| \le \varepsilon(r)$. We will exploit these asymptotics to get the asymptotics of $g_{u,i}^{-1}$, and thus of the radii $r_{N,u}$ and $R_{N,u}$.
We will denote by $f_u$ the asymptotic equivalent of $g_u$, that is:

$$f_u(r) := C_u e^{-\gamma_u r} r^{-\frac{d-1}{2}}, \quad r > 0.$$

Although the function $g_u$ is not necessarily invertible, the function $f_u$ is, and its inverse is easily expressed with the classical Lambert W function, which satisfies $W(y)e^{W(y)} = y$.

**Proposition 27.** Writing $p = \frac{2}{d-1} > 0$ and $f_u : r \mapsto C_u e^{-\gamma_u r} r^{-\frac{1}{p}}$, one has

$$f_u^{-1}(y) = \frac{1}{p\gamma_u} W(p\gamma_u C_u^p y^{-p}).$$

We remind the classical asymptotics of the function $W$ as $x \to +\infty$, see e.g. [11, Eq. (4.18)]:

$$W(x) = \log x - \log \log x + o(1). \qquad (15)$$

**Proposition 28.** We note $f_u^{-1}(y) = a_u W(b_u y^{-p})$, where $a_u$ and $b_u$ are positive and continuous functions of $u$. Then, when $y \to 0$,

$$f_u^{-1}(y) = a_u \left[ p \log \left( \frac{1}{y} \right) - \log \log \left( \frac{1}{y} \right) + \log \left( \frac{b_u}{p} \right) \right] + o(1)$$

and these asymptotics are uniform in $u$, that is, there exists a function $\varepsilon$ independent of $u$ such that $\varepsilon(y) \xrightarrow{y \to 0} 0$ and $|o(1)| \le \varepsilon(y)$.

*Proof.* Using (15), there exists a function $\varepsilon$ such that $\varepsilon(x) \xrightarrow{x \to +\infty} 0$ and such that $W(x) = \log x - \log \log x + \varepsilon(x)$. Therefore,

$$\begin{aligned} W(b_u y^{-p}) &= \log(b_u y^{-p}) - \log \log(b_u y^{-p}) + \varepsilon(b_u y^{-p}) \\ &= p \log(y^{-1}) + \log(b_u) - \log\left(p \log(y^{-1})\left(1 + \frac{\log(b_u)}{p \log(y^{-1})}\right)\right) + \varepsilon(b_u y^{-p}) \\ &= p \log(y^{-1}) + \log(b_u) - \log(p) - \log \log(y^{-1}) - \log\left(1 + \frac{\log(b_u)}{p \log(y^{-1})}\right) + \varepsilon(b_u y^{-p}). \end{aligned}$$

We now prove that $\log\left(1 + \frac{\log(b_u)}{p \log(y^{-1})}\right) + \varepsilon(b_u y^{-p})$ tends to 0, uniformly in $u$. Since $b_u$ is a positive continuous function of $u$, and $u$ lies in the compact set $S^{d-1}$, there exist constants $m, M > 0$ such that $\forall u \in S^{d-1}, m \le b_u \le M$. Therefore, $my^{-p} \le b_u y^{-p}$, hence

$$|\varepsilon(b_u y^{-p})| \le \sup_{x \ge my^{-p}} \varepsilon(x) \xrightarrow{y \to 0} 0,$$

where the bound $\sup_{x \ge my^{-p}} \varepsilon(x)$ no longer depends on $u$. Moreover, $\log(m) \le \log(b_u) \le \log(M)$, so

$$\left| \log \left( 1 + \frac{\log(b_u)}{p \log(y^{-1})} \right) \right| \le \max \left( \left| \log \left( 1 + \frac{\log(m)}{p \log(y^{-1})} \right) \right|, \left| \log \left( 1 + \frac{\log(M)}{p \log(y^{-1})} \right) \right| \right) \xrightarrow{y \to 0} 0.$$

In conclusion, if we define

$$\varepsilon(y) = \left( \max_{u \in S^{d-1}} a_u \right) \left[ \sup_{x \ge my^{-p}} \varepsilon(x) + \max \left( \left| \log \left( 1 + \frac{\log(m)}{p \log(y^{-1})} \right) \right|, \left| \log \left( 1 + \frac{\log(M)}{p \log(y^{-1})} \right) \right| \right) \right],$$

we have $\varepsilon(y) \xrightarrow{y \to 0} 0$ and

$$\left| f_u^{-1}(y) - a_u \left[ p \log \left( \frac{1}{y} \right) - \log \log \left( \frac{1}{y} \right) + \log \left( \frac{b_u}{p} \right) \right] \right| \le \varepsilon(y). \quad \square$$
**Lemma 29.** For $i \in \{1, 2\}$, $g_{i,u}^{-1}(y) \underset{y \to 0}{\longrightarrow} +\infty$ uniformly in $u$, that is:

$$ \forall A > 0, \quad \exists y_0 > 0, \quad \forall u \in \mathbb{S}^{d-1}, \quad \forall y \in (0, y_0), \quad g_{i,u}^{-1}(y) \ge A. \qquad (16) $$

*Proof.* Since $g_u = f_u \cdot (1 + o(1))$ at $+\infty$, uniformly in $u$, there exists $r_0 > 0$ such that for every $r \ge r_0$ and every $u \in \mathbb{S}^{d-1}$, $\frac{1}{2}f_u(r) \le g_u(r) \le 2f_u(r)$.

Let $A > r_0$. We define

$$ y_1 := \frac{1}{2} \min \{g_u(r) \mid 0 \le r \le r_0, u \in \mathbb{S}^{d-1}\}. $$

Note that $y_1$ is well defined by continuity of $(r, u) \mapsto g_u(r)$ and compactness of $\bar{B}(0, r_0) \times \mathbb{S}^{d-1}$. We also set

$$ y_2 := \frac{1}{2} \min_{u \in \mathbb{S}^{d-1}} f_u(A), $$

which is well defined because $u \mapsto f_u(A)$ is continuous and $\mathbb{S}^{d-1}$ is compact. We have $y_1, y_2 > 0$. We define $y_0 := \min(y_1, y_2)$. Let $u \in \mathbb{S}^{d-1}$ and $y \in (0, y_0)$. To prove that $g_{1,u}^{-1}(y) \ge A$, we need to prove that for every $r \le A$, $g_u(r) > y$. If $r \le r_0$, then $g_u(r) > y_1 \ge y$. If $r \in [r_0, A]$, then

$$ g_u(r) \ge \frac{1}{2} f_u(r) \ge \frac{1}{2} f_u(A) \ge y_2 \ge y_0 > y. $$

We thus have $g_{1,u}^{-1}(y) \ge A$. In conclusion, we proved (16) for $i = 1$. The result for $i = 2$ follows by observing that $g_{2,u}^{-1} \ge g_{1,u}^{-1}$. $\square$

**Proposition 30.** The functions $g_{i,u}^{-1}$ for $i \in \{1, 2\}$ have the same asymptotics as $f_u^{-1}$ in Proposition 28, that is, with $p = \frac{2}{d-1}$:

$$ g_{i,u}^{-1}(y) = a_u \left[ p \log \left( \frac{1}{y} \right) - \log \log \left( \frac{1}{y} \right) + \log \left( \frac{b_u}{p} \right) \right] + o(1) $$

uniformly in $u$. As a consequence, uniformly in $u$,

$$ r_{N,u} = \frac{1}{\gamma_u} \log N - \frac{d-1}{2\gamma_u} \log \log N + \frac{1}{\gamma_u} \log \left( \frac{\gamma_u^{\frac{d-1}{2}} C_u}{\alpha} \right) + o(1), $$

$$ R_{N,u} = \frac{1}{\gamma_u} \log N - \frac{d-1}{2\gamma_u} \log \log N + \frac{1}{\gamma_u} \log \left( \frac{\gamma_u^{\frac{d-1}{2}} C_u}{\beta} \right) + o(1). $$

*Proof.* Let $y > 0$ and $x = g_{i,u}^{-1}(y)$. By (14), we have $g_u(x) = y$. Besides, since $g_u \underset{+\infty}{\sim} f_u$, uniformly in $u$, there exists a function $q_u$ that tends to 1 at $+\infty$, uniformly in $u$, such that $g_u = q_u f_u$. Therefore

$$ g_{i,u}^{-1}(y) = x = f_u^{-1}(f_u(x)) = f_u^{-1}\left(\frac{g_u(x)}{q_u(x)}\right) = f_u^{-1}\left(\frac{y}{q_u(g_{i,u}^{-1}(y))}\right). $$

We denote $\tilde{q}_{i,u}(y) = q_u(g_{i,u}^{-1}(y))$. Since $g_{i,u}^{-1}(y) \underset{y \to 0}{\longrightarrow} +\infty$ uniformly in $u$ according to Lemma 29, and $q_u(r) \underset{r \to +\infty}{\longrightarrow} 1$ uniformly in $u$, we have $\tilde{q}_{i,u}(y) \underset{y \to 0}{\longrightarrow} 1$ uniformly in $u$. Therefore, Proposition 28
leads to the existence of a function $\varepsilon$ that tends to 0 and does not depend on $u$ such that

$$
\begin{align*}
g_{i,u}^{-1}(y) &= f_u^{-1}\left(\frac{y}{\tilde{q}_{i,u}(y)}\right) \\
&= a_u \left[ p \log(y^{-1}\tilde{q}_{i,u}(y)) - \log \log(y^{-1}\tilde{q}_{i,u}(y)) + \log\left(\frac{b_u}{p}\right) \right] + \varepsilon\left(\frac{y}{\tilde{q}_{i,u}(y)}\right) \\
&= a_u \left[ p \log(y^{-1}) + p \log(\tilde{q}_{i,u}(y)) - \log\left(\log(y^{-1}) + \log(\tilde{q}_{i,u}(y))\right) + \log\left(\frac{b_u}{p}\right) \right] + \varepsilon\left(\frac{y}{\tilde{q}_{i,u}(y)}\right) \\
&= a_u \left[ p \log(y^{-1}) + p \log(\tilde{q}_{i,u}(y)) - \log\left(\log(y^{-1})\left(1 + \frac{\log(\tilde{q}_{i,u}(y))}{\log(y^{-1})}\right)\right) + \log\left(\frac{b_u}{p}\right) \right] + \varepsilon\left(\frac{y}{\tilde{q}_{i,u}(y)}\right) \\
&= a_u \left[ p \log(y^{-1}) - \log \log(y^{-1}) + \log\left(\frac{b_u}{p}\right) \right] \\
&\qquad + a_u p \log(\tilde{q}_{i,u}(y)) - a_u \log\left(1 + \frac{\log(\tilde{q}_{i,u}(y))}{\log(y^{-1})}\right) + \varepsilon\left(\frac{y}{\tilde{q}_{i,u}(y)}\right). \tag{17}
\end{align*}
$$

We shall prove that the three terms appearing in (17) tend to 0 uniformly in $u$ as $y \to 0$. It is clear for $a_u p \log(\tilde{q}_{i,u}(y))$, because $a_u$ is bounded (by a continuity plus compactness argument) and $\tilde{q}_{i,u}(y)$ tends to 1 uniformly in $u$. Similarly, $\frac{\log(\tilde{q}_{i,u}(y))}{\log(y^{-1})}$ tends to 0 uniformly in $u$ as $y \to 0$, so the same conclusion holds for the second term in (17). Finally, $\frac{y}{\tilde{q}_{i,u}(y)}$ tends to 0 uniformly in $u$ as $y \to 0$. Therefore, $\varepsilon(\frac{y}{\tilde{q}_{i,u}(y)})$ tends to 0 uniformly in $u$. $\square$

We can summarize it into the following theorem.

**Theorem 31.** We have

$$ \lim_{N \to +\infty} \frac{r_{N,u}}{\log N} = \lim_{N \to +\infty} \frac{R_{N,u}}{\log N} = \frac{1}{\gamma_u} = \frac{1}{\Gamma^{-1}(u) \cdot u} $$

and these limits are uniform in $u$. This means that the limit shape of the LASM rescaled by $\log N$ is delimited by the curve

$$ C = \left\{ \frac{u}{\Gamma^{-1}(u) \cdot u} \mid u \in \mathbb{S}^{d-1} \right\}, \qquad (18) $$

and this limit is uniform in the sense of Definition 22.

Besides, since the asymptotics of $r_{N,u}$ and $R_{N,u}$ are the same until the constant term, there exists a constant $\tau$ such that

$$ \sup_{N \in \mathbb{N}, u \in \mathbb{S}^{d-1}} (R_{N,u} - r_{N,u}) \le \tau, $$

hence

$$ \sup_{u \in \mathbb{S}^{d-1}} \frac{1}{\log N} (R_{N,u} - r_{N,u}) \le \frac{\tau}{\log N}. $$

In the next result we characterise the limit shape (18) in terms of duality (see Appendix A for various reminders concerning duality in our context).

**Proposition 32.** The limit shape of the sandpile defined by (18) and the set $\partial T := \rho^{-1}(\{1\})$ introduced in Definition 7 are dual curves, in the sense of Definition 49.
Proof. Let $t_0 \in \partial \mathcal{T}$. The normalized gradient of $\rho$ at $t_0$ is, by definition, $\Gamma(t_0)$, so the normal cone to $\mathcal{T}$ at $t_0$ consists of the non-negative multiples of $\Gamma(t_0)$. By Definition 49, the point $t_0^*$ associated with $t_0$ in the dual curve of $\partial \mathcal{T}$ is $t_0^* = \frac{\Gamma(t_0)}{\Gamma(t_0) \cdot t_0}$. Writing $t_0^* = ru$ in spherical coordinates, we get $u = \Gamma(t_0)$ and $r = \frac{1}{\Gamma(t_0) \cdot t_0} = \frac{1}{u \cdot \Gamma^{-1}(u)}$, hence $t_0^* = \frac{1}{u \cdot \Gamma^{-1}(u)}u$ is in the limit shape defined in (18).

Conversely, if $\frac{1}{\Gamma^{-1}(u) \cdot u}u$ is in the limit shape defined by (18), then setting $t_0 := \Gamma^{-1}(u)$, the previous computations show that $\frac{1}{\Gamma^{-1}(u) \cdot u}u = t_0^*$ is in the dual of $\partial \mathcal{T}$. $\square$

Since duality preserves convexity (see Proposition 52), we obtain the following.

**Proposition 33.** The limit curve $\mathcal{C}$ in (18) is the boundary of a convex set.

# 6 Behavior of the limit curve when the leakiness parameter goes to zero or infinity

We study the limit curve $\mathcal{C}$ given in (18) in two special regimes, first when the leakiness parameter goes to infinity, and second when it goes to zero (the massless case). In both cases, the general strategy is to work on the level set $\partial \mathcal{T}$, which is the dual curve of the limit shape of the sandpile, according to Proposition 32, and to use duality to get the limit in the regimes studied. Useful properties of dual curves, especially Proposition 56 which links the convergence of convex sets with the convergence of their duals, are recalled in Appendix A.

## 6.1 Infinite leakiness parameter case

Throughout Section 6.1, in addition to the running Assumption 1, we will assume that the supports of the random walks are all bounded.

Let us first consider the uncolored case, i.e., $p = 1$. In this simpler case, the limit shape of the sandpile as $m$ goes to infinity is easily understood using the support of the random walk. The Laplace transform of the associated KRW is denoted by $L\mu(t) = \sum_{x \in \mathbb{Z}^d} e^{t \cdot x} \mu(x)$ ($\mu$ is $\mu_{1,1}$ in our notation (3)). Let $\mathcal{C}$ be the curve that defines the boundary shape of the LASM, see (18).

**Proposition 34.** We denote $\tilde{\mu} = m\mu$, which is a probability measure and does not depend on $m$, according to notation (3). The rescaled level set $\frac{1}{\log m}\{t \in \mathbb{R}^d | L\mu(t) \le 1\}$ is between the polytopes

$$ \bigcap_{x \in \text{supp } \mu} \{s \in \mathbb{R}^d \mid s \cdot x \le 1\} \qquad (19) $$

and

$$ \bigcap_{x \in \text{supp } \mu} \left\{ s \in \mathbb{R}^d \mid s \cdot x \le 1 - \frac{\log \tilde{\mu}(x)}{\log m} \right\}. $$

As a consequence, if the support of $\mu$ is finite, the rescaled level set converges uniformly to the polytope (19) in the sense of Definition 22 (and therefore also for the Hausdorff distance). In that case, the limit polytope (19) is the dual of the convex hull of the support of $\mu$.

See Figure 3 for an illustration of Proposition 34. Using general results on duality, the previous result implies that, when $p = 1$, the curve $(\log m) \cdot \mathcal{C}$ converges to the convex hull of the support of the random walk as $m$ goes to infinity. Moreover, Proposition 34 can be interpreted using first passage times, as will be done in Section 6.2.
Proof. Let $t \in \mathbb{R}^d$. Writing $s = \frac{t}{\log m}$, we have

$$L\mu(t) = \sum_{x \in \text{supp } \mu} e^{(\log m)(s \cdot x - 1)} \tilde{\mu}(x).$$

If $s \cdot x \le 1$ for every $x \in \text{supp } \mu$, then $L\mu(t) \le \sum_{x \in \text{supp } \mu} \tilde{\mu}(x) = 1$. Conversely, if $L\mu(t) \le 1$, then each term of the sum is no bigger than 1, that is $s \cdot x \le 1 - \frac{\log \tilde{\mu}(x)}{\log m}$.

Using the notations of Definition 22, writing $s = ru$, we can define

$$r_{N,u} = \frac{1}{u \cdot x} \quad \text{and} \quad R_{N,u} = \frac{1}{u \cdot x} + \frac{1}{\log m} \max_{\substack{y \in \text{supp } \mu \\ y \cdot u > 0}} \frac{-\log \tilde{\mu}(y)}{u \cdot y},$$

where $x$ is the element of the support that maximizes $u \cdot x$. It is then clear that $r_{N,u}$ and $R_{N,u}$ both converge uniformly to $\frac{1}{u \cdot x}$.

The fact that (19) is the dual of the convex hull of $\text{supp } \mu$ follows directly from Proposition 58. $\square$

Let us now consider the general case in $p$. Our main objective is to state and prove Theorem 39. We first introduce useful notation and preliminary results. For $q \in \{1, \dots, p\}$, define the polytope $Q_{\pm}^{(q)}$ as follows. First, for $q = 1$, we let

$$Q_{\pm}^{(1)} = \bigcap_{1 \le i \le p} \bigcap_{x \in \text{supp } \mu_{i,i}} \{s \in \mathbb{R}^d \mid s \cdot x \le 1 \pm c_m^{(i)}\},$$

where $c_m^{(i)}$ are some constants. For $q \ge 2$, we define

$$Q_{\pm}^{(q)} = \bigcap_{i_1, i_2, \dots, i_q} \bigcap_{\substack{x_{1,2} \in \text{supp } \mu_{i_1, i_2} \\ \dots \\ x_{q,1} \in \text{supp } \mu_{i_q, i_1}}} \{s \in \mathbb{R}^d \mid s \cdot (x_{1,2} + x_{2,3} + \dots + x_{q,1}) \le q \pm c_m^{(i_1, \dots, i_q)}\},$$

where $c_m^{(i_1, \dots, i_q)}$ are some constants. Then we introduce

$$P_{\pm} = \bigcap_{1 \le q \le p} Q_{\pm}^{(q)} \qquad (20)$$

The above polytopes obviously satisfy $P_{-} \subset P_{+}$. They both depend on a collection of constants

$$\{c_m^{(i_1, \dots, i_q)} \mid 1 \le i_1, \dots, i_q \le p, 1 \le q \le p\}.$$

For the sake of brevity, we will name these constants $c_m$. When they are all equal to zero, we have $P_{-} = P_{+}$, which we simply denote $P$.

**Proposition 35.** Let $p \ge 1$. We assume that $m_1 = \dots = m_p = m \to \infty$. There exist constants $c_m$ going to 0 such that the rescaled level set $\frac{1}{\log m}\{t \in \mathbb{R}^d \mid \rho(t) \le 1\}$ lies between the two polytopes $P_{-}$ and $P_{+}$, as defined in (20). In particular, it converges to the polytope $P$.

Proof. Performing the change of variable $s = \frac{t}{\log m}$, we first prove that the level set $\frac{1}{\log m}\{t \in \mathbb{R}^d \mid \rho(t) \le 1\}$ is included in some $P_{+}$, for suitable constants $c_m$ going to zero when $m \to \infty$.

If $\rho(t) \le 1$ then for all $q \in \{1, \dots, p\}$, one has $\text{Tr}((L\mu)^q) \le p$. For $q = 1$, we obtain that

$$\sum_{1 \le i \le p} \sum_{x \in \text{supp } \mu_{i,i}} e^{\log m(s \cdot x - 1)} \tilde{\mu}_{i,i}(x) \le p,$$
using our notation $\tilde{\mu}_{i,i} = m\mu_{i,i}$ ($\tilde{\mu}_{i,i}$ is a sub-probability measure). Accordingly, for all $1 \le i \le p$ and all $x \in \text{supp } \mu_{i,i}$, we should have

$$s \cdot x \le 1 + \frac{\log \frac{p}{\tilde{\mu}_{i,i}(x)}}{\log m}.$$

This immediately gives that $s$ should belong to $Q_+^{(1)}$, with some explicit constants

$$c_m^{(i)} = \sup_{x \in \text{supp } \mu_{i,i}} \frac{\log \frac{p}{\tilde{\mu}_{i,i}(x)}}{\log m},$$

which obviously go to zero when $m \to \infty$.

We now look at the condition $\text{Tr}((L\mu)^q) \le p$ for some $q \ge 2$. Using the definition of the trace, we get

$$\text{Tr}((L\mu)^q) = \sum_{i_1, i_2, \dots, i_q} \sum_{\substack{x_{1,2} \in \text{supp } \mu_{i_1, i_2} \\ \dots \\ x_{q,1} \in \text{supp } \mu_{i_q, i_1}}} e^{\log m(s \cdot (x_{1,2} + x_{2,3} + \dots + x_{q,1}) - q)} \tilde{\mu}_{i_1, i_2}(x_{1,2}) \dots \tilde{\mu}_{i_q, i_1}(x_{q,1}) \le p. \quad (21)$$

By the same argument as above, we deduce that $s \in Q_+^{(q)}$, with the constants

$$c_m^{(i_1, \dots, i_q)} = \sup_{\substack{x_{1,2} \in \text{supp } \mu_{i_1, i_2} \\ \dots \\ x_{q,1} \in \text{supp } \mu_{i_q, i_1}}} \frac{\log \frac{p}{\tilde{\mu}_{i_1, i_2}(x_{1,2}) \dots \tilde{\mu}_{i_q, i_1}(x_{q,1})}}{\log m}.$$

We immediately conclude that $s \in \mathcal{P}_+$, with constants $c_m$ all going to 0 (due to the boundedness of the supports).

We now assume that $s \in \mathcal{P}_-$ for some constants $c_m$ going to zero and prove that $\rho(t) \le 1$. Using (21) and the definition of $\mathcal{P}_-$, we obtain

$$\text{Tr}((L\mu)^q) \le \sum_{i_1, i_2, \dots, i_q} \sum_{\substack{x_{1,2} \in \text{supp } \mu_{i_1, i_2} \\ \dots \\ x_{q,1} \in \text{supp } \mu_{i_q, i_1}}} e^{-(\log m) c_m^{(i_1, \dots, i_q)}} \tilde{\mu}_{i_1, i_2}(x_{1,2}) \dots \tilde{\mu}_{i_q, i_1}(x_{q,1}) \le \sum_{i_1, i_2, \dots, i_q} e^{-(\log m) c_m^{(i_1, \dots, i_q)}}.$$

We see that adjusting the constants $c_m$ (for instance $\frac{1}{\sqrt{\log m}}$), all the traces $\text{Tr}((L\mu)^q)$ can be made as small as wanted. Using then Lemma 36, one can also control the spectrum $\rho(t)$ and make it less than one. The proof is complete. $\square$

**Lemma 36.** *One has*

$$\sup \{\rho(M) \mid M \in \mathcal{M}_p(\mathbb{R}), |\text{Tr}(M)| \le \alpha, \dots, |\text{Tr}(M^p)| \le \alpha\} \xrightarrow{\alpha \to 0} 0.$$

*Proof.* Let us write $a_i(M)$ for the coefficients of the characteristic polynomial of $M$, so that

$$\chi_M = X^p + a_{p-1}(M)X^{p-1} + \dots + a_1(M)X + a_0(M).$$

By continuity of the roots of monic polynomials with fixed degree, if we write $\rho(P)$ for the largest root of $P$ in modulus, the function

$$(a_{p-1}, \dots, a_1, a_0) \in \mathbb{R}^p \mapsto \rho(X^p + a_{p-1}X^{p-1} + \dots + a_1X + a_0)$$
<figure><img src="image_6.png" /><figcaption>Figure 3: Left picture: the step set $\{(1,0), (2,-2), (0,-1), (-1,0), (1,2)\}$. Middle picture: the limit level set. Right picture: the Newton polytope associated with the step set</figcaption></figure>

is continuous at $(0, \dots, 0)$, where it takes the value 0. Therefore, it is sufficient to prove that for all $i \in \{0, \dots, p-1\}$,

$$ \sup \{|a_i(M)| \mid M \in \mathcal{M}_p(\mathbb{R}), |\mathrm{Tr}(M)| \le \alpha, \dots, |\mathrm{Tr}(M^p)| \le \alpha\} \xrightarrow{\alpha \to 0} 0. \quad (22) $$

To do so, we use the Newton identities to write the elementary symmetric polynomials in terms of Newton sums. With Vieta's formulas, it leads to an expression of the non-leading coefficients of a monic polynomial $P$ of degree $p$ in terms of the $\sum_{i=1}^p \lambda_i^k$ for $k \in \{1, \dots, p\}$, where the $\lambda_i$ are the roots of $P$. Applying this to the characteristic polynomial of a matrix $M$, it means that all the $a_i(M)$ have a polynomial expression in terms of the $\mathrm{Tr}(M^k)$ for $k \in \{1, \dots, p\}$, and it is easy to see that this polynomial expression has no constant term. Therefore, the convergences mentioned in (22) become obvious. $\square$

We now use general results on duality to prove that the rescaled curve $\mathcal{C}$ converges to the dual polytope of $\mathcal{P}$. In order to apply Appendix A, we need $\mathcal{P}$ to be compact, hence the following lemma.

**Lemma 37.** *The limit polytope $\mathcal{P}$ is bounded.*

*Proof.* The polytope $\mathcal{P}$ is convex. Therefore, to prove the lemma, we only need to show that $\mathcal{P}$ does not contain a half-line starting from 0.

Let $u \in S^{d-1}$. We write $H(q)$ for the following assertion: there exist $i_1, \dots, i_q \in \{1, \dots, p\}$ and $x_{1,2} \in \mathrm{supp} \, \mu_{i_1, i_2}, \dots, x_{q,1} \in \mathrm{supp} \, \mu_{i_q, i_1}$ such that $u \cdot (x_{1,2} + \dots + x_{q,1}) > 0$. We want to prove that $H(q)$ is true for some $q \in \{1, \dots, p\}$, so that the condition $s \cdot (x_{1,2} + \dots + x_{q,1}) \le q$ prevents $\mathcal{P}$ from containing the half-line $[0, +\infty)u$.

By the assumption of irreducibility, $H(q)$ is true for some $q \in \mathbb{N}$. We choose the minimal $q$. Assume by contradiction that $q > p$. Then, in the path $i_1 \to \dots \to i_q \to i_1$, there is a cycle $i_k \to \dots \to i_\ell = i_k$ with $1 \le \ell - k \le p < q$. Since $H(\ell - k)$ is false, $u \cdot (x_{k,k+1} + \dots + x_{\ell-1,\ell}) \le 0$. Therefore, we have

$$ \begin{aligned} u \cdot (x_{1,2} + \dots + x_{k-1,k} + x_{\ell,\ell+1} + \dots + x_{q,1}) &= u \cdot (x_{1,2} + \dots + x_{q,1}) - u \cdot (x_{k,k+1} + \dots + x_{\ell-1,\ell}) \\ &\ge u \cdot (x_{1,2} + \dots + x_{q,1}) > 0, \end{aligned} $$

where

$$ x_{1,2} \in \mathrm{supp} \, \mu_{i_1, i_2}, \dots, x_{k-1,k} \in \mathrm{supp} \, \mu_{i_{k-1}, i_k}, x_{\ell, \ell+1} \in \mathrm{supp} \, \mu_{i_\ell, i_{\ell+1}} = \mathrm{supp} \, \mu_{i_k, i_{\ell+1}}, \dots, x_{q,1} \in \mathrm{supp} \, \mu_{i_q, i_1}. $$

This proves $H(q - (\ell - k))$, which contradicts the minimality of $q$. $\square$
**Proposition 38.** The limit polytope $\mathcal{P}$ is the dual of the convex hull of the set

$$ \mathcal{X} = \bigcup_{q \in \{1, \dots, p\}} \bigcup_{i_1, \dots, i_q \in \{1, \dots, p\}} \left\{ \frac{x_{1,2} + \dots + x_{q,1}}{q} \mid x_{1,2} \in \text{supp } \mu_{i_1, i_2}, \dots, x_{q,1} \in \text{supp } \mu_{i_q, i_1} \right\}. $$

*Proof.* It is exactly Proposition 58. □

**Theorem 39.** If the support of $\mu$ is finite, the rescaled shape of the sandpile, $(\log m) \cdot \mathcal{C}$, where $\mathcal{C}$ is defined in (18), converges to the convex hull of $\mathcal{X}$ as $m$ goes to $+\infty$. The convergence is uniform in the sense of Definition 22 and consequently holds for the Hausdorff distance.

In particular, in the uncolored case $p = 1$ the limit shape of the sandpile when $m \to \infty$ is the convex hull of the support of the measure $\mu$.

*Proof.* Using Propositions 32 and 55, $(\log m) \cdot \mathcal{C}$ is the dual curve of $\frac{1}{\log m} \{\rho(t) \le 1\}$. The inclusions of Proposition 35 imply that

$$ (1 - \varepsilon_m) \mathcal{P} \subset \frac{1}{\log m} \{\rho(t) \le 1\} \subset (1 + \varepsilon_m) \mathcal{P}, $$

where $\varepsilon_m = \max \frac{c_m^{(i_1, \dots, i_q)}}{q} \xrightarrow{m \to +\infty} 0$. Therefore, according to Proposition 56, $(\log m) \cdot \mathcal{C}$ converges to $\mathcal{P}^*$, which is the convex hull of $\mathcal{X}$. □

## 6.2 First passage time interpretation of the limiting polytope in the infinite leakiness parameter region

Our main result in this section is the following:

**Theorem 40.** If the support of $\mu$ is finite, as $m$ goes to $+\infty$ the rescaled shape of the sandpile converges to the region of points that can be reached by $n$ steps of the random walk, rescaled by $n$.

*Proof.* The theorem is immediate from Theorem 39 and Proposition 41 below. □

Let $\mathcal{X}$ be as in Proposition 38. Given $(x, i) \in \mathbb{Z}^d \times \{1, \dots, p\}$, let $n_{x,i}$ be the smallest positive integer such that $(x, i)$ can be reached in $n_{x,i}$ steps by our killed random walk. Denote by $A_n$ the set of points for which this first passage time is at most $n$, i.e.,

$$ A_n = \{(x, i) \in \mathbb{Z}^d \times \{1, \dots, p\} : n_{x,i} \le n\}. $$

**Proposition 41.** Suppose the support of $\mu$ is finite. As $n$ converges to $\infty$, the set $A_n$, scaled by $n$, converges to the convex hull of $\mathcal{X}$ in the Hausdorff distance, i.e.

$$ \lim_{n \to \infty} d_H \left( \text{Conv}(\mathcal{X}), \frac{1}{n} \cdot A_n \right) = 0. $$

*Proof.* First, consider the case $p = 1$. We will ignore the color index. Let $s = |\chi|$ and write $\mathcal{X} = \{x_1, x_2, \dots, x_s\}$ be the support of $\mu$, which is a finite subset of $\mathbb{Z}^d$ such that the $\mathbb{N}^d$-span of $\mathcal{X}$ is all of $\mathbb{Z}^d$. Without loss of generality we can assume $0 \in \mathcal{X}$.

If $x \in A_n$, then there exist points $x_{i_1}, \dots, x_{i_n} \in \mathcal{X}$ such that $x = \sum_{j=1}^n \frac{1}{n} x_{i_j}$. Thus

$$ \frac{1}{n}x = \sum_{j=1}^{n} \frac{1}{n}x_{i_j} \in \text{Conv}(\mathcal{X}), $$
so $\frac{1}{n} \cdot A_n \subset \text{Conv}(\mathcal{X})$.

To complete the argument we need to show that for any $y \in \text{Conv}(\mathcal{X})$, its distance from $\frac{1}{n} \cdot A_n$ converges to 0 uniformly in $y$ as $n \to \infty$. If $h$ is any fixed integer, then the Hausdorff distance between $\frac{1}{n} \cdot A_n$ and $\frac{1}{n} \cdot A_{n+h}$ converges to 0 as $n \to \infty$ since the support of $\mu$ is finite. Thus, it is enough to show that the distance of $y$ from $\frac{1}{n} \cdot A_{n+h}$ converges to 0 uniformly in $y$ as $n \to \infty$. For the constant $h$ we take the largest number of steps the random walk needs to reach a point in $s \cdot \text{Conv}(\mathcal{X})$, i.e. $h := \max\{n_x : x \in s \cdot \text{Conv}(\mathcal{X}) \cap \mathbb{Z}^d\}$.

Given $y \in \text{Conv}(\mathcal{X})$, there is a point $x \in n \cdot \text{Conv}(\mathcal{X}) \cap \mathbb{Z}^d$ such that $\|y - x/n\| \le \frac{1}{n}$. Since $x \in n \cdot \text{Conv}(\mathcal{X}) \cap \mathbb{Z}^d$, there exist non-negative real numbers $\alpha_1, \dots, \alpha_s$, which sum to 1, such that $x = \sum_{j=1}^s n\alpha_i x_i$. We can write

$$x = \sum_{i=1}^{s} \lfloor n\alpha_i \rfloor x_i + \sum_{i=1}^{s} \{n\alpha_i\} x_i, \qquad (23)$$

where $\lfloor t \rfloor$ stands for the integer part of $t$ and $\{t\}$ for $t - \lfloor t \rfloor$. Since both $x$ and the first sum on the right-hand side are from $\mathbb{Z}^d$, so is the second sum. Since each coefficient $\{n\alpha_i\} \in [0, 1)$, we have that $\sum_{i=1}^s \{n\alpha_i\} \in [0, s)$. Thus $\sum_{i=1}^s \{n\alpha_i\} x_i \in s \cdot \text{Conv}(\mathcal{X})$, so it can be written as the sum of at most $h$ elements of $\mathcal{X}$. Since $\sum_{i=1}^s \lfloor n\alpha_i \rfloor \le n$, we see that $x$ can be written as the sum of at most $n+h$ elements from $\mathcal{X}$, hence $x \in A_{n+h}$. It follows that $d(y, \frac{1}{n} \cdot A_{n+h}) \le \frac{1}{n}$, completing the argument in the case $p=1$.

We now move to the case of a general $p$. First, let us introduce some notation. Given colors $i$ and $j$, let $k_{i,j}$ be the minimum number of steps the random walk requires to move from color $i$ to color $j$, and let $\bar{x}_{i,j}$ be the total displacement of one such minimum walk. Let $k = \max\{k_{i,1} : i \in \{1, \dots, p\}\}$.

Now, suppose $(x, i) \in A_n \setminus A_{n-1}$. Thus, there are steps $x_{0,1}, \dots, x_{n-1,n}$, with $x_{\ell, \ell+1} \in \text{supp} \mu_{i_\ell, i_{\ell+1}}$ for some colors $i_0 = 1, i_1, \dots, i_{n-1}$, $i_n = i$, such that $x = \sum_{j=1}^n x_{j-1,j}$. By adding $k_{i,1}$ more steps $x_{n,n+1}, \dots, x_{n+k_{i-1},n+k_i}$ with $x_{n+\ell,n+\ell+1} \in \mu_{i_{n+\ell}, i_{n+\ell+1}}$, $i_{n+k_i} = 1$ and

$$x_{n,n+1} + \dots + x_{n+k_{i-1},n+k_i} = \bar{x}_{i,1},$$

we can ensure our walk starts and ends with the first color 1. We will show that $x + \bar{x}_{i,1} \in (n+k) \cdot \text{Conv}(\mathcal{X})$ by induction on $n$. If no colors in $\{i_0, \dots, i_{n+k}\}$ repeat, we are done. If there are repeating colors, there exist $0 \le j_1 < j_1 + t \le n+k$ such that the colors $i_{j_1}, i_{j_1+1}, \dots, i_{j_1+t}$ are all distinct except the first and the last one. Removing the portion $x_{i_{j_1},i_{j_1+1}}, \dots, x_{i_{j_1+t-1},i_{j_1+t}}$ from the walk, we get a walk of length $n - t + k_i$ which respects colors and starts and ends with color 1, so by induction,

$$x + \bar{x}_{i,1} - \sum_{\ell=1}^{t} x_{i_1+\ell-1, i_1+\ell} \in (n+k-t) \cdot \text{Conv}(\mathcal{X}).$$

Since $\sum_{\ell=1}^t x_{i_1+\ell-1, i_1+\ell} \in t \cdot \text{Conv}(\mathcal{X})$, we get $x + \bar{x}_{i,1} \in (n+k) \cdot \text{Conv}(\mathcal{X})$ which implies

$$\frac{1}{n+k}(x + \bar{x}_{i,1}) \in \text{Conv}(\mathcal{X}). \qquad (24)$$

Now, suppose $x \in n \cdot \text{Conv}(\mathcal{X})$. If we ignore colors, by the case $p=1$, there exists $h \in \mathbb{N}$ such that $x$ can be written as the sum of at most $n+h$ elements of $\mathcal{X}$. In the proof of the case $p=1$, in Equation (23), instead of separating the integer part, for each term of the form $\frac{x_{1,2}+\dots+x_{q,1}}{q} \in \mathcal{X}$, we can separate a multiple of $p!$, thus ensuring that the coefficient of $\frac{x_{1,2}+\dots+x_{q,1}}{q} \in \mathcal{X}$ will be a multiple of $q$. We can thus write $x$ as a non-negative integer linear combination of at most $n+h$ terms of the form $x_{1,2} + \dots + x_{q,1}$. Each of these terms corresponds to a walk that starts and ends at the
same color. We can arrange this sum, so that all the walks that start with the same color are next to each other. For example we can start with all our walks that start and end in color 1, followed by all the walks that start and end in color 2, etc. This will not correspond to a color respecting random walk from $A_{n+h}$, but by adding at most $p-1$ walks of the type $\bar{x}_{i,j}$ between these, each of which has length at most $k$, we get $x + \sum \bar{x}_{i,j} \in A_{n+h+k(p-1)}$, so

$$ \frac{1}{n} \left( x + \sum \bar{x}_{i,j} \right) \in \frac{1}{n} \cdot A_{n+h+k(p-1)}. $$

Combining this with (24) gives the desired result since $h, k$ and $\bar{x}_{i,j}$ are independent of $x$ and $n$. $\square$

## 6.3 Zero leakiness parameter case

In this section, we no longer assume that the supports of the random walks are bounded. We suppose that $m_1 = \dots = m_p = m \to 1$. We denote $\tilde{\mu}_{i,j} = m\mu_{i,j}$, which no longer depends on $m$, and $\tilde{\rho} = \rho(L\tilde{\mu})$. With these notations, by Proposition 32 the shape of the sandpile when $N \to +\infty$ for fixed $m$ is the dual of the set

$$ \{\tilde{\rho} = m\}. $$

We assume that $\nabla\tilde{\rho}(0) = 0$. It means that the process with jumps $\tilde{\mu}$ is centered (see [6, Prop. 2.16]). If the process is not centered, then the set $\{\tilde{\rho} = 1\}$ was studied in [6] and does not require any scaling.

Under these assumptions, $\tilde{\rho}$ is convex, has a strict, global minimum at 0, and we remind [6, Prop. 2.14] the asymptotics of $\tilde{\rho}$ at 0:

$$ \tilde{\rho}(t) = 1 + \frac{1}{2} t^{\top} \sigma t + o(\|t\|^2), \qquad (25) $$

where $\sigma$ is the energy matrix from [6]. When $p = 1$, the energy matrix is simply the covariance matrix of $\tilde{\mu}$.

**Proposition 42.** Let $\delta \in (0, 1)$. There exists $m_0 > 1$ such that for all $m \in (1, m_0)$, we have

$$ \{s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le 1 - \delta\} \subset \frac{1}{\sqrt{m-1}} \{t \in \mathbb{R}^d \mid \tilde{\rho}(t) \le m\} \subset \{s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le 1 + \delta\}. $$

As a consequence, the rescaled set $\frac{1}{\sqrt{m-1}}\{\tilde{\rho} \le m\}$ converges uniformly to the ellipsoid $\{s \in \mathbb{R}^d \mid \frac{1}{2}s^\top \sigma s \le 1\}$ as $m \to 1$.

*Proof.* Since $\sigma$ is positive-definite by [6, Prop. 2.8], $\alpha := \min_{u \in \mathbb{S}^{d-1}} u^\top \sigma u > 0$. Let $\varepsilon \in (0, \frac{\alpha}{2})$. The Taylor expansion (25) ensures that there exists a neighbourhood $V$ of 0 such that for $t \in V$,

$$ \left| \tilde{\rho}(t) - 1 - \frac{1}{2} t^{\top} \sigma t \right| \le \frac{\varepsilon}{2} \|t\|^2. \qquad (26) $$

Let $t \in \mathbb{R}^d$ such that $\tilde{\rho}(t) \le m$. Since $\tilde{\rho}$ is convex and has a global strict minimum at 0, if $m$ is close enough to 1, then $t \in V$, which allows to use (26). We set $s = \frac{t}{\sqrt{m-1}}$. Then

$$ \frac{1}{2} s^{\top} \sigma s = \frac{1}{2(m-1)} t^{\top} \sigma t \le \frac{1}{m-1} \left( \tilde{\rho}(t) - 1 + \frac{\varepsilon}{2} \|t\|^2 \right) \le 1 + \frac{\varepsilon}{2} \|s\|^2. \qquad (27) $$

Writing $s = \|s\|u$ with $u \in \mathbb{S}^{d-1}$, (27) becomes $\frac{1}{2}\|s\|^2 u^\top \sigma u \le 1 + \frac{\varepsilon}{2}\|s\|^2$, so

$$ \|s\|^2 \le \frac{2}{u^\top \sigma u - \varepsilon} \le \frac{4}{\alpha}, $$
and (27) finally becomes

$$ \frac{1}{2} s^{\top} \sigma s \le 1 + \frac{2}{\alpha} \varepsilon. $$

Conversely, let $\eta = \frac{1}{1+\frac{\varepsilon}{\alpha}}$. Assume that $\frac{1}{2} s^\top \sigma s \le \eta$. Then, writing $s = \|s\|u$ for $u \in \mathbb{S}^{d-1}$, we get $\|s\|^2 \le \frac{2\eta}{u^\top \sigma u} \le \frac{2\eta}{\alpha}$. Therefore, $\|t\| \le \sqrt{\frac{2\eta(m-1)}{\alpha}} \le \sqrt{\frac{2(m-1)}{\alpha}}$. This ensures that if $m$ is close enough to 1, we have $t \in V$, which allows us to use (26). It yields

$$ \tilde{\rho}(t) \le 1 + \frac{1}{2} t^{\top} \sigma t + \frac{\varepsilon}{2} \|t\|^2 \le 1 + (m-1)\eta + \varepsilon \frac{(m-1)\eta}{\alpha} = m. $$

To summarize, we have proved that for $m$ close enough to 1,

$$ \left\{ s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le \frac{1}{1 + \frac{\varepsilon}{\alpha}} \right\} \subset \frac{1}{\sqrt{m-1}} \left\{ t \in \mathbb{R}^d \mid \tilde{\rho}(t) \le m \right\} \subset \left\{ s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le 1 + \frac{2}{\alpha} \varepsilon \right\}. $$

Choosing $\varepsilon > 0$ small enough so that $\frac{2}{\alpha}\varepsilon \le \delta$ and $\frac{1}{1+\frac{\varepsilon}{\alpha}} \ge 1-\delta$, we obtain the announced inclusions. $\square$

**Theorem 43.** If $\nabla \tilde{\rho}(0) = 0$, then the rescaled shape of the sandpile, $\sqrt{m-1} \cdot C$, where $C$ is defined in (18), converges to the ellipsoid

$$ \{s \in \mathbb{R}^d \mid 2s^\top \sigma^{-1} s \le 1\} $$

as $m$ tends to 1. The convergence is uniform in the sense of Definition 22 and consequently holds for the Hausdorff distance.

*Proof.* Using Propositions 32 and 55, $\sqrt{m-1} \cdot C$ is the dual curve of $\frac{1}{\sqrt{m-1}} \{\tilde{\rho} \le m\}$. The inclusions of Proposition 42 can be written as

$$ \sqrt{1-\delta} \left\{ s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le 1 \right\} \subset \frac{1}{\sqrt{m-1}} \left\{ t \in \mathbb{R}^d \mid \tilde{\rho}(t) \le m \right\} \subset \sqrt{1+\delta} \left\{ s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le 1 \right\}. $$

Therefore, Proposition 56, shows that $\sqrt{m-1} \cdot C$ converges to the set

$$ \left\{ s \in \mathbb{R}^d \mid \frac{1}{2} s^\top \sigma s \le 1 \right\}^* = \left\{ s \in \mathbb{R}^d \mid 2s^\top \sigma^{-1} s \le 1 \right\}, $$

according to the computation of duals of ellipsoids in Proposition 57. $\square$

## A Dual curves and dual convex sets

In this appendix, we remind the definition of dual curves and give the properties that are used through the article. We denote by $K$ a compact, convex subset of $\mathbb{R}^d$ and $C = \partial K$. Keep in mind that the cases we are interested in are the case where $K$ is a polytope and the smooth case, that is the case when $K$ has the form

$$ K = \{x \in \mathbb{R}^d \mid f(x) \le 1\} $$

for some smooth, convex function $f$ that takes some values strictly larger than 1.

**Definition 44.** Let $x \in C$. The normal cone to $K$ at $x$ is the set

$$ N_K(x) = \{n \in \mathbb{R}^d \mid \forall y \in K, n \cdot (y - x) \le 0\}. $$
**Proposition 45.** In the smooth case, the normal cone at a point $x \in C$ consists of the non-negative multiples of $\nabla f(x)$.

*Proof.* We first prove that $\nabla f(x) \in N_K(x)$. Let $y \in K$. If $\nabla f(x) \cdot (y - x) > 0$, then a Taylor expansion shows that for $\varepsilon > 0$ small enough, $f(x+\varepsilon(y-x)) > f(x) = 1$, therefore $x+\varepsilon(y-x) \notin K$, which contradicts the convexity of $K$. Therefore, $\nabla f(x) \cdot (y - x) \le 0$, i.e., $\nabla f(x) \in N_K(x)$, so its non-negative multiples also belong to the cone.

Conversely, let $n \in N_K(x) \setminus \{0\}$. By definition of the normal cone and of $K = f^{-1}((-\infty, 1])$, we have

$$g(x) = \max_{f(y)=1} g(y),$$

where $g(y) = n \cdot y$. Besides, $\nabla f(x) \ne 0$ because $f$ is convex and we assumed that $1$ is not the minimum of $f$. Therefore, the Lagrange multiplier theorem ensures that $\nabla g(x) = n$ is a multiple of $\nabla f(x)$, which we write $n = \lambda \nabla f(x)$. It remains to prove that $\lambda > 0$. We know that $\nabla f(x) \in N_K(x)$, so if $\lambda \le 0$, we would have $n \cdot (y - x) = \lambda \nabla f(x) \cdot (y - x) \ge 0$ and $n \cdot (y - x) \le 0$ for all $y \in K$, so $K \subset x + n^{\perp}$. But we assumed that $f$ takes values strictly smaller than $1$, so $K$ has a non-empty interior and cannot be a subset of a hyperplane. $\square$

**Lemma 46.** 1. Normal cones are closed subsets of $\mathbb{R}^d$.

2. For every $x \in C$, $N_K(x) \ne \{0\}$.

3. Let $u \in \mathbb{R}^d \setminus \{0\}$. There exists $x \in C$ such that $u \in N_K(x)$.

4. For every $t > 0$ and $x \in C$, $N_{tK}(tx) = N_K(x)$.

5. For every $U \in O_d(\mathbb{R})$ and $x \in C$, $N_{UK}(Ux) = UN_K(x)$.

*Proof.* Let $x \in C$. Then

$$N_K(x) = \bigcap_{y \in K} g_y^{-1}((-\infty, 0]),$$

where $g_y : n \in \mathbb{R}^d \mapsto n \cdot (y - x)$ is continuous. This shows that $N_K(x)$ is the intersection of closed subsets of $\mathbb{R}^d$, and thus it is a closed subset of $\mathbb{R}^d$.

We denote by $p_K$ the projection onto the closed convex set $K$, which is well defined and continuous according to the Hilbert projection theorem. The point $x$ belongs to $C = \partial K$, so there exists a sequence $(z_k)_{k \ge 0}$ of $\mathbb{R}^d \setminus K$ that tends to $x$. Let $u_k = \frac{z_k - p_K(z_k)}{\|z_k - p_K(z_k)\|}$. The Hilbert projection theorem ensures that

$$u_k \cdot (y - p_K(z_k)) \le 0 \text{ for every } y \in K. \qquad (28)$$

By compactness of the unit sphere, there is a subsequence of $(u_k)_{k \ge 0}$ that converges to a unit vector $u$. Besides, $p_K(z_k)$ tends to $p_K(x) = x$ because $p_K$ is continuous. Therefore, taking the subsequential limit of (28), we get $u \cdot (y - x) \le 0$ for all $y \in K$, hence $u \in N_K(x) \setminus \{0\}$.

Let $x \in K$ such that $u \cdot x = \max_{y \in K} u \cdot y$. Then $x \in C = \partial K$, because if $x \in \bar{K}$, then for $t > 0$ small enough, $x + tu \in K$ and $u \cdot (x + tu) = u \cdot x + t > u \cdot x$, which contradicts the maximality of $u \cdot x$. By maximality of $u \cdot x$, we have $u \in N_K(x)$.

We have

$$\begin{align*}
n \in N_{tK}(tx) &\iff \forall w \in tK, \quad n \cdot (w - tx) \le 0 \\
&\iff \forall y \in K, \quad n \cdot (ty - tx) \le 0 \\
&\iff \forall y \in K, \quad n \cdot (y - x) \le 0 \\
&\iff n \in N_K(x).
\end{align*}$$
Similarly

$$
\begin{align*}
n \in N_{UK}(Ux) &\iff \forall w \in UK, \quad n \cdot (w - Ux) \le 0 \\
&\iff \forall y \in K, \quad n \cdot (Uy - Ux) \le 0 \\
&\iff \forall y \in K, \quad U^{-1}n \cdot (y - x) \le 0 \\
&\iff U^{-1}n \in N_K(x) \\
&\iff n \in UN_K(x).
\end{align*}
\hfill $\square$
$$

**Definition 47.** Let $x \in C$. The tangent hyperplanes to $K$ at $x$ are the hyperplanes

$$H_n := \{y \in \mathbb{R}^d \mid n \cdot y = n \cdot x\}$$

for $n \in N_K(x) \setminus \{0\}$. If $n \cdot x > 0$, the normalized equation of $H_n$ is $\frac{n}{n \cdot x} \cdot y = 1$.

In the smooth case, we recover the usual definition of the tangent space.

**Lemma 48.** If 0 lies in the interior of $K$, then every tangent hyperplane has a normalized equation, that is:

$$\forall x \in C, \quad \forall n \in N_K(x) \setminus \{0\}, \quad n \cdot x > 0.$$

*Proof.* Let $x \in C$ and $n \in N_K(x) \setminus \{0\}$. For $t > 0$ small enough, $tn \in K$ because $0 \in \mathring{K}$, so we have $n \cdot (tn - x) \le 0$, that is, $n \cdot x \ge t\|n\|^2 > 0$. $\square$

**Definition 49.** We assume that $0 \in \mathring{K}$. Let $x \in C$. The dual points to $K$ at $x$ are the elements of

$$x^* = \left\{ \frac{n}{n \cdot x} \mid n \in N_K(x) \setminus \{0\} \right\}. \qquad (29)$$

In other words, it is the set of $n$ such that the normalized equations of tangent hyperplanes to $K$ at $x$ are the $\{n \cdot y = 1\}$, that is

$$x^* = \{z \in N_K(x) \mid z \cdot x = 1\}.$$

The set $C^* = \bigcup_{x \in C} x^*$ is the dual curve of $C$.

Let us do two remarks on Definition 49. First, in the definition (29) of $x^*$, one can consider only unit vectors $n$ of the normal cone, which we will often do in the proofs. Second, one can get rid of the condition $0 \in \mathring{K}$ to define the dual curve, but the dual curve is then defined in a projective space, which we do not need here.

From now on, we will always assume that $0 \in \mathring{K}$.

**Lemma 50.** The dual curve $C^*$ is compact.

*Proof.* Let $r > 0$ such that $\bar{B}(0, r) \subset K$. Let $z = \frac{u}{u \cdot x} \in C^*$, where $x \in C$ and $u \in N_K(x) \cap \mathbb{S}^{d-1}$. Then

$$u \cdot x = u \cdot (ru + x - ru) = r + u \cdot (x - ru).$$

But $ru \in K$ and $u \in N_K(x)$, so $u \cdot (x - ru) \ge 0$, hence $u \cdot x \ge r$, i.e., $\|z\| \le \frac{1}{r}$. We proved that $C^* \subset \bar{B}(0, \frac{1}{r})$, so $C^*$ is bounded.

Let $z \in \bar{C}^*$ and $(z_k)_{k \ge 0}$ a sequence of $C^*$ that tends to $z$. For all $k \ge 0$, there exist $x_k \in C$ and $u_k \in N_K(x_k) \cap \mathbb{S}^{d-1}$ such that $z_k = \frac{u_k}{u_k \cdot x_k}$. The sequence $(x_k, u_k)_{k \ge 0}$ lies in the compact set $C \times \mathbb{S}^{d-1}$, so it has a subsequential limit $(x, u) \in C \times \mathbb{S}^{d-1}$ and $z = \frac{u}{u \cdot x}$. Finally, taking the subsequential limit of the inequalities $u_k \cdot (y - x_k) \le 0$ for $y \in K$ leads to $u \cdot (y - x) \le 0$, hence $u \in N_K(x)$, and thus $z \in x^*$. Therefore, $z \in C^*$, so $C^*$ is closed.

In conclusion, $C^*$ is a bounded, closed subset of $\mathbb{R}^d$, so it is compact. $\square$
**Lemma 51.** Let $K$ be a compact, convex subset of $\mathbb{R}^d$ such that $0 \in \mathring{K}$. Let $z \in \partial K$. Then for all $t > 1$, $tz \notin K$.

*Proof.* Let $r > 0$ such that $B(0, r) \subset K$. We claim that if there exists $t > 1$ such that $tz \in K$, then $B(z, (1 - \frac{1}{t})r) \subset K$. Indeed, let $u \in \mathbb{S}^{d-1}$ and $0 \le \varepsilon < (1 - \frac{1}{t})r$. We have $z + \varepsilon u = \frac{1}{t}tz + (1 - \frac{1}{t}) \frac{\varepsilon}{1 - \frac{1}{t}} u$, which is a convex combination of $tz$ and $\frac{\varepsilon}{1 - \frac{1}{t}}u$, with $tz \in K$ and $\frac{\varepsilon}{1 - \frac{1}{t}}u \in K$, so $z + \varepsilon u \in K$, which proves $B(z, (1 - \frac{1}{t})r) \subset K$. This contradicts $z \in \partial K$. $\square$

**Proposition 52.** *Duality preserves convexity, that is, there exists a compact convex subset $K^*$ of $\mathbb{R}^d$ such that $0 \in \mathring{K}^*$ and $C^* = \partial(K^*)$. The set $K^*$ is called the dual set of $K$.*

*Proof.* We define $K^*$ as the convex hull of $C^*$. Since $C^*$ is a compact subset of $\mathbb{R}^d$ according to Proposition 50, its convex hull $K^*$ is compact.

Let $u \in \mathbb{S}^{d-1}$. By item 3 of Lemma 46, there exists $x \in C$ such that $u \in N_K(x)$. Therefore, $C^*$ contains a positive multiple $r_u u$ of $u$. Applying this to the vectors $(\pm e_i)_{1 \le i \le d}$, where $(e_1, \dots, e_d)$ is the canonical basis of $\mathbb{R}^d$, we obtain that $K^*$ contains the convex hull of $\{\pm r_{\pm e_i} e_i\}$, which contains the ball centered in 0 of radius $\min_i r_{\pm e_i}$ for the norm $\| \cdot \|_1$. Therefore, 0 lies in the interior of $K^*$.

Let $z \in \partial(K^*)$. We prove that $z \in C^*$. The preceding result shows that $z \ne 0$. Therefore, according to item 3 of Lemma 46, there exists $x \in C$ such that $z \in N_K(x)$. To get $z \in C^*$, all that remains to be proved is that $z \cdot x = 1$. We start by proving that $z \cdot x \le 1$. By definition of $K^*$, we can write $z = \sum_{i=0}^d \lambda_i z_i$ where $\lambda_0, \dots, \lambda_d \ge 0$, $\sum_{i=0}^d \lambda_i = 1$ and $z_0, \dots, z_d \in C^*$. Writing each of the $z_i$ as $z_i = \frac{u_i}{u_i \cdot x_i}$ where $x_i \in C$ and $u_i \in N_K(x_i)$, we get

$$z = \sum_{i=0}^{d} \lambda_i \frac{u_i}{u_i \cdot x_i}.$$

But $u_i \in N_K(x_i)$ and $x \in K$, so $u_i \cdot (x - x_i) \le 0$, i.e., $\frac{u_i \cdot x}{u_i \cdot x_i} \le 1$. Summing up, we get

$$z \cdot x = \sum_{i=0}^{d} \lambda_i \frac{u_i \cdot x}{u_i \cdot x_i} \le \sum_{i=0}^{d} \lambda_i = 1.$$

Let us move to the second inequality, $z \cdot x \ge 1$. Assume by contradiction that $z \cdot x < 1$. Then $\frac{z}{z \cdot x} \in C^* \subset K^*$. This contradicts Lemma 51 with $K = K^*$ and $t = \frac{1}{z \cdot x}$. In conclusion, $z \cdot x = 1$ so $z \in C^*$, which proves $\partial(K^*) \subset C^*$.

Conversely, let $z \in C^*$. It means that there exists $x \in C$ such that $z \in N_K(x)$ and

$$z \cdot x = 1. \qquad (30)$$

We know that $z \in K^*$. If $z \notin \partial(K^*)$, then there exists $t > 1$ such that $tz \in \partial(K^*)$. Using the first inclusion, we get $tz \in C^*$, which means that there exists $\tilde{x} \in C$ such that $tz \in N_K(\tilde{x})$ and $tz \cdot \tilde{x} = 1$. Because $x \in K$ and $tz \in N_K(\tilde{x})$, we have $tz \cdot (x - \tilde{x}) \le 0$, that is $tz \cdot x \le tz \cdot \tilde{x} = 1$. Therefore, $z \cdot x \le \frac{1}{t} < 1$, which contradicts (30). In conclusion, $z \in \partial(K^*)$ and the proof is complete. $\square$

The dual set $K^*$ satisfies the same assumptions as $K$, so $C^*$ has a dual curve $(C^*)^*$, which we will simply write $C^{**}$.

**Proposition 53.** *Duality is an involution, i.e., $C^{**} = C$.*
**Proof.** Since $K$ is a compact, convex set such that $0 \in \check{K}$, Lemma 51 ensures that for every $u \in \mathbb{S}^{d-1}$, there is a unique $r_u > 0$ such that $r_u u \in \partial K = \mathcal{C}$, and $\mathcal{C} = \{r_u u \mid u \in \mathbb{S}^{d-1}\}$. Similarly, for every $u \in \mathbb{S}^{d-1}$, there is a unique $r_u^{**} > 0$ such that $r_u^{**} u \in \mathcal{C}^{**}$, and $\mathcal{C}^{**} = \{r_u^{**} u \mid u \in \mathbb{S}^{d-1}\}$. Therefore, if we prove an inclusion between $\mathcal{C}^{**}$ and $\mathcal{C}$, we will have $r_u = r_u^{**}$ for every $u \in \mathbb{S}^{d-1}$ and thus equality between the two curves.

Let $x \in \mathcal{C}$. Let $z \in \mathcal{C}^*$ such that $z \in N_K(x)$ and $z \cdot x = 1$ (such a $z$ exists: we can choose $w \in N_K(x) \setminus \{0\}$ thanks to Lemma 46 and since $w \cdot x > 0$ according to Lemma 48, we only need to set $z = \frac{w}{w \cdot x}$). To get $x \in \mathcal{C}^{**}$, it remains to show that $x \in N_{K^*}(z)$, that is

$$ \forall \tilde{z} \in K^*, \quad x \cdot \tilde{z} \le x \cdot z = 1. \qquad (31) $$

Every $\tilde{z} \in K^*$ is a convex combination of elements of $\mathcal{C}^*$, so it is enough to prove (31) when $\tilde{z} \in \mathcal{C}^*$. In that case, we can write $\tilde{z} = \frac{u}{u \cdot y}$ with $y \in \mathcal{C}$ and $u \in N_K(y) \setminus \{0\}$, so $x \cdot \tilde{z} = \frac{u \cdot x}{u \cdot y}$. But $x \in K$ and $u \in N_K(y)$, so $u \cdot x \le u \cdot y$, i.e., $\frac{u \cdot x}{u \cdot y} \le 1$, i.e., $x \cdot \tilde{z} \le 1$. In conclusion, we showed (31), hence $x \in \mathcal{C}^{**}$. As mentioned at the beginning of the proof, the inclusion $\mathcal{C} \subset \mathcal{C}^{**}$ is enough to conclude. $\square$

**Proposition 54.** Let $K_1, K_2$ be two compact, convex subsets of $\mathbb{R}^d$, with $0 \in \check{K}_1 \cap \check{K}_2$, such that $K_1 \subset K_2$. Then $K_2^* \subset K_1^*$.

**Proof.** It is enough to prove that $C_2^* \subset K_1^*$. Let $z_2 = \frac{u_2}{u_2 \cdot x_2} \in C_2^*$, where $x_2 \in \mathcal{C}_2$ and $u_2 \in N_{K_2}(x_2) \cap \mathbb{S}^{d-1}$. By item 3 of Lemma 46, there exists $x_1 \in \mathcal{C}_1$ such that $u_2 \in N_{K_1}(x_1)$, hence $\frac{u_2}{u_2 \cdot x_1} \in K_1^*$. But $x_1 \in K_1 \subset K_2$ and $u_2 \in N_{K_2}(x_2)$, so $u_2 \cdot (x_1 - x_2) \le 0$, i.e., $u_2 \cdot x_1 \le u_2 \cdot x_2$. Therefore

$$ z_2 = \frac{u_2}{u_2 \cdot x_2} \in \left[ 0, \frac{u_2}{u_2 \cdot x_1} \right] \subset K_1^*. \qquad \square $$

**Proposition 55.** Let $K$ be a compact, convex subset of $\mathbb{R}^d$ such that $0 \in \check{K}$, $t > 0$ and $U \in O_d(\mathbb{R})$. Then $(tK)^* = \frac{1}{t}K^*$ and $(UK)^* = U(K^*)$.

**Proof.** Let $z \in \mathbb{R}^d$. We have

$$ \begin{align*} z \in (tK)^* &\iff \exists x \in \mathcal{C}, \quad \exists u \in N_{tK}(tx), \quad z = \frac{u}{u \cdot (tx)} \\ &\iff \exists x \in \mathcal{C}, \quad \exists u \in N_K(x), \quad z = \frac{u}{u \cdot (tx)} = \frac{1}{t} \frac{u}{u \cdot x} \\ &\iff z \in \frac{1}{t} K^*, \end{align*} $$

where the second equivalence used the equality $N_{tK}(tx) = N_K(x)$ of Lemma 46. Similarly,

$$ \begin{align*} z \in (UK)^* &\iff \exists x \in \mathcal{C}, \quad \exists u \in N_{UK}(Ux), \quad z = \frac{u}{u \cdot Ux} \\ &\iff \exists x \in \mathcal{C}, \quad \exists v \in N_K(x), \quad z = \frac{Uv}{Uv \cdot Ux} = U \frac{v}{v \cdot x} \\ &\iff z \in U(K^*), \end{align*} $$

where the second equivalence uses the equality $N_{UK}(Ux) = UN_K(x)$ of Lemma 46. $\square$

**Proposition 56.** Let $K$ and $(K_k)_{k \ge 0}$ be compact, convex subsets of $\mathbb{R}^d$ such that $0 \in \check{K}$ and $0 \in \check{K}_k$ for all $k \ge 0$. We assume that there exist sequences $(\varepsilon_k^-)_{k \ge 0}$ and $(\varepsilon_k^+)_{k \ge 0}$ of positive numbers that tend to 0 such that for every $k \ge 0$,

$$ (1 - \varepsilon_k^-) K \subset K_k \subset (1 + \varepsilon_k^+) K. $$
Then, uniformly in the sense of Definition 22, and consequently for the Hausdorff distance,

$$K_k \xrightarrow{k \to +\infty} K \quad \text{and} \quad K_k^* \xrightarrow{k \to +\infty} K^*.$$

*Proof.* Let us first prove the result for $K_k$. The assumptions on the convex sets imply that there is a bounded family $(r_u)_{u \in \mathbb{S}^{d-1}}$ such that $K = \{ru | u \in \mathbb{S}^{d-1}, 0 \le r \le r_u\}$. Therefore, using the notations of Definition 22, the inclusions $(1 - \varepsilon_k^-) K \subset K_k \subset (1 + \varepsilon_k^+) K$ allow us to set

$$r_{u,k} = (1 - \varepsilon_k^-) r_u \quad \text{and} \quad R_{u,k} = (1 + \varepsilon_k^+) r_u.$$

Therefore,

$$\|r_{u,k} - r_u\|_{\infty} = \varepsilon_k^- \sup_{u \in \mathbb{S}^{d-1}} r_u \xrightarrow{k \to +\infty} 0 \quad \text{and} \quad \|R_{u,k} - r_u\|_{\infty} = \varepsilon_k^+ \sup_{u \in \mathbb{S}^{d-1}} r_u \xrightarrow{k \to +\infty} 0,$$

hence the uniform convergence of $K_k$.

For the duals, we use Propositions 54 and 55 to obtain

$$\frac{1}{1 + \varepsilon_k^+} K^* \subset K_k^* \subset \frac{1}{1 - \varepsilon_k^-} K^*$$

and use the first case with $\tilde{\varepsilon}_k^- = 1 - \frac{1}{1+\varepsilon_k^+}$ and $\tilde{\varepsilon}_k^+ = \frac{1}{1-\varepsilon_k^-} - 1$. $\square$

**Proposition 57 (Ellipsoid).** Let $\sigma$ be a symmetric, positive-definite matrix and

$$K_{\sigma} = \{x \in \mathbb{R}^d | x^{\top}\sigma x \le 1\}.$$

Then $K_{\sigma}^* = K_{\sigma^{-1}}$.

*Proof.* When $\sigma$ is diagonal, it is a simple computation. In the general case, we write $\sigma = U^{\top}DU$ where $D$ is diagonal and $U \in O_d(\mathbb{R})$. Then

$$K_{\sigma} = \{x \in \mathbb{R}^d | (Ux)^{\top} D (Ux) \le 1\} = U^{-1}K_D, \quad (32)$$

so Proposition 55 yields

$$K_{\sigma}^* = U^{-1}(K_D^*).$$

Using the diagonal case and (32) with $\sigma^{-1} = U^T D^{-1} U$ instead of $\sigma$, we finally get

$$K_{\sigma}^* = U^{-1}(K_D^*) = U^{-1}K_{D^{-1}} = K_{\sigma^{-1}}. \quad \square$$

**Proposition 58 (Polytope).** Let $x_1, \dots, x_n \in \mathbb{R}^d$ and

$$K = \bigcap_{i \in \{1, \dots, n\}} \{s \in \mathbb{R}^d | s \cdot x_i \le 1\}.$$

The set $K$ is compact if and only if 0 lies in the interior of the convex hull of $\{x_1, \dots, x_n\}$. In that case, this convex hull is the dual convex $K^*$ of $K$.

*Proof.* The convex $K$ is clearly a closed subset of $\mathbb{R}^d$, so it is compact if and only if it is bounded, which comes down to saying that for all $u \in \mathbb{S}^{d-1}$, there exists $i \in \{1, \dots, n\}$ such that $u \cdot x_i > 0$. If 0 lies in the interior of the convex hull of $\{x_1, \dots, x_n\}$, then for $\varepsilon > 0$ small enough, $\varepsilon u$ can be written as a convex combination of $\{x_1, \dots, x_n\}$. Since $u \cdot \varepsilon u > 0$, there has to be at least one $x_i$
such that $u \cdot x_i > 0$, so $K$ is bounded. Conversely, if $K$ is compact, we will obtain that 0 lies in the interior of the convex hull of $\{x_1, \dots, x_n\}$ as a consequence of the computation of $K^*$ and of Proposition 52, which ensures that 0 lies in the interior of $K^*$.

We assume that $K$ is compact. To compute $K^*$, we first note that

$$K = \bigcap_{x \in \text{extr}(\{x_1, \dots, x_n\})} \{s \in \mathbb{R}^d \mid s \cdot x \le 1\},$$

because if one of the $x_i$ is a convex combination of the others, then the condition $s \cdot x_i \le 1$ is redundant. The faces of the polytope $K$ are thus delimited by the hyperplanes $H_x = \{s \in \mathbb{R}^d \mid s \cdot x = 1\}$ for $x \in \text{extr}(\{x_1, \dots, x_n\})$. At a point $y$ that lies in the interior of the face delimited by $H_x$, the normal cone is $\mathbb{R}_{\ge 0}x$, so the dual point associated with $y$ is $x$. In particular, $K^*$ contains $\text{extr}(\{x_1, \dots, x_n\})$. Similarly, if $y$ is at the intersection of the hyperplanes $H_{x_i}$, $i \in I$, then the normal cone at $y$ consists of convex combinations of the $x_i$, $i \in I$, so $K^*$ is a subset of the convex hull of $\text{extr}(\{x_1, \dots, x_n\})$. In conclusion, $K^*$ is a subset of the convex hull of $\text{extr}(\{x_1, \dots, x_n\})$ and it contains $\text{extr}(\{x_1, \dots, x_n\})$. By convexity of $K^*$, it means that $K^*$ is the convex hull of $\text{extr}(\{x_1, \dots, x_n\})$, that is the convex hull of $\{x_1, \dots, x_n\}$. $\square$

## B An interpretation in terms of amoebas

In this appendix, we assume that $\mu$ has finite support. We give a geometric interpretation of the Doob transform and the set $\partial T$ with tools from real algebraic geometry.

A quantity directly related to the real Laplace transform $L\mu$ is the Fourier transform of $\mu$, denoted by $\mathcal{F}\mu: (\mathbb{C}^*)^d \to \mathbb{C}^*$, defined by:

$$\forall z = (z_1, \dots, z_d) \in (\mathbb{C}^*)^d, \quad \mathcal{F}\mu(z) := \sum_{x \in \mathbb{Z}^d} z^x \mu_{i,j}(x),$$

where $z^x$ stands for $\prod_{i=1}^d z_i^{x_i}$. With the convention that $e^t = (e^{t_1}, \dots, e^{t_d})$ for $t \in \mathbb{R}^d$, we have the following relation:

$$\mathcal{F}\mu(e^t) = L\mu(t).$$

The Fourier transform $\mathcal{F}\mu$ is a $d \times d$ matrix, whose entries are Laurent polynomials in $z_1, \dots, z_d$. The characteristic polynomial of the random walk is

$$P(z) = \det(I_d - \mathcal{F}\mu(z)),$$

and its spectral hypersurface is the algebraic variety in $(\mathbb{C}^*)^d$ consisting of the zero set of $P$:

$$\mathcal{C} = \{z \in (\mathbb{C}^*)^d \mid P(z) = 0\}.$$

In the special case when $d = 2$, the spectral hypersurface is an algebraic curve, and we rather use the term spectral curve. This terminology is directly borrowed from dimer models [21] and spanning trees [20], but has been used for much longer in exactly solvable models in statistical mechanics, and integrable systems.

An interesting object to apprehend the geometry of this spectral hypersurface is its so-colled amoeba $A(\mathcal{C})$, defined as:

$$A(\mathcal{C}) = \{\text{Log}(z) ; z \in \mathcal{C}\} \subset \mathbb{R}^d,$$

where $\text{Log}(z) = (\log |z_1|, \dots, \log |z_d|)$. See for example [14, Chap. 6] for more details.
Gershgorin's theorem guarantees that when $|z| = (|z_1|, \dots, |z_d|) = (1, \dots, 1)$, $\mathcal{F}\mu(z)$ has no eigenvalue near 1. This implies that the point $(0, \dots, 0)$ is in the complement of $A(\mathcal{C})$.

When $d=2$, and the graph of allowed transitions between vertices for the (killed) walk on $\mathbb{Z}^2 \times \{1, \dots, p\}$ can be realized as a planar, biperiodic graph, then the spectral curve is a Harnack curve [27, 20]. In particular, the boundary of the amoeba is given by the image by Log of the real locus of $\mathcal{C}$, that is, the points $z = (z_1, \dots, z_d)$ of $\mathcal{C}$ for which all $z_i$ are real. This is not true in general. However, the Perron-Frobenius theorem implies the following result:

**Lemma 59.** The boundary of the connected component of the complement of $A(\mathcal{C})$ containing the origin is the image by Log of a connected component of $\mathcal{C}^+$, the real positive locus of $\mathcal{C}$:

$$ \mathcal{C}^+ = \{z \in (\mathbb{R}_{>0})^d ; P(z) = 0\}. $$

*Proof.* By general results of real algebraic geometry [27], each connected component of the complement of the amoeba are convex. Let $r = (r_1, \dots, r_d)$ a point on the boundary of that connected component. Then the whole interval $[0, r)$ is contained in the connected component of the complement. Define $M(s) = \mathcal{F}\mu(e^{sr_1}, \dots, e^{sr_d})$ for $s \in [0, 1]$. This matrix has non-negative entries, so Perron-Frobenius theorem applies to $M(s)$, which should have spectral radius given by its largest positive eigenvalue $\lambda(s)$. The point $s=1$ is the first point where $\det(I_d - \mathcal{F}\mu(z)) = 0$, for $|z_j| = e^{sr_j}$, and the maximum principle implies that the only $z$ with this constraint satisfying this equation is $z = e^s$. The eigenvalue which is then equal to one must be the largest one, otherwise, by the intermediate value theorem, the Perron-Frobenius eigenvalue would have reached one for some $s < 1$, which would be in contradiction with the definition of $r$. $\square$

This allows one to make a connection between the geometry of the amoeba $A(\mathcal{C})$ and objects introduced earlier.

* The closure of the connected component of the complement of $A(\mathcal{C})$ containing the origin is the set $\mathcal{T}$ introduced in Definition 7.

* The fact that starting from the origin, every direction $u \in \mathbb{S}^{d-1}$ corresponds bijectively to a point of $\mathcal{T}$ is a reformulation of the Ney-Spitzer homeomorphism in this context (see [8] for related ideas).

* For every $t \in \partial\mathcal{T}$, the eigenvector $\varphi_t$ is a non-zero vector in the kernel of the matrix $I_d - \mathcal{F}\mu(z)$ for $z = (e^{t_1}, \dots, e^{t_d}) = e^t$.

* Performing the Doob transform with $\varphi_t$ to obtain $\mu_t$ implies that for any $w \in (\mathcal{C}^*)^d$,

$$ \mathcal{F}\mu_t(w) = \Phi_t^{-1} \cdot \mathcal{F}\mu_t(e^t w) \cdot \Phi_t, $$

where $\Phi_t$ is the diagonal matrix whose diagonal entries are the values of $\varphi_t$, and the product $e^t w$ is understood coordinatewise.

* In particular the last point implies that the amoeba $A(\mathcal{C}_t)$ for $\mu_t$ is obtained from the one for $\mu$ by a translation of vector $-t$, which means that the origin is now on the boundary of $A(\mathcal{C}_t)$.

* The normal to $\mathcal{T}$ at $t$, which is parallel to $\nabla\rho(t)$, can be interpreted as the direction of the drift for the conditioned walk, as noted in [6, Def. 1.3, Prop. 2.16].
## Acknowledgments

TB, CB and KR are supported by the project DIMERS (ANR-18-CE40-0033). TB and KR are supported by the project RAWABRANCH (ANR-23-CE40-0008), funded by the French National Research Agency, and by the France 2030 program Centre Henri Lebesgue (ANR-11-LABX-0020-01). SM was partly supported by Simons Foundation Grant No. 422190. Part of this research was performed while CB was visiting the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation (Grant No. DMS-1925919). TB and KR thank the VIASM (Hanoï, Vietnam) for their hospitality and wonderful working conditions. For part of this research SM was visiting Yerevan State University through the support of the Fulbright Scholar program, and would like to thank them for their hospitality. The authors would like to thank Ian Alevy and Sanjay Ramassamy for very interesting discussions at an early stage of the project.

## References

[1] Hayk Aleksanyan and Henrik Shahgholian. Discrete Balayage and boundary sandpile. *J. Anal. Math.*, 138(1):361–403, 2019.

[2] Ian Alevy and Sevak Mkrtchyan. The limit shape of the leaky abelian sandpile model. *Int. Math. Res. Not.*, 2022(16):12767–12802, 2022.

[3] Yacine Aoun, Sébastien Ott, and Yvan Velenik. Ornstein-Zernike behavior for Ising models with infinite-range interactions. *Ann. Inst. Henri Poincaré Probab. Stat.*, 60(1):167–207, 2024.

[4] Søren Asmussen. *Applied probability and queues.*, volume 51 of *Appl. Math. (N. Y.)*. New York, NY: Springer, 2nd revised and extended ed. edition, 2003.

[5] Per Bak, Chao Tang, and Kurt Wiesenfeld. Self-organized criticality: An explanation of the 1/f noise. *Phys. Rev. Lett.*, 59:381–384, Jul 1987.

[6] Théo Ballu. Asymptotics of Green functions for Markov-additive processes: an approach via dyadic splitting of integrals. Preprint, arXiv:2407.10685 [math.PR] (2024), 2024.

[7] A. A. Borovkov and A. A. Mogul'skij. The second rate function and the asymptotic problems of renewal and hitting the boundary for multidimensional random walks. *Sib. Math. J.*, 37(4):647–682, 1996.

[8] Cédric Boutillier and Kilian Raschel. Martin boundary of killed random walks on isoradial graphs. *Potential Anal.*, 57(2):201–226, 2022.

[9] Raphaël Cerf. *The Wulff crystal in Ising and percolation models. École d'Été de Probabilités de Saint-Flour XXXIV – 2004*, volume 1878 of *Lect. Notes Math*. Berlin: Springer, 2006.

[10] Qipin Chen, Shane Chern, and Lin Jiu. Multi-headed lattices and Green functions. *J. Phys. A, Math. Theor.*, 57(46):21, 2024. Id/No 465204.

[11] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth. On the Lambert $w$ function. *Adv. Comput. Math.*, 5(4):329–359, 1996.

[12] R. A. Doney. An analogue of the renewal theorem in higher dimensions. *Proc. Lond. Math. Soc. (3)*, 16:669–684, 1966.
[13] Matthieu Dussaule. The Martin boundary of a free product of abelian groups. *Ann. Inst. Fourier*, 70(1):313-373, 2020.

[14] Israel M. Gelfand, Mikhail M. Kapranov, and Andrei V. Zelevinsky. *Discriminants, Resultants, and Multidimensional Determinants*. Birkhäuser Boston, 1994.

[15] A. J. Guttmann and T. Prellberg. Staircase polygons, elliptic integrals, heun functions, and lattice green functions. *Phys. Rev. E*, 47:R2233-R2236, Apr 1993.

[16] Anthony J. Guttmann. Lattice Green's functions in all dimensions. *J. Phys. A, Math. Theor.*, 43(30):26, 2010. Id/No 305205.

[17] S. Hassani, Ch. Koutschan, J.-M. Maillard, and N. Zenine. Lattice Green functions: the $d$-dimensional face-centered cubic lattice, $d = 8, 9, 10, 11, 12$. *J. Phys. A, Math. Theor.*, 49(16):30, 2016. Id/No 164003.

[18] P-L. Hennequin. Processus de Markoff en cascade. *Ann. Inst. Henri Poincaré*, 18:109-193, 1963.

[19] G. S. Joyce. On the simple cubic lattice Green function. *Philos. Trans. R. Soc. Lond., Ser. A*, 273:583-610, 1973.

[20] Richard Kenyon. Spanning forests and the vector bundle Laplacian. *Ann. Probab.*, 39(5):1983-2017, 2011.

[21] Richard Kenyon, Andrei Okounkov, and Scott Sheffield. Dimers and amoebae. *Ann. of Math. (2)*, 163(3):1019-1056, 2006.

[22] J. F. C. Kingman. A convexity property of positive matrices. *Quart. J. Math. Oxford Ser. (2)*, 12:283-284, 1961.

[23] Lionel Levine, Wesley Pegden, and Charles K. Smart. Apollonian structure in the Abelian sandpile. *Geom. Funct. Anal.*, 26(1):306-336, 2016.

[24] Lionel Levine and Yuval Peres. Strong spherical asymptotics for rotor-router aggregation and the divisible sandpile. *Potential Anal.*, 30(1):1-27, 2009.

[25] R. J. Messikh. The surface tension near criticality of the 2d-Ising model. Preprint, arXiv:math/0610636 [math.PR] (2006), 2006.

[26] Emmanuel Michta and Gordon Slade. Asymptotic behaviour of the lattice Green function. *ALEA Lat. Am. J. Probab. Math. Stat.*, 19(1):957-981, 2022.

[27] G. Mikhalkin. Real algebraic curves, the moment map and amoebas. *Annals of Mathematics-Second Series*, 151(1):309-326, 2000.

[28] P. Ney and F. Spitzer. The Martin boundary for random walk. *Trans. Am. Math. Soc.*, 121:116-132, 1966.

[29] Wesley Pegden and Charles K. Smart. Convergence of the Abelian sandpile. *Duke Math. J.*, 162(4):627-642, 2013.

[30] Wesley Pegden and Charles K. Smart. Stability of patterns in the Abelian sandpile. *Ann. Henri Poincaré*, 21(4):1383-1399, 2020.
[31] Robin Pemantle and Mark C. Wilson. *Analytic combinatorics in several variables.*, volume 140 of *Camb. Stud. Adv. Math*. Cambridge: Cambridge University Press, 2013.

[32] Lucas Rey. The Doob transform and the tree behind the forest, with application to near-critical dimers. Preprint, arXiv:2401.13599 [math.PR] (2024), 2024.

[33] E. Seneta. *Non-negative matrices and Markov chains*. Springer Series in Statistics. Springer, New York, 2006. Revised reprint of the second (1981) edition [Springer-Verlag, New York; MR0719544].

[34] Wolfgang Woess. *Random walks on infinite graphs and groups*, volume 138 of *Camb. Tracts Math*. Cambridge: Cambridge University Press, 2000.