Part III: Phase Transitions | Chapter 2

The Ising Model

1D exact solution via transfer matrix, 2D Onsager solution, Monte Carlo simulation, and spontaneous magnetization

Historical Context

Ernst Ising solved the one-dimensional version of the model proposed by his advisor Wilhelm Lenz in his 1924 thesis. Finding no phase transition in 1D, Ising incorrectly concluded that the model would fail in all dimensions. Rudolf Peierls proved in 1936 that the 2D model does exhibit spontaneous magnetization, and Lars Onsager achieved the celebrated exact solution of the 2D Ising model in 1944 -- one of the greatest feats in mathematical physics.

The Ising model is the โ€œhydrogen atomโ€ of statistical mechanics: the simplest system exhibiting a genuine phase transition. It serves as the testing ground for every new theoretical and computational technique, from mean-field theory to renormalization group to Monte Carlo methods.

1. The Ising Hamiltonian

The Ising model consists of spins \(s_i = \pm 1\) on a lattice, interacting with nearest neighbors and an external field \(h\):

\[H = -J\sum_{\langle ij\rangle} s_i s_j - h\sum_i s_i\]

where \(J > 0\) for ferromagnetic coupling (parallel spins preferred) and the first sum runs over nearest-neighbor pairs. The partition function is:

\[Z = \sum_{\{s_i\}} e^{-\beta H} = \sum_{s_1=\pm 1}\cdots\sum_{s_N=\pm 1} \exp\left(\beta J\sum_{\langle ij\rangle}s_is_j + \beta h\sum_i s_i\right)\]
โ†‘โ†‘โ†‘โ†“โ†“โ†“= spin up (+1)= spin down (-1)domain wall
Figure 1. 2D Ising model (8ร—8 lattice) in a partially ordered state near Tแถœ. Gold circles represent up spins (+1) and dark circles represent down spins (-1). The dashed red line marks the domain wall separating two ordered regions.

2. Exact Solution in 1D: Transfer Matrix

Derivation 1: Transfer Matrix Method

For the 1D chain with periodic boundary conditions (\(s_{N+1} = s_1\)):

\[Z = \sum_{\{s_i\}} \prod_{i=1}^{N} \exp\left[K s_i s_{i+1} + \frac{h'}{2}(s_i + s_{i+1})\right]\]

where \(K = \beta J\) and \(h' = \beta h\). Define the transfer matrix:

\[T_{s_i,s_{i+1}} = \exp\left[K s_i s_{i+1} + \frac{h'}{2}(s_i + s_{i+1})\right]\]
\[\mathbf{T} = \begin{pmatrix} e^{K+h'} & e^{-K} \\ e^{-K} & e^{K-h'} \end{pmatrix}\]

The partition function becomes a trace: \(Z = \text{Tr}(\mathbf{T}^N) = \lambda_+^N + \lambda_-^N\), where the eigenvalues are:

\[\boxed{\lambda_{\pm} = e^K\cosh h' \pm \sqrt{e^{2K}\sinh^2 h' + e^{-2K}}}\]

In the thermodynamic limit (\(N \to \infty\)), \(Z \approx \lambda_+^N\)and the free energy per spin is:

\[f = -k_BT\ln\lambda_+ = -k_BT\ln\left[e^K\cosh h' + \sqrt{e^{2K}\sinh^2 h' + e^{-2K}}\right]\]

At \(h = 0\): \(\lambda_+ = 2\cosh K\), giving\(f = -k_BT\ln(2\cosh K)\). Since \(\lambda_+\) is analytic for all finite \(T > 0\), there is no phase transition in 1D.

3. Correlation Function in 1D

Derivation 2: Exponential Decay of Correlations

The spin-spin correlation function in 1D at \(h = 0\) can be computed exactly:

\[\langle s_i s_j \rangle = \left(\frac{\lambda_-}{\lambda_+}\right)^{|i-j|} = (\tanh K)^{|i-j|}\]

This has the Ornstein-Zernike form \(\langle s_i s_j\rangle = e^{-|i-j|/\xi}\) with correlation length:

\[\boxed{\xi = -\frac{1}{\ln(\tanh K)} = \frac{1}{\ln\coth(J/k_BT)}}\]

As \(T \to 0\), \(\xi \to \frac{1}{2}e^{2J/(k_BT)} \to \infty\): the correlation length diverges, but only at \(T = 0\). This is consistent with\(T_c = 0\) in 1D: a true phase transition exists only at absolute zero.

4. The 2D Ising Model: Onsager's Solution

Derivation 3: Critical Temperature and Free Energy

Onsager's 1944 solution of the 2D Ising model on a square lattice (at \(h = 0\)) is one of the monumental achievements of theoretical physics. The critical temperature satisfies:

\[\boxed{\sinh\left(\frac{2J}{k_BT_c}\right) = 1, \qquad \frac{k_BT_c}{J} = \frac{2}{\ln(1 + \sqrt{2})} \approx 2.269}\]

The exact free energy per spin is (Onsager, 1944):

\[-\beta f = \ln 2 + \frac{1}{2\pi}\int_0^{\pi}d\theta\,\ln\left[\cosh^2(2K) + \frac{1}{k_1}\sqrt{1 + k_1^2 - 2k_1\cos(2\theta)}\right]\]

where \(k_1 = 2\sinh(2K)/\cosh^2(2K)\). The specific heat diverges logarithmically:

\[C \sim -\frac{2k_B}{\pi}\left(\frac{2J}{k_BT_c}\right)^2\ln\left|1 - \frac{T}{T_c}\right| + \text{const}\]

giving \(\alpha = 0\) (logarithmic divergence, not a power law).

5. Spontaneous Magnetization

Derivation 4: Yang's Result

C.N. Yang derived the spontaneous magnetization in 1952:

\[\boxed{M = \left[1 - \sinh^{-4}\left(\frac{2J}{k_BT}\right)\right]^{1/8}}\]

Near \(T_c\): \(M \propto (T_c - T)^{1/8}\), giving the exact critical exponent \(\beta = 1/8\) (compared to the mean-field value \(\beta = 1/2\)). At \(T = 0\): \(M = 1\) (perfect alignment). The complete set of exact 2D Ising exponents is:

Exact 2D Ising Critical Exponents

\(\alpha = 0\) (logarithmic divergence of \(C\))

\(\beta = 1/8\) (spontaneous magnetization)

\(\gamma = 7/4\) (susceptibility)

\(\delta = 15\) (critical isotherm)

\(\nu = 1\) (correlation length)

\(\eta = 1/4\) (correlation function at \(T_c\))

6. Monte Carlo Simulation

Derivation 5: Metropolis Algorithm

The Metropolis algorithm generates spin configurations distributed according to the Boltzmann weight. At each step, propose flipping a single spin \(s_i \to -s_i\). The energy change is:

\[\Delta E = 2J s_i \sum_{j \in \text{nn}(i)} s_j + 2hs_i\]

Accept the flip with probability:

\[\boxed{P_{\text{accept}} = \min\left(1, e^{-\beta\Delta E}\right)}\]

This satisfies detailed balance: \(P(a \to b)/P(b \to a) = e^{-\beta(E_b - E_a)}\), ensuring convergence to the Boltzmann distribution. The algorithm can be optimized using the Wolff cluster algorithm, which flips correlated clusters of spins and dramatically reduces critical slowing down near \(T_c\).

Why No Phase Transition in 1D?

Peierls' argument: in 1D, a domain wall costs energy \(2J\) but gains entropy\(k_B\ln N\) (it can be placed at any of \(N\) positions). Since\(\Delta F = 2J - k_BT\ln N < 0\) for any \(T > 0\) in the thermodynamic limit, domain walls proliferate and destroy long-range order. In 2D, a domain wall is a line that costs energy proportional to its length \(L\), giving\(\Delta F = 2JL - k_BT \cdot S(L)\). For \(T < T_c\), the energy cost wins and order persists.

7. Computational Exploration

This simulation implements a Metropolis Monte Carlo simulation of the 2D Ising model, computing magnetization, energy, and specific heat across the phase transition.

2D Ising Model: Monte Carlo Simulation of Magnetization, Energy, and Specific Heat

Python
script.py180 lines

Click Run to execute the Python code

Code will be executed with Python 3 on the server

8. Summary and Key Results

Core Formulas

  • 1D: \(f = -k_BT\ln(2\cosh K)\), no transition
  • 2D \(T_c\): \(\sinh(2J/k_BT_c) = 1\)
  • Yang: \(M = [1 - \sinh^{-4}(2K)]^{1/8}\)
  • 1D correlation: \(\xi = -1/\ln(\tanh K)\)

Physical Insights

  • No phase transition in 1D at finite T (Peierls argument)
  • 2D Ising exponents differ from mean-field values
  • Metropolis algorithm: detailed balance ensures equilibrium
  • Critical slowing down near \(T_c\): cluster algorithms help

Practice Problems

Problem 1:Solve the 1D Ising model exactly using the transfer matrix method. Find $Z$, $\langle m\rangle$, and show no phase transition.

Solution:

1. For the 1D Ising chain with $N$ spins and periodic boundary conditions, the Hamiltonian is $H = -J\sum_{i=1}^N s_i s_{i+1} - h\sum_i s_i$. Define the $2\times 2$ transfer matrix:

$$\mathbf{T} = \begin{pmatrix} e^{K+h'} & e^{-K} \\ e^{-K} & e^{K-h'}\end{pmatrix}, \quad K = \beta J,\; h' = \beta h$$

2. The partition function is $Z = \text{Tr}(\mathbf{T}^N) = \lambda_+^N + \lambda_-^N$, where the eigenvalues are:

$$\lambda_\pm = e^K\cosh h' \pm \sqrt{e^{2K}\sinh^2 h' + e^{-2K}}$$

3. Since $\lambda_+ > \lambda_-$ for all finite $T$, in the thermodynamic limit $Z \to \lambda_+^N$. The free energy per spin is $f = -k_BT\ln\lambda_+$.

4. The magnetization per spin at $h = 0$ is:

$$\langle m\rangle = -\frac{\partial f}{\partial h}\bigg|_{h=0} = \frac{\sinh(\beta h)}{\sqrt{\sinh^2(\beta h) + e^{-4K}}}\bigg|_{h=0} = 0$$

5. Since $\lambda_+$ is analytic and non-degenerate for all $T > 0$, the free energy is analytic everywhere, so:

$$\boxed{\text{No phase transition in 1D at finite temperature}}$$

The eigenvalues only become degenerate at $T = 0$ (where $K\to\infty$), consistent with $T_c = 0$.

Problem 2:Derive the mean-field equation $m = \tanh(\beta Jzm + \beta h)$ and find $T_c = Jz/k_B$.

Solution:

1. In mean-field theory, replace the interaction of spin $s_i$ with its $z$ neighbors by an effective field. Write $s_j = m + (s_j - m)$ and neglect fluctuation terms:

$$H_{\text{MF}} = -\sum_i(Jzm + h)s_i + \frac{1}{2}NJzm^2$$

2. This is a system of independent spins in an effective field $h_{\text{eff}} = Jzm + h$. The partition function factorizes:

$$Z = e^{-\beta NJzm^2/2}\left[2\cosh\beta(Jzm + h)\right]^N$$

3. The self-consistency condition $m = \langle s_i\rangle$ gives:

$$\boxed{m = \tanh[\beta(Jzm + h)]}$$

4. At $h = 0$, expand for small $m$: $m \approx \beta Jzm - \tfrac{1}{3}(\beta Jzm)^3 + \cdots$. A non-trivial solution ($m \neq 0$) first appears when $\beta Jz = 1$:

$$\boxed{T_c = \frac{Jz}{k_B}}$$

5. For the square lattice ($z = 4$): $k_BT_c^{\text{MF}} = 4J$, which overestimates the exact Onsager result $k_BT_c \approx 2.269J$. Mean-field theory always overestimates $T_c$ because it neglects fluctuations.

Problem 3:Calculate the critical exponents $\beta$, $\gamma$, $\delta$ in mean-field theory.

Solution:

1. Exponent $\beta$ (magnetization near $T_c$, $h = 0$): Expand $m = \tanh(\beta Jzm)$ for small $m$:

$$m \approx \beta Jzm - \frac{1}{3}(\beta Jz)^3 m^3$$

Using $\beta Jz = T_c/T$ and $t = (T_c - T)/T_c$: $m^2 \approx 3t$, giving:

$$m \sim t^{1/2} \implies \boxed{\beta = 1/2}$$

2. Exponent $\gamma$ (susceptibility): Differentiate the mean-field equation with respect to $h$ at $h = 0$:

$$\chi = \frac{\partial m}{\partial h} = \frac{\beta(1-m^2)}{1 - \beta Jz(1-m^2)}$$

For $T > T_c$ ($m = 0$): $\chi = \beta/(1 - \beta Jz) \propto 1/(T - T_c)$, so:

$$\chi \sim |t|^{-1} \implies \boxed{\gamma = 1}$$

3. Exponent $\delta$ (critical isotherm, $T = T_c$): At $T = T_c$, the mean-field equation gives:

$$m = \tanh(m + \beta_c h) \approx m + \beta_c h - \frac{1}{3}(m + \beta_c h)^3$$

For small $m$ and $h$: $\beta_c h \approx \tfrac{1}{3}m^3$, so $h \propto m^3$:

$$h \sim m^{\delta} \implies \boxed{\delta = 3}$$

4. These mean-field exponents satisfy the Widom scaling relation $\gamma = \beta(\delta - 1) = \tfrac{1}{2}(3-1) = 1$, confirming internal consistency.

5. Mean-field exponents are exact above the upper critical dimension $d_c = 4$. For $d < 4$, fluctuations modify these values (e.g., the 2D Ising model has $\beta = 1/8$, $\gamma = 7/4$, $\delta = 15$).

Problem 4:Write the Metropolis algorithm for the 2D Ising model and estimate $T_c$.

Solution:

1. Initialize an $L \times L$ lattice with spins $s_{ij} = \pm 1$. For each Monte Carlo step:

2. Select a random spin $s_{ij}$. Compute the energy change for flipping it:

$$\Delta E = 2Js_{ij}(s_{i+1,j} + s_{i-1,j} + s_{i,j+1} + s_{i,j-1})$$

3. Accept the flip with probability $P = \min(1, e^{-\beta\Delta E})$:

- If $\Delta E \leq 0$: always accept (energy decreases).

- If $\Delta E > 0$: accept with probability $e^{-\beta\Delta E}$.

4. One "sweep" consists of $N = L^2$ attempted flips. After equilibration (typically 100-1000 sweeps), measure observables: $\langle |M|\rangle$, $\langle E\rangle$, $C_V = \beta^2(\langle E^2\rangle - \langle E\rangle^2)/N$.

5. The critical temperature is identified by the peak in $C_V(T)$ or $\chi(T)$. For $L = 20$, the peak occurs near $T \approx 2.3\,J/k_B$, close to the exact value:

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

Finite-size effects shift the apparent $T_c$ upward; extrapolation to $L\to\infty$ using finite-size scaling recovers the exact result.

Problem 5:Show that the correlation length $\xi$ diverges as $\xi \sim |T - T_c|^{-\nu}$ with $\nu_{\text{MF}} = 1/2$.

Solution:

1. In mean-field (Landau) theory, the free energy functional includes a gradient term:

$$F[m] = \int d^d r\left[\frac{1}{2}a\,t\,m^2 + \frac{1}{4}b\,m^4 + \frac{1}{2}c(\nabla m)^2 - hm\right]$$

where $t = (T - T_c)/T_c$ is the reduced temperature.

2. Consider a small perturbation $m(\mathbf{r}) = m_0 + \delta m(\mathbf{r})$. Linearizing the equation of state $\delta F/\delta m = 0$:

$$(at + 3bm_0^2)\,\delta m - c\,\nabla^2\delta m = 0$$

3. In Fourier space: $\delta m(\mathbf{q}) \propto 1/(at + 3bm_0^2 + cq^2)$. The correlation function is:

$$G(\mathbf{q}) = \frac{k_BT}{at + 3bm_0^2 + cq^2} = \frac{k_BT/c}{q^2 + \xi^{-2}}$$

4. The correlation length is $\xi^{-2} = (at + 3bm_0^2)/c$. For $T > T_c$ ($m_0 = 0$): $\xi^2 = c/(at) \propto t^{-1}$. For $T < T_c$ ($m_0^2 = -at/b$): $\xi^2 = c/(2a|t|) \propto |t|^{-1}$.

5. In both cases:

$$\boxed{\xi \sim |t|^{-1/2} \implies \nu_{\text{MF}} = \frac{1}{2}}$$

The exact 2D Ising value is $\nu = 1$, showing that fluctuations double the correlation length exponent. Mean-field theory underestimates divergences because it ignores the long-range correlations that dominate near criticality.

โ† Phase TransitionsMean-Field Theory โ†’
Rate this chapter: