A stochastic variant of the abelian sandpile model

TL;DR

This paper introduces a stochastic variant of the abelian sandpile model (SSP), exploring its mathematical properties, physical behavior, and potential as a simplified model for understanding the LLL algorithm in lattice cryptography.

Contribution

The paper develops a basic theoretical framework for SSP, a stochastic extension of ASM, and investigates its properties and potential applications in cryptography.

Findings

SSP shares mathematical properties with ASM

SSP exhibits different physical behavior and steady state shape

Numerical study of SSP's behavior supports its theoretical analysis

Abstract

We introduce a natural stochastic extension, called SSP, of the abelian sandpile model(ASM), which shares many mathematical properties with ASM, yet radically differs in its physical behavior, for example in terms of the shape of the steady state and of the avalanche size distribution. We establish a basic theory of SSP analogous to that of ASM, and present a brief numerical study of its behavior. Our original motivation for studying SSP stems from its connection to the LLL algorithm established in another work by the authors [5]. The importance of understanding how LLL works cannot be stressed more, especially from the point of view of lattice-based cryptography. We believe SSP serves as a tractable toy model of LLL that would help further our understanding of it.

Equations

$$\log \text{mbox}(RHF)=\frac{1}{(L+1{)}^2}\sum _{i=1}^L(L+1-i)r_i,$$

$$T(c):=T_{\mathcal{A}(c)}.$$

$$T(c):=\sum _{v\in V\backslash \{s\}}\mathcal{A}(c)(v)T_v(c).$$

$$c(v)-\sum _{i=1}^{k(v)}{\gamma }_{v,i}\cdot \mathrm{deg}(v)+\sum _{w\sim v}\sum _{i=1}^{k(w)}{\gamma }_{w,i}\mathrm{deg}(w,v)<T,$$

$$c(v)-\sum _{i=1}^{r(v)}{\gamma }_{v,i}\cdot \mathrm{deg}(v)+\sum _{w\sim v}\sum _{i=1}^{r(w)}{\gamma }_{w,i}\mathrm{deg}(v,w)<T.$$

$$(C^{\circ }+D{)}^{\circ }=(C+D{)}^{\circ },$$

$$C^{\circ }\oplus D=C\oplus D.$$

$$\mu (\{T_1:T_1(k(v))=S(k(v))\})\cdot \mu (\{T_2:T_2(l(v)-k(v))=S^{\prime }(l(v)-k(v))\})$$

$$=\mu (\{R:R(l(v))=S(l(v))\}),$$

$$C^{\circ }\oplus D=\sum p_i{q}_j(c_i^{\circ }\oplus d_j)=\sum p_i{q}_j(c_i\oplus d_j)=C\oplus D,$$

$$S\oplus C=S\mbox{ for all $C\in M$.}$$

$$A=\left[\begin{array}{cccccc}p(1)& 1& 0& 0& \dots & 0\\ p(2)& 0& 1& 0& \dots & 0\\ p(3)& 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ p(n)& 0& 0& \dots & \dots & 0\end{array}\right]$$

$$\left[\begin{array}{c}r(1)\\ r(2)\\ ⋮\\ r(n)\end{array}\right]$$

$${\chi }_A(\lambda )={\lambda }^n-p(1){\lambda }^{n-1}-\dots -p(n-1)\lambda -p(n).$$

$${\chi }_A^{\prime }(\lambda )=n{\lambda }^{n-1}-(n-1)p(1){\lambda }^{n-2}-\dots -p(n-1)$$

$$C={\left(\sum _{i=1}^n\sum _{j=i}^n\wp (j)\right)}^{-1}.$$

$$S(c):=\sum _{k_1=0}^{n-1}\cdots \sum _{k_L=0}^{n-1}\left(\prod _{j=1}^{L}r(k_j){\tau }_j^{k_j}\right)[c].$$

$$\sum _{k_1=0}^{n-1}\cdots \sum _{k_L=0}^{n-1}\left(\prod _{j=1}^{L}r(k_j){\tau }_j^{k_j}\right)T_{v_i}^{I(k_i=0)}[c]$$

$$=\sum _{k_1=0}^{n-1}\cdots \sum _{k_L=0}^{n-1}\left(\prod _{j=1}^{L}r(k_j){\tau }_j^{k_j}\right){\tau }_i[c]$$

$$=S({\tau }_{i}c).$$

$$S(R):=\frac{1}{|R|}\sum _{c\in R}S(c).$$

$$(S(R)+v_i{)}^{\circ }={\left(\frac{1}{|R|}\sum _{c^{\prime }\in R+v_i}S(c^{\prime })\right)}^{\circ }={\left(\frac{1}{|R|}\sum _{c^{\prime }\in R+v_i}S(c^{\prime }{)}^{\circ }\right)}^{\circ }={\left(\frac{1}{|R|}\sum _{c\in R}S(c{)}^{\circ }\right)}^{\circ }=S(R{)}^{\circ },$$

$$\left(\prod _{i=1}^L{T}_{v_i}\right)c=\sum _{k_1=0}^{n-1}\cdots \sum _{k_L=0}\left(\prod _{j=1}^L\wp (k_j){\tau }_j^{k_j}\right)[c].$$

$$a>\frac{1}{(L+1{)}^2}\sum _{i=1}^L(L+1-i)s_i.$$

$$\left(0,\dots ,\frac{a(L+1{)}^2}{L+1-i},\dots ,0\right)$$

$$\frac{a^L(L+1{)}^{2L}}{(L!{)}^2},$$

$$a<C^{-1}/e^2.$$

Full text

Our original motivation for studying SSP stems from its connection to the LLL algorithm established in another work by the authors [5]. The importance of understanding how LLL works cannot be stressed more, especially from the point of view of lattice-based cryptography. We believe SSP serves as a tractable toy model of LLL that would help further our understanding of it.

1. Introduction

1.1. Overview

Let us start by recalling ASM ([1], [3], also see [4]) on a one-dimensional lattice. Let $A_{L}$ be the cycle graph whose vertices $V(A_{L})$ consist of $L+1$ elements, say $V(A_{L})=\{v_{1},\ldots,v_{L+1}\}$, and the (undirected) edges are defined by $E(A_{L})=\{(v_{i},v_{i+1}):i=1,\ldots,L\}\cup\{(v_{L+1},v_{1})\}$. We designate $v_{L+1}$ as the sink. Fix two positive integers $T$ and $I$, preferably $T\geq 2I$.

We temporarily refer to $I$ as the “toppling strength” and $T$ as the threshold. That is, toppling at site $v_{i}$ is to subtract $2I$ grains from the pile on $v_{i}$ and add $I$ grains to each of its neighbors. A configuration on $A_{L}$ is called stable if and only if all non-sink vertices have $<T$ grains. Fix another positive integer $J$ coprime to $I$, preferably $J\approx I$. One then considers the following Markov chain on the space of the stable configurations of ASM: given a stable configuration on $A_{L}$, obtain the next stable configuration by first adding $J$ grains to any randomly chosen non-sink site, and then stabilizing the resulting configuration.111The only reason to impose the condition $(J,I)=1$ is to prevent the pile heights at each site from being concentrated on a select few congruence classes modulo $I$; it is not so much an essential condition as a cosmetic one.

On the other hand, consider the following simple stochastic variant of ASM, which we named SSP. SSP is exactly the same as ASM except for the toppling procedure: first one samples $\gamma$ uniformly from $\{1,2,\ldots,2I\}$ (so that $\gamma$ has mean $\approx I$) and then subtract $2\gamma$ grains from site $v_{i}$ and add $\gamma$ to each of the two neighboring sites. One can of course associate a Markov chain to SSP analogously to the above discussion: add $I$ grains to a random site and then stabilize.

It is rather surprising that this simple and natural extension of ASM has never been considered in the literature so far. In fact, our motivation for studying SSP does not come from physics, but from a study of the practical behavior of the LLL algorithm [8]. The LLL algorithm is one of the fundamental tools in computational mathematics used for lattice basis reduction, with numerous applications to number theory and cryptography — see the book [10] for an extensive treatment on the subject. Despite its celebrated status, much of its behavior in practice has been shrouded in complete mystery. For example, the most well-known problem concerns the quantity called the root Hermite factor(RHF), defined as

[TABLE]

where $r_{i}$ is the pile height at $v_{i}$, in sandpile language. It is a theorem from [8] that the RHF of LLL has a sharp upper bound by $(4/3)^{1/4}\approx 1.075$, but empirically one observes RHF $\approx 1.02$ most of the time — see [9] for details. This phenomenon has been well-known since the birth of LLL in 1982, yet there has been not even a vague heuristic as to why this must happen.

Recently, we argued in [5] that LLL behaves like a sandpile model, and that this idea explains much of its practical behavior all at once, in particular this RHF problem. The remaining difficulty is that LLL as a sandpile model is nonabelian, which is difficult to study analytically. Hence we invented SSP, the closest abelian version of LLL, to be used as a toy model.

Figure 1 compares the empirical average “steady state density” of LLL and SSP. The resemblance therein should come as surprising, since the two algorithms originate from very different fields of study, lattice reduction and statistical physics, respectively. We hope that the present paper helps build a bridge between them, to the benefit of both areas.

1.2. In this paper

In the next section, we develop a basic mathematical theory of SSP. As in the case of ASM, we impose a monoid structure on the set of stable configurations, which we show is abelian (Theorem 4), and we prove the existence of the unique steady state (Theorem 5). Many of our results can be proved using the operator algebra method e.g. as in [12], but here we introduce a different idea, which is useful for picturing the “shape” of the steady state. With its help we are able to prove that the average $\log(\mathrm{RHF})$ of SSP over the steady state is strictly bounded away from the worst-case by a constant independent of system size $L$ (Proposition 8). This is precisely what one wants to show with LLL, and hopefully our idea generalizes to this case, the only obstacle being its non-abelian toppling rule.

In Section 3, we provide some preliminary numerical studies on the statistical behavior of SSP, concerning its avalanche size distribution and its steady state. Though we still need a more extensive investigation, it seems that SSP has much in common with other one-dimensional stochastic models, such as the abelian Manna model and the Oslo rice-pile model. From the LLL perspective, it is especially interesting that both the Oslo and the Manna models exhibit the same boundary phenomenon as can be seen in Figure 1 (see [6] and [12]).

1.3. Acknowledgment

We thank Deepak Dhar, Phong Nguyen, and Su-Chan Park for helpful comments and suggestions.

2. Mathematical properties of SSP

2.1. Setup

In the literature on sandpiles, there seem to exist two sets of notations, one that is physics-oriented and one that is mathematics-oriented — heights versus configurations, sites versus vertices, relaxations versus stabilizations, and so on. This paper mostly employs the latter, and there are several words that we made up ourselves. In any case, we provide the definition of every term that we use below.

Like ASM, SSP is played on a finite undirected graph $G=(V,E)$, equipped with a designated vertex called sink, which we denote here by $s$. Throughout this paper, we assume that the graph obtained from $G$ by removing $s$ is connected. In addition, let $\wp$ be a probability mass function supported on a finite subset of $\mathbb{Z}_{>0}$, and let $T$ be an integer satisfying $T\geq\max\{i:i\in\mathrm{supp}\,\wp\}$, which we sometimes refer to as the threshold.

A pure configuration is a function $c:V\backslash\{s\}\rightarrow\mathbb{Z}$. A (mixed) configuration is a formal finite sum of pure configurations of form $\sum p_{i}[c_{i}]$, where $p_{i}$ are positive real numbers such that $\sum p_{i}=1$, and $c_{i}$ are pure configurations. The role the bracket notation $[.]$ here is to clarify that this is a formal sum. We define configurations this way, since toppling in SSP leads to multiple probabilistic results, as will be explained below.

A pure configuration $c$ is called stable if $c(v)<T$ for all $v\in V\backslash\{s\}$, and called nonnegative if $c(v)\geq 0$ for all $v\in V\backslash\{s\}$. A configuration is called stable or nonnegative if it is a formal sum of stable or nonnegative pure configurations, respectively.

Toppling a pure configuration $c$ at a non-sink vertex $v$ goes as follows: i) sample $\gamma$ from $I$ according to $\wp$, ii) subtract $\gamma\cdot\mathrm{deg}(v)$ from $c(v)$ iii) for each $w\sim v$ add $\gamma\cdot\mathrm{deg}(v,w)$ to $c(w)$, where $\mathrm{deg}(v,w)$ is the number of edges connecting $v$ and $w$. Denote the outcome by $T_{v}(c)$. If $c(v)\geq T$, we call this toppling legal. It is clear that, unless $\wp$ is supported on a singleton set, toppling a pure configuration results in a mixed configuration.

A choice algorithm is a map $\mathcal{A}:\{\mbox{unstable pure configurations}\}\rightarrow V\backslash\{s\}$, such that $c(\mathcal{A}(c))\geq T$. Fixing a choice algorithm $\mathcal{A}$, we can define what it means to topple a mixed configuration $C=\sum p_{i}[c_{i}]$. For a pure configuration $c$, define

[TABLE]

Then we define $T(C):=\sum p_{i}T(c_{i})$, the outcome of toppling $C$ once. The stabilization $C^{\circ}$ of $C$ is the outcome of toppling $C$ repeatedly until it becomes stable. Later we will show that $C^{\circ}$ is independent of $\mathcal{A}$.

Remark*.*

It is possible to extend the definition of a choice algorithm to a map into the set of probability mass functions on $V\backslash\{s\}$, so that $\mathcal{A}(c)(v)>0$ only if $c(v)\geq T$, and

[TABLE]

All the results in this paper carry over to this generalization.

As in the case of ASM, we can define an operation $\oplus$ on the set of nonnegative stable configurations. If $c$ and $d$ are pure stable configurations, we define $c\oplus d=(c+d)^{\circ}$, where $c+d$ is defined so that $(c+d)(v)=c(v)+d(v)$. In general, if $C=\sum p_{i}c_{i}$ and $D=q_{j}d_{j}$, then we define $C+D=\sum_{i,j}p_{i}q_{j}[c_{i}+d_{j}]$ and accordingly $C\oplus D=(C+D)^{\circ}=\sum_{i,j}p_{i}q_{j}(c_{i}\oplus d_{j})$.

Because we assumed that $T\geq\max\{i:i\in I\}$, for any nonnegative configurations $c$ and $d$, $c\oplus d=(c+d)^{\circ}$ is nonnegative as well. Later, we will show that $\oplus$ is abelian and associative, that the set of nonnegative stable configurations form a semigroup under this operation, whose minimal ideal consists of a single element which is the steady state.

2.2. Stabilization is independent of choice algorithms

We start by defining what we temporarily call an increment stack. An increment stack consists of $|V\backslash\{s\}|$ sequences, each of which is indexed by an element of $V\backslash\{s\}$, whose entries are chosen from the elements of $\mathrm{supp}\,\wp$; thus it has the form $S=(\gamma_{v,i})_{v\in V\backslash\{s\}\atop i\geq 1}$, with each $\gamma_{v,i}\in\mathrm{supp}\,\wp$. The set of all increment stacks is given a topology and a measure as follows: the cylinder sets, i.e. sets of form $Y(m,(v_{k},i_{k},p_{k})_{k=1}^{m})=\{(\gamma_{v,i})_{v\in V\backslash\{s\}\atop i\geq 1}:\gamma_{v_{k},i_{k}}=p_{k},k=1,\ldots,m\}$, form the base for the closed sets, and the measure $\mu$ is defined by $\mu(Y(m,(v_{k},i_{k},p_{k})_{k=1}^{m}))=\prod_{k=1}^{m}\wp(p_{k})$. This notion simulates the randomness in the amount being toppled during a stabilization process: SSP is equivalent to first sampling an increment stack $S$ according to $\mu$, and if one must topple at $v$ for the $i$-th time, remove $\gamma_{v,i}\cdot\mathrm{deg}(v)$ unit of sand from $v$ and distribute it equally among the neighbors of $v$ (respecting the multiplicities).

For a function $k:V\backslash\{s\}\rightarrow\mathbb{Z}_{\geq 0}$, a substack $S(k(v))$ of an increment stack $S=(\gamma_{v,i})_{v\in V\backslash\{s\}\atop i\geq 1}$ is a finite-length sequence of form $(\gamma_{v,i(v)})_{v\in V\backslash\{s\}\atop 1\leq i(v)\leq k(v)}$. We define the following partial order on the set of substacks: $S(k(v))\leq S(l(v))$ if $k(v)\leq l(v)$ for $v\in V\backslash\{s\}$. Given a pure configuration $c$, if for every $v\in V\backslash\{s\}$

[TABLE]

or equivalently, if toppling, legally or illegally, according to $S(k(v))$ stabilizes $c$ — then $S(k(v))$ is called a stabilizing substack of $c$.

Lemma 1**.**

For any given increment stack $S$ and pure configuration $c$, there exists a unique minimal stabilizing substack of $c$.

Proof.

Suppose $P:=S(k(v))$ and $Q:=S(l(v))$ are two distinct minimal stabilizing substacks. Define $r(v)=\min(k(v),l(v))$, and let $R:=S(r(v))$ be another substack; by assumption $R<P$ and $R<Q$. Then, for every non-sink $v$,

[TABLE]

Therefore $R$ is a stabilizing substack, contradicting the minimality of $P$ and $Q$. ∎

Given an increment stack $S$ and a pure configuration $c$, a choice algorithm $\mathcal{A}$ induces a stabilizing substack $S(k_{\mathcal{A}}(v))$, where $k_{\mathcal{A}}(v)$ is the number of times $\mathcal{A}$ would topple on $v$ until it stabilizes $c$.

Proposition 2**.**

Any choice algorithm induces the minimal stabilizing substack.

Proof.

Write $M:=S(m(v))$ for the minimal stabilizing substack, and $N:=S(k_{\mathcal{A}}(v))$ for the substack induced by a choice algorithm $\mathcal{A}$. Imagine $\mathcal{A}$ toppling an input configuration $c$, one vertex at a time, and consider a situation in which $\mathcal{A}$ has just toppled on $v_{1}\in V$ for $m(v_{1})$ times, $v_{1}$ being the first vertex at which this occurs. Observe that, at this point, the configuration is stabilized at $v_{1}$. Hence $\mathcal{A}$ cannot topple more on $v_{1}$, at least not until it topples on some neighborhood $w$ of $v_{1}$ more than $m(w)$ times.

Suppose $v_{2}$ is the second vertex at which $\mathcal{A}$ topples on $v_{2}$ for $m(v_{2})$ times. By the same argument as above, $\mathcal{A}$ cannot topple on $v_{2}$ more than $m(v_{2})$ times until it topples on some $w\sim v_{2}$ more than $m(w)$ times. Repeating, we see that on no vertices $v$ can $\mathcal{A}$ topple more than $m(v)$ times. This proves that $N\leq M$, but since $M$ is minimal we must have $M=N$. ∎

An immediate consequence of Proposition 2 is that, for any configuration $C$, pure or non-pure, its stabilization $C^{\circ}$ is independent of the choice algorithm $\mathcal{A}$. Indeed, $C^{\circ}$ is determined solely by the measure $\mu$ on the set of all increment stacks, and the minimal stabilizing substack of each increment stack. Therefore, the notion of a choice algorithm is somewhat superfluous from a theoretical viewpoint.

2.3. $\oplus$ is abelian and associative

It is clear by definition that $\oplus$ is abelian. We will now show that $\oplus$ is associative as well. Since $(F\oplus G)\oplus H=((F+G)^{\circ}+H)^{\circ}$ and $F\oplus(G\oplus H)=(F+(G+H)^{\circ})^{\circ}$, it suffices to show that

Proposition 3**.**

For nonnegative configurations $C$ and $D$,

[TABLE]

or equivalently,

[TABLE]

Proof.

Assume first that $C$ and $D$ are pure configurations. Fix an increment stack $S=(\gamma_{v,i})_{v\in V\backslash\{s\}\atop i\geq 1}$. Let $S(k(v))$ and $S(l(v))$ be the minimal stabilizing substacks for $C$ and $C+D$, respectively. Note $S(k(v))\leq S(l(v))$ — this is where we use the nonnegativity of $C$ and $D$.

Consider the left-hand side of (2), namely $(C^{\circ}+D)^{\circ}$. A priori, we need two increment stacks for each of the two stabilizations: $T_{1}$ for the inner term $C^{\circ}$ and $T_{2}$ for the outer one. Write $S^{\prime}=(\gamma_{v,i+k(v)})_{v\in V\backslash\{s\}\atop i\geq 1}$. It suffices to prove that

[TABLE]

or, in other words, that carrying out both stabilizations on the left-hand side of (2) using only one increment stack $S$ induces no distortion on the measure $\mu$. But this is clear from the definition of $\mu$.

It remains to consider the case in which $C$ and $D$ are not necessarily pure configurations. Write $C=\sum p_{i}[c_{i}]$ and $D=\sum q_{j}[d_{j}]$. Then

[TABLE]

as desired. ∎

Remark*.*

If the nonnegativity assumption is not present, Proposition 3 fails. In fact, it fails for ASM as well.

Thanks to our results so far, we have the following

Theorem 4**.**

The set of all nonnegative stable configurations of a SSP forms a commutative monoid under $\oplus$.

2.4. The existence of unique steady state

Denote the monoid in Theorem 4 by $M$. A steady state $S\in M$ is an element such that

[TABLE]

It is clear from this definition that a steady state, if exists, is unique.

Theorem 5**.**

Assume the greatest common divisor of the elements of $\mathrm{supp}\,\wp$ equals $1$. Then $M$ has a steady state.

There is no loss of generality incurred by the assumption in the theorem, since if we divide all elements of $\mathrm{supp}\,\wp$ by the g.c.d. then we obtain essentially the same model.

Our proof of Theorem 5 relies almost entirely on the following lemma.

Lemma 6**.**

Let $p:\{1,2,\ldots,n\}\rightarrow\mathbb{R}_{\geq 0}$ such that $\sum p(i)=1$, and such that the g.c.d. of the elements of $\mathrm{supp}\,p$ equals $1$. Consider the $n\times n$ matrix

[TABLE]

with the first column filled with $p(i)$, the superdiagonal with $1$, and [math] elsewhere. Then its (complex right) eigenvalue with the largest modulus is 1, which has multiplicity 1, with the eigenvector

[TABLE]

where $r(k)=\sum_{i\geq k}p(i)$.

Proof.

The characteristic equation of $A$ equals

[TABLE]

We first show that $1$ is an eigenvalue of multiplicity one. Certainly $\chi_{A}(1)=0$. Also since

[TABLE]

and $n>(n-1)p(1)+\ldots p(n-1)$, $\chi^{\prime}_{A}(1)\neq 0$. It is straightforward to check that (3) is a solution to $Ax=x$. This proves the latter two claims of the lemma.

If $|\lambda|>1$, then $|\lambda|^{n}>\sum p(i)|\lambda|^{n-i}$, so it cannot be a root of $\chi_{A}$. Therefore it suffices to show that $1$ is the only complex eigenvalue with modulus $1$. Suppose $|\lambda|=1$ and $\chi_{A}(\lambda)$ = 0. Then for all $i\in\mathrm{supp}\,p$, $\lambda^{n-i}=\lambda^{n}$ must hold. In other words, $\lambda$ is an $i$-th root of unity for all $i\in\mathrm{supp}\,p$. But since the g.c.d. of the elements of $\mathrm{supp}\,p$ equals 1 by assumption, $\lambda=1$ is forced. ∎

An immediate consequence of Lemma 6 is that, for any $x\in\mathbb{R}^{n}$, $A^{k}x$ approaches a constant multiple of (3) as $k\rightarrow\infty$. This is how we use Lemma 6 in the argument below.

Proof of Theorem 5.

Write $V\backslash\{s\}=\{v_{1},v_{2},\ldots,v_{L}\}$, and let $n$ be the smallest positive integer such that $\mathrm{supp}\,\wp\subseteq\{1,\ldots,n\}$. Define $\tau_{i}$ to be the toppling operator at $v_{i}$ with increment $1$ (hence $\tau_{i}$ coincides with the toppling operator of the standard ASM on $G$), so that $T_{v_{i}}=\sum_{j}\wp(j)\tau_{i}^{j},$ where $T_{v_{i}}$ is the toppling operator of SSP.

Define $r:\{0,\ldots,n-1\}\rightarrow\mathbb{R}$ by $r(k)=C\sum_{i>k}^{n}\wp(i)$, where the constant factor $C$ is set so that $\sum r(k)=1$, that is,

[TABLE]

For any pure configuration $c$, define

[TABLE]

One can think of $S(c)$ as a distribution supported on $\prod_{i=1}^{L}k_{i}$ points on the space of (pure) configurations, which turn out to form the shape of a parallelepiped; see Figure 2 for an illustration in case $L=2,\mathrm{supp}\,\wp=\{1,2,3,4,5\}$. The point on the upper-right corner corresponds to $[c]$, and other points are obtained by toppling $[c]$ in various directions; see again Figure 2 for an illustration.

Write $I(\cdot)$ for the indicator function of the statement inside the parenthesis, i.e. if it is satisfied, $I(\cdot)=1$, otherwise zero. From Lemma 6, it follows that

[TABLE]

This means that, if we “push” the parallelepiped-shaped distribution in the direction of $\tau_{i}$ by applying $T_{v_{i}}$ to the face of the parallelepiped corresponding to that direction, we obtain the same distribution on the parallepiped, except that the upper-right corner of the parallelepiped moves from $[c]$ to $[\tau_{i}c]$.

We choose $R$ to be the set of all recurrent configurations of ASM on $G$, whose threshold on each vertex is set uniformly to $T$, rather than the conventional $\deg v$. For convenience, let us temporarily write

[TABLE]

We will show that $S(R)^{\circ}$ is the steady state.

By Proposition 3, it suffices to show that $(S(R)+v_{i})^{\circ}=S(R)^{\circ}$ for all $i=1,\ldots,L$. By the theory of ASM, $R+v_{i}=\{c+v_{i}:c\in R\}$ contains exactly one representative of each of the equivalence classes of configurations i.e. for each $c^{\prime}\in R+v_{i}$, there exists exactly one $c\in R$ such that $c^{\prime}$ stabilizes to $c$ in ASM. By (4), $S(c^{\prime})^{\circ}=S(c)^{\circ}$, and therefore, by Proposition 3 again,

[TABLE]

as desired. ∎

The point of the following proposition is that, in a SSP, if we start from any configuration, and repeat sufficiently many times the Markov process of adding a sand grain at random — e.g. add $\frac{1}{L}\sum v_{i}$ — and re-stabilizing, then we obtain a configuration that is arbitrarily close to the steady state. In the context of LLL, in place of the Markov chain one takes a large input, which is why the proposition below is stated the way it is.

Proposition 7**.**

We continue with the notations of Theorem 5. In addition, below we interpret a pure configuration as an $L$-dimensional vector with integer entries, and a mixed configuration as a finitely supported function on the set of pure configurations. Then we impose the typical Euclidean metric on the space of pure or mixed configurations accordingly.

For any $\varepsilon>0$, there exists $D>0$ such that the stabilization $c^{\circ}$ of any non-negative pure configuration $c$ whose distance from origin is greater than $D$ is within an $\varepsilon$ distance of $S(c^{\prime})^{\circ}$ with respect to the uniform norm, where $c^{\prime}$ is the pure configuration that is recurrent and is equivalent to $c$ under the ASM on $G$.

Proof.

If $c$ is far enough from the origin, then it must be toppled at each vertex arbitrarily many times in order for it to become stabilized. So without loss of generality, assume that $c$ can be legally toppled at every vertex. Write $V\backslash\{s\}=\{v_{1},v_{2},\ldots,v_{L}\}$ as earlier, and consider the configuration

[TABLE]

(5), like $S(c)$, may be thought of as a distribution on a parallelepiped illustrated in Figure 2. By Lemma 6, if we “push” the parallelepiped — in the manner similar to (4) — from every direction sufficiently many times, the weight on each and every point will be arbitrarily close to $\prod_{j=1}^{L}r(k_{j})$. Hence, if $c$ is sufficiently far enough from the origin, the stabilization of (5) will be arbitrarily close to $S(c^{\prime})$. This completes the proof. ∎

The steady state alone forms the minimal ideal of the commutative monoid $M$ of the nonnegative stable configurations of SSP. One may think this is disappointing, since its counterpart in ASM has a very rich theory. However, the theory of the abelian sandpile group can be extended to SSP in a straightforward way, if phrased in terms of the cosets of the “translations” $<\tau_{i}:i=1,\ldots,L>$, since SSP respects this structure. Indeed, in our proof of Theorem 5, it can be seen that the underlying ASM theory plays a role. Our SSP theory is building up on top of the ASM theory, rather than replacing it.

2.5. The shape of the steady state on $A_{L}$

Ideally, we would like to prove the various quantitative properties of the steady state that we observe in Figure 1: the middle values are almost identical to one another, but near the boundaries there is a sharp drop in pile height, et cetera. This can be interpreted as the following combinatorial problem. Consider the parallelepiped in Figure 2, and suppose $[c]$ is a stable configuration; it may help to think of $[c]$ as the maximal stable configuration, i.e. has height $T-1$ on all sites. Then some configurations on the parallelepiped are stable e.g. $[\tau_{1}^{3}\tau_{2}^{2}c]$, but others may not be. Because $[c]$ is stable, we cannot legally push the entire face of the parallelepiped as we have done in the proofs of Theorem 5 and Proposition 7. Instead, we must “fold” the parts of the parallelepiped that are outside the set of stable configurations. Understanding the steady state then amounts to understanding the resulting distribution on the “folded” parallelepiped — they are identical. In some sense, this parallelepiped-folding replaces the role of the burning algorithm (see [3], also [2]) in the study of the steady state.

The difficulty lies in describing this folded distribution in a mathematically concise manner. The best we can prove at the moment is that there exists a gap, whose size is independent of the system size $L$, between the average-case and the worst-case RHF (1) of SSP on $A_{L}$. This is in analogy with the folklore observation that there exists a gap between the average-case and the worst-case RHF of the LLL algorithm.

Proposition 8**.**

Let $C$ be as in the proof of Theorem 5. Then for any $\varepsilon>0$ and $L$ sufficiently large, the steady state of SSP on $A_{L}$ has the average RHF bounded from above by $T/2-C^{-1}/e^{2}+\varepsilon$.

Remark*.*

Note that the maximum $\log(\mathrm{RHF})$ equals $(T-1)L/2(L+1)\approx T/2$. For the SSP in the introduction i.e. $\wp=$ (uniform distribution on $\{1,\ldots,2I-1\}$), $C^{-1}\approx I$ so the proposition yields the bound $T/2-I/e^{2}$, whereas experimentally we have $\approx T/2-I/4$.

Proof of Proposition 8.

As in the remark above, the greatest possible $\log(\mathrm{RHF})$ equals $M:=(T-1)L/2(L+1)$. For $a>0$, we will estimate the number of pure stable configurations whose RHF is greater than $M-a$. If $c$ satisfies such a condition, then writing $s_{i}=T-c(i)$, from (1) it follows that

[TABLE]

Moreover, this is an if-and-only-if condition. Hence it comes down to measuring the volume of the set of vectors $(s_{1},\ldots,s_{L})$, $s_{i}>0$, such that (6) holds. This set forms a simplex, with one vertex at origin, and $L$ other vertices given by

[TABLE]

for each $i=1,\ldots,L$, whose only non-zero entry is the $i$-th entry. The volume of this simplex equals

[TABLE]

which, by Stirling’s formula, is strictly bounded by $a^{L}e^{2L}$ for all sufficiently large $L$. Furthermore, the maximum possible probability mass of the steady state is given by $C^{L}$. Hence, for $L$ large, the portion of the steady state lying on the set $(\mathrm{RHF})>M-a$ is arbitrary small if

[TABLE]

This completes the proof. ∎

3. Behavior of SSP on one-dimensional lattices

In this section, we present a pilot numerical study on the behavior of SSP pertaining to quantitative properties of its steady state, whose existence is established in the previous section. We restrict our attention to SSP on one-dimensional grids $A_{L}$ for $L=400,4000,40000$, with $T=400$ and $\wp$ the uniform distribution on $\{1,\ldots,100\}$, adding $50$ grains to a random site for each step in the associated Markov chain. We hope to be able to carry out more extensive experiments at a later time.

3.1. Avalanche size distribution

Figure 3 presents the log-log scale graph of the avalanche size — i.e. the number of sites on which at least one toppling occurred during a stabilization — distribution of SSP. For each system size $L=400,4000,40000$, we made one million trials. The $x$-axis represents $\log_{2}$ of avalanche sizes, and the $y$-axis represents $\log_{2}$ of the number of occurences.

Figure 3 clearly suggests a power law, plus a delta distribution near $\log_{2}L$. The data points for $L=40000$, excluding those at the tail, form a line of slope $\approx-0.98$. Of course, more experiments are needed to precisely determine the exponent of SSP.

It may be amusing to recall the behavior of one-dimensional ASM, which is well-known to behave in the totally opposite way, i.e. (avalanche size) $\propto$ (frequency). On the other hand, several stochastic models in dimension 1 are known to demonstrate a power law ([6], [7], [12]).

3.2. Average heights of the piles in the steady state

Figure 4 presents the (empirical) average steady state density of SSP for system sizes 400, 4000, 40000, respectively. The values in the middle seem to stabilize at $\approx 176.6$, and in the extremes at $\approx 152.5$.

The obvious curiosity is the fall-off shape at the boundaries. Figure 5 puts together the three graphs in Figure 4 for a comparison. Here we can see that the boundary areas overlap, which indicate that the effect of a boundary is independent of system size. Also observe that, in Figure 5, the first half of the average output shape of $L=400$ overlaps with the other shapes; the latter half, of course, must be the mirror image of the first half. Indeed, the first half of the average output shape of $L=4000$ is also found to overlap with that of $L=40000$.

At this point, we are able to only speculate as to why the average height diminishes near the boundaries, as follows. When a toppling occurs at one of the middle vertices, it is expected that some of the toppled sand comes back when the neighboring vertex topples. Some grains of sand may travel several steps away from where it started and then come back. However, when a toppling occurs at the left extreme, for example, the sand that moved to the left is lost forever.

Other stochastic sandpile models in dimension one, such as the abelian Manna model ([12]) and the Oslo rice-pile model ([6]), have also been observed to exhibit similar-looking shapes. [6] attributes it to the finite-size scaling theory — see equation (7) in [6]. However, in the SSP case, the same explanation seems to account for only the first 10 points or so away from the boundary.

It is also natural to ask whether the general features of the average steady state density are preserved if $\wp$ is replaced with something other than a uniform distribution. To this end, we ran several brief experiments with the Poisson distribution, and also the family of fat-tailed distributions of index $\alpha$ i.e. $\mathrm{Pr}(X>x)\sim x^{-\alpha}$, which is often used to produce unusual outcomes.222These distributions have infinite support, and thus they are not strictly covered by the theory developed here. However, the parallelepiped-pushing idea still goes through. In all cases, we still observe both the flatness of the middle values and the diminishing amount of sand near the boundaries, suggesting that these are general properties of one-dimensional stochastic sandpiles. For instance, see Figure 6, the outcome for the case of fat-tailed distributions of indices $3$ and $4$, with mean $\approx 15,13$ and standard deviation $\approx 17,14$ respectively. It may be interesting to study how the average output shape of SSP and the representative values of $\wp$ are related; for example, it seems reasonable to conjecture that the mean of $\wp$ correlates with the pile height in the middle, and the standard deviation with the boundary effect.