5 Phase Transitions

How do we figure out if this indeed happens? Just because both liquid and gas states can exist, doesn’t mean that they can cohabit. It might be that one is preferred over the other. We already saw some conditions that must be satisfied in order for two systems to sit in equilibrium back in Section 1. Mechanical and thermal equilibrium are guaranteed if two systems have the same pressure and temperature respectively. But both of these are already guaranteed by construction for our two liquid and gas solutions: the two solutions sit on the same isotherm and at the same value of $p$. We’re left with only one further requirement that we must satisfy which arises because the two systems can exchange particles. This is the requirement of chemical equilibrium,

${\mu }_{\mathrm{liquid}}={\mu }_{\mathrm{gas}}$

(5.182)

Because of the relationship (4.175) between the chemical potential and the Gibbs free energy, this is often expressed as

$g_{\mathrm{liquid}}=g_{\mathrm{gas}}$

(5.183)

where $g=G/N$ is the Gibbs free energy per particle.

Notice that all the equilibrium conditions involve only intensive quantities: $p$, $T$ and $\mu$. This means that if we have a situation where liquid and gas are in equilibrium, then we can have any number $N_{\mathrm{liquid}}$ of atoms in the liquid state and any number $N_{\mathrm{gas}}$ in the gas state. But how can we make sure that chemical equilibrium (5.182) is satisfied?

I should confess that there’s something slightly dodgy about the Maxwell construction. We already argued that the part of the isotherm with $dp/dv>0$ suffers an instability and is unphysical. But we needed to trek along that part of the curve to derive our result. There are more rigorous arguments that give the same answer.

For each isotherm, we can determine the pressure at which the liquid and gas states are in equilibrium. The gives us the co-existence curve, shown by the dotted line in Figure 37. Inside this region, liquid and gas can both exist at the same temperature and pressure. But there is nothing that tells us how much gas there should be and how much liquid: atoms can happily move from the liquid state to the gas state. This means that while the density of gas and liquid is fixed, the average density of the system is not. It can vary between the gas density and the liquid density simply by changing the amount of liquid. The upshot of this argument is that inside the co-existence curves, the isotherms simply become flat lines, reflecting the fact that the density can take any value. This is shown in graph on the right of Figure 37.

To illustrate the physics of this situation, suppose that we sit at some fixed density $\rho =1/v$ and cool the system down from a high temperature to at a point inside the co-existence curve so that we’re now sitting on one of the flat lines. Here, the system is neither entirely liquid, nor entirely gas. Instead it will split into gas, with density $1/v_{\mathrm{gas}}$, and liquid, with density $1/v_{\mathrm{liquid}}$ so that the average density remains $1/v$. The system undergoes phase separation. The minimum energy configuration will typically be a single phase of liquid and one of gas because the interface between the two costs energy. (We will derive an expression for this energy in Section 5.5). The end result is shown on the right. In the presence of gravity, the higher density liquid will indeed sink to the bottom.

The van der Waals states which lie between the spinodal curve and the co-existence curve are good states. But they are meta-stable. One can show that their Gibbs free energy is higher than that of the liquid-gas equilibrium at the same $p$ and $T$. However, if we compress the gas very slowly we can coax the system into this state. It is known as a supercooled vapour. It is delicate. Any small disturbance will cause some amount of the gas to condense into the liquid. Similarly, expanding a liquid beyond the co-existence curve results in an meta-stable, superheated liquid.

Either side of the line, all particles are either in the gas or liquid phase. We know from (5.183) that the Gibbs free energies (per particle) of these two states are equal,

$g_{\mathrm{liquid}}=g_{\mathrm{gas}}$

So $G$ is continuous as we move across the line of phase transitions. Suppose that we sit on the line itself and move along it. How does $g$ change? We can easily compute this from (4.174),

$\mathrm{ }d{G}_{\mathrm{liquid}}=-S_{\mathrm{liquid}}dT+V_{\mathrm{liquid}}dp\mathit{\hspace{1em}\hspace{1em} }\Rightarrow \mathit{\hspace{1em}\hspace{1em} }d{g}_{\mathrm{liquid}}=-s_{\mathrm{liquid}}dT+v_{\mathrm{liquid}}dp$

where, just as $g=G/N$ an $v=V/N$, the entropy density is $s=S/N$. Equating this with the free energy in the gaseous phase gives

$d{g}_{\mathrm{liquid}}=-s_{\mathrm{liquid}}dT+v_{\mathrm{liquid}}dp=d{g}_{\mathrm{gas}}=-s_{\mathrm{gas}}dT+v_{\mathrm{gas}}dp$

This can be rearranged to gives us a nice expression for the slope of the line of phase transitions in the $p-T$ plane. It is

${\displaystyle \frac{dp}{dT}}={\displaystyle \frac{s_{\mathrm{gas}}-s_{\mathrm{liquid}}}{v_{\mathrm{gas}}-v_{\mathrm{liquid}}}}$

We usually define the specific latent heat

$L=T(s_{\mathrm{gas}}-s_{\mathrm{liquid}})$

This is the energy released per particle as we pass through the phase transition. We see that the slope of the line in the $p-T$ plane is determined by the ratio of latent heat released in the phase transition and the discontinuity in volume. The result is known as the Clausius-Clapeyron equation,

${\displaystyle \frac{dp}{dT}}={\displaystyle \frac{L}{T(v_{\mathrm{gas}}-v_{\mathrm{liquid}})}}$

(5.185)

There is a classification of phase transitions, due originally to Ehrenfest. When the $n^{\mathrm{th}}$ derivative of a thermodynamic potential (either $F$ or $G$ usually) is discontinuous, we say we have an $n^{\mathrm{th}}$ order phase transition. In practice, we nearly always deal with first, second and (very rarely) third order transitions. The liquid-gas transition releases latent heat, which means that $S=-\partial F/\partial T$ is discontinuous. Alternatively, we can say that $V=\partial G/\partial p$ is discontinuous. Either way, it is a first order phase transition. The Clausius-Clapeyron equation (5.185) applies to any first order transition.

As we approach $T\to T_c$, the discontinuity diminishes and $S_{\mathrm{liquid}}\to S_{\mathrm{gas}}$. At the critical point $T=T_c$ we have a second order phase transition. Above the critical point, there is no sharp distinction between the gas phase and liquid phase.

For most simple materials, the phase diagram above is part of a larger phase diagram which includes solids at smaller temperatures or higher pressures. A generic version of such a phase diagram is shown to the right. The van der Waals equation is missing the physics of solidification and includes only the liquid-gas line.

Dubious as its theoretical foundation is, the law of corresponding states is the first suggestion that something remarkable happens if we describe a gas in terms of its reduced variables. More importantly, there is striking experimental evidence to back this up! Figure 42 shows the Guggenheim plot, constructed in 1945. The co-existence curve for 8 different gases in plotted in reduced variables: $\overline{T}$ along the vertical axis; $\overline{\rho }=1/\overline{v}$ along the horizontal. The gases vary in complexity from the simple monatomic gas $Ne$ to the molecule $C{H}_4$. As you can see, the co-existence curve for all gases is essentially the same, with the chemical make-up largely forgotten. There is clearly something interesting going on. How to understand it?

First, we can ask what happens to the difference in (inverse) densities $v_{\mathrm{gas}}-v_{\mathrm{liquid}}$ as we approach the critical point along the co-existence curve. For , or equivalently , the reduced van der Waals equation (5.187) has two stable solutions,

$\overline{p}={\displaystyle \frac{8\overline{T}}{3{\overline{v}}_{\mathrm{liquid}}-1}}-{\displaystyle \frac{3}{{\overline{v}}_{\mathrm{liquid}}^2}}={\displaystyle \frac{8\overline{T}}{3{\overline{v}}_{\mathrm{gas}}-1}}-{\displaystyle \frac{3}{{\overline{v}}_{\mathrm{gas}}^2}}$

If we solve this for $\overline{T}$, we have

$\overline{T}={\displaystyle \frac{(3{\overline{v}}_{\mathrm{liquid}}-1)(3{\overline{v}}_{\mathrm{gas}}-1)({\overline{v}}_{\mathrm{liquid}}+{\overline{v}}_{\mathrm{gas}})}{8{\overline{v}}_{\mathrm{gas}}^2{\overline{v}}_{\mathrm{liquid}}^2}}$

Notice that as we approach the critical point, ${\overline{v}}_{\mathrm{gas}},{\overline{v}}_{\mathrm{liquid}}\to 1$ and the equation above tells us that $\overline{T}\to 1$ as expected. We can see exactly how we approach $\overline{T}=1$ by expanding the right right-hand side for small $ϵ\equiv {\overline{v}}_{\mathrm{gas}}-{\overline{v}}_{\mathrm{liquid}}$. To do this quickly, it’s best to notice that the equation is symmetric in ${\overline{v}}_{\mathrm{gas}}$ and ${\overline{v}}_{\mathrm{liquid}}$, so close to the critical point we can write ${\overline{v}}_{\mathrm{gas}}=1+ϵ/2$ and ${\overline{v}}_{\mathrm{liquid}}=1-ϵ/2$. Substituting this into the equation above and keeping just the leading order term, we find

$\overline{T}\approx 1-{\displaystyle \frac{1}{16}}{({\overline{v}}_{\mathrm{gas}}-{\overline{v}}_{\mathrm{liquid}})}^2$

Or, re-arranging, as we approach $T_c$ along the co-existence curve,

$v_{\mathrm{gas}}-v_{\mathrm{liquid}}\sim {(T_c-T)}^{1/2}$

(5.188)

This is the answer to our first question.

Our second variant of the question is: how does the volume change with pressure as we move along the critical isotherm. It turns out that we can answer this question without doing any work. Notice that at $T=T_c$, there is a unique pressure for a given volume $p(v,T_c)$. But we know that $\partial p/\partial v={\partial }^{2}p/\partial v^2=0$ at the critical point. So a Taylor expansion around the critical point must start with the cubic term,

$p-p_c\sim {(v-v_c)}^3$

(5.189)

This is the answer to our second question.

Our third and final variant of the question concerns the compressibility, defined as

$\kappa =-{{\displaystyle \frac{1}{v}}{\displaystyle \frac{\partial v}{\partial p}}|}_T$

(5.190)

We want to understand how $\kappa$ changes as we approach $T\to T_c$ from above. In fact, we met the compressibility before: it was the feature that first made us nervous about the van der Waals equation since $\kappa$ is negative in the unstable region. We already know that at the critical point ${\partial p/\partial v|}_{T_c}=0$. So expanding for temperatures close to $T_c$, we expect

${{\displaystyle \frac{\partial p}{\partial v}}|}_{T;v=v_c}=-a(T-T_c)+\mathrm{\dots }$

This tells us that the compressibility should diverge at the critical point, scaling as

$\kappa \sim {(T-T_c)}^{-1}$

(5.191)

We now have three answers to three questions: (5.188), (5.189) and (5.191). Are they right?! By which I mean: do they agree with experiment? Remember that we’re not sure that we can trust the van der Waals equation at the critical point so we should be nervous. However, there is also reason for some confidence. Notice, in particular, that in order to compute (5.189) and (5.191), we didn’t actually need any details of the van der Waals equation. We simply needed to assume the existence of the critical point and an analytic Taylor expansion of various quantities in the neighbourhood. Given that the answers follow from such general grounds, one may hope that they provide the correct answers for a gas in the neighbourhood of the critical point even though we know that the approximations that went into the van der Waals equation aren’t valid there. Fortunately, that isn’t the case: physics is much more interesting than that!

The experimental results for a gas in the neighbourhood of the critical point do share one feature in common with the discussion above: they are completely independent of the atomic make-up of the gas. However, the scaling that we computed using the van der Waals equation is not fully accurate. The correct results are as follows. As we approach the critical point along the co-existence curve, the densities scale as

$v_{\mathrm{gas}}-v_{\mathrm{liquid}}\sim {(T_c-T)}^{\beta }\mathit{\hspace{1em}\hspace{1em} }\mathrm{with}\mathit{\hspace{1em}\hspace{1em} }\beta \approx 0.32$

(Note that the exponent $\beta$ has nothing to do with inverse temperature. We’re just near the end of the course and running out of letters and $\beta$ is the canonical name for this exponent). As we approach along an isotherm,

$p-p_c\sim {(v-v_c)}^{\delta }\mathit{\hspace{1em}\hspace{1em} }\mathrm{with}\mathit{\hspace{1em}\hspace{1em} }\delta \approx 4.8$

Finally, as we approach $T_c$ from above, the compressibility scales as

$\kappa \sim {(T-T_c)}^{-\gamma }\mathit{\hspace{1em}\hspace{1em} }\mathrm{with}\mathit{\hspace{1em}\hspace{1em} }\gamma \approx 1.2$

The quantities $\beta$, $\gamma$ and $\delta$ are examples of critical exponents. We will see more of them shortly. The van der Waals equation provides only a crude first approximation to the critical exponents.

This is simplest to see in the grand canonical ensemble. Recall that back in Section 1 that we argued that $\mathrm{\Delta }N/N\sim 1/\sqrt{N}$, which allowed us to happily work in the grand canonical ensemble even when we actually had fixed particle number. In the context of the liquid-gas transition, fluctuating particle number is the same thing as fluctuating density $\rho =N/V$. Let’s revisit the calculation of $\mathrm{\Delta }N$ near the critical point. Using (1.45) and (1.48), the grand canonical partition function can be written as $\log 𝒵=\beta Vp(T,\mu )$, so the average particle number (1.42) is

$⟨N⟩={V{\displaystyle \frac{\partial p}{\partial \mu }}|}_{T,V}$

We already have an expression for the variance in the particle number in (1.43),

$\mathrm{\Delta }{N}^2={{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ⟨N⟩}{\partial \mu }}|}_{T,V}$

Dividing these two expressions, we have

${\displaystyle \frac{\mathrm{\Delta }{N}^2}{N}}={{{\displaystyle \frac{1}{V\beta }}{\displaystyle \frac{\partial ⟨N⟩}{\partial \mu }}|}_{T,V}{\displaystyle \frac{\partial \mu }{\partial p}}|}_{T,V}={{\displaystyle \frac{1}{V\beta }}{\displaystyle \frac{\partial ⟨N⟩}{\partial p}}|}_{T,V}$

But we can re-write this expression using the general relationship between partial derivatives ${{{\partial x/\partial y|}_z\partial y/\partial z|}_x\partial z/\partial x|}_y=-1$. We then have

${\displaystyle \frac{\mathrm{\Delta }{N}^2}{N}}=-{{{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ⟨N⟩}{\partial V}}|}_{p,T}{\displaystyle \frac{1}{V}}{\displaystyle \frac{\partial V}{\partial p}}|}_{N,T}$

This final expression relates the fluctuations in the particle number to the compressibility (5.190). But the compressibility is diverging at the critical point and this means that there are large fluctuations in the density of the fluid at this point. The result is that any simple equation of state, like the van der Waals equation, which works only with the average volume, pressure and density will miss this key aspect of the physics.

Understanding how to correctly account for these fluctuations is the subject of critical phenomena. It has close links with the renormalization group and conformal field theory which also arise in particle physics and string theory. You will meet some of these ideas in next year’s Statistical Field Theory course. Here we will turn to a different phase transition which will allow us to highlight some of the key ideas.

To see this, consider the same $d$-dimensional lattice as before, but now with particles hopping between lattice sites. These particles have hard cores, so no more than one can sit on a single lattice site. We introduce the variable $n_i\in \{0,1\}$ to specify whether a given lattice site, labelled by $i$, is empty ($n_i=0$) or filled ($n_i=1)$. We can also introduce an attractive force between atoms by offering them an energetic reward if they sit on neighbouring sites. The Hamiltonian of such a lattice gas is given by

$E=-4J{\displaystyle \sum _{⟨ij⟩}}{n}_i{n}_j-\mu {\displaystyle \sum _i}{n}_i$

where $\mu$ is the chemical potential which determines the overall particle number. But this Hamiltonian is trivially the same as the Ising model (5.192) if we make the identification

$s_i=2n_i-1\in \{-1,1\}$

The chemical potenial $\mu$ in the lattice gas plays the role of magnetic field in the spin system while the magnetization of the system (5.194) measures the average density of particles away from half-filling.

Here we’ll develop an approximate method to evaluate $Z$ known as mean field theory. We write the interactions between neighbouring spins in term of their deviation from the average spin $m$,

	$s_i{s}_j$	$=$	$[(s_i-m)+m][(s_j-m)+m]$
		$=$	$(s_i-m)(s_j-m)+m(s_j-m)+m(s_i-m)+m^2$

The mean field approximation means that we assume that the fluctuations of spins away from the average are small which allows us to neglect the first term above. Notice that this isn’t the statement that the variance of an individual spin is small; that can never be true because $s_i$ takes values $+1$ or $-1$ so $⟨s_i^2⟩=1$ and the variance $⟨{(s_i-m)}^2⟩$ is always large. Instead, the mean field approximation is a statement about fluctuations between spins on neighbouring sites, so the first term above can be neglected when summing over ${\sum }_{⟨ij⟩}$. We can then write the energy (5.192) as

	$E_{\mathrm{mf}}$	$=$	$-J{\displaystyle \sum _{⟨ij⟩}}[m(s_i+s_j)-m^2]-B{\displaystyle \sum _i}{s}_i$		(5.195)
		$=$	${\displaystyle \frac{1}{2}}JNq{m}^2-(Jqm+B){\displaystyle \sum _i}{s}_i$

where the factor of $Nq/2$ in the first term is simply the number of nearest neighbour pairs ${\sum }_{⟨ij⟩}$. The factor or $1/2$ is there because ${\sum }_{⟨ij⟩}$ is a sum over pairs rather than a sum of individual sites. (If you’re worried about this formula, you should check it for a simple square lattice in $d=1$ and $d=2$ dimensions). A similar factor in the second term cancelled the factor of $2$ due to $(s_i+s_j)$.

We see that the mean field approximation has removed the interactions. The Ising model reduces to the two state system that we saw way back in Section 1. The result of the interactions is that a given spin feels the average effect of its neighbour’s spins through a contribution to the effective magnetic field,

$B_{\mathrm{eff}}=B+Jqm$

In fact the low temperature behaviour requires slightly more explanation. For small values of $B$ and $T$, there are again three solutions to (5.197). This follows simply from continuity: there are three solutions for and $B=0$ shown in Figure 44 and these must survive in some neighbourhood of $B=0$. One of these solutions is again unstable. However, of the remaining two only one is now stable: that with $\mathrm{sign}(m)=\mathrm{sign}(B)$. The other is meta-stable. We will see why this is the case shortly when we come to discuss the free energy.

The net result of our discussion is depicted in the figures above. When $B=0$ there is a phase transition at $T=T_c$. For , the system can sit in one of two magnetized states with $m=\pm m_0$. In contrast, for $B\ne 0$, there is no phase transition as we vary temperature and the system has at all times a preferred magnetization whose sign is determined by that of $B$. Notice however, we do have a phase transition if we fix temperature at and vary $B$ from negative to positive. Then the magnetization jumps discontinuously from a negative value to a positive value. Since the magnetization is a first derivative of the free energy (5.194), this is a first order phase transition. In contrast, moving along the temperature axis at $B=0$ results in a second order phase transition at $T=T_c$.

We can calculate critical exponents for the Ising model. To compare with our discussion for the liquid-gas critical point, we will compute three quantities. First, consider the magnetization at $B=0$. We can ask how this magnetization decreases as we tend towards the critical point. Just below $T=T_c$, $m$ is small and we can Taylor expand (5.197) to get

$m\approx \beta Jqm-{\displaystyle \frac{1}{3}}{(\beta Jqm)}^3+\mathrm{\dots }$

The magnetization therefore scales as

$m_0\sim \pm {(T_c-T)}^{1/2}$

(5.199)

This is to be compared with the analogous result (5.188) from the van der Waals equation. We see that the values of the exponents are the same in both cases. Notice that the derivative $dm/dT$ becomes infinite as we approach the critical point. In fact, we had already anticipated this when we drew the plot of the magnetization in Figure 45.

Secondly, we can sit at $T=T_c$ and ask how the magnetization changes as we approach $B=0$. We can read this off from (5.197). At $T=T_c$ we have $\beta Jq=1$ and the consistency condition becomes $m=\tanh (B/Jq+m)$. Expanding for small $B$ gives

$m\approx {\displaystyle \frac{B}{Jq}}+m-{\displaystyle \frac{1}{3}}{\left({\displaystyle \frac{B}{Jq}}+m\right)}^3+\mathrm{\dots }\approx {\displaystyle \frac{B}{Jq}}+m-{\displaystyle \frac{1}{3}}{m}^3+𝒪(B^2)$

So we find that the magnetization scales as

$m\sim B^{1/3}$

(5.200)

Notice that this power of $1/3$ is again familiar from the liquid-gas transition (5.189) where the van der Waals equation gave $v_{\mathrm{gas}}-v_{\mathrm{liquid}}\sim {(p-p_c)}^{1/3}$.

Finally, we can look at the magnetic susceptibility $\chi$, defined as

$\chi ={N{\displaystyle \frac{\partial m}{\partial B}}|}_T$

This is analogous to the compressibility $\kappa$ of the gas. We will ask how $\chi$ changes as we approach $T\to T_c$ from above at $B=0$. We differentiate (5.197) with respect to $B$ to get

$\chi ={\displaystyle \frac{N\beta }{{\cosh }^2\beta Jqm}}\left(1+{\displaystyle \frac{Jq}{N}}\chi \right)$

We now evaluate this at $B=0$. Since we want to approach $T\to T_c$ from above, we can also set $m=0$ in the above expression. Evaluating this at $B=0$ gives us the scaling

$\chi ={\displaystyle \frac{N\beta }{1-Jq\beta }}\sim {(T-T_c)}^{-1}$

(5.201)

Once again, we see that same critical exponent that the van der Waals equation gave us for the gas (5.191).

There is actually a version of the Ising model for which the mean field theory is exact: it is the $d=\mathrm{\infty }$ dimensional lattice. This is unphysical (even for a string theorist). Roughly speaking, mean field theory works for large $d$ because each spin has a large number of neighbours and so indeed sees something close to the average spin.

But what about dimensions of interest? Mean field theory gets things most dramatically wrong in $d=1$. In that case, no phase transition occurs. We will derive this result below where we briefly describe the exact solution to the $d=1$ Ising model. There is a general lesson here: in low dimensions, both thermal and quantum fluctuations are more important and invariably stop systems forming ordered phases.

In higher dimensions, $d\ge 2$, the crude features of the phase diagram, including the existence of a phase transition, given by mean field theory are essentially correct. In fact, the very existence of a phase transition is already worthy of comment. The defining feature of a phase transition is behaviour that jumps discontinuously as we vary $\beta$ or $B$. Mathematically, the functions must be non-analytic. Yet all properties of the theory can be extracted from the partition function $Z$ which is a sum of smooth, analytic functions (5.193). How can we get a phase transition? The loophole is that $Z$ is only necessarily analytic if the sum is finite. But there is no such guarantee when the number of lattice sites $N\to \mathrm{\infty }$. We reach a similar conclusion to that of Bose-Einstein condensation: phase transitions only strictly happen in the thermodynamic limit. There are no phase transitions in finite systems.

What about the critical exponents that we computed in (5.199), (5.200) and (5.201)? It turns out that these are correct for the Ising model defined in $d\ge 4$. (We will briefly sketch why this is true at the end of this Chapter). But for $d=2$ and $d=3$, the critical exponents predicted by mean field theory are only first approximations to the true answers.

For $d=2$, the exact solution (which goes quite substantially past this course) gives the critical exponents to be,

	$m_0\sim {(T_c-T)}^{\beta }$	$\mathrm{ }\mathit{\hspace{1em}}\mathrm{with}$	$\beta ={\displaystyle \frac{1}{8}}$
	$m\sim B^{1/\delta }$	$\mathrm{with}$	$\delta =15$
	$\chi \sim {(T-T_c)}^{-\gamma }$	$\mathrm{with}$	$\gamma ={\displaystyle \frac{7}{4}}$

The biggest surprise is in $d=3$ dimensions. Here the critical exponents are not known exactly. However, there has been a great deal of numerical work to determine them. They are given by

$\beta \approx 0.32\mathit{\hspace{1em}\hspace{1em} },\delta \approx 4.8\mathit{\hspace{1em}\hspace{1em} },\gamma \approx 1.2$

But these are exactly the same critical exponents that are seen in the liquid-gas phase transition. That’s remarkable! We saw above that the mean field approach to the Ising model gave the same critical exponents as the van der Waals equation. But they are both wrong. And they are both wrong in the same, complicated, way! Why on earth would a system of spins on a lattice have anything to do with the phase transition between a liquid and gas? It is as if all memory of the microscopic physics — the type of particles, the nature of the interactions — has been lost at the critical point. And that’s exactly what happens.

What we’re seeing here is evidence for universality. There is a single theory which describes the physics at the critical point of the liquid gas transition, the 3d Ising model and many other systems. This is a theoretical physicist’s dream! We spend a great deal of time trying to throw away the messy details of a system to focus on the elegant essentials. But, at a critical point, Nature does this for us! Although critical points in two dimensions are well understood, there is still much that we don’t know about critical points in three dimensions. This, however, is a story that will have to wait for another day.

The energy of the system (5.192) can be trivially rewritten as

$E=-J{\displaystyle \sum _{i=1}^N}{s}_i{s}_{i+1}-{\displaystyle \frac{B}{2}}{\displaystyle \sum _{i=1}^N}(s_i+s_{i+1})$

We will impose periodic boundary conditions, so the spins live on a circular lattice with $s_{N+1}\equiv s_1$. The partition function is then

$Z={\displaystyle \sum _{s_1=\pm 1}}\mathrm{\dots }{\displaystyle \sum _{s_N=\pm 1}}{\displaystyle \prod _{i=1}^N}\exp \left(\beta J{s}_i{s}_{i+1}+{\displaystyle \frac{\beta B}{2}}(s_i+s_{i+1})\right)$

(5.202)

The crucial observation that allows us to solve the problem is that this partition function can be written as a product of matrices. We adopt notation from quantum mechanics and define the $2\times 2$ matrix,

$⟨s_i|T|s_{i+1}⟩\equiv \exp \left(\beta J{s}_i{s}_{i+1}+{\displaystyle \frac{\beta B}{2}}(s_i+s_{i+1})\right)$

(5.203)

The row of the matrix is specified by the value of $s_i=\pm 1$ and the column by $s_{i+1}=\pm 1$. $T$ is known as the transfer matrix and, in more conventional notation, is given by

The sums over the spins ${\sum }_{s_i}$ and product over lattice sites ${\prod }_i$ in (5.202) simply tell us to multiply the matrices defined in (5.203) and the partition function becomes

$Z=\mathrm{Tr}\left(⟨s_1|T|s_2⟩⟨s_2|T|s_3⟩\mathrm{\dots }⟨s_N|T|s_1⟩\right)=\mathrm{Tr}{T}^N$

(5.204)

where the trace arises because we have imposed periodic boundary conditions. To complete the story, we need only compute the eigenvalues of $T$ to determine the partition function. A quick calculation shows that the two eigenvalues of $T$ are

${\lambda }_{\pm }=e^{\beta J}\cosh \beta B\pm \sqrt{e^{2\beta J}{\cosh }^2\beta B-2\sinh 2\beta J}$

(5.205)

where, clearly, . The partition function is then

$Z={\lambda }_{+}^N+{\lambda }_{-}^N={\lambda }_{+}^N\left(1+{\displaystyle \frac{{\lambda }_{-}^N}{{\lambda }_{+}^N}}\right)\approx {\lambda }_{+}^N$

(5.206)

where, in the last step, we’ve used the simple fact that if ${\lambda }_{+}$ is the largest eigenvalue then ${\lambda }_{-}^N/{\lambda }_{+}^N\approx 0$ for very large $N$.

The partition function $Z$ contains many quantities of interest. In particular, we can use it to compute the magnetisation as a function of temperature when $B=0$. This, recall, is the quantity which is predicted to undergo a phase transition in the mean field approximation, going abruptly to zero at some critical temperature. In the $d=1$ Ising model, the magnetisation is given by

$m={{\displaystyle \frac{1}{N\beta }}{\displaystyle \frac{\partial \log Z}{\partial B}}|}_{B=0}={{\displaystyle \frac{1}{{\lambda }_{+}\beta }}{\displaystyle \frac{\partial {\lambda }_{+}}{\partial B}}|}_{B=0}=0$

We see that the true physics for $d=1$ is very different than that suggested by the mean field approximation. When $B=0$, there is no magnetisation! While the $J$ term in the energy encourages the spins to align, this is completely overwhelmed by thermal fluctuations for any value of the temperature.

There is a general lesson in this calculation: thermal fluctuations always win in one dimensional systems. They never exhibit ordered phases and, for this reason, never exhibit phase transitions. The mean field approximation is bad in one dimension.

The partition function is given by the sum over all states, weighted by $e^{-\beta E}$. At low temperatures, this is always dominated by the lowest lying states. For the Ising model, we have

$Z={\displaystyle \sum _{\{s_i\}}}\exp \left(\beta J{\displaystyle \sum _{⟨ij⟩}}{s}_i{s}_j\right)$

The low temperature limit is $\beta J\to \mathrm{\infty }$, where the partition function can be approximated by the sum over the first few lowest energy states. All we need to do is list these states.

The ground states are easy. There are two of them: spins all up or spins all down. For example, the ground state with spins all up looks like

Each of these ground states has energy $E=E_0=-2NJ$.

The first excited states arise by flipping a single spin. Each spin has $q=4$ nearest neighbours – denoted by red lines in the example below – each of which leads to an energy cost of $2J$. The energy of each first excited state is therefore $E_1=E_0+8J$.

There are, of course, $N$ different spins that we we can flip and, correspondingly, the first energy level has a degeneracy of $N$.

To proceed, we introduce a diagrammatic method to list the different states. We draw only the “broken” bonds which connect two spins with opposite orientation and, as in the diagram above, denote these by red lines. We further draw the flipped spins as red dots, the unflipped spins as blue dots. The energy of the state is determined simply by the number of red lines in the diagram. Pictorially, we write the first excited state as

$\mathit{\hspace{1em}\hspace{1em}\hspace{1em} }\begin{array}{c}{E}_1=E_0+8J\hfill \\ \mathrm{Degeneracy}=N\hfill \end{array}$

The next lowest state has six broken bonds. It takes the form

$\mathit{\hspace{1em}\hspace{1em}\hspace{1em} }\begin{array}{c}{E}_2=E_0+12J\hfill \\ \mathrm{Degeneracy}=2N\hfill \end{array}$

where the extra factor of 2 in the degeneracy comes from the two possible orientations (vertical and horizontal) of the graph.

Things are more interesting for the states which sit at the third excited level. These have 8 broken bonds. The simplest configuration consists of two, disconnected, flipped spins

$\mathit{\hspace{1em}\hspace{1em}\hspace{1em} }\begin{array}{c}{E}_3=E_0+16J\hfill \\ \mathrm{Degeneracy}=\frac{1}{2}N(N-5)\hfill \end{array}$

(5.207)

The factor of $N$ in the degeneracy comes from placing the first graph; the factor of $N-5$ arises because the flipped spin in the second graph can sit anywhere apart from on the five vertices used in the first graph. Finally, the factor of $1/2$ arises from the interchange of the two graphs.

There are also three further graphs with the same energy $E_3$. These are

$\mathit{\hspace{1em}\hspace{1em}\hspace{1em} }\begin{array}{c}{E}_3=E_0+16J\hfill \\ \mathrm{Degeneracy}=N\hfill \end{array}$

and

$\mathit{\hspace{1em}\hspace{1em}\hspace{1em} }\begin{array}{c}{E}_3=E_0+16J\hfill \\ \mathrm{Degeneracy}=2N\hfill \end{array}$

where the degeneracy comes from the two orientations (vertical and horizontal). And, finally,

$\mathit{\hspace{1em}\hspace{1em}\hspace{1em} }\begin{array}{c}{E}_3=E_0+16J\hfill \\ \mathrm{Degeneracy}=4N\hfill \end{array}$

where the degeneracy comes from the four orientations (rotating the graph by ${90}^{\circ }$).

Adding all the graphs above together gives us an expansion of the partition function in power of $e^{-\beta J}\ll 1$. This is

$Z=2e^{2N\beta J}\left(1+N{e}^{-8\beta J}+2N{e}^{-12\beta J}+{\displaystyle \frac{1}{2}}(N^2+9N)e^{-16\beta J}+\mathrm{\dots }\right)$

(5.208)

where the overall factor of 2 originates from the two ground states of the system. We’ll make use of the specific coefficients in this expansion in Section 5.3.4. Before we focus on the physics hiding in the low temperature expansion, it’s worth making a quick comment that something quite nice happens if we take the log of the partition function,

$\log Z=\log 2+2N\beta J+N{e}^{-8\beta J}+2N{e}^{-12\beta J}+{\displaystyle \frac{9}{2}}N{e}^{-16\beta J}+\mathrm{\dots }$

The thing to notice is that the $N^2$ term in the partition function (5.208) has cancelled out and $\log Z$ is proportional to $N$, which is to be expected since the free energy of the system is extensive. Looking back, we see that the $N^2$ term was associated to the disconnected diagrams in (5.207). There is actually a general lesson hiding here: the partition function can be written as the exponential of the sum of connected diagrams. We saw exactly the same issue arise in the cluster expansion in (2.85).

The fact that both energy and entropy scale with $L$ means that there is an interesting competition between them. At temperatures where the droplets are important, the partition function is schematically of the form

$Z\sim {\displaystyle \sum _L}{e}^{\alpha L}{e}^{-2\beta JL}$

For large $\beta$ (i.e. low temperature) the partition function converges. However, as the temperature increases, one reaches the critical temperature

$k_B{T}_c\approx {\displaystyle \frac{2J}{\alpha }}$

(5.209)

where the partition function no longer converges. At this point, the entropy wins over the energy cost and it is favourable to populate the system with droplets of arbitrary sizes. This is the how one sees the phase transition in the partition function. For temperature above $T_c$, the low-temperature expansion breaks down and the ordered magnetic phase is destroyed.

We can also use the droplet argument to see why phase transitions don’t occur in $d=1$ dimension. On a line, the boundary of any droplet always consists of just two points. This means that the energy cost to forming a droplet is always $E=2J$, regardless of the size of the droplet. But, since the droplet can exist anywhere along the line, its degeneracy is $N$. The net result is that the free energy associated to creating a droplet scales as

$F\sim 2J-k_{B}T\log N$

and, as $N\to \mathrm{\infty }$, the free energy is negative for any $T>0$. This means that the system will prefer to create droplets of arbitrary length, randomizing the spins. This is the intuitive reason why there is no magnetic ordered phase in the $d=1$ Ising model.

We again work with zero magnetic field, $B=0$ and write the partition function as

$Z={\displaystyle \sum _{\{s_i\}}}\exp \left(\beta J{\displaystyle \sum _{⟨ij⟩}}{s}_i{s}_j\right)={\displaystyle \sum _{\{s_i\}}}{\displaystyle \prod _{⟨ij⟩}}{e}^{\beta J{s}_i{s}_j}$

There is a useful way to rewrite $e^{\beta J{s}_i{s}_j}$ which relies on the fact that the product $s_i{s}_j$ only takes $\pm 1$. It doesn’t take long to check the following identity:

	$e^{\beta J{s}_i{s}_j}$	$=$	$\cosh \beta J+s_i{s}_j\sinh \beta J$
		$=$	$\cosh \beta J\left(1+s_i{s}_j\tanh \beta J\right)$

Using this, the partition function becomes

	$Z$	$=$	${\displaystyle \sum _{\{s_i\}}}{\displaystyle \prod _{⟨ij⟩}}\cosh \beta J\left(1+s_i{s}_j\tanh \beta J\right)$		(5.210)
		$=$	${(\cosh \beta J)}^{qN/2}{\displaystyle \sum _{\{s_i\}}}{\displaystyle \prod _{⟨ij⟩}}\left(1+s_i{s}_j\tanh \beta J\right)$

where the number of nearest neighbours is $q=4$ for the 2d square lattice.

With the partition function in this form, there is a natural expansion which suggests itself. At high temperatures $\beta J\ll 1$ which, of course, means that $\tanh \beta J\ll 1$. But the partition function is now naturally a product of powers of $\tanh \beta J$. This is somewhat analogous to the cluster expansion for the interacting gas that we met in Section 2.5.3. As in the cluster expansion, we will represent the expansion graphically.

We need no graphics for the leading order term. It has no factors of $\tanh \beta J$ and is simply

$Z\approx {(\cosh \beta J)}^{2N}{\displaystyle \sum _{\{s_i\}}}1=2^N{(\cosh \beta J)}^{2N}$

That’s simple.

Let’s now turn to the leading correction. Expanding the partition function (5.210), each power of $\tanh \beta J$ is associated to a nearest neighbour pair $⟨ij⟩$. We’ll represent this by drawing a line on the lattice:

$=s_i{s}_j\tanh \beta J$

But there’s a problem: each factor of $\tanh \beta J$ in (5.210) also comes with a sum over all spins $s_i$ and $s_j$. And these are $+1$ and $-1$ which means that they simply sum to zero,

${\displaystyle \sum _{s_i,s_j}}{s}_i{s}_j=+1-1-1+1=0$

How can we avoid this? The only way is to make sure that we’re summing over an even number of spins on each site, since then we get factors of $s_i^2=1$ and no cancellations. Graphically, this means that every site must have an even number of lines attached to it. The first correction is then of the form

$\mathrm{ }={(\tanh \beta J)}^4{\displaystyle \sum _{\{s_i\}}}{s}_1{s}_2{s}_2{s}_3{s}_3{s}_4{s}_4{s}_1=2^4{(\tanh \beta J)}^4$

There are $N$ such terms since the upper left corner of the square can be on any one of the $N$ lattice sites. (Assuming periodic boundary conditions for the lattice). So including the leading term and first correction, we have

$Z=2^N{(\cosh \beta J)}^{2N}\left(1+N{(\tanh \beta J)}^4+\mathrm{\dots }\right)$

We can go further. The next terms arise from graphs of length 6 and the only possibilities are rectangles, oriented as either landscape or portrait. Each of them can sit on one of $N$ sites, giving a contribution

$\mathrm{ }\mathit{\hspace{1em}}+\mathit{\hspace{1em}}=\mathrm{\hspace{0.25em}2}N{(\tanh \beta J)}^4$

Things get more interesting when we look at graphs of length 8. We have four different types of graphs. Firstly, there are the trivial, disconnected pair of squares

$\mathrm{ }\mathit{\hspace{1em}}=\mathit{\hspace{1em}}{\displaystyle \frac{1}{2}}N(N-5){(\tanh \beta J)}^8$

Here the first factor of $N$ is the possible positions of the first square; the factor of $N-5$ arises because the possible location of the upper corner of the second square can’t be on any of the vertices of the first, but nor can it be on the square one to the left of the upper corner of the first since that would give a graph that looks like which has three lines coming off the middle site and therefore vanishes when we sum over spins. Finally, the factor of $1/2$ comes because the two squares are identical.

The other graphs of length 8 are a large square, a rectangle and a corner. The large square gives a contribution

$\mathrm{ }\mathit{\hspace{1em}}=\mathit{\hspace{1em}}N{(\tanh \beta J)}^8$

There are two orientations for the rectangle. Including these gives a factor of 2,

$\mathrm{ }\mathit{\hspace{1em}}=\mathrm{\hspace{0.25em}\hspace{0.25em}2}N{(\tanh \beta J)}^8$

Finally, the corner graph has four orientations, giving

$\mathrm{ }\mathit{\hspace{1em}}=\mathrm{\hspace{0.25em}\hspace{0.25em}4}N{(\tanh \beta J)}^8$

Adding all contributions together gives us the first few terms in high temperature expansion of the partition function

	$Z=2^N{(\cosh \beta J)}^{2N}(1$	$+$	$N{(\tanh \beta J)}^4+2N{(\tanh \beta J)}^6$		(5.211)
			$+{\displaystyle \frac{1}{2}}(N^2+9N){(\tanh \beta J)}^8+\mathrm{\dots })$

There’s some magic hiding in this expansion which we’ll turn to in Section 5.3.4. First, let’s just see how the high energy expansion plays out in the $d=1$ dimensional Ising model.

We can express the statement of the duality more clearly. The Ising model at temperature $\beta$ is related to the same model at temperature $\tilde{\beta }$, defined as

$e^{-2\tilde{\beta }J}=\tanh \beta J$

(5.213)

This way of writing things hides the symmetry of the transformation. A little algebra shows that this is equivalent to

$\sinh 2\tilde{\beta }J={\displaystyle \frac{1}{\sinh 2\beta J}}$

Notice that this is a hot/cold duality. When $\beta J$ is large, $\tilde{\beta }J$ is small. Kramers-Wannier duality is the statement that, when $B=0$, the partition functions of the Ising model at two temperatures are related by

	$Z[\beta ]$	$=$	${\displaystyle \frac{2^N{(\cosh \beta J)}^{2N}}{2e^{2N\tilde{\beta }J}}}Z[\tilde{\beta }]$		(5.214)
		$=$	$2^{N-1}{(\cosh \beta J\sinh \beta J)}^{N}Z[\tilde{\beta }]$

This means that if you know the thermodynamics of the Ising model at one temperature, then you also know the thermodynamics at the other temperature. Notice however, that it does not say that all the physics of the two models is equivalent. In particular, when one system is in the ordered phase, the other typically lies in the disordered phase.

One immediate consequence of the duality is that we can use it to compute the exact critical temperature $T_c$. This is the temperature at which the partition function in singular in the $N\to \mathrm{\infty }$ limit. (We’ll discuss a more refined criterion in Section 5.4.3). If we further assume that there is just a single phase transition as we vary the temperature, then it must happen at the special self-dual point $\beta =\tilde{\beta }$. This is

$k_{B}T={\displaystyle \frac{2J}{\log (\sqrt{2}+1)}}\approx 2.269J$

The exact solution of Onsager confirms that this is indeed the transition temperature. It’s also worth noting that it’s fully consistent with the more heuristic Peierls droplet argument (5.209) since .