Dropping grains one at a time drives a sandpile to a critical state, with avalanches of every size and no parameter to tune. It is the simplest model of self-organized criticality.
Drop grains of sand on a single spot. At first the pile only grows; once it gets steep, grains start to slide, and a single grain can trigger anything from a trickle to an avalanche that rearranges the pile.
The regrowing forest of the previous essay produced fires of every size without anyone setting its density. Bak, Tang and Wiesenfeld (1987) turned the sandpile into a grid model that evolves to a critical state on its own, and called this self-organized criticality. Below, the order of toppling turns out not to matter, avalanche sizes are measured on four grids and compared with a random variant of the rule, a test sorts the piles the drive can reach from the ones it cannot, a single spot of sand grows a fractal, and the model's original motivation, that it explains 1/f noise, gets checked.
Each cell holds 0 to 3 grains. A cell that reaches 4 topples: it loses 4 grains and sends one to each of its four neighbours, which may topple in turn. Grains pushed off the edge are lost.
Click to drop a grain, or auto-drop. Early drops do nothing; cascades begin once most cells hold 2 or 3 grains. On this 100×100 grid it takes roughly 21,000 grains before the mean height stops climbing and settles near 2.11. On an infinite grid it would be exactly 17/8 = 2.125; grains lost over the edges keep a finite pile a little lower. Skip to critical state drops 30,000 grains at once.
Figure 1. The sandpile on a 100×100 grid, shaded by grain count from 0 (faintest) to 3 (strongest). Click the grid to drop a grain where you click; auto-drop adds grains at random sites (seeded, so a reload repeats the run).
When several cells are over the threshold at once, something has to decide which topples first. Figure 1 keeps a queue; a parallel update would topple them all at once. For this rule the choice makes no difference: the final pile, and the number of times each cell topples, are the same in every order (Dhar 1990). Figure 2 relaxes one overloaded pile twice.
Figure 2. One 9×9 pile with three overloaded cells, relaxed twice. Left: always topple the first unstable cell in reading order. Right: topple an unstable cell chosen at random. Numbers are grain counts, amber at 4 or more; cyan shading counts how often each cell has toppled (darkest at 6). Partway through, the two copies disagree at up to 48 cells (default seed). When both stop, they agree everywhere, in grains and in toppling counts. New pile draws another seed.
That commutativity is why the model is called Abelian. It also means a stable pile can be added to another cell by cell and relaxed with a well-defined answer, which a later section uses.
In the stationary state just over half the grains topple nothing, and a few topple tens of thousands of cells. The probability of an avalanche of size $s$ (topplings per grain) falls on a straight line on log-log axes, $P(s) \sim s^{-\tau}$, so there is no typical size. On a finite grid the line ends where avalanches run into the edge, and that cutoff grows as a power of the grid width $L$. If both hold,
and plotting $s^{\tau} P(s)$ against $s / L^{D}$ puts every grid size on one curve. Figure 3 runs four grids from 16 to 128 cells wide, for this rule and for the stochastic rule of Manna (1991), in which a cell holding 2 or more grains empties and throws each grain to a neighbour chosen at random.
Figure 3. Left: avalanche-size distributions at $L$ = 16, 32, 64 and 128 (faint to strong), log-binned, as $s^{\tau}P(s)$ against $s/L^{D}$. With $D = 0$ the bumps where avalanches meet the edge sit at different sizes; raising $D$ slides them onto each other. Right: if $\langle s^q \rangle \sim L^{\sigma(q)}$ and the scaling form holds, $\sigma(q)$ is a straight line of slope $D$ once $q$ is past $\tau - 1$, so its local slope, plotted here, levels off at $D$. Each $\sigma(q)$ is a least-squares fit over the four grids. Moment estimate sets $D$ to the slope at $q = 3.5$ and $\tau = 2 - \sigma(1)/D$. Every grid runs a warm-up of $3L^2$ grains (BTW) or $2L^2$ (Manna) first, with fixed seeds, while the figure is on screen.
For the BTW rule the value of $\tau$ has been hard to pin down: about 1.22 (Manna 1991), exactly 6/5 by a scaling argument (Priezzhev, Ktitarev and Ivashkevich 1996), 1.29 on grids up to 4,096 wide (Lübeck and Usadel 1997) and 1.27 (Chessa and colleagues 1999). The right panel shows one reason. The Manna slope is near 2.7 by $q = 1.5$ and moves by only 0.06 from there to $q = 3.75$, so a single $D$ describes every moment. The BTW slope climbs from 2.35 at $q = 1$, rises by 0.16 over the same stretch, and has not levelled off at the last point. Lübeck (2000) saw the same bend on grids up to 4,096 wide and concluded that the two rules belong to different universality classes; Tebaldi, De Menech and Stella (1999) had already found the BTW toppling numbers multifractal, with no single exponent. Chessa and colleagues had read their moment analysis the other way, as one class with $D \approx 2.7$ and $\tau = 1.27$ for both. On these small grids the moment estimate gives $D = 2.62$, $\tau = 1.27$ for BTW and $D = 2.73$, $\tau = 1.31$ for Manna, whose $\tau$ Lübeck and Usadel put at 1.275. Both rules give $\sigma(1) \approx 1.9$, a little short of the exact value 2: a grain wanders about $L^2$ steps before it falls off the edge, so the mean avalanche grows as $L^2$.
Some stable piles can never come back once the drive has started. Two neighbouring cells that both hold 0 are an example: a cell drops to 0 only by toppling, and toppling either cell gives the other a grain. The piles the drive keeps returning to are called recurrent; the rest are transient. Majumdar and Dhar (1992) gave a test. Light a fire at the edge and burn any cell holding at least as many grains as it has unburnt neighbours; repeat. If every cell burns the pile is recurrent, and whatever refuses to burn is a forbidden pattern like the pair of 0s (Dhar 1999).
Recurrent piles form a group under adding cell by cell and relaxing. Its identity $e$ satisfies $(e + c)^{\circ} = c$ for every recurrent $c$, where $^{\circ}$ means relax, and takes two relaxations to compute: put 6 grains on every cell and relax, subtract the result from 6, and relax again.
Figure 4. Left: a 48×48 pile, shaded by grain count as in Figure 1. Right: the burning test, each cell coloured by the round in which it burns (faint early, strong late, on a fixed square-root scale up to 180 rounds); red cells never burn. Driven pile is a pile after 6,912 random grains (seed 7). Add identity adds $e$ and relaxes: a recurrent pile comes back unchanged, and any other pile comes back recurrent, since adding grains to a recurrent pile keeps it recurrent.
The identity is mostly 2s and 3s (85% of cells on this grid), a central square of 2s ringed by 3s and threaded with lines of 0s and 1s. A random stable pile is never recurrent in practice. Dhar (1990) counted the recurrent piles: their number is the determinant of the toppling matrix, about $10^{1178}$ for this grid against $4^{2304} \approx 10^{1387}$ stable piles. A recurrent pile also needs at least one grain per bond between neighbouring cells (Dhar 1999), 4,512 here, so its mean height is at least $2(L-1)/L$: 1.96 on this grid, close to 2 on large ones. The pile in Figure 1 settles a little above that floor, at 2.11.
Drop the grains on a single cell instead. Figure 5 piles $N = 2^n$ grains on the centre of an unbounded grid and relaxes them. The pile spreads into a disc whose radius grows as $\sqrt{N}$, and the inside fills with patches of constant height at many scales. Pegden and Smart (2013) proved that the picture, shrunk by $\sqrt{N}$, converges to a single limit as $N$ grows, which they characterised through an obstacle problem from the theory of elliptic equations.
Figure 5. Left: $2^n$ grains relaxed from the centre cell in parallel sweeps (every cell with $h \ge 4$ sends $\lfloor h/4 \rfloor$ grains to each neighbour), drawn as it spreads and scaled to $\sqrt{N}$, so a larger $n$ adds detail without growing. Right: radius (the farthest cell holding a grain) against $N$ for each $n$ computed so far, with a dashed guide of slope 1/2. $n = 17$ needs about 300 million topplings and takes a few seconds.
The toppling counts are a check on the Abelian property at scale: relaxing $2^{16}$ grains in sweeps takes 77,107,818 topplings, and so does toppling cell by cell from a stack.
Bak, Tang and Wiesenfeld introduced self-organized criticality as an explanation for 1/f noise, the flicker whose power falls off inversely with frequency and turns up in the flow of the Nile, the light from quasars and highway traffic. The sandpile has two time series whose spectra you can check.
The first is the avalanche size after each grain. Each grain lands at a random site on a stationary pile, so the series has almost no memory and its spectrum is flat, apart from a dip at the lowest frequencies: a big avalanche drains the pile, and the drops after it are a little quieter. The second counts topplings per step of the avalanche's own clock, with avalanches laid end to end. Jensen, Christensen and Fogedby (1989) showed that this flow falls as 1/f², since a sequence of isolated bursts with a broad spread of durations gives 1/f² above the inverse of the longest ones. Laurson, Alava and Zapperi (2005) later found the exponent significantly below 2, 1.59 for this model, equal to the exponent that relates an avalanche's size to its duration.
Time series. Avalanche size for 1,000 successive grains.
Power spectrum. Fitted slope: –. Gray guides have slopes 0, −1 and −2.
Figure 6. A fresh 60×60 pile (run 1), brought to its stationary state with 10,800 unrecorded grains, then measured. Size per grain records the avalanche size for each of 16,384 grains. Topplings per step runs 8,000 grains with parallel updates and records how many cells topple at each step, with one quiet step for every grain, so grains that topple nothing add quiet steps between avalanches. Spectra are averaged over segments and binned logarithmically; the red dashed line is a least-squares fit over the range it spans.
The size per grain comes out flat, with a fitted slope near 0. The flow of topplings has a slope near −1.6, matching Laurson and colleagues, and flattens at the lowest frequencies, on time scales longer than the longest avalanches. Neither signal is 1/f.
In How Nature Works (1996) Per Bak went on to argue that earthquakes, stock market crashes and mass extinctions follow sandpile-like power laws, and that self-organized criticality is how much of nature operates. Whether or not that holds, the sandpile is the clearest example of the idea: no temperature or coupling is set to a critical value, and the pile fills up to its stationary state and stays there. The next essay turns from avalanches to growth, with clusters built one random walker at a time.