The BTW Model: How the Sandpile Works

The Rules

The Bak-Tang-Wiesenfeld (BTW) sandpile model operates on a two-dimensional square lattice of L x L cells. Each cell i holds a non-negative integer z_i representing its “height” — the number of sand grains at that position. The dynamics are governed by two operations: addition and toppling.

Addition. A single grain is added to a randomly selected cell: z_i <- z_i + 1. The cell is chosen uniformly at random from the entire grid. This is the external driving.

Toppling. If any cell has z_i >= z_c (the critical threshold, equal to 4 for a square lattice — the number of orthogonal neighbors), it topples:

z_i <- z_i - 4

z_neighbor <- z_neighbor + 1 for each of the four orthogonal neighbors.

If a neighbor is outside the grid boundary, the grain is lost from the system. The boundary is open — it acts as a drain that permanently removes grains. For interior cells, the toppling is conservative: the cell loses 4 grains and distributes exactly 4 grains to neighbors. For boundary cells, some grains leave the system, making the boundary the system’s dissipation mechanism.

The critical rule governing timescale: toppling is instantaneous relative to addition. After a grain is added, all topplings triggered by that addition are resolved — the entire avalanche completes — before the next grain is added. This separation of timescales is not a computational convenience. It is a necessary condition for the model to exhibit self-organized criticality. If grains were added continuously during avalanche propagation, the system would not self-organize to a critical state.

The update order within an avalanche can be parallel (all cells above threshold topple simultaneously at each sub-step) or sequential (cells are toppled one at a time in some order). For the BTW model, Deepak Dhar proved in 1990 that the final configuration after an avalanche is independent of the toppling order — the abelian property. Whether you topple cells left-to-right, right-to-left, randomly, or all at once, the same final state results. The avalanche’s trajectory (which cells topple in which order) may differ, but the endpoint is unique.

The threshold value z_c = 4 is not arbitrary — it equals the coordination number of the square lattice (the number of neighbors). On a triangular lattice, the threshold would be 6; on a hexagonal lattice, 3. The key structural requirement is that each toppling redistributes exactly one grain to each neighbor. This conservation property is what makes the model abelian and analytically tractable.

What Happens When It Runs

Starting from an empty grid (all z_i = 0), the system passes through three qualitative stages.

Transient filling. Grains are added one at a time. Initially, no cell reaches the threshold, and every addition simply increments a random cell. No topplings occur. The average height increases linearly with the number of additions. This phase is subcritical: the system is far from the threshold everywhere, and perturbations (grain additions) have only local effects.

Approach to criticality. As grains accumulate, cells begin reaching the threshold. The first topplings occur — small, localized events that redistribute a few grains. Each toppling may push a neighbor past threshold, creating a small chain reaction. Avalanches are still predominantly small, but their average size grows. The average height of the grid approaches a value near (but below) the threshold. The system is approaching its critical state.

Steady state (criticality). After sufficient additions (on the order of L^2 for an L x L grid, though the exact transient time depends on initial conditions), the system reaches a dynamic equilibrium. The average number of grains entering the system per addition (one grain) equals the average number leaving through the boundary per addition. Large avalanches periodically sweep sections of the grid, dissipating grains through the boundary and keeping the average height near the critical value. The system has self-organized to criticality.

The steady-state configuration is not uniform. It is a complex, heterogeneous arrangement with long-range spatial correlations. Cells near the boundary tend to have lower heights (they lose grains to the boundary more easily). Cells in the interior tend to have higher heights. But the height field is not a smooth gradient — it has structure at all scales, reflecting the history of avalanches that have shaped it. Visually, the steady state looks like a rough landscape, not a smooth hill.

A crucial feature: the steady state is not a single fixed configuration. The system is dynamically critical — it fluctuates around the critical state, never settling to a single configuration but maintaining statistical stationarity. The distribution of heights, the spatial correlations, and the avalanche statistics all become stationary (time-independent) in the steady state, even though the specific configuration changes with every grain addition.

Avalanche Statistics

An avalanche is defined as the set of all topplings triggered by a single grain addition. Three quantities characterize each avalanche:

Size (s). The total number of topplings in the avalanche. A cell that topples multiple times (receiving grains from multiple neighbors that push it past threshold repeatedly) contributes each toppling to the count.

Area (a). The number of distinct cells that topple at least once during the avalanche. Area measures the spatial extent of the event.

Duration (T). The number of parallel update steps in the avalanche. In parallel updating, all cells above threshold topple simultaneously at each step; the duration counts the number of such steps until all cells are below threshold.

In the steady state, the probability distributions of these quantities follow power laws:

P(s) ~ s^(-tau_s) with tau_s approximately 1.1-1.2

P(a) ~ a^(-tau_a) with tau_a approximately 1.3-1.4

P(T) ~ T^(-tau_T) with tau_T approximately 1.5

These exponents are measured numerically through large-scale simulations. Their exact values have been debated — finite-size effects, logarithmic corrections, and the question of whether the BTW model produces exact power laws or merely approximate ones have generated a substantial literature. Priezzhev (1996) computed exact results for specific avalanche quantities in one dimension. In two dimensions, analytical calculation of the exponents remains an open problem.

The power-law distribution means there is no characteristic avalanche size. The system produces avalanches of all sizes, from single-cell adjustments to system-spanning cascades, with a probability that decreases as a power of size. This is the signature of criticality — the absence of a preferred scale. It is also what distinguishes the sandpile’s behavior from simpler threshold systems that produce exponential size distributions (where large events are exponentially rare rather than power-law rare).

The scaling relations between the exponents are predicted by finite-size scaling theory: a ~ s^(d_a) (area scales as a power of size) and T ~ s^(d_T) (duration scales as a power of size), where the scaling exponents d_a and d_T relate the avalanche exponents through standard scaling relations. These relationships have been verified numerically and provide internal consistency checks on the measured exponents.

The Abelian Property

Deepak Dhar’s 1990 paper, “Self-Organized Critical State of Sandpile Automaton Models” in Physical Review Letters, proved a remarkable property of the BTW model: the final configuration after an avalanche is independent of the order in which cells topple. This is the abelian property, and it has profound consequences for the model’s mathematical analysis.

The name comes from algebra: the toppling operators commute. If cells A and B are both above threshold, toppling A then B produces the same final state as toppling B then A. This extends to any number of cells and any order of toppling. The abelian property is not generic — most cellular automata do not have it. It depends on the specific structure of the BTW toppling rule: each toppling adds exactly one grain to each neighbor, and the threshold is exactly the number of neighbors.

The mathematical consequences are substantial. Dhar showed that the set of configurations reachable from any initial state through the addition-and-toppling dynamics — the recurrent configurations — forms an abelian group under the operation of grain addition. The group structure is isomorphic to the cokernel of the graph Laplacian, connecting the sandpile to algebraic graph theory. The number of recurrent configurations equals the determinant of the reduced Laplacian (the Laplacian matrix with one row and column removed), which equals the number of spanning trees of the grid by Kirchhoff’s matrix tree theorem.

This connection to graph theory has made the abelian sandpile model a rich object of mathematical study beyond its physical interpretation. The “identity” element of the sandpile group — the recurrent configuration that acts as a neutral element under grain addition — produces visually striking fractal-like patterns when computed for large grids. The theory of chip-firing games on graphs (Biggs, 1999) is essentially the same mathematical object viewed from a combinatorial perspective.

For the physical understanding of SOC, the abelian property has a specific implication: it guarantees that the steady state is well-defined and does not depend on simulation details (the order in which you process above-threshold cells). This eliminates a class of potential artifacts and makes numerical results for the BTW model more reliable than for non-abelian cellular automata, where the update rule (parallel, sequential, random) can change the dynamics qualitatively.

Further Reading