Your Phase-Transition Simulation Must Pass a 1944 Exam First

A square-lattice Ising code that cannot find $T_c = 2.269\ J/k_B$ (Onsager's exact value [1]) and an order-parameter exponent of 1/8 is not ready to report a transition in any harder model. I hold this firmly. The evidence level here is derivation plus published results. I ran no new simulation for this essay, and I say below which numbers come from my earlier Lab run.

The answer key exists

The exact critical temperature of the nearest-neighbour square-lattice Ising model follows from Onsager's 1944 solution, published as Physical Review 65, 117 [1]:

$$\frac{k_B T_c}{J} = \frac{2}{\ln(1+\sqrt{2})}$$

Hand work: $\ln(1+\sqrt{2}) = 0.881374$, so $k_B T_c / J = 2.26919$. Units: $J$ is the coupling energy, so the ratio has no units. This is the value that Onsager's solution gives [1]. I derived the same value in an earlier post on duality in four lines, and this essay extends that post: there I found the number, here I use it as an exam.

Onsager announced the spontaneous magnetisation without proof. Yang derived it in 1952 (Phys. Rev. 85, 808) [2]:

$$m = \left(1 - \sinh^{-4}(2J/k_BT)\right)^{1/8}, \quad T < T_c$$

So $\beta = 1/8$ exactly, against the mean-field value 1/2 [2]. A code that returns something near 0.5 has a mean-field bug or a lattice too small to see anything. A code that returns 0.125 within error has passed a hard test, because the exponent is small and the window where it holds is narrow.

One caveat on the key itself. Later authors note that Yang's method makes an assumption that is not fully rigorous, and that the magnetisation has never been solved with complete mathematical rigour [2]. The free energy is rigorous [1]. I treat the magnetisation formula as a very well-confirmed answer, not as a theorem. That is still a far better reference than any simulation of a model with no exact solution.

Where this fails: the key covers only nearest-neighbour, zero-field, square-lattice spins. Change the lattice, add a field or add next-nearest couplings, and the exact numbers are gone.

What the exam asks, with order-of-magnitude costs

Start with an estimate. A Metropolis code near $T_c$ has autocorrelation time growing as $L^z$ with $z \approx 2$. At $L = 128$ that is roughly $10^4$ sweeps per independent sample. Wolff's cluster algorithm removes most of that. His 1989 paper reports no visible critical slowing down on lattices up to $128^2$, with autocorrelation times of 1 to 2 steps per spin for long-range estimators [3]. The Metropolis figure is my own order-of-magnitude recollection of the standard result, not a number from [3], so treat it as a factor of 3 estimate.

So the exam is cheap. A few minutes of a cluster code on a laptop gives statistical errors small enough to test the three results below. I give what a passing code must show.

1. The critical temperature

Measure the Binder cumulant

$$U_L = 1 - \frac{\langle m^4 \rangle}{3\langle m^2\rangle^2}$$

for several $L$. The curves cross near $T_c$. The crossing value is universal, so there is a second test beside the temperature: for the square lattice with periodic boundaries, $U^* = 0.61069$ [4]. That value depends on boundary conditions and system shape [4], so compare only like with like. A square periodic box is the standard case.

My own Lab run gave $T_c = 2.26914 \pm 0.00018\ J/k_B$ against the exact $2.269185$ [1], a deviation of $-0.25\sigma$ (that post). It also gave $\nu = 0.9997 \pm 0.0075$ against the exact 1. I quote those from my earlier work and did not rerun them here. Crossings from small $L$ carry corrections to scaling, so the quoted error depends on the pairs of sizes used. That is a model-form term and it belongs next to the number.

2. The magnetisation exponent

Fit $m(T)$ below $T_c$ in a window $t = 1 - T/T_c$ that you state. The exact curve above lets you do better than fit: compute the exact $m$ at each simulated temperature and compare directly. Finite-size rounding makes the simulated $m$ larger than exact close to $T_c$, so the comparison is clean only for $L \gg \xi$, where $\xi$ is the correlation length. At $t = 0.05$ the exact $m$ is of order $t^{1/8} \approx 0.69$ (hand work, amplitude near 1). A mean-field value would give $t^{1/2} \approx 0.22$. The two predictions differ by a factor of 3 at $t=0.05$, so the exam separates the physics cleanly. The harder part is the fit, because the exponent shifts with the window, and the window is a free parameter.

Where this fails: if your fit window includes $t \gtrsim 0.3$, the pure power law no longer describes the exact curve, and a fitted $\beta$ will differ from 1/8 for a reason that has nothing to do with your code.

3. The specific heat

The specific heat of the infinite lattice diverges logarithmically, not as a power law. Replacing $|t|$ by $1/L$ in the Onsager singularity [1] gives, at $T_c$,

$$\frac{C}{k_B} \approx \frac{2}{\pi}\,(2K_c)^2 \ln L + \text{const}, \quad K_c = J/k_BT_c = 0.4407$$

Hand work: $(2K_c)^2 = 0.7769$, so the slope is $0.495$ per unit of $\ln L$. This is a heuristic scaling argument, and the exact finite-size treatment is the classic Ferdinand and Fisher result, which I have not reread for this essay. Treat 0.495 as my derived estimate, to be checked against the exact slope, not as a quoted value. The point stands: a code whose peak $C$ grows like a power of $L$ has the wrong universality class, and it takes only four lattice sizes to see that a log is not a power. Beware that $\alpha = 0$ here, so a fit of $C \sim L^{\alpha/\nu}$ will return a small number that looks like noise. Fit the log.

Where this fails: with only $L$ from 8 to 32, a log and a small power law differ by less than the typical statistical error. State the range of $L$ and the error next to the slope.

Why a pass means something

Each test targets a different kind of bug.

Test	Catches	Misses
$T_c$ from crossings	wrong Boltzmann weights, wrong cluster probability, bad RNG	errors that cancel in a ratio
$U^*$	wrong boundary handling, non-ergodic updates	slow equilibration at one size only
$\beta$ against exact $m(T)$	magnetisation estimator bugs, too-small lattices	bugs only present at other couplings
slope of $C$ in $\ln L$	wrong energy estimator, wrong universality class	very small errors in energy offsets

The failures are not exotic. A Wolff bond probability of $1 - e^{-2K}$ typed as $1 - e^{-K}$ gives a sharp-looking crossing at the wrong temperature. A magnetisation taken as the mean of $m$ instead of $|m|$ gives zero on a finite lattice below $T_c$. Neither bug looks wrong on a plot with no reference line. Both fail the exam at once. This is why I dislike any figure that has no convergence check or reference curve. The reader cannot tell a good code from one that merely produces a smooth curve.

The strongest objection

The strongest objection is this. Passing the Ising exam proves little about a harder model. The sign problem, frustration, a first-order transition, glassy dynamics and long-range forces are absent in the square-lattice Ising model. A code can pass all three tests and still be wrong on the model you care about. I accept the premise. The exam is necessary, not sufficient.

I made the same structural claim about energy drift in my post on energy conservation: a small drift is necessary for a good trajectory and does not make it right. Exact conservation checks the integrator's setup. It does not check that the equations are the ones you meant. The Ising exam is the same kind of check. It tests the sampler, the estimators, the RNG and the error analysis. It does not test the physics of the new model.

That is not a small thing. In a harder model, no exact answer exists to flag a bug. If the sampler is untested, then every disagreement with a paper or an experiment has two possible causes, the model and the code, and you cannot separate them. After a pass on Onsager's model, one cause is much less likely. A failure on the harder model then points to the model, the size of the lattice, or the equilibration time. I would not call that proof. I would call it a reduction of the search space, and it costs a few minutes of CPU.

A second form of the objection says the benchmark is gamed. People tune a code until it hits $2.269$. That can happen, and the defence is the same as for any test: fix the code before the run, publish the seeds, the sizes and the error method, and let someone else rerun it. A pass with a stated jackknife or bootstrap error and a stated fit window is evidence. A pass quoted as "agrees with Onsager" is not.

What follows if I am right

If the exam is the minimum, three things change. A paper that reports a phase transition from Monte Carlo should state that its code reproduced $T_c$ and $U^*$ on the square-lattice Ising model, with errors and sizes. A referee should treat that line as a reproducibility item like a data-availability statement. And a code that fails should be treated as a bug report, not as a "model dependence".

What would change my mind: a documented case where a code passed the Ising exam at the stated tolerance and still had a sampler bug that the exam could have caught. I know of none, but I have not searched the literature for one, and I put no number on the claim. My own computation for the missing piece, the exact finite-size specific-heat slope against a Wolff run at $L = 16$ to $256$ with a convergence check, is still pending.