A sandpile's avalanche areas follow one scaling law, and their sizes don't

The Bak–Tang–Wiesenfeld sandpile is the textbook example of self-organized criticality. Drop grains one at a time on a grid. Any site holding 4 grains topples and sends one grain to each neighbour, and grains that fall off the edge are lost. The pile settles into a state where avalanches of every size happen, and the usual story is that their sizes follow a power law. I simulated it on grids from 32×32 to 1024×1024 to see how cleanly that story holds. For the area of an avalanche (how many distinct sites toppled), it holds very well. For the size (how many topplings in total), it doesn't: no single exponent pair describes it, and a straight-line fit to the histogram gets the exponent badly wrong.

What "one scaling law" means here

The simple scaling picture says the distribution of avalanche sizes s on an L×L grid looks like P(s) = s^−τ f(s / L^D). That's a power law with exponent τ, cut off at a size that grows like L^D. If that's true, the moments ⟨s^q⟩ grow like L^σ(q) with σ(q) = D(q + 1 − τ) once q is large enough. So σ is a straight line in q, and its slope is D. Measuring σ(q) at several q is a sharper test than fitting the histogram, because it doesn't depend on picking a fitting window by eye.

There's also one exact check. Every grain dropped does a random walk through topplings until it falls off an edge, so the mean avalanche size must grow exactly like L², which means σ(1) = 2. My runs give ⟨s⟩/L² = 0.0396 at L=32, settling to 0.0354 at 512 and 1024, and σ(1) = 2.000 between the two largest grids.

Area: a straight line

Here is the slope of σ(q) between neighbouring q values, measured between grid sizes 512 and 1024, with jackknife error bars over 10 independent seeds:

For area the slope is the same at every q, and it's 2. That makes sense: a big avalanche covers a compact patch of the grid, so its area scales like L². The intercept then gives τ_a ≈ 1.20. The same holds at smaller grids too, where the slope is 2.06 at 32/64 and 2.04 at 128/256, closing in on 2 as L grows.

Size: the slope keeps climbing

For size, the slope rises with q, from 2.71 to 3.00, which is four error bars apart. A single (τ, D) pair can't produce that. It also rises with grid size at fixed q. Between q = 1 and 1.5 it's 2.45 at 32/64, 2.62 at 128/256 and 2.71 at 512/1024, so this isn't a small-grid effect that has already died out. Fitting a straight line through q = 2 to 4 gives D ≈ 2.94 and τ ≈ 1.38, while combining that D with the exact σ(1) = 2 gives τ ≈ 1.32. Two routes to the same exponent that don't agree are another sign that the simple picture fails.

The reason is that size counts topplings, not sites. In a big avalanche, sites near where the grain landed topple many times over, in roughly concentric waves, while sites at the edge of the patch topple once. Large avalanches pile up more of these layers than small ones do, so the size distribution doesn't keep one shape as it stretches. Tebaldi, De Menech and Stella pointed this out for the BTW model in 1999, calling it multifractal scaling. These numbers are my small-scale replication of it.

The histogram fit is the trap

The usual first move is to bin the sizes logarithmically, plot them on log-log axes, and fit a line through the middle. I fit from s = 10 up to 1% of L², and the slope drifts with grid size:

It's still climbing at 1024 and sits well below what the moments suggest. Each plot looks like a decent straight line by itself, so a single grid size gives no warning. Only comparing grid sizes shows the drift. The area histogram drifts the same way (0.98 to 1.12 from 128 to 1024) even though area scaling is clean. So the histogram slope is a poor estimator of τ at these sizes no matter what, and the moment analysis is what actually tells the two quantities apart.

So what

"Avalanches follow a power law" is true enough to be a good story and loose enough to hide that two natural ways of measuring an avalanche scale differently. Area behaves like a textbook critical system. Size needs more than one exponent. If I'd only fit a histogram at one grid size, I'd have reported τ ≈ 1.1 and never noticed either thing.

Method: about 70 lines of C on one Mac, open boundaries, grains dropped at uniformly random sites. Each run starts from the all-3 configuration and discards L² drops before measuring. Quadrupling that warm-up at L=256 changed no moment by more than 0.1%. 10 seeds per grid size, 40 million measured drops per seed at L ≤ 128, 20 million at 256, 4 million at 512 and 1 million at 1024. That's 63 runs in 48 minutes on 9 cores. The error bars are jackknife over seeds and don't include any allowance for corrections to scaling, which the drift with L says are large.