Percolation Simulator

Picture a porous stone lowered into a bucket of water. Inside the stone is a maze of tiny passages, some open and some blocked, and the question is which part of it gets wet. That picture is how one survey of the subject introduces percolation, and it is close to the form in which Simon Broadbent and John Hammersley posed the problem in a paper of 1957.

Their model is almost insultingly simple. Take a regular lattice. Declare each link open with probability p and closed otherwise, every link deciding independently of all the others. Then ask whether the open links join up into something large. Nothing in the rules mentions distance, flow or pressure. There is only a coin, tossed once for each link.

What makes the model worth sixty years of work is that the answer changes abruptly. Below a certain value of p the open links form small islands. Above it there is a connected piece that runs across everything. The dividing value is the percolation threshold, and the argument of this essay is that the threshold is easy to state, hard to compute, and impossible to see directly, because it belongs to an infinite lattice and nobody has ever simulated one.

Two ways to toss the coin

There are two standard versions. In bond percolation every site of the lattice is present and the coin is tossed for each link between neighbours. In site percolation the coin is tossed for each site, and two open sites are joined when they are neighbours. The lattice matters too. On the square lattice a site has four neighbours; on the triangular lattice it has six.

These choices change the threshold, and they change how much is known about it. For bond percolation on the square lattice the threshold is exactly 1/2. The lower half of that statement was proved by Harris in 1960, who showed that no infinite cluster exists at 1/2, and the harder upper half by Harry Kesten in 1980. The square lattice is its own dual, which is what makes 1/2 the natural candidate, but turning that symmetry into a proof took twenty years.

Site percolation on the triangular lattice also has threshold 1/2, and it is the case where the deepest results are proved. Site percolation on the square lattice has no such luck. Its threshold is known only from computation. A large simulation reported in 2001 gave 0.59274621 with an uncertainty of 0.00000013, and a later calculation by a different method gives fourteen digits, beginning 0.5927. No formula for the number is known.

So of the three lattices offered by a typical simulator, two have thresholds that are theorems and one has a threshold that is a measurement. It is worth keeping the two kinds apart, because a simulation can check a theorem but can only add to a measurement.

Why a finite grid has no threshold

Now take a grid you can actually hold in memory, say 128 sites on a side, and ask a concrete question: is there a path of open sites from the top row to the bottom row? Call that spanning. For any p strictly between 0 and 1 the answer is sometimes yes and sometimes no, depending on the coins.

The smallest case shows why. On a grid of 2 by 2 sites, a path from top to bottom needs both sites of one column open. Counting the sixteen configurations gives a spanning probability of exactly 2 p squared minus p to the fourth. At p equal to 1/2 that is 0.4375. It is a smooth polynomial, rising from 0 to 1 with no jump anywhere.

The same is true at any finite size: the spanning probability is a polynomial in p, of higher degree, and a polynomial cannot jump. What happens instead is that the rise gets steeper as the grid grows. Near the threshold the typical cluster radius grows like a power of the distance from the threshold, with an exponent of 4/3 in the plane, so the width of the rise shrinks roughly as the side of the grid to the power of minus 3/4. Double the side and the rise narrows by a factor of about 1.7.

Only in the limit of an infinite grid does the curve become a step. The threshold is the position of that step. Everything a simulation reports is therefore a statement about a finite grid, and the honest quantity to quote is not the threshold but something like the p at which half of all grids of a given size span.

That halfway point is not fixed either. It drifts as the size changes and it depends on what you call spanning. Count a crossing in either direction instead of top to bottom, and the halfway point starts well below the threshold and creeps up as the grid grows.

A few shapes are special. On a rhombus cut from the triangular lattice exactly one of two things happens: open sites cross from top to bottom, or closed sites cross from side to side. At p equal to 1/2 the two are equally likely, so the crossing probability is exactly 1/2 at every size, a fact that can be confirmed by counting every configuration of a small rhombus. Cardy predicted the crossing probabilities of rectangles at the threshold from a conformal symmetry, and Smirnov proved his formula for this lattice in 2001.

A worked example with its error bars

Here is a run you can repeat. Choose site percolation on the square lattice, open sides, spanning from top to bottom, a grid of 128 by 128 and seed number 2026. Raise p slowly. This grid first spans at p equal to 0.601. At that moment 9,827 of its 16,384 sites are open, they form 446 clusters, and the spanning cluster holds 5,503 sites. One grid tells you about one seed. The number 0.601 is exact for that seed and says very little about the lattice. To learn about the lattice you need many seeds, and then you need to say how sure you are.

Run 500 seeded grids at each of four sizes. The p at which half of them span comes out at about 0.585 for 16 by 16, 0.593 for 32 by 32, 0.594 for 64 by 64 and 0.593 for 128 by 128. The width of the rise, from a quarter of grids spanning to three quarters, falls from 0.081 to 0.047, 0.033 and 0.020. The curves steepen exactly as the argument above says they should.

Each of those halfway points carries an interval. At a fixed p, each grid either spans or does not, so the count of spanning grids is binomial, and a proportion from a binomial count has a well-studied confidence interval. The one used here is the Wilson score interval, which a government statistics handbook recommends because its worth depends little on the sample size or the proportion and its lower limit cannot go negative. For 128 by 128 with 500 trials the interval for the halfway point runs from 0.591 to 0.594.

Three digits is all that interval supports, so three digits is all that should be written. With 10,000 trials at 256 by 256 the same procedure gives about 0.5928 with an interval from 0.5926 to 0.5930. That agrees with the published 0.5927, and it is still a statement about a 256 by 256 grid.

The cluster at the edge of spanning

Stop the example grid at 0.601 and look at the spanning cluster. It is not a solid block. It is a stringy thing full of holes, with dangling branches and narrow necks where a single closed site would cut it in two. It occupies about a third of the grid and yet reaches across all of it.

This raggedness is not an accident of the seed. At the threshold the large clusters are fractal: the number of sites within a distance r grows as r to the power 91/48, a little under 2, so the cluster fills a vanishing fraction of the plane as it grows. The sizes of the other clusters follow a power law as well, with no typical scale, which is why a chart of cluster sizes on logarithmic axes comes out close to a straight line near the threshold and bends sharply away from it elsewhere.

Move p up by a few hundredths and the picture changes character. The spanning cluster swallows its neighbours, the holes close, and the remaining clusters are small. Move p down and the big cluster falls apart into pieces of ordinary size. The interesting geometry lives in a narrow window, and the window gets narrower as the grid gets larger.

What the model is used for, and what it cannot do

The original application was fluid in a porous material, and percolation is still used for the flow of fluids through rock. The same sums have been borrowed widely. Forest fires and epidemics are natural fits, because both spread from neighbour to neighbour with some probability. One standard model of an epidemic on a network turns out to be exactly bond percolation, with the chance of transmission playing the part of p and the final size of an outbreak equal to the size of a cluster.

Networks are the other large use. Remove routers from a network at random and ask when it falls into disconnected pieces: that is site percolation on a graph that is not a lattice. For networks whose connections follow a broad power law the answer is striking. One calculation for the structure of the Internet found that more than 99 percent of nodes would have to fail at random before it came apart, although the same networks are fragile when the best connected nodes are removed on purpose.

A lattice simulator says nothing reliable about any of these real systems. Rock is not a square grid, pores are not independent, and a forest does not burn with one fixed probability. What the simulator teaches is the shape of the phenomenon: a smooth curve that sharpens into a step, a threshold that depends on the lattice and on the kind of coin, and a critical cluster with a geometry of its own.

It also cannot settle the fourteenth digit of anything. An estimate from a few thousand small grids is good to three or four digits, and its error bar is as much a part of the result as the number. The useful habit is the one the experiment forces on you: say which lattice, which rule, which size, how many trials, and how wide the interval is.