MCMC in Physics

Imagine a magnet made of a trillion spinning electrons, each able to point up or down. The number of possible configurations is $2^{10^{12}}$ — a number so large it dwarfs the count of atoms in the observable universe. Yet physicists routinely compute the magnet's average energy, heat capacity and magnetic moment to exquisite precision. How?

The answer is Markov Chain Monte Carlo (MCMC). Instead of visiting every configuration — an impossibility — MCMC constructs a random walk that, over time, visits each configuration with a probability proportional to its Boltzmann weight $e^{-E/kT}$. Once the walk has run long enough, simple averages over the visited states give exact thermodynamic quantities.

The key insight, discovered by Nicholas Metropolis and colleagues in 1953, is that you do not need to know the partition function Z (the sum over all states). A simple accept/reject rule based on energy differences is enough to guarantee the walk converges to the correct distribution. That idea — local decisions, global correctness — is one of the most influential algorithms ever devised.

Below is a 16×16 Ising model — a grid of spins, each +1 (white) or −1 (black). The Metropolis algorithm proposes flipping one random spin per step and accepts the flip with probability $\min(1,\, e^{-\Delta E/T})$. At low temperature the system orders into domains; at high temperature it disorders.

<div class="controls">
  <label>Temperature: <span id="tval">2.5</span>
    <input type="range" id="temp" min="0.5" max="5.0" step="0.05" value="2.5">
  </label>
  <div class="btns">
    <button id="runBtn" type="button">▶ Run</button>
    <button id="stepBtn" type="button">Step ×100</button>
    <button id="resetBtn" type="button" class="ghost">Reset</button>
  </div>
</div>
<canvas id="canvas" width="320" height="320"></canvas>
<div class="stats">
  <span>Steps: <b id="stepCount">0</b></span>
  <span>|M|: <b id="mag">—</b></span>
  <span>⟨E⟩/spin: <b id="eng">—</b></span>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; }
.controls { display: flex; flex-wrap: wrap; gap: .6rem; align-items: center; margin-bottom: .5rem; }
label { font-size: .85rem; display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 130px; accent-color: #1d3557; }
.btns { display: flex; gap: .4rem; }
button { font: 600 13px system-ui; padding: .35rem .75rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
canvas { display: block; border: 1px solid #cdd9e3; border-radius: 6px; image-rendering: pixelated; max-width: 100%; }
.stats { display: flex; gap: 1rem; font-size: .82rem; margin-top: .4rem; color: #444; }
.stats b { color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Metropolis algorithm on a 2-D Ising model (periodic boundaries).
// Spin +1 (white) or -1 (black). Hamiltonian H = -J * sum_{<i,j>} s_i * s_j (J=1).
// Accept flip with probability min(1, exp(-dE / T)).
// The Boltzmann distribution is the unique stationary distribution (detailed balance).
const N = 32;         // grid size NxN
const CELL = 10;      // pixels per spin
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
canvas.width = N * CELL;
canvas.height = N * CELL;

let spin = new Int8Array(N * N); // +1 or -1
let steps = 0;
let running = false;
let rafId = null;

function idx(r, c) { return ((r + N) % N) * N + ((c + N) % N); }

function initRandom() {
  for (let i = 0; i < N * N; i++) spin[i] = Math.random() < 0.5 ? 1 : -1;
  steps = 0;
}

function siteEnergy(r, c) {
  const s = spin[idx(r, c)];
  const nb = spin[idx(r-1,c)] + spin[idx(r+1,c)] + spin[idx(r,c-1)] + spin[idx(r,c+1)];
  return -s * nb; // J=1, divide by 2 later for total energy
}

function metropolisStep(T) {
  // N*N attempted flips = 1 Monte Carlo sweep
  for (let k = 0; k < N * N; k++) {
    const i = (Math.random() * N * N) | 0;
    const r = (i / N) | 0;
    const c = i % N;
    const dE = -2 * siteEnergy(r, c); // energy change if we flip
    if (dE <= 0 || Math.random() < Math.exp(-dE / T)) {
      spin[i] = -spin[i];
    }
  }
  steps += N * N;
}

function draw() {
  const img = ctx.createImageData(N * CELL, N * CELL);
  for (let r = 0; r < N; r++) {
    for (let c = 0; c < N; c++) {
      const v = spin[idx(r, c)] === 1 ? 240 : 30;
      for (let dr = 0; dr < CELL; dr++) {
        for (let dc = 0; dc < CELL; dc++) {
          const p = ((r * CELL + dr) * N * CELL + c * CELL + dc) * 4;
          img.data[p] = v; img.data[p+1] = v; img.data[p+2] = v; img.data[p+3] = 255;
        }
      }
    }
  }
  ctx.putImageData(img, 0, 0);
}

function updateStats() {
  let m = 0, e = 0;
  for (let r = 0; r < N; r++) {
    for (let c = 0; c < N; c++) {
      m += spin[idx(r, c)];
      const nb = spin[idx(r-1,c)] + spin[idx(r+1,c)] + spin[idx(r,c-1)] + spin[idx(r,c+1)];
      e += -spin[idx(r,c)] * nb;
    }
  }
  e /= 2; // each bond counted twice
  const nn = N * N;
  document.getElementById('mag').textContent = Math.abs(m / nn).toFixed(3);
  document.getElementById('eng').textContent = (e / nn).toFixed(3);
  document.getElementById('stepCount').textContent = fmt('%v% sweeps', { v: (steps / (N * N)).toFixed(0) });
}

function getT() { return parseFloat(document.getElementById('temp').value); }

function loop() {
  if (!running) return;
  metropolisStep(getT());
  draw();
  updateStats();
  rafId = requestAnimationFrame(loop);
}

document.getElementById('runBtn').addEventListener('click', () => {
  running = !running;
  document.getElementById('runBtn').textContent = running ? '⏸ Pause' : '▶ Run';
  if (running) loop();
});

document.getElementById('stepBtn').addEventListener('click', () => {
  const T = getT();
  for (let i = 0; i < 100; i++) metropolisStep(T);
  draw(); updateStats();
});

document.getElementById('resetBtn').addEventListener('click', () => {
  running = false;
  if (rafId) cancelAnimationFrame(rafId);
  document.getElementById('runBtn').textContent = '▶ Run';
  initRandom(); draw(); updateStats();
});

document.getElementById('temp').addEventListener('input', () => {
  document.getElementById('tval').textContent = getT().toFixed(2);
});

// Boot
initRandom();
draw();
updateStats();

<div class="controls">
  <label>Temperature: <span id="tval">2.5</span>
    <input type="range" id="temp" min="0.5" max="5.0" step="0.05" value="2.5">
  </label>
  <div class="btns">
    <button id="runBtn" type="button">▶ Run</button>
    <button id="stepBtn" type="button">Step ×100</button>
    <button id="resetBtn" type="button" class="ghost">Reset</button>
  </div>
</div>
<canvas id="canvas" width="320" height="320"></canvas>
<div class="stats">
  <span>Steps: <b id="stepCount">0</b></span>
  <span>|M|: <b id="mag">—</b></span>
  <span>⟨E⟩/spin: <b id="eng">—</b></span>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; }
.controls { display: flex; flex-wrap: wrap; gap: .6rem; align-items: center; margin-bottom: .5rem; }
label { font-size: .85rem; display: flex; align-items: center; gap: .4rem; }
input[type=range] { width: 130px; accent-color: #1d3557; }
.btns { display: flex; gap: .4rem; }
button { font: 600 13px system-ui; padding: .35rem .75rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
canvas { display: block; border: 1px solid #cdd9e3; border-radius: 6px; image-rendering: pixelated; max-width: 100%; }
.stats { display: flex; gap: 1rem; font-size: .82rem; margin-top: .4rem; color: #444; }
.stats b { color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Metropolis algorithm on a 2-D Ising model (periodic boundaries).
// Spin +1 (white) or -1 (black). Hamiltonian H = -J * sum_{<i,j>} s_i * s_j (J=1).
// Accept flip with probability min(1, exp(-dE / T)).
// The Boltzmann distribution is the unique stationary distribution (detailed balance).
const N = 32;         // grid size NxN
const CELL = 10;      // pixels per spin
const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
canvas.width = N * CELL;
canvas.height = N * CELL;

let spin = new Int8Array(N * N); // +1 or -1
let steps = 0;
let running = false;
let rafId = null;

function idx(r, c) { return ((r + N) % N) * N + ((c + N) % N); }

function initRandom() {
  for (let i = 0; i < N * N; i++) spin[i] = Math.random() < 0.5 ? 1 : -1;
  steps = 0;
}

function siteEnergy(r, c) {
  const s = spin[idx(r, c)];
  const nb = spin[idx(r-1,c)] + spin[idx(r+1,c)] + spin[idx(r,c-1)] + spin[idx(r,c+1)];
  return -s * nb; // J=1, divide by 2 later for total energy
}

function metropolisStep(T) {
  // N*N attempted flips = 1 Monte Carlo sweep
  for (let k = 0; k < N * N; k++) {
    const i = (Math.random() * N * N) | 0;
    const r = (i / N) | 0;
    const c = i % N;
    const dE = -2 * siteEnergy(r, c); // energy change if we flip
    if (dE <= 0 || Math.random() < Math.exp(-dE / T)) {
      spin[i] = -spin[i];
    }
  }
  steps += N * N;
}

function draw() {
  const img = ctx.createImageData(N * CELL, N * CELL);
  for (let r = 0; r < N; r++) {
    for (let c = 0; c < N; c++) {
      const v = spin[idx(r, c)] === 1 ? 240 : 30;
      for (let dr = 0; dr < CELL; dr++) {
        for (let dc = 0; dc < CELL; dc++) {
          const p = ((r * CELL + dr) * N * CELL + c * CELL + dc) * 4;
          img.data[p] = v; img.data[p+1] = v; img.data[p+2] = v; img.data[p+3] = 255;
        }
      }
    }
  }
  ctx.putImageData(img, 0, 0);
}

function updateStats() {
  let m = 0, e = 0;
  for (let r = 0; r < N; r++) {
    for (let c = 0; c < N; c++) {
      m += spin[idx(r, c)];
      const nb = spin[idx(r-1,c)] + spin[idx(r+1,c)] + spin[idx(r,c-1)] + spin[idx(r,c+1)];
      e += -spin[idx(r,c)] * nb;
    }
  }
  e /= 2; // each bond counted twice
  const nn = N * N;
  document.getElementById('mag').textContent = Math.abs(m / nn).toFixed(3);
  document.getElementById('eng').textContent = (e / nn).toFixed(3);
  document.getElementById('stepCount').textContent = fmt('%v% sweeps', { v: (steps / (N * N)).toFixed(0) });
}

function getT() { return parseFloat(document.getElementById('temp').value); }

function loop() {
  if (!running) return;
  metropolisStep(getT());
  draw();
  updateStats();
  rafId = requestAnimationFrame(loop);
}

document.getElementById('runBtn').addEventListener('click', () => {
  running = !running;
  document.getElementById('runBtn').textContent = running ? '⏸ Pause' : '▶ Run';
  if (running) loop();
});

document.getElementById('stepBtn').addEventListener('click', () => {
  const T = getT();
  for (let i = 0; i < 100; i++) metropolisStep(T);
  draw(); updateStats();
});

document.getElementById('resetBtn').addEventListener('click', () => {
  running = false;
  if (rafId) cancelAnimationFrame(rafId);
  document.getElementById('runBtn').textContent = '▶ Run';
  initRandom(); draw(); updateStats();
});

document.getElementById('temp').addEventListener('input', () => {
  document.getElementById('tval').textContent = getT().toFixed(2);
});

// Boot
initRandom();
draw();
updateStats();

Try dragging the Temperature slider. Below the critical temperature (~2.27 in these units) large ordered domains form — the system "magnetizes." Above it, thermal noise wins and spins randomize. Watch the Average magnetization track the phase transition in real time.

MCMC is provably correct — in the limit of infinite steps. The hard question is how many steps are needed to get close to the true distribution, a quantity called the mixing time.

• Detailed balance: the Metropolis rule satisfies $\pi(x)P(x\to y) = \pi(y)P(y\to x)$, guaranteeing the Boltzmann distribution is a fixed point of the walk. This is the proof of correctness.
• Ergodicity: the walk must be able to reach any configuration from any other, or it will be stuck in a local pocket of state space and never converge.
• Mixing time can be exponential. For systems with a rugged energy landscape — many local minima separated by high barriers — the chain may take exponentially many steps to escape a valley. Simulated annealing (slowly cooling T) is one heuristic fix.
• The phase transition itself is the enemy. Near the critical temperature, correlation lengths diverge and the mixing time grows as a power law — a phenomenon called critical slowing down. Specialized algorithms (Wolff cluster updates, parallel tempering) were invented precisely to fight it.
• #P-hardness lurks beneath: exactly computing the partition function Z is #P-hard, which is why MCMC approximation is so valuable — and why proving tight mixing-time bounds is an active research frontier.

The Metropolis idea escaped physics decades ago and now powers an enormous slice of modern science:

• Statistical mechanics: Ising models, spin glasses, liquid simulations, and phase-transition studies all rely on MCMC. It is the workhorse of condensed-matter and materials physics.
• Lattice QCD: the strong nuclear force is simulated on a discrete spacetime lattice. MCMC is the only practical way to compute quark masses and hadron spectra from first principles.
• Bayesian inference: posterior distributions over model parameters are high-dimensional Boltzmann-like distributions. MCMC (especially Hamiltonian Monte Carlo) is the standard tool in probabilistic programming and data analysis.
• Protein folding and drug design: conformational energy landscapes are explored via MCMC to predict how molecules fold and bind — directly connected to the protein folding problem.
• Gravitational-wave astronomy: LIGO uses MCMC to infer source parameters (masses, spins, sky location) from noisy detector output.
• Machine learning: training energy-based models (Boltzmann machines, diffusion models' score estimators) draws on the same ideas.

Anywhere there is a probability distribution too complex to sum analytically, MCMC turns the problem into a random walk.

The Metropolis algorithm is a beautiful piece of engineering: by never needing to compute the total weight of all configurations, it turns an exponential enumeration into a walk you can afford to run. The price is that correctness comes only in the limit, and the mixing time — how long you actually need to run — can itself be exponential for hard landscapes.

That tension between the algorithm's elegance and the harshness of rugged energy surfaces is the engine of an entire research programme. Wolff updates, parallel tempering, Hamiltonian Monte Carlo — each is a new answer to the same question: how do you make the walk mix fast enough to be useful?

In the end, MCMC is a lesson in the power of local rules to produce global order — the same theme behind Bayesian inference, neural network training, and the statistical physics it was born to serve.