$\newcommand{\I}{\mathrm{i}} \newcommand{\E}{\mathrm{e}} \newcommand{\D}{\mathop{}\!\mathrm{d}} \DeclareMathOperator{\Tr}{Tr} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\braket}[1]{\langle{#1}\rangle} \newcommand{\bbraket}[1]{\langle\!\langle{#1}\rangle\!\rangle} \newcommand{\bm}{\boldsymbol}$

Lectures on statistical mechanics.

Onsager phase transition: the square lattice ferromagnet

Following an original proposition of Pierre Weiss, Lenz around 1920 thought that a magnetic material was formed by elementary interacting magnets, however, at variance to Langevin and later Weiss, he assumed discrete and not continuous values of the magnetic moment or “magneton” (Weiss, 1911)¹. Ising introduced in his thesis (1924)² the idea of nearest neighbor interaction between “spins” taking values $s=\pm1$. The ising hamiltonian in a twodimensional lattice reads,

$$H = -J_x \sum_{x,y} s_{xy}s_{x+1y} - J_y \sum_s s_{xy}s_{xy+1} - B \sum_{x,y} s_{xy}\,,$$

where $x=0,1, \ldots, L-1$ and $y=0,1, \ldots, M-1$ are the coordinates of the lattice sites (the length unit is the lattice step $a=1$). The exchange energies in the $x$ and $y$ directions are $J_x$ and $J_y$, $B$ is the external applied field (in energy units), and $N=LM$ is the number of sites. The onedimensional case corresponds to $L=1$ (or $M=1$). The configuration space is generated by the $2^N$ vectors of the form,

$$s \in S = \{(s_0, \ldots, s_n, \ldots, s_{N-1}), \; s_n = \pm 1, \; n=0,\ldots,N-1\}\,.$$

Before considering the planar case, we rework the linear ising model in a way such that we can generalize the transfer matrix method. We write the onedimensional hamiltonian as,

\begin{align*}H_{1D} &= - J \sum_x (s_x s_{x+1} - 1) - B \sum_x s_x \\ &= \frac{J}{2} \sum_x (s_x - s_{x+1})^2 - \frac{B}{2} \sum_x (s_x + s_{x+1}) \end{align*}

($J$ is the coupling constant in 1D) which gives a ground state at zero energy (for $B=0$). The second form is symmetric in $s_x, s_{x+1}$. Using the first equation we can write the partition function in terms of a product of $2 \times 2$ matrices,

$$Z(K,h) = \Tr \E^{-H_{1D}/T} = \Tr (V_1 V_3)^N$$

or using the cyclic property of the trace, we can obtain the symmetric form,

$$Z(K,h) = \Tr \big(V_3^{1/2} V_1 V_3^{1/2} \big)^N = \Tr V^N \,,$$

where $V= V_3^{1/2} V_1 V_3^{1/2}$ is the so called transfer matrix,³ and

$$V_1(K) = \begin{pmatrix} 1 & \E^{-2K} \\ \E^{-2K} & 1 \end{pmatrix} \,,$$

and

$$V_3(h) = \begin{pmatrix} \E^{h} & 0\\ 0 & \E^{-h} \end{pmatrix} \,,$$

with $K=J/T$ and $h=B/T$. Using instead the symmetric form (the transfer matrix $V$),

$$V = \begin{pmatrix} \E^{h} & \E^{-2K} \\ \E^{-2K} & \E^{-h} \end{pmatrix}\,,$$

we get a simpler expression to compute $Z$. Indeed, for the symmetric form of the transfer matrix the partition function can be evaluated using its eigenvalues $\lambda_\pm$,

$$Z(K,h) = \Tr V^N = \lambda_+^N + \lambda_-^N \approx \lambda_+^N\,, \quad \lambda_+ > \lambda_-$$

where the last approximation is valid in the thermodynamic limit $N\rightarrow \infty$. This is not possible in the other case, where $V$ is expressed as a product of two noncommuting matrices.

Instead of directly computing the eigenvalues of $V$, we observe that for twodimensional matrices we have the identity EX,

$$\E^{\varepsilon \hat{\bm n} \cdot \bm \sigma} = \cosh(\varepsilon)\, 1_2 + \sinh(\varepsilon)\, (\hat{\bm n} \cdot \bm \sigma)$$

which will allow us to transform $V_1, V_3$ and $V$ into the exponentials of a matrix; in the previous formula $\varepsilon$ is a constant, $\hat{\bm n}$ a unit vector, $1_2$ the unit matrix, and $\bm \sigma$ the vector of pauli matrices. Identifying the exponential form of $V_1$ and $V_3$ is straightforward:

$$V_3 = \E^{h \sigma_3}$$

and

$$V_1 = 1 + \E^{-2K}\sigma_1 = \frac{\E^{K^\star \sigma_1}}{\cosh(K^\star)}$$

where we defined the new constant,

$$\tanh(K^\star) = \E^{-2K}\,.$$

In the case of the matrix

$$V = \sqrt{\det V} \E^{-H}, \quad H = \varepsilon \hat{\bm n} \cdot \sigma\,.$$

we should solve the system of equations in $\varepsilon$ and $\hat{\bm n}$, to obtain EX,

\begin{align*} \cosh \varepsilon & = \cosh K^\star \cosh h \\ \hat{\bm n} &= \frac{(\sinh K^\star, 0, \cosh K^\star \sinh h)}{\sqrt{\cosh K^\star \cosh h - 1}}\,. \end{align*}

The eigenvalues of $H=\varepsilon \hat{\bm n} \cdot \bm \sigma$ are readily obtained,

$$\varepsilon = \pm \ln\left[ \cosh K^\star \cosh h + \sqrt{\cosh^2 K^\star \cosh^2 h -1} \right]$$

(the eigenvalues of $\hat{\bm n} \cdot \bm \sigma$ are obviously $\pm1$).

It is interesting to note that the relation between $K$ and $K^\star$ is dual, in the sense that EX,

$$\tanh(K^\star) = \E^{-2K}\,,\quad \tanh(K) = \E^{-2K^\star}\,.$$

EX Compute the free energy per spin and show that the mean magnetization is given by
$$\braket{s} = \frac{\sinh h}{\sqrt{\cosh^2 K^\star \cosh^2 h - 1}}$$
and show that this result is equivalent to the one of section Lattices.

The alternative formulation in terms of the product $V_1V_3$, can be handled in a particular limit, which will allow us to conveniently combine the two exponentials. The interest of the method is that it can be generalized to our case of interest, the twodimensional lattice. The Baker-Campbell-Hausdorff formula says, for two operators operators $A,B$,

$$\E^{\tau A}\E^{\tau B} = \E^{\tau (A+B) + (\tau^2/2) [A,B] + \cdots}$$

where $\tau$ is a real constant, and the omitted terms are higher order in $\tau$. In the limit $\tau \rightarrow 0$ one may replace the product of exponentials by the exponential of the sum. We define new parameters $\lambda$ and $\tau$, such that,

$$K^\star = \lambda \tau, \quad h = \tau$$

and write,

$$V_1 V_3 \approx \frac{1}{\cosh K^\star} \E^{-\tau H}, \quad H_\lambda = -\lambda \sigma_1 - \sigma_3\,.$$

(Note that the prefactor in this expression is independent of $h$, thus irrelevant for the computation of the magnetic properties of the system.) Therefore, in the low energy limit the onedimensional ising model can be described by a hamiltonian $H_\lambda$ corresponding to a single quantum spin. In this limit $\tau \sim 1/T$ becomes a continuous imaginary time of a quantum system: $\tau = \I t$.⁴ One may establish the analogy:

Quantum system	Statistical system
$D$ dimensions	$D+1$ dimensions
ground state	equilibrium state
expected value	ensemble average
fundamental energy	free energy

The eigenvalues of $\tau H_\lambda$ are,

$$\varepsilon_\pm = \pm \tau \sqrt{\lambda^2 + 1}$$

which gives a free energy (per spin)

$$f(\lambda) = -T \varepsilon_- = -T \tau \sqrt{\lambda^2 +1},$$

remark that only the lowest eigenvalue (fundamental level of $H_\lambda$) contributes to the free energy in the thermodynamic limit, which gives in turn, the mean magnetization per spin,

$$\braket{s} = - \frac{\partial f}{\partial h} = \frac{h}{\sqrt{K^{\star 2} + h^2} }.$$

Now we may compare the result of the exact method with the continuum limit result. In the limit of small $K^\star$ and $h$, the exact eigenvalue $\varepsilon$ can be approximated by:

$$\varepsilon = \ln\left[ \cosh K^\star \cosh h + \sqrt{\cosh^2 K^\star \cosh^2 h -1} \right] \approx \sqrt{K^{\star 2} + h^2}, \quad (K^*,h \rightarrow 0),$$

in accordance with the $H_\lambda$ result EX.

We apply now the continuum limit method used to obtain $H_\lambda$, to the twodimensional lattice case. We consider the $B=0$ case, which is enough to exhibit the existence of a phase transition, our main goal, through the existence of a singularity in the free energy for some value of the parameters $T, J_x, J_y$. Our first task is to find the transfer matrix. We then plan to transform the transfer matrix into a unique exponential of an effective hamiltonian easily amenable to a diagonal form.

Twodimensional square lattice in the $(x,\tau)$ coordinates; the $x$ axis remains discrete while the $\tau$ axis goes to the continuum limit.

In analogy with the onedimensional case we will treat the $y$ direction as a (imaginary) time, and the other direction as a normal discrete linear lattice; the transfer matrix will then represent the evolution operator of the discrete system in $\tau$ (see the above figure). The physical properties of the system depend on two parameters:

$$K_x = J_x/T,\quad K_\tau = J_y/T\,.$$

With this notation the partition function is,

$$Z = \sum_s \exp\left[K_x \sum_\tau \sum_{x=0}^{L-1} s_{x}(\tau) s_{x+1}(\tau) \right] \exp\left[ K_\tau \sum_\tau \sum_{x=0}^{L-1} \big( s_x(\tau) s_x(\tau + 1) -1\big) \right]\,,$$

where we substracted a constat energy term in $H$. The first exponential is diagonal in the variable $\tau$, and can be compared with the $h$ term (and the matrix $V_3$) in the onedimensional case, while the second exponential represents the interaction between successive rows, and can be treated in the variable $\tau$ as the onedimensional system in the variable $x$ (and the matrix $V_1$): instead of summing on the two orientations of the spin, we sum over the whole configuration of the row. Therefore, we can write:

$$Z = \Tr (V_1 V_3)^M$$

with

$$V_1 = \frac{\E^{K_\tau^\star \sum_x \sigma_1(x)}}{(\cosh K_\tau^\star)^L}$$

and

$$V_3 = \E^{K_x \sum_x \sigma_3(x) \sigma_3(x+1)}\,.$$

At variance to the onedimensional lattice, the $V$ matrices are $2^L \times 2^L$ which makes the search for the eigenvalues much more difficult.

The exact computation of the square lattice system partition function was obtained in 1944 by Onsager.⁵ Later, in 1964, Schultz, Mattis and Lieb found a mapping to a fermion problem that allowed them to greatly simplify the original Onsager’s computation.⁶ Here, we are interested in the thermodynamic limit, in which the system is translationally invariant and independent of the boundary conditions. In this limit the computation of the free energy using the fermion method is straightforward. However, before discussing the continuous limit, let us explore the parameter phase space $(K_x,K_\tau)$.

Parameter phase space of the anisotropic 2D ising model. The line $\tanh K_x = \E^{-2K_\tau}$ separates the disorderd (paramagnetic), when both $K_x,K_\tau$ are small, and ordered (ferromagnatic) phases, for large coupling. The isotropic case corresponds to the line $K_x = K_\tau$; the continuous limit is on the small $K_x = \tau$ side: it contains a transition point (black dot). A stretching of the vertical axis allows to recover the isotropic case (black arrow) $K_c = -0.5\ln(\sqrt{2}-1)$.

In the isotropic case we have $K_x = K_\tau$ and the duality relation:

$$\tanh K^\star = \E^{-2K},\, \text{or } \sinh 2K^\star \sinh 2K = 1\,.$$

This last relation shows that when $K$ is large (at fixed $J$, this means low temperature), then $K^\star$ must be small (high temperature), and vice versa. The common region where both low and high temperature properties overlap, is the transition region, corresponding then to the fixed point

$$K=K^\star=K_c,\quad \tanh K_c = \E^{-2K_c}$$

which defines the critical temperature $T_c=-2/\ln(\sqrt{2} -1) = 2.269\,J$: below $T_c$ (small $K^\star$) the system is ordered and disordered above it (large $K^\star$).

In the anistropic case we have the duality relations:

$$\tanh K_\tau^\star = \E^{-2K_\tau}\,, \quad \tanh K_\tau^\star = \E^{-2K_\tau}$$

and, a point in the plane $(K_x, K_\tau)$, corresponds to a point in the dual plane $(K^\star_\tau, K^\star_x)$ (note the inversion of indices in the dual plane); the self-dual points are then on the line

$$(K_x, K_\tau) = (K_\tau^\star, K_x^\star)$$

contains the critical points:

$$\sinh 2K_x \sinh 2K_\tau = 1$$

represented in the figure to separate the two phases.

Therefore, using the parametrization,

$$K_\tau^\star = \lambda \tau,\quad K_x = \tau$$

we know that for some value of $\lambda$, in fact $\lambda_c = 1$ as we shall demonstrate, there is a phase transition, we can combine the two exponentials in the small $\tau$ limit.This parametrization rescale the $\tau$ axis to recover, in the continuum limit, the isotropic transition (corresponding to the arrow mapping the anisotropic critical point to the isotropic one, in the figure above). Consequently, from the expressions of $V_1$ and $V_3$, in the small $\tau$ limit we write:

\begin{equation} Z(\lambda) = \Tr \frac{\E^{-\tau H(\lambda)}}{(\cosh K_\tau^\star)^L}, \quad H(\lambda) = - \lambda \sum_x \sigma_1(x) - \sum_x \sigma_3(x) \sigma_3(x+1)\,. \label{e:HQ} \end{equation}

When $\lambda \rightarrow 0$ we can drop the first term and deal with the interaction term; because it only depends on $\sigma_3$, the full spin up or full spin down states are ground energy eigenstates; in this case the ground state is degenerate and corresponds to an ordered phase. We deduce that small $\lambda$ corresponds to the low temperature regime. Instead, when $\lambda \rightarrow \infty$ we are in a disordered state. Indeed, in this limit we keep the firs term, corresponding to independent spins, and any state which is a product of eigenstates of $\sigma_1$ is an eigenstate. The lowest energy state is the product of the eigenvectors of $\sigma_1$ with eigenvalue 1, and therefore it is non-degenerated. In the first case ($\lambda \ll 1$) the expected value (magnetization per site) $\braket{\sigma_3} = \pm1$, in the second case ($\lambda \gg 1$) the mean magnetization vanishes $\braket{\sigma_3} = 0$; for some intermediate value we have a transition between the ferromagnetic and the paramagnetic states.

To diagonalize $H(\lambda)$ we map the spin operators (pauli matrices) into a pair of majorana fermions $c$ and $d$, using the Jordan-Wigner transformation. A fermion operator $f_x$, which annihilates a fermion at site $x$, can be represented by two real fermion operators (majoranas):

$$f_x = \frac{c_x + \I d_x}{\sqrt{2}}, \quad f^\dagger_x = \frac{c_x - \I d_x}{\sqrt{2}} \,,$$

which satisfy the anticommutation rules,

$$\{c(x),d(x)\} = \{c(x), d(y)\} = 0, \quad \{c(x),c(y)\} = \{d(x),d(y)\} = \delta_{xy} \,.$$

The analogous spin operators are,

$$\sigma^- = \frac{\sigma_1 - \I \sigma_2}{2},\quad \sigma^+ = \frac{\sigma_1 + \I \sigma_2}{2},$$

which satisfy the mixed (boson-fermion) commutation relations,

$$[\sigma^+(x),\sigma^-(x)]=0 ,\quad \{\sigma^+(x),\sigma^-(y)\}=0 \; x \ne y$$

therefore the relation between the two set of operators is not trivial, even if both have the same dimension, the two orientations of the spin, and the two values $0$ or $1$ of the occupation number. The relation which ensures the correct fermion commutation rules is in fact nonlocal EX:

\begin{align*} c(x) &= \frac{1}{\sqrt{2}}\left[\prod_{y<x} \sigma_1(y)\right] \sigma_2(x), \\ d(x) &= \frac{1}{\sqrt{2}} \left[\prod_{y<x} \sigma_1(y)\right] \sigma_3(x) \end{align*}

These formulas can easily be inverted EX:

\begin{align*} \sigma_1(x) &= -2\I c(x)d(x) \\ \sigma_2(x) &= \sqrt{2} \left[\prod_{y<x} \big( -2\I c(y)d(y)) \big) \right] c(x) \\ \sigma_3(x) &= \sqrt{2} \left[\prod_{y<x} \big( -2\I c(y)d(y)) \big) \right] d(x)\,. \end{align*}

The product of $\sigma_3$ is transformed as,

$$\sigma_3(x) \sigma_3(x+1) = 2\I c(x) d(x+1)\,.$$

Combining the above expresions we obtain the representation of $H(\lambda)$ in terms of the majoranaas,

$$H(\lambda) = 2\I \lambda \sum_x c(x) d(x) - 2\I \sum_x c(x) d(x+1)\,.$$

The first term is diagonal in $x$ but not the second; however, because of the translation invariance we can perform a fourier transformation to $k$ space:

$$c(x) = \int_{-\pi}^\pi \frac{\D k}{2\pi} \E^{\I kx} c_k\,,$$

and an analogous formula for $d(x) \rightarrow d_k$. Replacing this expresion into $H(\lambda)$, we obtain EX,

$$H(\lambda) = \int_{0}^\pi \frac{\D k}{2\pi} H_k(\lambda),\quad H_k(\lambda) = 2\I \lambda (c^\dagger_k d_k + c_k d^\dagger_k) - \big( \E^{\I k} c^\dagger_k d_k + \E^{-\I k} c_k d^\dagger_k \big)\,.$$

Note that the integral goes over $k > 0$, to explicitly exhibit the hermitian structure of the hamiltonian. We used the formula, valid in the limit $L \gg 1$:

$$\lim_{L \rightarrow \infty} \sum_{x=L/2+1}^{L/2} \E^{\I k x} = 2\pi \delta(k)\,.$$

We can conveniently put $H_k(\lambda)$ in matrix form:

$$H_k(\lambda) = \begin{pmatrix} c^\dagger_k & d^\dagger_k \end{pmatrix} \begin{pmatrix} 0 & 2\I (\lambda - \E^{\I k}) \\ -2\I (\lambda - \E^{-\I k}) & 0 \end{pmatrix} \begin{pmatrix} c_k \\ d_k \end{pmatrix} \equiv C^\dagger \hat{H} C\,,$$

more suitable for diagonalization. The $2 \times 2$ matrix $\hat{H}$ has the eigenvalues $\pm \varepsilon_\lambda$:

\begin{equation} \label{eigen} \varepsilon_\lambda(k) = 2 \sqrt{(\lambda - 1)^2 + 2\lambda(1 - \cos k)} \end{equation}

we can thus define a unitary transformation $U$ such that

$$H_k = C^\dagger \hat{H} C = A^\dagger (U^\dagger \hat{H} U) A, \quad C = U A, \quad A = \begin{pmatrix} a_k \\ b_k \end{pmatrix}$$

in the new basis $A$ the hamiltonian becomes,

$$H_k(\lambda) = \varepsilon_\lambda(k) (a^\dagger_k a_k - b^\dagger_k b_k)\,,$$

whose lowest energy eigenvalue is $-\varepsilon$. Note moreover that $\varepsilon$ is an even function of $k$. Is just this eigenvalue that will contribute to the partition function, giving us an exact solution of the phase transition as we see now. The explicit expression of $U$ is EX,

$$U_k(\lambda) = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & -\E^{\I \theta_k} \\ \E^{-\I \theta_k} & 1 \end{pmatrix}, \quad \E^{\I \theta_k(\lambda)} = \I \frac{\lambda - \E^{\I k}}{|\lambda - \E^{\I k}|} \,.$$

We can now compute the partition function (for $L \rightarrow \infty$):

\begin{align*} Z(\lambda) &= \Tr (V_1V_3)^M = \frac{1}{(\cosh K_\tau^\star)^L} \Tr \exp\left[ - M\tau H(\lambda) \right] \\ &= \frac{1}{(\cosh K_\tau^\star)^L} \Tr \exp \left[- M\tau \int_{-\pi}^\pi \frac{\D k}{2\pi} H_k(\lambda) \right]\,, \end{align*}

then, taking the trace, we obtain the explicit formula,

$$Z(\lambda) = \exp\left[ -M\tau \int_{-\pi}^\pi \frac{\D k}{2\pi} \varepsilon_\lambda(k) \right]\,,$$

from which we compute the free energy per spin $f = F/M = - (T/M) \ln Z$,

$$f(\lambda) = f_0(K_\tau) - J_x \int_{-\pi}^\pi \frac{\D k}{2\pi} \sqrt{(\lambda - 1)^2 + 2\lambda(1 - \cos k)}\,,$$

we took the lowest eigenvalue, recall also that $\tau = K_x = J_x/T$, in addition, we included in the first term the prefactor in $Z$, which do not contribute to the free energy near the singularity. Although the integral is of the elliptic type, we are only interested in finding its singular behavior, which must occur when

$$t \equiv \frac{T-T_c}{T_c} = \lambda - 1$$

is small ($t=0$ is then the critical point). To see this we compute the specific heat near $t=1$; it is given by the second derivative

$$c_V = -T \frac{\partial^2 f}{\partial T^2} = - \frac{T}{T_c^2} \frac{\partial^2 f}{\partial t^2}.$$

In summary, in the neighborhood of $t\rightarrow0$ we may write EX,

$$c_V \sim \int_0^\pi \frac{\D k}{\sqrt{t^2 + 2(1-\cos k)}} \,,$$

(the symbol $\sim$ means that we keep only the dominant terms, without irrelevant factors or regular terms) changing the variable $k/t \rightarrow k$, we obtain,

$$c_V \sim \int_0^{\pi/t} \frac{\D k}{\sqrt{1 + (2/t^2)(1-\cos tk)}} \sim \int_0^{\pi/t} \frac{\D k}{\sqrt{1 + k^2}}\,,$$

($\cos x \approx 1 - x^2/2$) or,

$$c_V \sim - \ln|t|\,.$$

Therefore, the free energy is singular leading to a logarithmically diverging specific heat. This is the signature of the ferromagnetic-paramagnetic transition at $T = T_c$.

The specific heat near the transition.

Exercise

Onsager found, from the computation of the spin-spin correlation at large distance, the exact formula of the mean magnetization:

$$\braket{s} = \left( 1 - \frac{1}{\sinh^2 K_x \sinh^2 2 K_\tau} \right)^{ \frac{1}{8} }\,,$$

here written in the general anisotropic case. Show that the reparametrization

$$K_x = \tau,\quad K^\star_\tau = \lambda \tau$$

where $\tanh K^\star_\tau = \E^{-2 K_\tau}$, leads in the $\tau \rightarrow 0$ limit, to the singular behavior ($\lambda \le 1$, in the ferromagnetic region)

$$\braket{s} \sim (1-\lambda)^{\frac{1}{8}},$$

near the critical point $\lambda = 1$. Therefore, the critical exponent is $\beta = 1/8$ (different to the mean field value $1/2$).

Notes

Exercises, generally intended to fill the omitted steps of a calculation, are signaled by EX.

We followed in this lecture the book: R. Shankar, Quantum Field Theory and Condensed Matter, Cambridge, 2017.

R.Voltz, V. Pierron-Bohnes, M.-C. Cadeville, Pierre Weiss (1865-1940) (.pdf) ↩
M. Niss, History of the Lenz-Ising Model 1920–1950: From Ferromagnetic to Cooperative Phenomena, Archive for History of Exact Sciences 59, 267 (2005) (.pdf) ↩
Kramers, H. A., and Wannier, G. H., Statistics of the Two-Dimensional Ferromagnet, Phys. Rev. 60, 252 (1941) ↩
E. Fradkin and L. Susskind, Order and disorder in gauge systems and magnets, Phys. Rev. D 17, 2637 (1978) ↩
L. Onsager, A two-dimensional model with an order-disorder transition, Phys. Rev. 65, 117 (1944) ↩
T. Schultz, D. Mattis and E. Lieb, Two-dimensional ising model as a soluble problem of many fermions, Rev. Mod. Phys. 36, 856 (1964) ↩

Quantum system	Statistical system
\(D\) dimensions	\(D+1\) dimensions
ground state	equilibrium state
expected value	ensemble average
fundamental energy	free energy