Percolation Theory, in depth: site and bond thresholds, critical scaling, and the Newman-Ziff union-find algorithm

The model

Take a graph, usually a lattice. In site percolation each vertex is open with probability p; in bond percolation each edge is open with probability p. A cluster is a connected component of the open subgraph. On a finite L x L box the usual question is whether an open path crosses from the top row to the bottom row; on an infinite lattice it is whether an infinite cluster exists.

Let theta(p) be the probability that a given site belongs to the infinite cluster. It is zero for p below the threshold p_c and positive above it. On a finite box the spanning probability R_L(p) rises from 0 to 1 as an S-curve that sharpens as L grows; in the limit it becomes a step at p_c. Everything interesting about the model is in how fast that S-curve sharpens and what happens right at the step.

Known thresholds

Three patterns are worth remembering. Higher coordination lowers the threshold, because each open site has more chances to connect. For a given lattice the site threshold is at least the bond threshold, because opening a site is a coarser event than opening an edge. And on tree-like graphs the threshold is exactly where the expected number of new branches per step reaches one: a path entering a site of degree z can leave by z - 1 edges, so the cluster grows forever when (z - 1)p exceeds 1. The same branching argument gives the Erdos-Renyi giant component at mean degree 1 and, with degree heterogeneity, the Molloy-Reed condition used in epidemic models on networks.

Critical behaviour and finite-size scaling

Near p_c the model is described by power laws with universal exponents, identical for every two-dimensional lattice and bond or site variant. In two dimensions the infinite cluster density grows as theta(p) ~ (p - p_c)^beta with beta = 5/36; the correlation length, the typical size of finite clusters, diverges as xi ~ |p - p_c|^(-nu) with nu = 4/3; the mean finite cluster size diverges with gamma = 43/18; and the spanning cluster at p_c is a fractal of dimension 91/48.

The exponent that matters most for simulation is nu. A box of side L cannot tell p from p_c once xi exceeds L, so the transition window on a finite box has width proportional to L^(-1/nu) = L^(-3/4). Two consequences follow. Estimates of p_c from finite boxes are shifted by an amount of that order, so you extrapolate in L rather than trusting the largest box. And the sample-to-sample spread of the point where spanning first appears shrinks like L^(-3/4), which is a useful self-check on any implementation.

Right at p_c the crossing probability of a large square tends to 1/2 for both site and bond percolation, a consequence of Cardy's formula for crossing probabilities in rectangles. That gives a second estimator: the p at which R_L(p) = 1/2.

The Newman-Ziff algorithm

The naive method fixes p, samples a grid, runs BFS or DFS to test spanning, and repeats for many p and many samples. Each test costs O(N) for N = L^2 sites, and you pay it for every point on the curve. Hoshen and Kopelman's 1976 algorithm improved cluster labelling with a one-pass union-find scan, but still per p.

Newman and Ziff (2000, 2001) observed that you can instead add sites one at a time in a random order and keep clusters in union-find. After n additions you have a uniform random configuration with exactly n open sites, so one sweep produces the whole curve Q_n, the observable at every occupation count, for nearly linear total cost. Spanning is monotone in n, so it is enough to record the first n at which top and bottom join.

import random

def newman_ziff_site(L, rng):
    """Open sites in random order; return how many were open when the
    top row first connected to the bottom row (site percolation, L x L)."""
    N = L * L
    TOP, BOT = N, N + 1
    parent = list(range(N + 2))
    size = [1] * (N + 2)
    occupied = bytearray(N)

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]      # path halving
            x = parent[x]
        return x

    def union(a, b):
        ra, rb = find(a), find(b)
        if ra == rb:
            return
        if size[ra] < size[rb]:
            ra, rb = rb, ra
        parent[rb] = ra                        # union by size
        size[ra] += size[rb]

    order = list(range(N))
    rng.shuffle(order)
    for k, site in enumerate(order, start=1):
        occupied[site] = 1
        r, c = divmod(site, L)
        if r == 0:
            union(site, TOP)
        if r == L - 1:
            union(site, BOT)
        for nr, nc in ((r - 1, c), (r + 1, c), (r, c - 1), (r, c + 1)):
            if 0 <= nr < L and 0 <= nc < L and occupied[nr * L + nc]:
                union(site, nr * L + nc)
        if find(TOP) == find(BOT):
            return k
    return N

To turn counts into a curve in p, convolve with the binomial distribution: the probability of spanning at p is the sum over n of C(N, n) p^n (1 - p)^(N - n) times the fraction of samples that had spanned by n. Compute the binomial weights in log space; for large N they underflow otherwise.

from math import lgamma, log, exp

def spanning_curve(first_counts, N, ps):
    """first_counts: list of k values from newman_ziff_site, one per sample."""
    hist = [0] * (N + 1)
    for k in first_counts:
        hist[k] += 1
    cum, Q = 0, []
    for n in range(N + 1):
        cum += hist[n]
        Q.append(cum / len(first_counts))       # fraction spanned with n sites open
    curve = []
    for p in ps:
        lp, lq = log(p), log(1 - p)
        total = sum(exp(lgamma(N + 1) - lgamma(n + 1) - lgamma(N - n + 1)
                        + n * lp + (N - n) * lq) * Q[n] for n in range(N + 1))
        curve.append(total)
    return curve

Worked example: estimating the square-lattice threshold

Running newman_ziff_site 200 times per size, with one seeded generator (seed 42) shared across the sizes in order, gave the following fractions of open sites at first spanning on the authoring machine:

Read it the way the critical theory says to. The means bracket the accepted 0.5927 and move toward it as L grows, with the small box visibly biased. The spread falls by a factor of about 4.1 from L = 16 to L = 128; the L^(-3/4) prediction for an eightfold change in L is about 4.8, and with only 200 samples the standard deviation itself is uncertain by roughly 5 percent per size, so the ratio is roughly consistent, with the remaining gap plausibly due to corrections to scaling at small L. Before trusting the code, it was checked on 300 small grids against a brute-force BFS that replays the same random order and tests spanning after every addition; the first-spanning counts matched exactly. That replay test is the one to keep in your own suite.

For a real estimate you would push to L of 512 or more in a compiled language, take thousands of samples, and fit the shift in the mean against L^(-3/4) to extrapolate.

Fullness and the backwash bug

A common variant asks not whether the system percolates but which sites are full, meaning connected to the top. With both virtual nodes in one union-find, once the system percolates every cluster touching the bottom row joins the bottom node, which is joined to the top, so bottom clusters that never reach the top are wrongly reported full. This is the backwash bug familiar from introductory algorithms courses. The fix is a second union-find with only the virtual top, used for fullness queries, or a single structure that records per root whether it touches the bottom, avoiding the virtual bottom node entirely. The spanning-only code above is immune because it never asks about fullness.

Where percolation shows up

Network robustness. Removing random nodes from a network is site percolation on it. Networks with heavy-tailed degree distributions have a vanishing threshold for random removal in the infinite-size limit, so they shrug off random failures, but removing the highest-degree nodes first fragments them quickly, a result reported by Albert, Jeong and Barabasi and by Cohen and colleagues in 2000. Run both experiments, random and targeted, before trusting a topology.

Epidemics. The final size of an SIR outbreak maps to bond percolation with the transmissibility as the bond probability; the epidemic article develops this.

Grids and images. Connected-component labelling of thresholded images, reachability in procedurally generated maps, and porous-media flow are all percolation questions; the BFS grid techniques in BFS grid patterns answer a single instance, and Newman-Ziff answers the whole family over p.

Failure modes

• Trusting one box size. Finite-size shift is of order L^(-3/4); extrapolate.
• Mixing definitions. Top-to-bottom crossing, crossing in either direction, wrapping on a torus and largest-cluster fraction all have different finite-size curves.
• Recursive DFS on large grids. Clusters near p_c are huge and stack depth explodes; use iterative BFS or union-find.
• Union-find without balancing. Without union by size or rank, chains form and the sweep slows from near linear to far worse.
• Index overflow. In C or CUDA, L*L exceeds 32-bit range once L passes 46,340.
• Correlated seeds. Seeding parallel workers with consecutive integers from a weak generator correlates samples; use independent streams.
• Backwash. Fullness queries with a shared virtual bottom report false positives.

Trade-offs

Newman-Ziff gives the entire curve per sample and is the default for threshold and curve estimation. Fixed-p sampling with BFS is simpler and better when you need one p and rich per-cluster statistics. Hoshen-Kopelman scans row by row with memory proportional to one row of labels, which matters for very tall systems. Exact enumeration of all configurations is possible only for tiny boxes but gives exact polynomials useful for testing. Union-find variants with rollback, described in union-find variants, let you remove sites as well, at the cost of losing path compression.

Choose the estimator to match the question. The mean fraction at first spanning is cheap and needs no convolution, but its finite-size shift depends on the boundary conditions. The p at which R_L(p) crosses 1/2 uses the convolved curve and tends to converge faster for open square boxes. Wrapping on a torus removes edge effects altogether and is the usual choice in high-precision work, at the price of a slightly more involved wrap-detection step: each union must track the displacement between a site and its root so that a cluster meeting itself across the boundary is recognised.

What to do next

1. Implement newman_ziff_site and its BFS replay test; make the test pass before running anything large.
2. Measure the mean and spread of the first-spanning fraction for L = 16 to 256 and check the spread scales like L^(-3/4).
3. Add the binomial convolution and plot R_L(p) for three sizes on one chart.
4. Switch to bond percolation on the square lattice and confirm the estimate approaches 1/2.
5. Apply the same sweep to a real network: remove nodes randomly and by degree, and record the largest-component fraction.
6. Port the inner loop to a compiled language before going past L = 512.