Ising Model Phase Transitions Explained
Ising Model Phase Transitions Explained
H = − J ∑ Si S j − Bµ ∑ Si
hiji i
119
Figure 6.1: Ising lattice examples in two dimensions: cubic (# neighbours z = 4, left),
triangular (z = 6, center) and hexagonal (honeycomb) (z = 3, right).
Here hiji indicates summation over nearest neighbours and µ is the magnetic moment
of a spin. In non-dimensional units we have
β H = − K ∑ Si S j − H ∑ Si
hiji i
where now both the coupling constant K = βJ and the external field H = βBµ depend
on temperature.
For J > 0 the symmetric configurations ↑↑ and ↓↓ are favorable and ↑↓ and ↓↑ are
unfavorable. Thus the system wants to avoid grain boundaries between regions with
up and down spins, at least at low temperature. At high temperature, grain boundaries
will proliferate because the correspond to a lot of entropy.
For B = 0 the system is invariant under Si ⇒ −Si . If B > 0, ↑-spins are favored. Using
the canonical formalism the partition sum for N spins reads:
ZN (K, H ) = ∑ ∑ ... ∑ e− βH
S1 =±1 S2 =±1 S N =±1
| {z }
2 N states
In practice one often uses periodic boundary conditions for lattice models to avoid
boundary effects, or finite size scaling to get rid of boundary effects by making the
system larger and larger.
Usually the Ising model is treated in magnetic language: J > 0 represents a ferromag-
netic interaction and B is a magnetic field. However, it can be used in many other
ways, e.g. in socioeconomic physics to present the spread of a opinion in a population
(spin up = believer, spin down = non-believer, you start to believe if your neighbors
are believers) or in biophysics to present the spread of a conformation in a molecular
ensembles (a molecule goes into another conformation if the neighbors have switched,
too).
In order to decide if the microscopic rules lead to a macroscopic change, one has to
introduce an order parameter. For a magnetic model like the Ising model, the natural
choice is the magnetisation: * +
N
M (K, H ) = µ ∑ Si
i =1
120
which is a measure for the averaged spin orientation. B > 0 will lead to M > 0. If for
B = 0 we find M 6= 0, then the system has spontaneously polarized itself (an example of
spontaneous symmetry breaking). In the following we will discuss two important results:
0
T
0
T Tc
If M changes in a smooth way at the transition (no jumps), we talk about a phase transi-
tion of second order or continuous phase transition. The 2D Ising model is the paradig-
matic case for such a transition at the critical temperature Tc . In the region around the
critical point, the system has very unusual properties (large fluctuations, universality,
critical slowing down).
121
6.2 The 1D Ising model
In one dimension the Ising model is an ‘Ising chain’ of spins (Figure 6.4). With periodic
boundary conditions this chain becomes a ring.
i
1 2 3 ... N
Without periodic boundary conditions and considering the external field to vanish,
hence H = 0, ZN becomes:
1
∑ ( Si Si + j ) e − β H
⇒ Si Si + j =
ZN { Si }
1 N −1
=
ZN ∑ ( Si Si + j ) e − ∑ i =1 K i Si Si + 1
{ Si } | {z }
= (Si Si+1 )(Si+1 Si+2 ) ... (Si+ j−1 Si+ j )
= Si Si + 1 Si + 1 Si + 2 . . . Si + j − 1 Si + j
| {z } | {z } | {z }
=1 =1 =1
= ∂Ki ∂Ki+1 ... ∂Ki+ j−1
122
ZN then can be calculated iteratively as above:
N −1
ZN = 2 N ∏ cosh Ki
i =1
j
cosh K1 ... sinh Ki ... sinh Ki+ j−1 ... cosh K N −1
= ∏ tanh Ki+k−1
⇒ Si Si + j =
cosh K1 ... cosh Ki ... cosh Ki+ j−1 ... cosh K N −1 k =1
Si Si+ j = (tanh K ) j
∀ i : Ki = K ⇒
The resulting spin correlations are shown in Figure 6.6. Despite the short-ranged inter-
action - we only consider nearest neighbour interactions - a longer ranged correlation
emerge from the statistical average, which decays exponentially with distance. Most
importantly, for T = 0 we have tanh K and the correlations do not decay at all.
1
0.8
0.6
0.4
tanh(x)
0.2
0
−0.2
−0.4
−0.6
−0.8
−1
−5 −3 −1 1 3 5
x
x →±∞
Figure 6.5: tanh( x ) as a function of x. For x > 0 tanh( x ) > 0; tanh( x ) → ±1.
∀i : h Si i = h S i
⇒ M = µN hSi
j→∞
Si Si + j → h Si i Si + j = h S i 2
(
µ2 N 2 T=0
M2 = µ2 N 2 lim Si Si+ j =
j→∞ 0 T>0
At finite T no spontaneous magnetisation occurs. At T = 0 we have a phase transition
(compare Figure 6.2). For T → 0 we have first made this limit and then the thermody-
namic limit N → ∞.
123
hSi , Si+j i = (tanh(x))j
1
0
0 1 2 3 4 5
j
j
Figure 6.6: Si , Si+ j = (tanh( x )) as a function of j. As can be seen in Figure 6.5 ,
tanh( x ) > 0 for x > 0. For the plot tanh( x ) was taken to be 0.5. For T = 0,
this curve would not decay.
β H = − K ∑ Si S j − H ∑ Si
hiji i
1
N 2
Figure 6.7: With periodic boundary conditions the one-dimensional Ising chain be-
comes a ring.
124
We define a ‘transfer function’ :
1
Ti,i+1 := eKSi Si+1 + 2 H (Si +Si+1 )
⇒ e− βH = T1,2 T2,3 ... TN,1
Each transfer function has four possible values which define a symmetric ‘transfer ma-
trix’:
e−K
K+ H
e
T=
e−K eK − H
In quantum mechanical notation:
1 0
| Si = + 1 i = | Si = − 1 i =
0 1
⇒ ZN = ∑ e− βH
{ Si }
= ∑ h S1 | T | S2 i h S2 | T | S3 i ... hS N | T |S1 i
{ Si }
= ∑ h S1 | T N | S1 i
S1 =±1
= TN + TN = tr T N
11 22
= λ1N + λ2N
We note that solving the Ising model amounts to an eigenvalue problem with λi being
the eigenvalues of T. This implies:
eK + H − λ e−K
det =0
e−K eK − H − λ
eK+ H − λ eK− H − λ − e−2K = 0
λ2 − 2eK cosh H λ + e2K − e−2K = 0
q
⇒ λ1,2 = eK cosh H ± e2K cosh2 H − 2 sinh 2K
q
= eK cosh H ± cosh2 H − 2e−2K sinh 2K
Thus we have arrived at an exact solution for the one dimension Ising model with
external field:
ZN = λ1N + λ2N
125
In the thermodynamic limit, only the larger eigenvalue λ1 is relevant:
N !
λ2 N →∞ N
ZN = λ1N 1 + → λ1
λ1
For H = 0 we get:
q
λ1 = e K + e2K − (e2K − e−2K ) = eK + e−K = 2 cosh K
⇒ ZN = (2 cosh K ) N for N 1
like before from the solution by direct summation (but different boundary conditions).
With the full solution we now can calculate any thermodynamic quantity of interest.
The thermal equation of state describes the magnetisation:
!
1
Z {∑
M( T, B) = µ ∑ Si e − β H
S}i i
µN
= µ∂ H ln ZN = ∂ H λ1
λ1
µN sinh H
=p
cosh2 H − 2e−2K sinh 2K
µN
T2 > T1
T1
M
−µN
0 B
Figure 6.8: The magnetisation M as a function of magnetic field B plotted for different
temperatures.
126
Next we calculate the entropy for B = 0:
F = − Nk B T ln (2 cosh K )
∂F
⇒ S=−= Nk B [ln (2 cosh K ) − K tanh K ] (Figure 6.9)
∂T
Considering the low and high temperature limits:
T →∞, K →0
S → Nk B ln 2
T →0, K →∞
S → Nk B (K − K ) = 0
where we recovered the third law of thermodynamics.
N kB ln(2)
S
0 T
Figure 6.9: The entropy S as a function of temperature T. For high temperatures S ap-
proaches S0 = Nk B ln 2 asymptotically.
βµ2
χT =
(1 − tanh K )
127
B
c
0 T
Figure 6.10: The heat capacity c B as a function of temperature shows a similar shape as
the one for the two state model (compare Fig. ??).
T →∞ 1
χT → law of Curie
T
In Figure 6.11 χ T is plotted as a function of temperature.
∝ 1/T
0
T
128
( ! )
1 ∂M βµ 1
Z {∑
χT = = ∂H µ ∑ S i e K ∑ Si S j + H ∑ Si
N ∂B T N Si } i
!
N N
βµ2 1
∑ ∑
N Z {S } i=1 j∑
= Si S j e − β H
i =1
!
βµ2 N N
For the one-dimensional Ising model and the limit N → ∞ the result becomes:
βµ2 ∞ 1
χT = N ∑ (tanh K ) j = βµ2
N j =0
1 − tanh K
1 N → ∞ (thermodynamical limit)
2 The range of the correlations must diverge, such that infinitively many terms give
a non-finite contribution. Therefore phase transitions are related to a divergence
of the ‘correlation length’. This implies that microscopic details become irrelevant
because the system becomes correlated on a macroscopic scale (‘critical fluctua-
tions’).
∆E = 2 · 2J
129
because there are two defects, each with an energy penalty 2J. The change in entropy
corresponds to the number of ways to choose the positions of the two defects:
N ( N − 1)
∆S = k B ln ≈ 2k B ln N
2
where we assume the number of lattice sites N 1 in the thermodynamic limit. Thus
the change in free energy reads
∆F = 4J − 2k B T ln N < 0
for any temperature T in the thermodynamic limit. This means that it is always favor-
able to create grain boundaries due to entropic reasons and a phase transition to order
cannot occur at finite temperature.
F
= 2Jx + k B T ( x ln x + (1 − x ) ln(1 − x ))
N
where x = M/N is the domain wall density. If we minimize F for x we get
1
xeq =
e2J/k B T + 1
thus at finite T there is always a finite domain wall density and correlations decay over a
finite distance. Moreover the system will not feel the effect of the boundary conditions.
Only at T = 0 we have xeq = 0, because then entropy does not matter.
130
where L is the contour length of the domain. A crude estimate for the number of pos-
sible shapes is 3 L , assuming a random walk on a 2D cubic lattice and neglecting inter-
sections and the fact that it has to close onto itself (at each lattice site, there are three
possibilities to proceed). Thus for entropy we have
∆S = k B ln 3 L .
Together we get
∆F = L(2J − k B T ln 3)
and thus ∆F < 0 only for T > Tc = 2J/(ln 3k B ) even in the thermodynamic limit
L → ∞. Thus this simple argument predicts that in 2D a phase transition can take place
at finite T, and the reason is a feature that is only present in two and higher dimensions,
namely shape fluctuations of the domain walls.
1 1 1
m=
Z ∑ e− βH − Z ∑ e− βH = Z ∑ e − β H (1 − Σ )
Σ+ Σ− Σ+
The first and second terms are sums over all configurations with a positive and nega-
tive central spin, respectively. The basic idea of the newly defined quantity Σ is that
each configuration with a positive central spin can be turned into one with a negative
central spin by flipping all spins in the surrounding positive domain. Importantly, the
difference in energy is simply 2Jl, where l is the length of the domain wall surrounding
this domain. Therefore one can write
∞
Σ= ∑ e−2Jβl = ∑ g(l )e−2Jβl
l =4
where the sum is now over all configurations which have been obtained by flipping. In
the second step we have rewritten the sum in terms of the length of the boundary. Here
g(l ) is the number of domains with length l. We note that the minimum l is 4 (one spin
flipped) and that one only will have even values (l = 4, 6, . . . ), because adding spins
one by one to the domain increases l by 2.
131
In order to prove the polarization, we have to show that Σ can be smaller than 1. We do
this by establishing an upper bound for g(l ):
l 1 l
g ( l ) < ( )2 · 4 · 3l −1 · = 3l
4 2l 24
The first term is the maximal area corresponding to the contour length l. The second
term is the number of possible paths starting from each point within this area: 4 for the
first step and 3 for each additional step (on a 2D simple cubic lattice). The last term
corrects for the fact that a path can go in two directions and can start at any point along
the contour of a boundary. We now transfer this into an upper bound for Σ:
∞ ∞
l 1 w4 (2 − w2 )
Σ< ∑ 24 wl = 24 ∑ (2n)w(2n) = 12(1 − w2 )2
l =4 n =2
where w = 3e−2βJ . We thus obtain Σ < 1 for w < wc = 0.87. This in turn translates into
a critical temperature
2J
Tc = = 1.6J/k B
k B ln(3/wc )
The exact result for the 2D Ising model is Tc = 2.269J/k B (see below). Thus the Peierls
argument does not only prove the transition, but even gives a reasonable first estimate
for its value. Note that here we have established only an upper bond for Σ. This does
not mean that Σ will be different from 1 above the critical temperature, we only showed
that it will certainly become smaller than this value at sufficiently low temperature. Our
argument is obviously very crude because we neglect interactions between boundary
loops, which will strongly bring down the number of possible paths.
This double integral cannot be reduced anymore. A phase transition occurs when the
logarithm diverges:
⇒ sinh 2Kc = 1
1 √
⇒ Kc = ln(1 + 2) ≈ 0.4407
2
√
⇒ Tc = 2J/ ln(1 + 2) ≈ 2.269J/k B
132
We define the ‘reduced temperature’:
T − Tc
e :=
Tc
and ‘critical exponents’ for the divergences (for B = 0) around Tc :
0
(
(−e)−α T < Tc
cB =
e−α T > Tc
(
(−e) β T < Tc
M=
0 T > Tc
From the exact solution one finds:
1 c B has a logarithmic divergence (Figure 6.12).
⇒ α = α0 = 0
18
2 M = 1 − sinh−4 2K (Figure 6.13)
⇒ β = 18
This result was announced by Onsager in 1948 at a conference, but never pub-
lished by himself.
B
c
Tc T
From the result for the magnetisation (which is the order parameter of the phase transi-
tion) one can construct the phase diagram. Figure 6.14 (left) shows the phase diagram
in the T-M-plane. Values for the magnetisation in the grey area (two-phase region) can-
not be realized in one system, because a self-polarized system jumps to the upper or
lower values of M. However, such a magnetisation can be realized by two systems, so
the system has to split into two. For example, M = 0 can be realized by two equally
133
M
Tc T
large systems with up and down magnetisation, respectively. Using the lever rule, each
desired value of M can be realized. Figure 6.14 (right) shows the phase diagram in the
T-B-plane. Now the two-phase region reduce to a line because any small external field
will immediately bias the system to up or down. Only for B = 0 phase coexistence can
occur.
M B
two-phase region
Tc T Tc T
Figure 6.14: Left: Phase diagram with excluded ‘two-phase region’ where the system
splits into two parts. Right: The two-phase region becomes a line in the
B( T ) diagram.
In order to understand the mechanisms underlying this phase transition, we now con-
sider the ‘mean field theory’ for the Ising model. This theory approximates a system of
interacting particles by a system of non-interacting particles. It can be made rigorous
by the ‘Gibbs-Bogoliubov-Feynman inequality’ and as such is a ‘perturbation theory’ (simi-
lar to the ‘Hartree-Fock approximation’ in quantum mechanics). In general, it is important
to have as many exactly solvable models in Statistical Physics as possible, even if they
might be physically not so realistic because they are built around some mathematical
134
trick to solve them. Nevertheless they can be very useful as starting points for pertur-
bative analyses.
Perturbation theory
We start from a model Hamiltonian H0 for which an exact solution is known:
H(λ) = H0 + λH1
F (0) = F0
F (1) = F result of interest
n o
tr H e − β(H0 +λH1 )
dF 1
⇒ = = hH1 i (λ)
tr e− β(H0 +λH1 )
dλ
n o n o 2
2
tr H1 e − β (H 0 + λ H 1 ) tr H1 e − β ( H0 + λH1 )
d2 F
= − β −
dλ2 tr e− β(H0 +λH1 ) tr e− β( H0 +λH1 )
D E
= − β H12 − hH1 i2 = − β (H1 − hH1 i)2 ≤ 0
dF
⇒ F ( λ ) ≤ F (0) + λ
dλ λ=0
λ =1
⇒ F ≤ Fu = F0 + hH1 i0 Bogoliubov inequality
A visualisation of the Bogoliubov inequality is sketched in Figure 6.15. Note that the
real F is everywhere concave, not only at λ = 0, so we can use λ = 1 without problems.
In order to optimize the approximation one minimizes the upper bound with respect to
the free model parameters. The modern master of this type of perturbation theory was
Richard Feynman.
135
λ
Figure 6.15: Sketch visualising the Bogoliubov inequality: F (λ) (solid line) ≤ F (0) +
dF
λ dλ (0) (dashed line).
However, we note that a spontaneous magnetization looks like there was an effective
magnetic field. We therefore choose as our unperturbed reference Hamiltonian
H 0 = − B ∑ Si
i
where we set µ = 1 for convenience and have introduced an effective magnetic field B.
For H = H0 we know the free energy expression:
F0 = − Nk B T ln e βB + e− βB
| {z }
=2 cosh βB
F ≤ F0 + hH − H0 i0
∑ + B ∑ hSi i0 = Fu
= − Nk B T ln (2 cosh( βB)) − J Si S j 0
hi,ji i
| {z } | {z }
N h S i0
N (z/2)hSi20
Here z is the number of nearest neighbours and we have to correct with a factor of 2 so
that we count each bond only once (compare Figure 6.1).
e βB − e− βB
h S i0 = = tanh βB
e βB + e− βB
136
We now fix B such that the upper bound Fu becomes minimal:
1 dFu d h S i0 d h S i0
0= = − hSi0 − Jz hSi0 + h S i0 + B
N dB dB dB
⇒ B = Jz hSi0 = Jz tanh βB
Note that a factor of 2 has canceled here. We note that our central result is a self-
consistent relation for the effective field B. We could have obtained this result directly
from a mean field reasoning, but it is more rigorous to derive it from the Bogoliubov
inequality.
1.2
0.8
0.3
tanh( x )
x
−0.2 Kz < 1
Kz
x
Kz > 1
Kz
−1
−0.7 +1
−1.2
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
X
Figure 6.16: M( x ) = tanh( x ) (blue) as a function of x. For Kz < 1 there is only one
x
intersection with g( x ) = Kz (red) at x = 0. For Kz > 1 there is also an
x
intersection with g( x ) = Kz (green) at finite x.
x
We define x = βB and have a look at the intersection of f ( x ) = tanh( x ) and g( x ) = Kz
(Figure 6.16). We note:
1 zJ
⇒ Kc = ⇒ Tc =
z kB
For the two-dimensional Ising model with cubic arrangement:
1
z=4 ⇒ Kc = = 0.25
4
137
Compare exact solution: Kc = 0.4407. Obviously the mean field theory is just a crude
approach because it predicts a phase transition in any dimension d. It becomes exact
for d → ∞.
How does magnetisation behave below Tc ? One can show that hSi = hSi0 using an
external field. Assuming a small magnetisation just below Tc , we can do a Taylor ex-
pansion:
1
hSi = tanh βB ≈ βB − ( βB)3
3
βB = zK hSi
1
⇒ hSi ≈ zK hSi − (zK hSi)3
3
12
1 − zK Tc − T
2 Kc :=1/z 1
⇒ hSi = (−3) ⇒ hSi = 3 2
zK Tc
1
We see that our approximative calculation yields a critical exponent β = 2 (compare
exact solution β = 81 ).
2. Averaging over all configurations in the Markov chain amount to doing the aver-
age with exact enumeration.
For the Ising model this means flipping one spin at random. We compare two configu-
rations i and j with:
pi
= e− β(Ei −Ej )
pj
We define pi→ j to be the ‘transition probability’ for one spin to go from state i to j.
⇒ ∑ pi → j = 1
j
138
We now require that locally we have detailed balance (should follow from time reversal
invariance):
pi → j pj
= = e− β(Ej −Ei )
p j →i pi
!
⇒ pi = ∑ pi → j pi = ∑ p j →i p j
j j
We note that pi is an eigenvector of the transition probability matrix. Thus a rule that
obeys detailed balance brings us to a steady state distribution { pi }. The simplest im-
plementation of this is the Metropolis algorithm:
i
Figure 6.17: Sketch visualising the Metropolis algorithm and how it recovers from local
minima.
139
Some applications of the Ising model
1 ferromagnetism:
The Ising model is the scalar version of the three-dimensional ‘Heisenberg model’ :
3 spin glasses:
now each bond is assigned an individual coupling constant J and they are drawn
from a random distribution. E.g. one can mix ferromagnetic and anti-ferromagnetic
couplings. This is an example for a structurally disordered system, on top of
which we can have a thermal order-disorder transition.
4 conformations in biomolecules:
a famous example is the helix-coil transition from biophysics. Si = 1 a hydrogen
bond in a DNA-molecules is closed; Si = −1 the bond is open. The phase tran-
sition is between a straight DNA-molecule (helix) and a coiled DNA-molecule.
Other examples are the oxygen-binding sites in hemoglobin, chemotactic recep-
tors in the receptor fields of bacteria, or the molecules building the bacterial flag-
ellum, which undergos a conformational switch if the flagellum is rotated in the
other direction (switch from run to tumble phases).
140
6.5 Real gases
We consider the Hamiltonian of an ideal gas with N particles:
N p2
H= ∑ 2mi
i =1
There is no dependence on
1. momenta (only positions)
2. relative orientations
of the particles.
An example for which the second assumption does not hold are liquid crystals (Fig-
ure 6.18).
141
higher
density
Figure 6.18: Liquid crystals: For increased density orientational, but not positional or-
der is established. This is the ‘isotropic-nematic transition’ of liquid crystals
that has been calculated by Lars Onsager in 1949.
1 a universal attraction between neutral atoms and molecules (‘van der Waals inter-
action’) proportional to 1/r6
σ 1.12σ r
-ε
For computer simulations one typically shifts and truncates the potential to achieve a
finite interaction range. These simulations can be done based on ‘Monte Carlo’ or ‘Molec-
ular Dynamics’ procedures.
solid
liquid
gas
Tc
T
Figure 6.20: A generic phase diagram typical for a simple one-component system. Tc in-
dicates the temperature of the ‘critical point’ where phase boundaries cease
to exist.
142
Figure 6.20 shows a phase diagram which is typical for a simple one-component sys-
tem (as for example described by the Lenard-Jones potential). We now return to the
analytical description .
N
1 1 1
Z Z
− βp2 /(2m)
Z= d~p e VN · N d N~q e− β ∑i< j U (rij )
N! h3 V
| {z } | {z }
= Zid := Zint
⇒ F = −k B T ln Z = Fid + Fint
∂F
p=− = pid + pint
∂V T,N
The interaction part does not factorise into single particle properties. Hence one needs
approximations. Because we understand the dilute case, we now introduce the ‘virial
expansion’, which is an expansion in low density around the ideal gas as a reference
system. We note that corrections to the ideal case pressure have to be of order ρ2 or
higher, because they arise if two particles or more collide.
∞
⇒ pint = k B T ∑ Bi ( T )ρi = k B TB2 ( T )ρ2 + O(ρ3 )
i =2
Here the Bi ( T ) are termed ‘virial coefficients’. In lowest order of the correction we thus
have
⇒ F = Nk B T ln ρλ3 − 1 + B2 ρ
p = ρk B T [1 + B2 ρ]
Calculation of B2 ( T )
In the following we use the grand canonical formalism to calculate B2 ( T ) from U (r ):
N
∞
ZG ( T, V, µ) = ∑ e βµ
Z ( T, V, N ) |{z}
N =0
| {z }
:= ZN fugacity z
⇒ ZG = Z0 + Z1 z + Z2 z2 + O z3
V
Z0 = 1, Z1 =
λ3
1 V4π
Z Z Z
− βU (|~
r1 −~
r2 |)
Z2 = d~
r1 d~
r2 e = dr r2 e− βU (r)
2!λ6 2λ6
143
Next we use the Euler relation for the grand canonical potential:
J = −k B T ln ZG = − pV
pV z 1
= ln ZG ≈ ln Z0 + Z1 z + Z2 z2
⇒
kB T
z 1 Z2
≈ Z1 z + Z2 z2 − 1 z2
2
x2
Were we used the approximation ln (1 + x ) ≈ x − 2 for x 1.
pV
= V ρ + B2 ρ2 + O ρ3
kB T
To make a comparison of coefficients we need the relation between z and ρ.
hNi ln z
ρ= , z = e βµ ⇒ µ =
V β
⇒ ∂µ = βz∂z
1
hNi = ∂µ ln ZG = z∂z ln ZG
β
≈ Z1 z + 2Z2 − Z12 z2
2Z2 − Z12
hNi
= z+ z2
Z1 Z1
|{z} | {z }
:=c := a
√
−1 + 1 + 4ac −1 + 1 + 12 (4ac) − 81 (4ac)2
⇒ z= ≈
2a 2a
= c(1 − ac)
h N i 2Z2 − Z12
hNi
= 1−
Z1 Z1 Z1
144
√
Here we used 1 + x ≈ 1 + 12 x − 18 x2 for x 1 in the first step.
pV
⇒ = ln ZG
kB T
= V ρ + B2 ρ2 + O ρ3 = h N i 1 + B2 ρ + O ρ2
1 2 2
= Z1 z + Z2 − Z1 z
2
h N i 2Z2 − Z12 Z12 h N i2
hNi
+ O ρ3
= Z1 1− + Z2 − 2
Z1 Z1 Z1 2 Z1
hNi 1
= h N i 1 + 2 −2Z2 − Z12 + Z2 − Z12 + O ρ2
Z1 2
2
Z hNi
= h N i 1 − Z2 − 1
2 Z12
Z2 1
⇒ B2 ( T ) = −V −
Z12 2
1
Z
=− d~r e− βU (r) − 1
2
Z
B2 ( T ) = −2π r2 dr e− βU (r) − 1
Higher virial coefficients follow in a systematic manner from the (graphical) cluster
expansion.
Examples
1 hard spheres
Spheres of radius d/2 which cannot penetrate each other. This yields an excluded
region of radius d (Figure 6.21).
Z d
2π 3 1
⇒ B2 ( T ) = −2π r2 (−1) dr = d = Vexcl = 4Vsphere > 0
0 3 2
Due to B2 being positive, a finite, excluded volume increases the pressure. B2 does
not depend on temperature, because there is no finite interaction energy.
2 square well
We consider a potential well of depth e between d and d + δ (Figure 6.22).
145
U
d r
Figure 6.21: Potential for the hard spheres with excluded region r < d.
r
-ε
δ
2π 3 a
⇒ B2 ( T ) = d − 2πd2 δβe = b −
3 kB T
with constants a, b > 0.
N
⇒ pV = Nk B T 1 + B2 ( T )
V
N2 N2
N Nk B T
= Nk B T 1 + b − a≈ N
− a
V V 1 − bV V
146
U
r
-ε
δ
Figure 6.23: Combination of the potentials for the hard spheres and the square well.
The resulting form is similar to the Lennard-Jones potential.
B2
B2(T) vanishes at T
`Boyle temperature´
V 1
Introducing the specific volume v = N = ρ this yields
kB T a
p= − van der Waals equation of state
v − b v2
The excluded volume (b) reduces the accessible volume for the particles while an
attractive interaction (a/v2 ) reduces pressure.
8a
For T < Tc = 27bk B , p(v) will have a minimum and maximum (see Figure 6.25).
dp
In the region between the minimum and maximum we have dv > 0. This implies
a local fluctuation to higher density (smaller v) to cause an increase in pressure
which then itself leads to a further increase in density. Due to this instability a part
of the system collapses and becomes a liquid. Hence we obtain a phase transition.
The details of the phase transition follow from the ‘Maxwell construction’. We consider
147
p
vmin vmax v
Figure 6.25: Pressure isotherm for a van der Waals gas below the critical temperature.
In the region between vmin and vmax the system is unstable.
G = E − TS + pV = µN
| {z }
:= F
F
⇒ µ= + p |{z}
v
N
|{z} =V/N
:= f
For two coexisting phases L and G in equilibrium the intensive parameters T, p and µ
have to be the same:
µ L ( T, p) = µG ( T, p)
⇒ f G − f L = pt (v L − v G )
Here pt is the transition pressure.
p
superheated liquid
undercooled gas
instable
vL vG v
Figure 6.26: Van der Waals isotherm with Maxwell construction based on the equality
of areas 1 and 2.
The left hand side can be calculated by integration along the isotherm:
Z vG Z vG
∂ f ( T, v)
fG − fL = dv =− dv p( T, v)
vL ∂v T vL
148
Z vL
⇒ p T (v L − vG ) = dv p( T, V )
vG
Geometrically this means that in Figure 6.26 the dotted area has to equal the one below
the solid line. Hence pt can be determined based on the equality of areas 1 and 2.
μ(T,p) v
μG
vG
μL
vL
pt
p pt
p
Figure 6.27: Left: The chemical potential µ as a function of pressure for phases G and L.
At the transition point µG = µ L , but the slopes have jumps.
∂µ
Right: The specific volume v = ∂p as a function of pressure has a jump
T
at the transition pressure pt .
In order to bring the fluid from liquid to gas, we need the ‘heat of evaporation’ or ‘latent
heat’ Q: Z Tt+ Z Tt+
Q= T dS = dH = HG − HL
Tt− Tt−
where we used dH = TdS + Vdp and p = pt = const.
H E + pV vdW eq a
⇒ h= = = e( T ) − + pv
N N |{z} v
kinetic energy contr.
Q a a a
q= = hG − h L = − + p (vG − v L ) ≈ + pvG (vG v L )
N vL vG vL
a
vL is the energy required to overcome attraction while pvG is the energy required for
expansion.
149
μ(T,p) S
G
Q/Tt
μL
L
μG
T T
Tt Tt
Figure 6.28: Left: The chemical potential µ as a function of temperature for phases G
and L. At the transition point µG = µ L , but the slopes have jumps.
Right: The entropy as a function of temperature jumps at the transition
point.
∂µ ∂µ
We conclude that both v = ∂p T and s = − ∂T p jump at the transition (compare Fig-
ures 6.27 and 6.28 respectively). Therefore this phase transition is called to be of ‘first
order’ or ‘discontinuous’. Both jumps disappear at the critical point, where the isotherm
becomes horizontal at the transition. From
∂2 p
∂p
= =0
∂v T ∂v2 T
one calculates the critical values:
8a a
vc = 3b, Tc = , pc =
27bk B 27b2
for water: pc = 217 atm, Tc = 647 K
pc vc 3
⇒ = = 0.375 independent of a and b
k B Tc 8
Experimental values are similar, but slightly smaller (around 0.3).
Figure 6.29 shows van der Waals isotherms for different temperatures with respect to
Tc .
150
3
2.5
T = 1.2 Tc
binodal
T = Tc
2
T = 0.95 Tc
1.5 T = 0.85 Tc
p in units of pc
spinodal
T = 0.75 Tc
1
0.5
−0.5
−1
−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4
log(v) (v in units of vc)
Figure 6.29: Van der Waals isotherms for different temperatures. The ‘spinodal’ is the
boundary between metastable and unstable states. The ‘binodal’ separates
metastable and absolutely stable states. The latter curve was calculated
numerically based on Maxwell’s construction for different temperatures.
This
reduced
equation leads to the ‘law of corresponding states’: Two fluids with the same
pe, ve, T
e are in equivalent states. Indeed experimental curves show surprisingly good
data collapse. Even more surprisingly, their behaviour becomes almost identical at the
critical point - large fluctuations render microscopic details irrelevant.
The ‘van der Waals’ equation of states predicts the fluid-fluid phase transition caused
by attractive interactions. The fluid-solid phase transition can be predicted by a simple
entropic argument. Recall the van der Waals theory for a hard sphere fluid:
Nλ3
F = Nk B T ln −1
V − Nb
∂F Nk B T
⇒ p=− =
∂V T,N V − Nb
b = 4Vs ⇒ V − Nb = V (1 − ρb) = αV
with α = 1 − ρ/ρ0 and ρ0 = 1/b. αV is the free volume in the fluid.
Based on Figure 6.30 and L = V 1/3 we define the free volume of a solid as:
3
" 13 #3
1 ρ
αV ≈ V − d 3 = 1− V
ρ0
151
L
Figure 6.30: Unit cell with hard spheres of diameter d. The grey shaded region indicates
the free volume.
α2
⇒ F1 − F2 = Nk B T ln
α1
Hence the phase with larger α is favored. For that reason the fluid F and the solid S
are stable at low and high densities, respectively. Figure 6.31 shows how the Maxwell
construction looks like in this case.
F
F S
Maxwell
construction
2p
ρ
Figure 6.31: Free energy for liquid (F) and solid phase(S) as a function of density ρ. The
tangent represents the Maxwell construction.
We now can understand the complete phase diagram of a simple one-component fluid:
Distribution functions
In contrast to thermodynamics, statistical physics not only predict phase behaviour, but
also structure. The central quantity in both theory and experiments is the ‘radial distri-
bution function’.
152
a) b)
T T
F S F S
Tc
G L
Tt
ρ ρ
c) d)
ρ p
S
L F
Tt Tc T Tt Tc T
Figure 6.32: Combining the two transitions in (a), one gets the complete phase diagram
in (b). In (c) we swap T and ρ axes. By replacing ρ by p, we get the final
phase diagram in (d). Two-phase coexistence regions become lines in this
representation.
e− β ∑i< j U (rij )
r1 , ..., r~N ) = R
p(~
r1 ...dr~N e− β ∑i< j U (rij )
d~
By defining W := ∑i< j U (rij ), the probability that any particle is at position ~r can be
written as:
N
n1 (~x ) = ∑ hδ (~x − ~rk )i
k =1
r1 ...dr~N ∑kN=1 δ (~x − ~ rk ) e− βW
R
d~
=
r1 ...dr~N e− βW
R
d~
r2 ...dr~N e− βW (~x,r~2 ,...,r~N )
d~
=N
r1 ...dr~N e− βW
R
d~
153
ideal gas N
= } =ρ
n1 (~x ) | {z
V
W =0
The probability that some particle is at x~1 and another at x~2 is:
x1 , x~2 ) = ∑ δ (~
n2 (~ x1 − ~ri ) δ x~2 − ~r j
i6= j
ideal gas N ( N − 1) N → ∞ 2
x1 , x~2 ) | {z
n2 (~ = } → ρ
V2
W =0
For the pairwise additive potential everything follows from n1 and n2 . Eg the averaged
interaction energy:
hW i = ∑ U (rij )
i< j
1
Z
2 i∑
= x1 − x~2 ) δ (~
dx~1 dx~2 U (~ x1 − ~ri ) δ x~2 − ~r j
6= j
1
Z
= x1 − x~2 ) n2 (~
dx~1 dx~2 U (~ x1 , x~2 )
2
In a homogeneous system:
x1 − x~2 |)
x1 , x~2 ) = n2 (|~
n2 (~
x1 − x~2 |) = ρ2 g (|~
n2 (|~ x1 − x~2 |)
N2
Z
⇒ hW i = d~r U (r ) g(r )
2V
ρg(r )4πr2 dr is the average number of particles in a spherical shell of width dr at a
distance r from any particle.
While for the ideal gas g = 1, g(r ) has damped oscillations for a real gas (compare
Figure 6.33).
The pair correlation function g(r ) can be measured in scattering experiments (x-rays,
neutrons, electrons, light):
154
g(r)
r
Figure 6.33: The correlation function g for a real gas as a function of distance r. The
oscillations damp away with distance.
Figure 6.34: Schematic sketch of a scattering experiment. If~k denotes the wave vector of
the incoming wave and ~k0 for the outgoing wave, the wave-fluid interaction
results in a momentum transfer q = ~k0 −~k with |~k0 | = |~k |.
The ‘form factor’ describes the interaction between probe and particles and depends on
the experiment.
f (~q) = |U (~q)|2
Z
U (~q) = d~r U (~r ) e−i~q·~r
The ‘structure factor’ represents the internal structure and is independent of the type of
probe. * +
1
Z
N i∑
i~q(~ri −~r j )
S (~q) = e = 1 + ρ d~r ( g(r ) − 1) ei~q·~r
6= j
We note that S (q) is essentially the Fourier transform of g(r ) (qualitative shape similar
to g(r ) - compare Figure 6.33).
155