Sandpile Model and Self-Organized Criticality — Running the BTW Model in the Browser

#What Is Self-Organized Criticality?

Drop grains of sand one at a time onto a sandpile. At first, nothing happens. Sand accumulates, slopes steepen, and at some point an avalanche occurs. Some avalanches are small, others engulf the entire pile.

This phenomenon where small inputs trigger unpredictably scaled chain reactions is called Self-Organized Criticality (SOC). Proposed in 1987 by Per Bak, Chao Tang, and Kurt Wiesenfeld, it drew attention as the mechanism behind the power laws found everywhere in nature -- earthquakes, forest fires, stock price fluctuations, and more.

The key is the "self-organized" part. Without external parameter tuning, the system naturally moves toward a critical state. The sandpile model (BTW model) is the simplest mathematical model expressing this phenomenon.

#BTW Model Rules

The rules are strikingly simple:

1. Each cell in a 2D grid holds a "grain count"
2. Drop one grain onto the center
3. Any cell with a grain count at or above the threshold (4) topples: subtract 4 from itself and distribute 1 to each of the 4 Von Neumann neighbors
4. Repeat until no more toppling occurs
5. Grains that fall off the edge are lost (open boundary conditions)

Before topple:    After topple:
  0 0 0             0 1 0
  0 4 0     →       1 0 1
  0 0 0             0 1 0

From just four rules, avalanches following a power-law distribution emerge.

#Core Logic: The topple Function

The heart of the implementation is the toppling process. The grid is managed as a 1D flat array using Uint8Array:

const TOPPLE_THRESHOLD = 4

function topple(grid: Uint8Array, cols: number, rows: number): number {
  let totalTopples = 0
  let unstable = true

  while (unstable) {
    unstable = false
    for (let y = 0; y < rows; y++) {
      for (let x = 0; x < cols; x++) {
        const idx = y * cols + x
        if (grid[idx] >= TOPPLE_THRESHOLD) {
          unstable = true
          const distribute = Math.floor(grid[idx] / 4)
          grid[idx] -= distribute * 4
          totalTopples += distribute

          // Von Neumann neighborhood (up/down/left/right)
          if (y > 0) grid[(y - 1) * cols + x] += distribute
          if (y < rows - 1) grid[(y + 1) * cols + x] += distribute
          if (x > 0) grid[y * cols + (x - 1)] += distribute
          if (x < cols - 1) grid[y * cols + (x + 1)] += distribute
        }
      }
    }
  }

  return totalTopples
}

Why loop with while (unstable)? Because a single scan may not settle all toppling. When a cell topples and distributes sand to neighbors, those neighbors may also exceed the threshold. This chain is the avalanche. The loop continues until the chain completely subsides.

Math.floor(grid[idx] / 4) handles cases where the threshold is greatly exceeded. When dropping 64 grains at once, a single cell ends up with a large amount of sand. Using grid[idx] -= 4 would require many loop iterations, but processing floor(n/4) topples at once speeds up convergence.

#The Significance of Open Boundary Conditions

How boundaries are handled carries important physical meaning:

// When an edge cell topples, distribution to the outside is simply lost
if (y > 0) grid[(y - 1) * cols + x] += distribute
if (y < rows - 1) grid[(y + 1) * cols + x] += distribute
if (x > 0) grid[y * cols + (x - 1)] += distribute
if (x < cols - 1) grid[y * cols + (x + 1)] += distribute

Just if guards. Grains destined for directions that fail the condition are simply discarded. This is the open boundary condition. Think of sand falling off the edge of the pile.

If the boundary were a torus (periodic boundary conditions), grains would never leave the system. Without energy dissipation, the system wouldn't reach a critical state, and SOC couldn't be observed. Open boundaries express that the system is dissipative, which is a necessary condition for SOC.

#High-Performance Rendering with ImageData

To render the grid every frame, rather than calling fillRect per cell, we directly manipulate Canvas API's ImageData:

const CELL_SIZE = 6

const GRAIN_COLORS: [number, number, number][] = [
  [17, 24, 39],    // 0 grains - dark background
  [59, 130, 246],  // 1 grain - blue
  [34, 197, 94],   // 2 grains - green
  [245, 158, 11],  // 3 grains - amber
]

function render(
  ctx: CanvasRenderingContext2D,
  grid: Uint8Array,
  cols: number,
  rows: number,
  ratio: number,
) {
  const w = cols * CELL_SIZE * ratio
  const h = rows * CELL_SIZE * ratio
  const imgData = ctx.createImageData(w, h)
  const data = imgData.data
  const scaledCell = CELL_SIZE * ratio

  for (let y = 0; y < rows; y++) {
    for (let x = 0; x < cols; x++) {
      const val = grid[y * cols + x]
      const color = GRAIN_COLORS[Math.min(val, 3)]

      const pxStartX = Math.floor(x * scaledCell)
      const pxStartY = Math.floor(y * scaledCell)
      const pxEndX = Math.floor((x + 1) * scaledCell)
      const pxEndY = Math.floor((y + 1) * scaledCell)

      for (let py = pxStartY; py < pxEndY; py++) {
        for (let px = pxStartX; px < pxEndX; px++) {
          const idx = (py * w + px) * 4
          data[idx] = color[0]
          data[idx + 1] = color[1]
          data[idx + 2] = color[2]
          data[idx + 3] = 255
        }
      }
    }
  }

  ctx.putImageData(imgData, 0, 0)
}

Colors are limited to 4 based on grain count. Math.min(val, 3) clamps values of 3 or above to the same color. In the stable state, each cell holds 0-3 grains, so 4 colors suffice. Values of 4+ occur only briefly during toppling and are virtually invisible at render time.

ratio is the devicePixelRatio. On Retina displays, double the pixels need to be drawn. Separating the Canvas's internal resolution from its CSS size is a fundamental Canvas rendering technique.

#Observing Power Laws

The highlight of the sandpile model is that the avalanche size distribution follows a power law:

P(s) ~ s^(-tau)

Small avalanches occur frequently; large avalanches occur rarely. But the distribution is power-law, not exponential. This means "impossibly large avalanches" occur at far higher probabilities than an exponential distribution would predict.

In the implementation, the return value of topple (totalTopples) is the avalanche size. Recording this in a histogram and plotting on a log scale reveals a linear distribution. The SandPile component simplifies this by displaying only the avalanche count and maximum size:

const topples = topple()

if (topples > 0) {
  avalanchesRef.current++
  if (topples > maxAvalancheRef.current) {
    maxAvalancheRef.current = topples
  }
}

When continuously dropping grains in Auto Drop mode, after the system reaches the critical state, the maximum size jumps sporadically. This is the essence of SOC: the system spontaneously maintains a "barely unstable" state.

#Auto Drop Timing Control

In Auto Drop mode, timing is controlled using the requestAnimationFrame frame count:

let frameCount = 0

const loop = () => {
  frameCount++

  if (autoDropRef.current && frameCount % Math.max(1, 60 - speed * 5) === 0) {
    const cx = Math.floor(cols / 2)
    const cy = Math.floor(rows / 2)
    dropGrains(cx, cy, dropAmount)
  }

  render()
  rafRef.current = requestAnimationFrame(loop)
}

Higher speed means a smaller modulo value, increasing drop frequency. Math.max(1, ...) prevents division by zero. Using rAF rather than setInterval synchronizes with rendering. The side effect of rAF pausing when the tab is inactive is also convenient.

#Why the Sandpile Model Is Fascinating

The fundamental insight of the BTW model is that complexity doesn't arise from complex rules. Four rules, one threshold, one neighborhood definition. That's all it takes to reproduce the power laws found in nature.

As is true for cellular automata in general, local rules generate global behavior. The Game of Life produces "life-likeness," SandFall produces "matter-likeness," and the sandpile model produces "criticality-likeness." All are built on the same skeleton of a grid and neighborhood rules.

#Tools Used

For grasping the big picture of self-organized criticality and complex systems, the books in the toolshelf serve as a reference. The mechanism by which complex behavior emerges from simple rules is closely connected to the software design philosophy of "loosely coupled modules composing complex systems."