Metastability and Nucleation in the 2D Ising Model

Metastability in the Ising model¶

Start from all spins up at $T < T_c$ and switch on a small negative field $h < 0$ . The fully-magnetized +1 state is no longer the thermodynamic ground state — the stable phase is now -1 — but it looks stable for a long time. Flipping the whole lattice requires creating a droplet of the new phase, and a droplet of radius $R$ costs surface free energy $\sim 2\pi R \sigma$ while gaining bulk free energy $\sim -2|h| \pi R^2 / (\text{unit cell})$ . The droplet only grows spontaneously past the critical radius

R_c \sim \frac{\sigma}{|h|}, \qquad \Delta F_c \sim \frac{\pi \sigma^2}{|h|}.

(1)

Waiting for a thermal fluctuation to produce such a droplet gives a mean first-passage time

\langle \tau \rangle \sim \exp\!\left(\frac{\Delta F_c}{k_B T}\right),

(2)

which diverges as $|h| \to 0$ — the hallmark of a rare event driven by a free-energy barrier.

import numpy as np
import matplotlib.pyplot as plt
from numba import njit

rng = np.random.default_rng(0)

Metropolis dynamics with a field¶

We add $-h \sum_i s_i$ to the Hamiltonian. Flipping $s_i$ changes the energy by $\Delta E = 2 J s_i \sum_{\langle j \rangle} s_j + 2 h s_i$ .

@njit
def metropolis_h_trace(spins, T=2.0, h=-0.08, J=1.0, n_sweeps=6000):
    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 + 2.0 * h * s
            if dE <= 0 or np.random.rand() < np.exp(-beta * dE):
                spins[i, j] = -s
        M_trace[sweep] = spins.mean()
    return spins, M_trace


@njit
def first_passage(spins, T=2.0, h=-0.08, J=1.0, max_sweeps=40000):
    N = spins.shape[0]
    beta = 1.0 / T
    for sweep in range(max_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 + 2.0 * h * s
            if dE <= 0 or np.random.rand() < np.exp(-beta * dE):
                spins[i, j] = -s
        if spins.sum() < 0:
            return sweep + 1
    return -1

A single nucleation event¶

A long plateau near $M = +1$ is broken by a sudden drop — that is a critical droplet getting lucky, growing past $R_c$ , and sweeping through the lattice.

N = 32
spins = np.ones((N, N), dtype=np.int64)
spins, M = metropolis_h_trace(spins, T=2.0, h=-0.10, n_sweeps=3000)

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.plot(M, lw=1, color='steelblue')
ax.axhline(0, color='k', lw=0.5)
ax.set(xlabel='MC sweep', ylabel=r'$M(t)$',
       title=f'Nucleation trajectory ($L={N}$, $T=2.0$, $h=-0.10$)')
ax.grid(True, ls=':', alpha=0.5)
fig.tight_layout()

A distribution of nucleation times¶

Nucleation is stochastic, so a single run is not representative. Running many independent realizations gives the first-passage-time distribution — it should be approximately exponential (Poisson events), with mean $\langle \tau \rangle$ .

n_runs = 60
taus = []
for k in range(n_runs):
    spins = np.ones((N, N), dtype=np.int64)
    tau = first_passage(spins, T=2.0, h=-0.10, max_sweeps=20000)
    if tau > 0:
        taus.append(tau)

taus = np.array(taus)
fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(taus, bins=20, color='steelblue', edgecolor='k', alpha=0.7)
ax.axvline(taus.mean(), color='crimson', lw=2, label=fr'$\langle\tau\rangle={taus.mean():.0f}$')
ax.set(xlabel='first-passage time $\\tau$ (sweeps)',
       ylabel='count',
       title=f'Nucleation time distribution ({len(taus)}/{n_runs} runs)')
ax.legend()
ax.grid(True, ls=':', alpha=0.5)
fig.tight_layout()

Problems¶

Field dependence of $\langle \tau \rangle$ . Measure $\langle \tau \rangle$ at $T = 2.0$ for $h \in \{-0.08, -0.10, -0.12, -0.15\}$ using $\sim 30$ runs each. Plot $\ln \langle \tau \rangle$ vs. $1/|h|$ . CNT predicts a linear relation with slope $\pi \sigma^2/(k_B T)$ — extract the effective surface tension $\sigma$ .
Temperature dependence. Fix $h = -0.10$ and measure $\langle \tau \rangle$ at $T = 1.8, 2.0, 2.15$ . Higher $T$ helps cross the barrier but also lowers the surface tension — does $\langle \tau \rangle$ decrease or increase, and why?
Critical droplet shape. Take a trajectory that nucleates and save the configuration a few sweeps before $M$ crosses zero. Identify the largest connected cluster of -1 spins (a flood-fill is enough) and record its size. Are nucleating droplets compact or ramified, and how does the typical size compare to $R_c = \sigma/|h|$ ?

Reference¶

K. Binder, Theory of first-order phase transitions, Rep. Prog. Phys. 50, 783 (1987). Binder (1987)