Measuring a Universal Constant on a Laptop: The Site Percolation Threshold

A Constant With No Formula

Some constants can be derived: the bond percolation threshold of the square lattice is exactly one half. Kesten proved it in 1980 by a self-duality argument, complete and beautiful [1]. Site percolation on the same lattice has no such argument, and nobody has found a closed form. Nobody has proved that none exists either, so the value is known only because people keep measuring it, over and over, with better and better computers.

Newman and Ziff reported 0.59274621(13) in 2000 [2]. Lee reported 0.59274598(4) in 2008, after an extended study of which pseudo-random generators can be trusted at that precision [3], and in 2014 Jacobsen obtained 0.59274605(3) by a graph-polynomial method, not Monte Carlo at all [4]. We took the last of those: newest, most precise, and reached by a route that shares none of our machinery.

Our question was narrower than theirs. More interesting for a school club, we think. Not what \(p_c\) is; eight places settle that. How close can a laptop get in one evening, and can we say honestly how close we got?

The second number is the one we care about. An estimate without an honest error bar is a guess with decimals attached. Getting 0.5928 is easy. The work is knowing we get four decimal places, and that the fourth is shaky.

What the Lattice Is, and What It Is Not

The club has no laboratory. No porous rock here, no filter paper, no forest and no epidemic. We have a room with computers in it, so our research is computation, and the computation is the experiment. Every number below came out of a program we wrote, running on an array in memory, so when this article says "measured", it means measured from our own simulated lattices. Flip the coin ten thousand times. Reasoning about coins is a different activity.

The model is site percolation, the thing Broadbent and Hammersley coined the word for in 1957, while thinking about fluid spreading through a random medium [5]. Ours is the plainest version of it:

• An \(L \times L\) square grid. Each cell is occupied with probability \(p\), independently of every other cell.
• Two occupied cells connect if they share an edge. Diagonal neighbours do not.
• Boundaries are free. Nothing wraps around.
• The lattice spans when a single cluster touches both end rows. Left-to-right crossings do not count.

Each of those four choices is a modelling decision. Each one changes the numbers: counting diagonal neighbours as connected would produce a different constant entirely. Periodic boundaries would converge toward the same limit faster, and counting a crossing in either direction would raise the crossing probability at every finite \(L\). We come back to the damage in §9.

The model leaves out most of what makes real spreading interesting. Sites here are independent. No clustering, no correlated disorder. The model holds no time in it: the lattice is static, and the only question is connectivity. Nothing flows along the path once it exists, so this says nothing about conductivity or flux, and real porous media are three-dimensional anyway, where the exponents differ and no exact values are known. We have a clean mathematical object that some real systems resemble, and the resemblance is a claim about universality classes, not a claim that our grid is a rock.

The Sweep

The naive way to measure a crossing probability. Pick a \(p\), generate a lattice, flood-fill it, record the verdict, repeat. Forty values of \(p\) cost forty times the work. Newman and Ziff found the better way [2][7], and the whole feasibility of this study rests on it.

Take a uniformly random permutation of all \(N = L^2\) sites and add them one at a time. After each addition, union the new site with any neighbours already present. The structure is a disjoint-set forest with union by size and path compression, the one Tarjan analysed in 1975 [8] and the ancestor of the cluster-labelling method Hoshen and Kopelman published a year later [9]. Each root carries two bits: top row touched, bottom row touched. Every union ORs the bits together. The first root to carry both marks the moment of spanning.

Record \(n_{\mathrm{span}}\), the number of occupied sites at that instant. Spanning is monotone. Adding sites can never destroy a crossing, so one sweep hands you the whole microcanonical curve

Converting that to the fixed-\(p\) curve anybody actually wants is an identity, not an approximation. If sites are independent with probability \(p\), the number occupied is Binomial\((N, p)\), so

One sweep therefore contributes to every point of the whole curve, hence 320 sweeps at \(L = 1024\), less absurd than it sounds. The same fact makes the sweeps our sample size, since points along a curve are not independent of one another. They are the same sweeps viewed through different binomial weights.

The convolution window in \(n\) holds essentially all the binomial mass. The narrowest carried 0.999999999939 of it, at \(L = 512\). We print that number in the output for a reason: a silently truncated convolution would bias every curve the same way, and would be close to impossible to spot afterwards.

Checking the Code Against Itself

Union-find code is short and fast. Path compression makes it shorter still, and easy to break in ways nobody sees. The failure mode that frightened us is not a crash. The dangerous version merges clusters correctly 99.99% of the time and shifts the threshold estimate by two parts in ten thousand, which is exactly the size of the effect we were trying to measure. So the code checks itself four ways before reporting anything, and the run halts if the first two fail.

Against an independent flood fill. 11,400 random configurations, six lattice sizes, \(p\) drawn uniformly from 0.30 to 0.80. Both the spanning verdict and the size of the largest cluster were computed twice: once by the union-find, once by a breadth-first search written separately, sharing no code with it. Mismatches: zero. Both quantities, every time.

Against the sweep's own claim. The sweep asserts that the lattice first spans at site \(n_{\mathrm{span}}\). So rebuild from the first \(n_{\mathrm{span}}\) sites of the permutation. Ask the flood fill whether it spans. It must. Remove the last site and ask again. It must not. Over 850 sweeps at three sizes, both conditions held every time.

Against exact enumeration. A 3×3 lattice has \(2^9 = 512\) configurations and a 4×4 has \(2^{16} = 65{,}536\), so brute force enumerated all of them, which gives the crossing probability as an exact polynomial in \(p\). Then the whole pipeline ran on those tiny lattices, 200,000 sweeps each. At \(p = 0.59274605\) the exact 3×3 answer is 0.566703307. We got 0.566421794. A gap of 0.25 standard errors. Worst of the ten comparisons: 0.50 standard errors. This check tests everything at once. A bug anywhere in the chain surfaces here.

Against conformal field theory. Cardy's 1992 formula gives the probability of a crossing at the critical point in any conformally mapped shape [10], and for a square of aspect ratio one that probability is exactly 1/2, independent of size. Smirnov later proved it rigorously for triangular-lattice site percolation [11]. So \(R_L(p_c)\) should converge to 1/2 as \(L\) grows. Panel A of Figure 3 shows what we got: at \(L = 16\) we read 0.5165, or 6.60 standard errors above a half. Not a bug. The finite-size correction does precisely what theory says it should. By \(L = 256\): 0.5005, or 0.08 standard errors out.

The Arithmetic

The calculation, with nothing hidden in it.

Seven lattice sizes: 16, 32, 64, 128, 256, 512, 1024. Sweep counts in the same order: 40,000, 40,000, 20,000, 8,000, 5,000, 1,200 and 320. Total 114,520 sweeps. Run every sweep to the last site: 1,241,989,120 site placements. The campaign took 399.9 seconds. Validation and the convergence study took the rest of the 526.5.

For each size, bisection finds the \(p\) where \(R_L(p) = 1/2\). Call it \(p^*(L)\). Finite-size scaling says

so \(L^{-1/\nu} = L^{-3/4}\). Plot \(p^*(L)\) against \(L^{-3/4}\). The intercept is \(p_c\). We did not fit the exponent \(\nu = 4/3\). It descends from den Nijs's 1979 conjecture for the thermal exponent of the \(q\)-state Potts model [12], confirmed by Nienhuis's Coulomb-gas work [13] and made rigorous for the triangular lattice by Smirnov and Werner [14].

Six values of \(L^{-3/4}\), 0.074325 down to 0.005524. Six matching values of \(p^*\). Weighted least squares, weights \(1/\sigma^2\). Intercept 0.5928171. Slope \(a = -0.007456\). Chi square 2.41 on 4 degrees of freedom, probability 0.66. The line fits about as well as lines ever fit six points.

The error bar comes from the sweeps and nothing else. Four hundred bootstrap replicates. Resample the sweeps at every size, independently and with replacement. Redo the convolution and the bisection at all six sizes. Redo the fit. The standard deviation of those 400 intercepts is 0.0001110. Subtract: \(0.5928171 - 0.59274605 = +0.0000710\). Divide by the error bar: \(+0.64\).

The two shaded rows are weak, for opposite reasons. \(L = 16\) has plenty of sweeps and a badly broken scaling law. \(L = 1024\) has an excellent scaling law and 320 sweeps. Hence its \(R_L(p_c)\) of 0.4798. Furthest from a half, still only 0.73 standard errors out.

Working Notes From the Club Table

The Answer, With Its Error Bar

0.5928171, plus or minus 0.0001110.

The accepted value is 0.59274605, and the difference is \(+7.105\times10^{-5}\), or 0.64 of our own standard error. Half the time, a one-sigma interval lands about that far out, which is neither suspicious nor impressive. A difference of 0.02 sigma would worry us more than one of 2.0.

The relative precision is \(1.87\times10^{-4}\), or 187 parts per million, and our one-sigma interval runs from 0.5927061 to 0.5929281. It holds both 0.5927 and 0.5928, so the fourth decimal is not safe. Every digit after that is noise we have no right to. Newman and Ziff got 0.59274621(13) in 2000, using a C implementation of the same algorithm and several orders of magnitude more sweeps [2], and our estimate differs from theirs by \(7.1\times10^{-5}\), which is 0.64 of our error bar and 546 of theirs. When two measurements disagree, the honestly large error bar is not the wrong one.

Panel B matters more than it looks. A bootstrap that came back skewed, or bimodal, or with a suspiciously small spread, would mean the bisection was landing on a flat part of the curve, or the resampling was failing to capture the real variability. It came back looking like a well-behaved error bar, and the 95% interval [0.5925836, 0.5930203] holds the accepted value with room on both sides.

Watching It Settle

An estimate right only at the end is poor comfort. So we reran the pipeline on the first \(m\) sweeps only, for ten values of \(m\) from 50 up to 40,000. Each reduced set bootstrapped 100 times.

That curve is the honest story of the run. The first fifty sweeps produced 0.5939574, standard error 0.0006601, an estimate 1.8 standard errors high and essentially worthless. At \(m = 100\) the central value landed at 0.5927641. Closer to the accepted value than our final answer, by 0.0000181. Call that luck. The error bar of 0.0004677 says so. From \(m = 1000\) onward the estimate stays within 0.0001 of where it finally lands, while the error bar shrinks by only a factor of 1.8 across a fortyfold increase in \(m\), because \(L = 512\) and \(L = 1024\) have 1,200 and 320 sweeps between them and no amount of extra work at \(L = 32\) can fix that.

That flattening is the real limit of this study. We are not short of Monte Carlo samples. We are short of large lattices.

The Strongest Objection We Can Make

Now the case against our own number.

The scaling form \(p^*(L) = p_c + aL^{-1/\nu}\) is the leading term of an expansion rather than a law, and the next term decays faster, with an exponent Ziff and Newman measured at about 0.90 for one class of estimators [15], which we simply dropped. Our own table shows the consequence. The \(L = 16\) point sits 0.0026 below \(p_c\) while the single power law fitted to the larger sizes predicts about 0.0018 below, and including it drags the intercept up to 0.5930250 with a claimed error of 0.0000955, which is 2.92 standard errors from the accepted value. Without the chi square to warn us we would have published a number wrong by three of its own error bars, and the error bar would have been the smaller of the two, because adding a point always tightens a least-squares intercept. That is the trap, and it deserves saying out loud: a tighter error bar from a wrong model is worse than a wider one from a right model.

Second objection. We fixed \(1/\nu = 3/4\) rather than fitting it. Let it float and the chi square is minimised at \(1/\nu = 0.300\), an absurd value implying \(\nu = 3.33\), where the chi square is 2.065 against 2.408 at the correct exponent. A difference of 0.34 in chi square across that whole range means our data carry essentially no information about \(\nu\), and they constrain \(p_c\) only because \(p_c\) is the intercept and the points sit close to it. So we import a theoretical result to get our headline number, and a reader who does not accept \(\nu = 4/3\) should not accept our \(p_c\) either.

How much does the import cost? Refit at \(1/\nu = 0.60\): 0.5928692. At \(1/\nu = 1.00\): 0.5927658. That range spans 0.0001034, or 0.93 of our error bar. Notice the direction. The wrong exponent 1.00 lands closer than the correct 0.75. We call that coincidence, and we will not adopt it.

Third objection, and the one we think is largest. Free boundaries. Most published codes use periodic boundaries in the transverse direction, which removes the edge effect that dominates our small lattices, and we chose free boundaries only because they are the easiest thing to define and the easiest thing to explain. That same choice forced us to drop \(L = 16\). A different boundary condition, with everything else identical, would have given a different set of \(p^*(L)\), a different slope, probably a usable \(L = 16\) point, and a different final digit. We do not know what that digit would have been, because we did not run it. Saying so beats guessing.

Fourth. Our criterion is \(R_L(p^*) = 1/2\). Use the mean occupancy at first spanning instead, a different estimator of the same limit, and the same arithmetic yields 0.5928214, a difference of \(+0.0000754\) from the accepted value rather than \(+0.0000710\). The two routes agree to within a twentieth of an error bar, reassuring though not independent: both come from the same sweeps.

Where We Fall Short

The threshold came out well. The fractal dimension did not. We would redo that part first.

At the critical point the largest cluster is a fractal. Its mass scales as \(S_{\max} \sim L^{D}\) with \(D = 91/48 = 1.8958333\), a value that follows from the same Coulomb-gas machinery that gives \(\nu = 4/3\) [13] and is rigorous on the triangular lattice [14]. We measured \(S_{\max}\) at fixed occupancy \(n = \mathrm{round}(p_c N)\), all seven sizes. Then a straight line through the log-log plot.

Over all seven sizes we get \(D = 1.86068\), 0.035 low, with chi square 160.4 on 5 degrees of freedom. That chi square is a confession. A fit that bad means the model is not a pure power law over that range, and forcing one through anyway produces a slope that means nothing. Drop \(L = 16\) and \(L = 32\): \(D = 1.88231\), chi square 0.8 on 3. Still 0.0135 below the exact value. But now the fit describes the points.

So our fractal dimension disagrees with the exact value. A real disagreement, not a statistical fluctuation. The reason is not mysterious. The compensated plot has not flattened even at \(L = 1024\), so corrections to scaling are still contributing at our largest size, and a two-parameter power-law fit over a range that short will always come out low. Fixing it needs either lattices several times larger, which the laptop cannot reach in a night, or a fit that includes a correction term, which seven points cannot constrain. So we report 1.8823 and the 0.0135 shortfall, rather than restricting the fit quietly until the number looked right.

The contrast with the threshold is useful. \(p_c\) is an intercept, and intercepts forgive: the corrections we ignored push \(p^*(L)\) around, but they vanish as \(L\) grows and the extrapolation absorbs most of what is left. \(D\) is a slope, and slopes do not forgive, because a slow drift in the amplitude looks exactly like a change in the exponent, which is why one answer here is good to four decimals and the other to two. Same data and same code, one evening of computing.

Reproducing This

Four decimal places of a constant that has no formula. Under nine minutes, on a machine a student already owns. Getting the same numbers takes two lines.

Python 3.12.3 and numpy 2.4.2. The script sets the BLAS thread-count variables to 1 before importing numpy. Then it runs eight worker processes of its own. Expect about nine minutes. Ours printed total runtime: 526.5 s on the last line. A machine that does not throttle will be faster, and one with fewer than eight cores slower, roughly in proportion.

The master seed is 20251215, and every worker stream is a numpy.random.SeedSequence with entropy [20251215, L, chunk_index]. So the output does not depend on chunk scheduling, or on the number of workers. If your output differs from ours in any digit, either numpy changed its generator or something is wrong, and we would like to hear about it.

The last section of the output file is a machine-readable block: lines beginning CURVE, SIZE, CONV, FIT, BOOTHIST and NSPAN, and every figure in this article and every number in the table came from those lines and from nothing else.

Three things we would change next time. Periodic boundaries in the transverse direction, for the reasons in §9. A cell-to-cell threshold estimator, which converges without needing the correction exponent [15]. And an inner loop in a compiled language, because our bottleneck is not the algorithm but CPython's pointer-chasing, and a factor of fifty there would put \(L = 4096\) within reach and produce a fractal dimension worth printing. The rest we would keep. Especially the four validations, which cost about a fifth of the runtime and are the only reason we are willing to put a number in a headline.

One last thought about why this matters outside a physics course. Percolation is the standard model for spreading on a spatial network, and Grassberger showed in 1983 that the general epidemic process maps onto exactly this transition [16]. The threshold is what makes such models interesting. Below it nothing spreads. Above it everything does. The crossover is as sharp as our Figure 1 says. Jacobsen's value is 3,700 times more precise than ours [4]. Ours took nine minutes and a laptop, and we can tell you exactly how much to trust it.

References