% ------- LaTeX file --------
BODY
%-----------------------------
% LaTeX file
%
\documentstyle[11pt]{article}
\topmargin -0.5in
\textwidth 158mm
\oddsidemargin 8mm
\textheight 9.1in
\pretolerance=2000
\parskip 5pt plus 1pt minus 0.5pt
\title{Exact solutions for a mean-field Abelian sandpile\thanks{
Supported in part by NSF Grant DMR89-18903}}
\author{Steven A. Janowsky\thanks{Supported in part by NSF Mathematical
Sciences Postdoctoral Research Fellowship DMS
90-07206}\hspace{.25em}\thanks{Address after August 1993:
Department of Mathematics, University of Texas,
Austin, TX 78712}\\
Department of Mathematics\\
Rutgers University\\New Brunswick, NJ 08903
\and
Claude A. Laberge\\
Department of Physics\\
Rutgers University\\Piscataway, NJ 08855-0849}
\date{June 1993}
\begin{document}
\maketitle
\renewcommand{\baselinestretch}{1.57}\large\normalsize
\begin{abstract}
We introduce a model for a sandpile, with $N$ sites, critical height $N$
and each site connected to every other site. It is thus a mean-field
model in the spin-glass sense. We find an exact solution for the
steady state probability distribution of avalanche sizes, and discuss
its asymptotics for large $N$.
\end{abstract}
%\subsection*{Introduction}
The publication of a series of articles by P. Bak and collaborators
\cite{Bak1,Bak2} generated a growing interest in the study of certain
class of cellular automata models now commonly known as ``sandpiles''
because of a crude analogy between the dynamical rules and the way sand
topples when building a sandpile. The basic motive was not to model the
way real sandpiles behave, although it did broaden the audience for
granular flow (see e.g.~\cite{Nag,Mehta,Puhl}), but to study a
phenomenon they coined ``self-organized criticality'' in which a
system reaches a critical
state, i.e.\ a state with no intrinsic length or time scales, without
the tuning of an external parameter. This generated a big influx of
theoretical, computational and even experimental studies of
self-organized criticality and of sandpiles in general. Unfortunately,
even though the number of models (i.e.\ sets of dynamical rules)
exhibiting self-organized criticality or falling into the more general
category of sandpiles has increased dramatically, very few exact results
are known at the present (see e.g.\ \cite{DM1,DM2,DM3,DR,Lee}).
In this letter we present a model in which
many quantities can be calculated exactly in a rigorous way.
%\subsection*{The model}
A sandpile model is basically a set of dynamical rules describing the
way that grains of sand are added to a system, the conditions under
which those grains can be redistributed inside the system, and the way
they are removed from the system. Here we consider a system of $N$ sites
and define $h(i)$ as the (integer) height of the column of sand at site
$i$, $i\in \{1\ldots N\}$. We drop a grain of sand on a site $i$ chosen
at random, thereby increasing its height by one: $h(i) \rightarrow
h(i)+1$. If this new height exceeds the maximum stable value $h_c$ then
that column topples and gives 1 grain of sand to each of the $N-1$ other
sites while one grains drops out of the system. (We take $h_c \ge N$ so
that $h(i)\ge 0$; in fact we are primarily interested in $h_c=N$.) We
then examine the system to see if any site has a column exceeding $h_c$
in which case we topple that column also. We keep toppling until all
the sites are stable (this characterizes an avalanche). We then repeat
the procedure of adding a grain at a randomly chosen site.
This model falls into the category of abelian sandpiles since we always
obtain the same stable configuration from an unstable one irrespective
of the order in which we executed the topplings. Dhar~\cite{Dhar} used
this fact to obtain several properties of those models which we use
extensively in our analysis. Furthermore, since all sites are connected
to all other sites, we can say that we have a mean field theory of
abelian sandpiles. It is this combination that permits the model to be
solved. It should be pointed out that our model is quite different from
previous attempts to study mean-field sandpiles \cite{Tang,Gaveau}.
%\subsection*{Results}
It is shown in~\cite{Dhar} that not all configurations are allowed in the
asymptotic regime, but that the allowed (recurrent) configurations all
have equal weight. In general, the number of recurrent configuration is
given by the determinant of the toppling matrix $\Delta$ which is an $N
\times N$ matrix where $\Delta_{ii}=h_c$ and $\Delta_{ij} = -1$ for each
site $j$ connected to $i$, so that row $i$ of $\Delta$ represents the
amount of sand lost by every site when site $i$ topples. Thus in the
model considered here this matrix takes the form
\begin{equation}
\Delta_{ij}= -1 + (h_c + 1)\delta_{ij},
\end{equation}
where $\delta_{ij}$ is the Kronecker delta.
Because of the highly symmetric nature of $\Delta$ in our model
it is straightforward to calculate the determinant and determine that
the number of recurrent configurations is
\begin{equation}
Z(N, h_c) = \det\Delta = (h_c + 1 - N)(h_c + 1)^{N-1}.
\end{equation}
For the case $h_c=N$ we have $Z(N,N)=(N+1)^{N-1}$. Incidentally this is
the number of spanning trees on a fully connected graph with $N+1$
sites. Such a result was expected from the general observation that
$\det\Delta$ is precisely the cofactor of the matrix tree (see e.g.\
\cite{Graph}) for graphs defined on a superset of the sandpile lattice
(there is an additional site, the ground, connected to every site). In
fact there is a one-to-one correspondence between the recurrent
configurations and the spanning trees, in the sense that a configuration
is recurrent if and only if it passes the ``burning algorithm'' test
described in \cite{DM1,DM3}. Actually this holds only for symmetric abelian
sandpiles, i.e.\ models where $\Delta$ is symmetric; for the asymmetric
case see \cite{Speer}.
Clearly the recurrent configurations fall into equivalent classes, where
two configurations are equivalent if one is a permutation of the other.
Because of the symmetry of the model we need only examine one member of
each equivalence class; we find it convenient to restrict ourselves to
the case $h(i)\le h(i+1)$, $i=1\ldots N-1$, where the amount of sand
increases as the site label increases.
Our analysis will focus on a particular subset of the recurrent
configurations, namely the minimal configurations. A given
configuration is minimal if there is no recurrent configuration that can
be obtained by removing sand from the given configuration. Some brief
computations show (one can use the so-called ``burning algorithm''
of~\cite{DM1}) that the configuration $h(i) = i$ is a minimal
configuration. All minimal configurations fall into this equivalence
class and thus there are exactly $N!$ minimal configurations, the
permutations of $h(i) = i$. If we now consider equivalence classes for
general configurations, we see that the restriction that a recurrent
configuration be greater than or equal to the minimal configuration
yields
\begin{equation}\label{stability}
h(i) \ge i, \qquad
h(i) \ge h(i-1)
\end{equation}
(set $h(0) = 0$).
Now consider how avalanches of various sizes come about in our model; we
limit the discussion to the case $h_c=N$. Size here means either the
number of grains that fall out of the system or the number of sites that
topple; in the model considered here there is no ambiguity because they
are equal. An avalanche begins when a grain of sand is dropped on a
site of height $N$. It will continue if the second-highest site is also
at height $N$. It will be at least of size three if the third-highest
site is at height $N-1$, and it will be at least of size four if the
fourth-highest site is at height $N-2$, etc. Therefore the size will be
(using the equivalence class notation)
\begin{equation}\label{N-aval}
N_{\rm aval}= N- \min \{j: h(i) > i {\rm\ for\ all\ } j\le i0$ fixed and $N\rightarrow\infty$, (\ref{exactav}) yields
\begin{equation}
P({\rm avalanche\ of\ size\ }k) \sim \frac1N\frac{k^{k-2}e^{-k}}{(k-1)!},
\end{equation}
and if k is also large ($1\ll k\ll N$) then this reduces to
\begin{equation}
P({\rm avalanche\ of\ size\ }k) \sim \frac{k^{-3/2}}{\sqrt{2\pi}N}.
\end{equation}
Thus we reproduce the exponent $-3/2$ previously derived on
trees~\cite{DM2} and via numerical and nonrigorous
arguments~\cite{Obukov,KNWZ}.
%\subsubsection*{Height distribution}
Another quantity of interest is the single-site probability distribution
for the heights, {\em i.e.}\ the probability that a given site contains
$H$ grains of sand in the stationary state. First look at
configurations where exactly $K$ sites have height $H$. (In the
equivalence class notation these sites must be consecutive.) Remove
these sites and consider the subsystem consisting of the remaining $N-K$
sites. The recurrent configurations of this subsystem are in one-to-one
correspondence with those of a system of size $N-K$ where the recurrent
configurations are determined by, in addition to the usual occupancy
restrictions (see (\ref{stability})), the conditions $j \le h(j) \le
N-1$ if $j \le H-K$ and $j+K-1 \le h(j)\le N-1$ if $j > H-K$.
The direct evaluation of the number of allowed configurations satisfying
the above criteria is possible but somewhat complicated; appropriately
summing over $K$ would yield the desired probabilities. However, it
turns out to be simpler instead to compare the number of allowed
configurations with $K$ sites having $H$ grains with the number of
allowed configurations with $K$ sites having $H-1$ grains, {\em e.g.}
by subtracting the latter from the prior. Which configurations are left?
In both cases $j \le h(j) \le N-1$ for $j < H-K$ and $j+K-1 \le h(j)\le
N-1$ if $j > H-K$. The only difference occurs at site $H-K$ which can
be occupied by between $H-1$ and $N-1$ grains in the first case, and
between $H-K$ and $N-1$ grains in the second. If $K=1$ these conditions
are identical and no allowed configurations are left. For $K \geq 2$,
however, site $H-K$ is now restricted so that $h(H-K) \le H-2$. Since
we need $h(j+1) \geq h(j)$ (for the configuration to be a recurrent
one), this implies that the first $H-K$ sites can have a maximum of
$H-2$ grains each and thus contribute a factor $Z(H-K,H-2)$ to the
partition function, excluding (global) symmetry factors. On the other
hand, (\ref{stability}) imposes no additional conditions on the last
$N-H$ sites, and it is easily seen that they contribute a factor
$Z(N-H,N-H)$. Defining $P(H)\;$= the probability that a given site will
have $H$ grains, putting in the appropriate symmetry factors and summing
over the index $K$ we get
\begin{equation}
P(H) - P(H-1) = Z(N,N)^{-1}Z(N-H, N-H) \sum_{K=2}^H {N \choose K, H-K, N-H}
\frac{K}{N} Z(H-K,H-2),
\end{equation}
which is exactly (\ref{avalanchesum})! Since $P(0)=0$ we can collapse
the telescoping sum and obtain
\begin{equation}\label{height-dist}
P(\mbox{site has $H$ grains}) = \sum_{k=1}^H P(\mbox{avalanche of size k}).
\end{equation}
Because of the symmetry of the avalanche distribution, we have here the
symmetry
\begin{equation}
P(\mbox{site has $H$ grains}) + P(\mbox{site has $N-H$ grains}) =
\frac{2}{N+1}
\end{equation}
for $0\le H \le N$.
The asymptotics of (\ref{height-dist}) are relatively easy to compute using our
previous results. For $N\gg1$ we have
\begin{equation}
P(\mbox{site has $H$ grains}) \sim \frac1N\sum_{k=1}^H
\frac{k^{k-2}e^{-k}}{(k-1)!}
\end{equation}
which can be used directly to compute the behavior for small $H$.
For $N\gg H\gg 1$ we have
\begin{equation}
P(\mbox{site has $H$ grains}) \sim \frac1{N+1}\left(1-\sqrt{\frac2{\pi
H}}\right).
\end{equation}
%\subsubsection*{Mass distribution in the stationary state}
We can also examine the total amount of sand in the system. The
distribution of masses (total number of grains of sand) of recurrent
configurations is most easily computed for mass $M$ near the maximum
value $N^2$. For example, the number of recurrent configurations with
mass $N^2-2$ is simply $N(N+1)/2$. This is easily generalized to
arbitrary mass (using inclusion-exclusion arguments) and we find that
the number of recurrent configurations with total mass $N^2-K$,
$Z_{N^2-K}(N,N)$, is
\begin{eqnarray}
\lefteqn{Z_{N^2-K}(N,N) = {N-1+K \choose N-1} +}\nonumber\\
&&\sum_{l=1}^{l(N,K)} \!(-1)^l {N \choose l}
\sum_{j=lN+l(l-1)/2}^K {N-l-1+K-j \choose N-l-1}
{j - lN - l(l+1)/2 -1 \choose l-1},
\label{massdist}\end{eqnarray}
where
\begin{equation}
l(N,K) = \left\lceil N - (N(N-5) + 2K + 17/4)^{1/2} - 3/2
\right\rceil
\end{equation}
simply indicates ``how far'' one must go in the inclusion-exclusion;
$\lceil A \rceil$ is the least integer $\ge A$. Then the probability of
a configuration having mass $M$ is $Z_M(N,N) / (N+1)^{N-1}$.
Unfortunately we are only able to find the asymptotic behavior of
(\ref{massdist}) for the tails of the distribution, i.e.\ $M$ near
$N(N+1)/2$ or $M$ near $N^2$, which is not where most of the (recurrent)
configurations reside. Expressions similar to that which were used
to find the avalanche distribution, e.g.\ (\ref{avalanchesum}), are
available:
\begin{equation}
Z_M(N,h_c) = \sum_{j=1}^N Z_{M-jh_c} (N-j, h_c-1) {N \choose j},
\end{equation}
and one could consider a generating function approach to eliminate the
restrictions on $j$ in the sum. Unfortunately the resulting expression
will only converge in the regions equivalent to very large or very small
mass, and thus provide no new information.
%\subsubsection*{Equivalent particle model; time distribution of avalanches}
The reformulation of the $N$-dimensional model in terms of a one
dimensional sandpile-like problem (with the required symmetry factor
attached to each configuration) greatly facilitated the computations of
``static'' quantities like the avalanche size, single-site height, and
total mass distributions. In order to study the ``dynamical'' quantities
associated to the evolution of an avalanche, another reformulation in
term of a one dimensional {\em particle} model proves to be useful; we
will use it to calculate the distribution of avalanche durations. This
reformulation also provides a much more efficient method for numerical
study than a straightforward implementation of the sandpile process.
The duration of an avalanche is the number of sweeps that are required
in order to reach a stable configuration once a grain of sand is dropped
on the system. A sweep consists of two steps: we first go through all
the sites and mark those with a number of grains exceeding the critical
value and then topple all those sites simultaneously. The process is
repeated until no site can topple. Note that this definition of
duration is not the only one; others are possible because of the abelian
nature of the model.
The particle model consists of $N$ particles moving on a one dimensional
lattice with sites labeled from left to right by $k \in [1,N+1]$ (where
N is the number of site in the original model). The correspondence
between sandpile configurations and particle configurations is as
follows: given a sandpile configuration $C_0$ place a particle at
position k for every site in $C_0$ containing k grains. If $n_k$ is the
resulting number of particles at site $k$, the stability condition
(\ref{stability}) translates to $n_k \leq k$ and $n_N >1$; since the
critical height is $N$ we also have $n_{N+1}=0$.
The boundary conditions are ``almost'' periodic as will be explained
below; we also introduce a marker (barrier) between sites $N$ and $N+1$
that will play an important role in the dynamics.
The dynamics of the particle model is as follow: pick a particle at
random and move it one site to the right (this is equivalent of dropping
a grain of a sand on $C_0$). If the barrier was not crossed then we have
a new stable configuration and can repeat the step. Thus, neglecting
the boundary, our model is a zero-range process; however, the boundary
plays a very important role.
If the chosen particle jumps through the barrier it means that one site
in $C_0$ is unstable and will topple. In the particle model such a
toppling will mean that all particles simultaneously jump one site to
the right except for the one at site $N+1$ which should move back to
site 1. An equivalent formulation which we find more convenient is to
move the barrier to the left by one unit and relabel the sites
accordingly (i.e.\ the barrier defines the position of site N and N+1 on
the periodic lattice).
As we moved the barrier, a number of particles (say $p_2$) might have gone
through it. In the sandpile model this means that the initial toppling
caused $p_2$ more sites to become unstable. To duplicate their toppling
the barrier should now be moved $p_2$ more units to the left. The process
continues until no particle crosses the barrier as it jumps. The total
number of barrier jumps required is the duration of the avalanche.
The probability that an avalanche has duration $T$ is equivalent to the
fraction of all particle configurations that will make the barrier jump
$T$ times if a particle at site $N$ is picked. Of course $T=0$
corresponds to the no avalanche case so that the probability that the
avalanche duration is zero is $(N-1)/(N+1)$ (see (\ref{noavalancheresult})).
If there are $T$ barrier jumps, let $p_1, p_2, \ldots p_T$ be the number of
particles crossing the boundary at each jump and let $P = p_1 + p_2 +
\cdots + p_T$; note that $p_i>0$ for $i 1$ for $T>1$. Provided $T>0$ we have
\begin{equation}
P(\mbox{duration $T$})=Z(N,N)^{-1}\sum_{p_1,\ldots,p_T\rm~restricted}
\frac{p_1}{N}W_N(1; p_1, p_2, \ldots, p_T)
\end{equation}
where $W_N(1; p_1, p_2, \ldots, p_T)$ is the number of configurations
that produce the specified barrier jumps. We can compute this in an
iterative fashion, in terms of $W_{M}(s; p_i, p_{i+1}, \ldots, p_T)$, the
number of (restricted) configurations on the subsystem consisting of
sites $1\ldots M$; $s$ is size of the current jump, and $p_i, p_{i+1},
\ldots, p_T$ are the remaining jump sizes. We have
\begin{eqnarray}
W_N(1; p_1, p_2, \ldots, p_T) &=& {N \choose p_1} w(p_1,1)
W_{N-1}(p_1-1; p_2, p_3, \ldots, p_T),\\
W_{N-1}(p_1-1; p_2, p_3, \ldots, p_T) &=& {N-p_1\choose p_2} w(p_2,p_1-1)
W_{N-p_1}(p_2; p_3, p_4, \ldots, p_T),\\
&&\mbox{and}\nonumber\\
W_{M}(p_{i-1}; p_i, p_{i+1}, \ldots, p_T) &=& {M \choose p_i} w(p_i,p_{i-1})
W_{M-p_{i-1}}(p_i; p_{i+1}, p_{i+2}, \ldots, p_T)
\end{eqnarray}
for $2*0$):
\begin{equation}
P(\mbox{duration $T$})=\frac{1}{(N+1)^{N-1}}\!\! \sum_{S=T}^N \frac{1}{N}
{N \choose S}(N-S+1)^{N-S-1}
\!\!\!\!\!\!\!\sum_{s_1 +s_2 +\cdots+ s_T=S} \frac{S!}{s_1!..s_T!}
s_1^{s_2} \ldots s_{T-1}^{s_T}
\end{equation}
where we have the restriction $s_1=1$, $s_i > 0$.
We identify $s_i$ with $p_{i-1}$ except that because of the need to
initiate the avalanche $s_1=1$ and $s_2 = p_1-1$. Of course $s_i$ is
nothing more than the size of the $i^{\rm th}$ jump.
%A plot of $P(T)$ reveals that it is not a monotonically decreasing
%function of $T$ and that the asymptotic ($T \sim N$) region does not
%exhibit any type of power law.
\subsubsection*{Acknowledgments:}
We would like to thank J. L. Lebowitz, S. Majumdar and E. R. Speer for
their help with this work.
\renewcommand{\baselinestretch}{1.1}\large\normalsize
\begin{thebibliography}{99}
\bibitem{Bak1} Bak P, Tang C and Wiesenfeld K 1987 Self-organized
criticality:An explanation of $1/f$ noise
{\em Phys. Rev. Lett.} {\bf 59} 381--384
\bibitem{Bak2} Bak P, Tang C and Wiesenfeld K 1988 Self-organized
criticality {\em Phys. Rev.} A {\bf 38} 364--374
\bibitem{Nag} Nagel S R 1992 Instabilities in a sandpile
{\em Rev. Mod. Phys.} {\bf 64} 321--325
\bibitem{Mehta} Mehta A 1992 Real sandpiles: dilatancy, hysteresis and
cooperative dynamics {\em Physica} A {\bf 186} 121--153
\bibitem{Puhl} Puhl H 1992 On the modelling of real sand piles
{\em Physica} A {\bf 182} 295--319
\bibitem{DM1} Majumdar S N and Dhar D 1991 Height correlations in the
Abelian sandpile model {\em J. Phys.} A {\bf 24} L357--L362
\bibitem{DM2} Dhar D and Majumdar S N 1990 Abelian sandpile model on
the Bethe lattice {\em J. Phys.} A {\bf 23} 4333--4350
\bibitem{DM3} Majumdar S N and Dhar D 1992 Equivalence of the Abelian
Sandpile Model and the $q \rightarrow 0$ Limit of the Potts Model
{\em Physica} A {\bf 185} 129--145
\bibitem{DR} Dhar D and Ramaswamy R 1989 Exactly Solved Model of
Self-Organized Critical Phenomena {\em Phys. Rev. Lett.} 63
1659--1662
\bibitem{Lee} Lee S-C, Liang N Y and Tzeng W-J 1991 Exact Solution of a
Deterministic Sandpile Model in One Dimension {\em Phys. Rev.
Lett.} {\bf 67} 1479--1481
\bibitem{Dhar} Dhar D 1990 Self-Organized Critical State of Sandpile
Automaton Models {\em Phys. Rev. Lett.} {\bf 64} 1613--1616
\bibitem{Tang} Tang C and Bak P 1988 Mean Field Theory of Self-Organized
Critical Phenomena {\em J. Stat. Phys.} {\bf 51} 797--802
\bibitem{Gaveau} Gaveau B and Schulman L S 1991 Mean-field self-organized
criticality {\em J. Phys.} A {\bf 24} L475--L480
\bibitem{Graph} Harary F 1990 {\em Graph Theory} (Reading, MA: Addison-Wesley)
\bibitem{Speer} Speer E R 1993 Asymmetric Abelian Sandpile Models
{\em J. Stat. Phys.} {\bf 71} 61--74
\bibitem{Obukov} Obukov S P 1988 in {\em Random Fluctuations and Pattern
Growth: Experiments and Models}, ed. H E Stanley and N
Ostrowsky, NATO ASI Series C {\bf 157} (Dordrecht:Kluwer)
\bibitem{KNWZ} Kadanoff L P, Nagel S R, Wu L and Zhou S-M 1989
Scaling and universality in avalanches
{\em Phys. Rev.} A {\bf 39} 6524--6537
\end{thebibliography}
\end{document}
ENDBODY
*