Glauber Dynamics: Time Evolution of Ising Spins

Metropolis samples. Glauber evolves.¶

Metropolis Monte Carlo samples the Boltzmann distribution — its “time” is arbitrary, a bookkeeping index. Glauber dynamics picks a specific, physical rate law for how a single spin flips in contact with a heat bath:

w(s_i \to -s_i) = \frac{1}{1 + e^{\beta \Delta E}},

(1)

where $\Delta E = 2 J\, s_i \sum_{\langle i,j\rangle} s_j$ is the energy cost of flipping $s_i$ . This heat-bath rate satisfies detailed balance and has a physical interpretation: each spin instantaneously equilibrates with its local mean field at unit rate. As a consequence $M(t)$ , time correlations, and relaxation times from a Glauber trajectory carry dynamical — not just statistical — content.

The cleanest experiment is a quench. Initialize all spins up and set $T > T_c \approx 2.269\, J/k_B$ , where the equilibrium magnetization is zero. The order parameter then decays from 1 toward 0 approximately exponentially, $M(t) \sim e^{-t/\tau(T)}$ , with $\tau(T)$ growing sharply as $T$ approaches $T_c$ from above — the onset of critical slowing down.

import numpy as np
import matplotlib.pyplot as plt
from numba import njit
from scipy.optimize import curve_fit

rng = np.random.default_rng(0)

A Glauber sweep¶

One sweep = $N^2$ attempted single-spin updates, each accepted with the heat-bath probability above. We store $M$ once per sweep, so the trace has one point per unit of Glauber time.

@njit
def glauber_ising(spins, T=3.0, J=1.0, n_sweeps=400):
    N = spins.shape[0]
    beta = 1.0 / T
    M_trace = np.empty(n_sweeps)
    for sweep in range(n_sweeps):
        for _ in range(N * N):
            i = np.random.randint(N)
            j = np.random.randint(N)
            s = spins[i, j]
            nb = (spins[(i+1) % N, j] + spins[(i-1) % N, j]
                  + spins[i, (j+1) % N] + spins[i, (j-1) % N])
            dE = 2.0 * J * s * nb
            if np.random.rand() < 1.0 / (1.0 + np.exp(beta * dE)):
                spins[i, j] = -s
        M_trace[sweep] = spins.mean()
    return spins, M_trace

Quenching from all-up at $T > T_c$ ¶

We start from a fully magnetized state, $s_i = +1$ everywhere, and let Glauber dynamics bring the system to equilibrium at four temperatures above $T_c$ . Hotter temperatures relax faster.

N = 32
n_sweeps = 400
Ts = [2.5, 3.0, 3.5, 4.0]

traces = {}
for T in Ts:
    spins = np.ones((N, N), dtype=np.int64)
    _, M = glauber_ising(spins, T=T, n_sweeps=n_sweeps)
    traces[T] = M

fig, ax = plt.subplots(figsize=(7, 4))
for T, M in traces.items():
    ax.plot(M, lw=1.2, label=f'$T = {T}$')
ax.axhline(0, color='k', lw=0.5)
ax.set(xlabel='MC sweep', ylabel=r'$M(t)$',
       title=f'Quench from all-up, $L={N}$')
ax.legend()
ax.grid(True, ls=':', alpha=0.5)
fig.tight_layout()

Extracting $\tau(T)$ ¶

Fit the early portion of each trajectory — while $M$ is still well above zero — to $M(t) = e^{-t/\tau}$ . The fitted $\tau$ grows rapidly as $T$ drops toward $T_c$ , the fingerprint of critical slowing.

def exp_decay(t, tau):
    return np.exp(-t / tau)

taus = {}
fig, ax = plt.subplots(figsize=(7, 4))
for T, M in traces.items():
    mask = M > 0.05
    t = np.arange(len(M))[mask]
    popt, _ = curve_fit(exp_decay, t, M[mask], p0=[20.0], maxfev=10000)
    taus[T] = popt[0]
    ax.plot(t, M[mask], lw=1.2, label=f'$T={T}$, $\\tau={popt[0]:.1f}$')
    ax.plot(t, exp_decay(t, *popt), 'k--', lw=0.8)
ax.set(xlabel='MC sweep', ylabel=r'$M(t)$', yscale='log',
       title='Exponential fits to the relaxation')
ax.legend()
ax.grid(True, ls=':', alpha=0.5, which='both')
fig.tight_layout()

print({T: f'{tau:.2f}' for T, tau in taus.items()})

Equilibrium autocorrelation¶

Once the system reaches equilibrium the magnetization fluctuates around zero. Its autocorrelation $C(t) = \langle M(0) M(t)\rangle / \langle M^2\rangle$ should decay on the same time scale $\tau$ — this is the equilibrium counterpart of the quench experiment, and is what the fluctuation–dissipation relation connects to linear response.

T = 3.0
spins = np.ones((N, N), dtype=np.int64)
spins, _ = glauber_ising(spins, T=T, n_sweeps=500)          # burn-in
_, M_eq = glauber_ising(spins, T=T, n_sweeps=4000)          # production

def autocorr(x, maxlag=200):
    x = np.asarray(x) - np.mean(x)
    var = np.dot(x, x) / len(x)
    lags = np.arange(maxlag)
    C = np.array([np.mean(x[:len(x)-t] * x[t:]) for t in lags]) / var
    return lags, C

lags, C = autocorr(M_eq, maxlag=200)

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(lags, C, color='crimson')
ax.axhline(0, color='k', lw=0.5)
ax.set(xlabel='lag (sweeps)', ylabel=r'$C(t)/C(0)$',
       title=f'Equilibrium magnetization autocorrelation at $T={T}$')
ax.grid(True, ls=':', alpha=0.5)
fig.tight_layout()

Problems¶

Critical slowing down. Measure $\tau(T)$ at $T = 2.5,\ 2.4,\ 2.35,\ 2.30$ (approaching $T_c \approx 2.269$ from above). Plot $\tau$ vs. $T - T_c$ on log axes. Does a power law $\tau \sim (T-T_c)^{-\nu z}$ describe the data? For the 2D Ising universality class $\nu = 1$ and $z \approx 2.17$ .
Glauber vs. Metropolis. Implement Metropolis spin-flip dynamics (accept with $\min(1, e^{-\beta \Delta E})$ ) and re-run the quench at $T = 3.0$ . Plot both $M(t)$ curves. Both must approach the same equilibrium — what differs is the time scale of relaxation.
Response to a small field. Equilibrate the system at $T = 3.0$ with $h = 0$ , then at $t = 0$ switch on a weak field $h = 0.1$ (add $-h s_i$ to the single-spin energy, i.e. add $2 h s_i$ to dE). Fit $M(t) \sim M_\infty (1 - e^{-t/\tau_h})$ and compare $\tau_h$ to the equilibrium autocorrelation time at the same $T$ — this is the fluctuation–dissipation relation at work.

Reference¶

R. J. Glauber, Time-Dependent Statistics of the Ising Model, J. Math. Phys. 4, 294 (1963). Glauber (1963)