The equilibrium average and its empty territory

Picture 224 rigid disks in a square. A configuration is not one photograph but a point in a space with hundreds of coordinates: an x and y position for every disk, constrained so that no two overlap. Even if every coordinate were allowed only one thousand values, the catalogue would have far more entries than any machine could enumerate. Thermodynamics nevertheless asks for averages over that space. Pressure, energy and other equilibrium quantities depend not on the most attractive arrangement alone, but on all arrangements weighted by how probable they are at a chosen temperature.

For a configuration x with potential energy E(x), the canonical equilibrium weight is exp(−E(x)/kT). The normalising constant, conventionally called Z, is the integral of that weight over all allowed configurations. An observable A therefore has equilibrium mean equal to the integral of A(x)exp(−E(x)/kT) divided by Z. The formula is compact; its domain is ferocious. Direct quadrature needs a grid in hundreds of dimensions. Uniform random sampling is legal in principle but profligate in practice, because most randomly generated arrangements carry negligible thermal weight or violate the hard-core constraint.

The pre-1953 Monte Carlo idea was already clear: replace an integral by an average over random points. Metropolis and Ulam had described that numerical outlook in 1949. The new paper changes where the random points come from. Instead of drawing configurations uniformly and attaching wildly unequal weights afterwards, it constructs a dependent sequence whose own visiting frequencies carry the Boltzmann weight. Once that sequence has settled, an unweighted average of A along the chain estimates the canonical expectation. The intractable Z cancels because the chain needs only ratios of weights.

This is importance sampling joined to a Markov process. “Importance” means that computational effort is concentrated where the integral receives its important contribution. “Markov” means that the next configuration is generated from the present one, not independently from the full space. Local dependence makes proposals cheap, but it also creates correlation and the possibility of slow travel. The 1953 paper is strongest on the sampling construction and its physical demonstration. It does not contain the modern diagnostic machinery that later disciplines developed for judging whether a finite chain has travelled far enough.

A proposal, an energy difference and a gate

The 1953 construction begins from one configuration. Select a particle, draw independent random displacements for its coordinates within a prescribed range, and calculate how the system’s energy would change if that particle moved. Because the potential is a sum of pair interactions, terms unrelated to the selected particle cancel. Local proposals are therefore computationally economical. In the hard-disk example, a proposed overlap has effectively infinite positive energy and must be rejected; any non-overlapping move has no energy penalty under the rigid-core potential and is accepted.

For a general potential, define ΔE = E(new) − E(old). If ΔE is negative or zero, accept the new configuration. If ΔE is positive, generate a uniform random number u between zero and one and accept only when u is smaller than exp(−ΔE/kT). A large uphill move is possible but rare. A tiny uphill move is commonly accepted. This asymmetric treatment of energy change is precisely calibrated: for a symmetric proposal mechanism, the ratio of forward and reverse accepted transitions matches the ratio of Boltzmann weights.

Three pieces must remain distinct. The proposal chooses what change to consider. The acceptance gate decides whether that proposed change becomes the next state. The estimator records the state that actually results. Alter one piece and the stationary frequencies may change. In particular, the simple exponential gate presumes proposal symmetry: proposing y from x is as likely as proposing x from y. The authors’ random displacement has that property when the displacement distribution is centred and used identically in both directions.

The chain is artificial. A proposed particle jump is not integrated from Newton’s laws, and one cycle of attempted moves is not a measured interval of molecular time. The sequence is designed to integrate an equilibrium distribution. It may reproduce static equilibrium expectations even though its path has no faithful kinetic meaning. Confusing sampling steps with physical seconds would turn a numerical device into an unsupported dynamical model.

Why a rejected move still counts

The algorithm’s most easily omitted instruction appears mundane: when the candidate is rejected, the previous configuration is retained for the next step and included again in the average. Suppose a low-energy configuration is surrounded by much higher-energy neighbours. Most outward proposals fail. Repeating the current state for each failure makes the chain dwell there. That duration is how the sequence represents high equilibrium probability. If an analyst logs only accepted moves, the dwelling time disappears and the resulting catalogue answers a different question.

The error is easiest to see with two states. Let L have weight 1 and H have weight e⁻². Suppose the proposal always offers the other state. Every H→L proposal is accepted, but only e⁻² of L→H proposals pass. The complete chain spends long runs at L: a rejected uphill attempt produces L again. If only accepted changes are kept, the log alternates L,H,L,H and falsely reports equal frequencies. The rejection is not missing data. It is the event that gives L its longer residence.

Self-transitions also help break periodicity. A chain that must change state on every step can fall into a rigid rhythm, revisiting classes of states on alternating times. Rejection gives it a positive chance of staying put. Modern Markov-chain theory separates this property from irreducibility, which asks whether the chain can eventually reach every relevant region. The 1953 authors discuss the need to reach configurations through successive displacements, but the mature vocabulary and general convergence theorems came later.

Pairwise flow and the stationary Boltzmann distribution

Let π(x) be proportional to e⁻ᵝᴱ⁽ˣ⁾, and let q(x,y) be the probability of proposing y from x. The local displacement rule is symmetric, so q(x,y)=q(y,x). Assume y has higher energy. The accepted flow from x to y is π(x)q(x,y)e⁻ᵝ⁽ᴱ⁽ʸ⁾⁻ᴱ⁽ˣ⁾⁾. Multiplying the exponentials reduces this to a constant times q(x,y)e⁻ᵝᴱ⁽ʸ⁾. The reverse move is downhill and accepted with probability one, giving π(y)q(y,x), the same quantity. Probability flow balances between the pair.

When this equality holds for every pair, the total flow into each state equals its total flow out. The target distribution is stationary: if the chain is distributed according to π before a step, it remains so afterwards. Contemporary texts call the equality detailed balance or reversibility. The 1953 paper establishes the essential transition-frequency argument without presenting today’s general measure-theoretic framework. Source fidelity therefore requires two sentences, not one: the pairwise proof is present; the standard label belongs to the later vocabulary used to organise it.

Stationarity is a property of a transition mechanism and a distribution. It is not evidence that an arbitrary starting state has already reached that distribution after twenty steps. To obtain convergence from a chosen start, the chain also needs access to the relevant configuration space and freedom from a deterministic cycle. Even then, the rate may be painfully slow. A theorem that identifies π as stationary does not specify how many cycles the MANIAC, a laptop or a supercomputer must run before estimates are dependable.

Two hundred and twenty-four disks on the MANIAC

The primary paper does not introduce its algorithm as abstract statistics. It asks for the equation of state of interacting molecules and then reports a preliminary two-dimensional rigid-disk calculation. The authors place N particles in a square and imagine that square repeated throughout the plane. A particle interacts through the minimum distance to another particle or one of its periodic images. This periodic construction reduces the artificial surface that would dominate a small finite box.

The reported system contains 224 disks, initially arranged on a 14 by 16 triangular lattice. A cycle consists of attempting one displacement for each particle. The step range is tuned so that moves are neither almost always rejected nor so timid that the configuration changes imperceptibly. On the Los Alamos MANIAC, one cycle took about three minutes. The resulting runs were tiny by present standards, which makes the paper’s agreement with analytic comparisons historically striking but also makes its cautious language important.

For hard disks, potential energy has a stark form: zero when disks do not overlap and infinite when they do. The Boltzmann gate therefore becomes a geometric gate. A non-overlapping trial is accepted; an overlapping trial is rejected and the old configuration is counted again. Pressure is inferred from how frequently small virtual changes in scale would create overlaps, connecting sampled geometry to the equation of state. The article compares the numerical results with a four-term virial expansion at lower densities and with a free-volume expression at high density.

The calculation does not report a Lennard-Jones fluid. The introduction says work with a two-dimensional Lennard-Jones potential was in progress and that three-dimensional work was being investigated. Later histories often compress these developments into a single origin story. The governing source does not permit that compression. Canon 0014 is a general acceptance method demonstrated on rigid disks, with later potentials and dimensions left for subsequent work.

From sampled geometry to pressure

Sampling configurations is only the first half of an equation-of-state calculation. The second half converts those configurations into a thermodynamic observable. For smooth pair potentials, pressure can be related to a virial involving interparticle forces and separations. Hard disks are less accommodating: their energy changes discontinuously at contact, so an ordinary force is zero away from collision and singular at collision. The authors instead ask how the accessible configuration volume changes under a minute compression. A compressed copy is forbidden when the shrinkage makes a pair overlap. The frequency of near-contact pairs therefore reveals the derivative needed for pressure.

In practice they accumulate a distribution of the smallest fractional compression that would cause an overlap. Near zero compression, the slope of that distribution supplies the contact information. This is a second importance-sampling lesson. The chain supplies equilibrium configurations; a carefully chosen geometric statistic extracts a derivative that cannot be read from the hard-core energy values themselves. The article’s headline algorithm and its pressure estimator solve different problems, and both are needed to turn random configurations into an equation-of-state point.

The low-density virial expansion and the high-density free-volume expression offer external structure against which the numerical points can be compared. A virial series organises pressure as powers of density and works best where particles interact relatively rarely. A free-volume argument imagines each disk confined by neighbours and becomes more plausible in a crowded, ordered regime. Agreement or disagreement with either approximation is scientifically informative, but neither comparison alone certifies sampling convergence. An approximate theory can be wrong where a chain is right, and a short chain can accidentally resemble an approximate curve.

Random numbers were part of the apparatus

The machine also needed repeatable pseudorandom numbers. Early computers did not provide a trusted random service, so generation had to be programmed and checked as part of the experiment. The 1953 implementation used a middle-square style generator and monitored it for short repetitions. Such a generator would not satisfy modern statistical standards, but its presence makes a general methodological point: a stochastic algorithm inherits the defects of its random source. A perfect acceptance equation fed by a cycling or correlated number stream may visit configurations with unintended regularity.

This apparatus-level view keeps the paper from becoming a disembodied formula. The calculation joins a representation of a physical material, a periodic geometry, a proposal scale, an energy or overlap test, a random-number generator, a state ledger, an observable estimator and a finite computer budget. The algorithm’s validity argument concerns the idealised transition probabilities. The credibility of a numerical result also depends on whether the implemented apparatus approximates those probabilities closely enough.

Initial relaxation and the choice of what to average

The initial triangular lattice is convenient, not drawn from the equilibrium distribution at every density. Early cycles allow the artificial order and chosen coordinates to relax. Modern readers call discarded initial iterations warm-up or burn-in, but the underlying decision is older: do not let an arbitrary starting construction dominate the reported average. Discarding a fixed prefix cannot guarantee equilibrium, especially near a slow transition. It merely limits one visible source of bias.

After relaxation, every resulting configuration contributes, including duplicates after rejection. The resulting average has two uncertainties. Thermal variability reflects genuine variation of the observable under the target distribution. Monte Carlo error reflects having observed only a correlated finite chain. Running longer reduces the second under suitable conditions but not the first, which belongs to the physical ensemble. Confusing them can produce either exaggerated confidence in a simulation or a mistaken demand that equilibrium fluctuations disappear.

Correlation, tuning and the gap between stationary and sampled

Every state depends on the previous state, so successive observations are correlated. A chain that records nearly the same configuration a hundred times has one hundred rows but less information than one hundred independent configurations. Modern practice summarises the loss through autocorrelation and effective sample size. These are later tools, useful for reading the 1953 mechanism but not part of its reported apparatus.

Proposal scale creates a trade-off. If displacements are huge in a dense system, almost every candidate overlaps another disk and the chain repeats itself. If displacements are microscopic, acceptance is high but positions diffuse through configuration space glacially. An intermediate scale can move efficiently. No universal acceptance percentage follows from the primary paper, and modern rules of thumb depend on dimension, proposal family and target geometry.

A harder failure occurs when the target has separated basins. Local moves may explore one basin beautifully while almost never crossing the barrier to another. Trace plots inside the visited basin can look stable. Repeated runs from dispersed starts, tempering, collective moves or problem-specific proposals can expose or reduce the problem, but no finite diagnostic proves that an undiscovered region does not exist. The acceptance formula preserves a target distribution; it does not guarantee convenient travel through every target.

Executable Metropolis laboratory

The bench below uses seven discrete positions with fixed energies. A proposal moves one position left or right with wraparound, chosen symmetrically. The exact target is therefore proportional to e⁻ᵝᴱ. Each press of Step produces exactly one recorded state. When a proposal fails, the trajectory advances horizontally while its vertical position remains unchanged. Run adds the selected number of steps using a deterministic seeded generator, making results reproducible for inspection and testing.

The landscape is deliberately small enough that exact probabilities can be calculated. Finite runs will fluctuate around them, especially at high β, where the two deep wells communicate through energetic barriers. Increase the run length and compare seeds. A correct implementation does not force every run to look exact; it makes long-run discrepancies shrink in a probabilistic way while preserving the target as stationary.

When the proposal is not symmetric

Hastings’s 1970 paper generalises the 1953 method. If q(x,y) and q(y,x) differ, the acceptance gate must account for both the target ratio and the proposal ratio. A common form accepts with probability min{1, π(y)q(y,x)/[π(x)q(x,y)]}. The correction prevents a convenient but directionally biased proposal from tilting the stationary distribution. When q is symmetric, the q terms cancel and the original Metropolis rule returns.

A five-name paper and a complicated history

The conventional name “Metropolis algorithm” is too established to avoid, but it should not erase the byline. Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller and Edward Teller all appear on the 1953 publication. Later recollections complicate what that order implies. Marshall Rosenbluth’s 2003 conference account credits Edward Teller with the initial statistical-mechanics suggestion, describes Marshall as developing the algorithm, and identifies Arianna as writing the MANIAC code. A later historical reconstruction by J. E. Gubernatis records those recollections and the wider Los Alamos setting.

Arianna Rosenbluth’s role matters especially because programming an unprecedented stochastic algorithm on an early computer was intellectual experimental work, not transcription. The American Physical Society’s later history reports that she programmed the calculation and ran the machine during difficult overnight sessions. This evidence corrects an omission in many compressed narratives. It remains retrospective evidence, gathered decades after the article, and should be described with that status rather than inserted into the primary text as if it were a contemporaneous authorship statement.

Hastings’s name also marks a genuine mathematical extension, not merely a renaming. The symmetric proposal assumption is hidden by the elegance of the original rule. Once proposals become state-dependent or directional, its cancellation no longer occurs. Hastings supplies a family of acceptance constructions and discusses estimation from Markov-chain samples. The combined name “Metropolis–Hastings” appropriately signals a lineage while still leaving room to name all five 1953 authors when teaching the origin.

Memory studio

A durable memory of the paper needs more than “accept downhill, sometimes accept uphill”. The four devices below preserve the governing integral, the self-loop, the symmetry condition and the separation between equilibrium sampling and physical dynamics.

What the paper establishes, and what it leaves open

The paper establishes a constructive way to make configuration visits carry canonical weight, verifies the transition logic for a symmetric local proposal and demonstrates the method on an early electronic computer. It also makes two engineering choices that became enduring: local energy updates and periodic boundary conditions. These are not decorative details around the acceptance formula. They make the numerical experiment feasible and connect a finite machine state to a bulk physical model.

It leaves open how best to diagnose convergence, how to quantify Monte Carlo uncertainty under strong correlation, how to design proposals for separated modes and how to scale the method across modern high-dimensional models. Later work developed answers, but none removes the need to inspect a sampler as a coupled system. A correct target, a lawful acceptance gate and a runnable program can coexist with an unusably slow chain.

For the wider canon, the enduring move is architectural: transform an impossible weighted sum into a process whose ordinary occupancy embodies the weights. That move recurs in probabilistic inference, statistical physics and stochastic computation. It is relevant to how intelligent systems may approximate uncertain expectations. It does not follow that a system which samples a distribution understands the samples, has goals, or is conscious. Computational mechanism and phenomenal interpretation remain separate questions.