0% found this document useful (0 votes)
27 views37 pages

Ising Model Phase Transitions Explained

The document summarizes key aspects of the Ising model in statistical physics. It discusses: 1) The Ising model is one of the most important and analytically solvable models with interactions and a phase transition. It was originally developed to model ferromagnetism. 2) The 1D Ising model only exhibits a phase transition at 0 Kelvin, with the magnetization vanishing above 0 Kelvin. 3) The 2D Ising model exhibits a continuous phase transition at a finite critical temperature Tc, above which the magnetization smoothly goes to zero.

Uploaded by

amit
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
27 views37 pages

Ising Model Phase Transitions Explained

The document summarizes key aspects of the Ising model in statistical physics. It discusses: 1) The Ising model is one of the most important and analytically solvable models with interactions and a phase transition. It was originally developed to model ferromagnetism. 2) The 1D Ising model only exhibits a phase transition at 0 Kelvin, with the magnetization vanishing above 0 Kelvin. 3) The 2D Ising model exhibits a continuous phase transition at a finite critical temperature Tc, above which the magnetization smoothly goes to zero.

Uploaded by

amit
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

6 Phase transitions

6.1 Ising model


As the two most important examples for classical phase transitions, we now discuss
the Ising model and the real gas. In fact the Ising model is the most important model
in statistical physics, both historically and conceptually. It is one of the few analytically
solvable models with interactions and a phase transition. In 1920 it was invented by
Wilhelm Lenz as an exercise in ferromagnetism and given to his PhD student Ernst
Ising, who solved the one-dimensional case (publication Z. Phys. 31, 253Ð258, 1925),
which has a phase transition only at T = 0 (some textbooks therefore state that it has
no phase transition at all, which is true as a statement on finite temperature). In 1933
Rudolf Peierls used scaling arguments to show that the 2D version must have a phase
transition at finite temperature. In 1944 Lars Onsager solved the two-dimensional Ising
model in vanishing magnetic field with algebraic methods (publication Phys. Rev. 65,
117, 1944). Generations of theoretical physicists worked on this and related models,
leading to the developments of concepts such as universality, renormalization group
theory and symmetry breaking. The two-dimensional Ising model with magnetic field
and the three-dimensional Ising model are still not solved and this is one of the greatest
challenges in theoretical physics (like solving the Navier-Stokes equation). However,
it has been extensively studied with numerical methods (again like the Navier-Stokes
equation), so that we can say that in principle, we know everything about it. Thus the
Ising model is for statistical physics what the harmonic oscillator is for mechanics, the
hydrogen atom for quantum physics and the fruit fly for biology.
The Ising model is defined on a lattice of given connectivity and dimension, thus in
principle, there are infinitely many Ising models. Figure 6.1 shows some examples in
two dimensions. Usually however one studies it on a simple cubic lattice. The Ising
model is an example for a lattice spin model. More complicated lattice models of this
kind are the n-vector model, the XY-model, the Heisenberg model and the Potts-model.
In each lattice site i, we consider a spin Si with two possible states: Si = ±1 (cor-
responding to |↑i and |↓i). Nearest neighbours interact with an energy scale J. In
addition there might be an external magnetic field B giving a preference for one direc-
tion. Note that although we talk about spins and the cartoon show vectors, we really
have only a scalar degree of freedom, which is binary; this makes it the simplest model
possible, like the flip of a coin in probability theory.
The Hamiltonian of the Ising model reads

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:

1 The one-dimensional Ising model shows a phase transition only at T = 0 (Fig-


ure 6.2).
2 The two-dimensional Ising model shows a phase transition at finite temperature
Tc (Figure 6.3).
M

0
T

Figure 6.2: 1D Ising model: Magnetisation M as a function of temperature T. For T 6= 0


M vanishes and shows a jumping behaviour at T = 0.
M

0
T Tc

Figure 6.3: 2D Ising model: Magnetisation M as a function of temperature T. For values


T ≤ Tc M has a finite value.

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

Figure 6.4: 1D Ising model: ‘Ising chain’ of N spins.

Without periodic boundary conditions and considering the external field to vanish,
hence H = 0, ZN becomes:

ZN = ∑ ∑ ... ∑ eK(S1 S2 +S2 S3 +...+SN −1 SN )


S1 =±1 S2 =±1 S N =±1

= ∑ ... ∑ eK(S1 S2 +...+SN −2 SN −1 ) ∑ eKSN −1 SN


S1 =±1 S N −1 =±1 S N =±1
| {z }
=eK +e−K =2 cosh K
= ZN −1 2 cosh K
N 1
= (2 cosh K ) N −1 Z1 ≈ (2 cosh K ) N = ZN
|{z}
=2
Hence the free energy expression becomes, remembering K = βJ:
 
J
F = −k B TN ln 2 cosh
kB T
Because this is an analytical function for finite temperature, one already expects that no
phase transition takes place at finite T. We show this by considering spin correlations:
For each spin pair we introduce a different coupling constant Ki :
N −1 N −1
βH = −K ∑ Si Si + 1 → − ∑ K i Si Si + 1
i =1 i =1

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.

Because the system is homogeneous:

∀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.

Transfer matrix method


Next we investigate the case H 6= 0 using the transfer matrix method introduced by
Onsager for the 2D Ising model:

β H = − K ∑ Si S j − H ∑ Si
hiji i

We now use periodic boundary conditions: S N +1 = S1 (compare Figure 6.7). In the


thermodynamic limit N → ∞, boundaries become irrelevant.

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

⇒ Ti,i+1 = hSi | T |Si+1 i

⇒ 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.

We note that M( T 6= 0, B = 0) = 0 and no spontaneous magnetisation at finite T occurs.


For strong fields, hence
M( T 6= 0, B → ±∞) → ±µN
magnetisation saturates. For T → 0 M turns into a step function (compare Figure 6.8).

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.

From this we calculate the heat capacity in absence of a field:


K2

∂S
cB = T = k B (Figure 6.10)
∂T B=0 cosh2 K
Finally we consider the isothermal susceptibility (needs H):

1 ∂M βµ
χT = = ∂H M
N ∂B T N
βµ2 2
= ∂ ln ZN = βµ2 ∂2H ln λ1
N H
In the special case B = 0, χ T becomes:

βµ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.

divergence stronger than ∝ 1/T


χT

∝ 1/T

0
T

Figure 6.11: The susceptibility χ T as a function of temperature. χ T diverges for T → 0,


hence for T approaching the phase transition. In general this is a typical
signature of a phase transition. For large temperatures: χ T ∝ 1/T (law of
Curie).

We also note an interesting relation between susceptibility (a response function like


viscosity) and spin correlations (describing thermally activated ‘fluctuations’):

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

∑ ∑ Si Sj = χT ‘fluctuation dissipation theorem’ for the Ising model



=
N i =1 j =1

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

which is the same as above.




As −1 ≤ Si S j ≤ 1 this implies χ T can only diverge when

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’).

6.3 The Peierls argument


Starting around 1933, Peierls published scaling arguments why a phase transition should
occur in 2D as opposed to 1D. Here we report a few of these kinds of arguments to
demonstrate their spirit. Note that their validity also comes from the fact that thanks to
Onsager, we have an exact solution and thus can check back if they describe the core of
the problem or not.

Simple argument for 1D


We consider an Ising chain with all spins up and then select a few neighboring spins
and flip the whole island over. This creates two domain walls (also called grain boundaries
or defects) in the chain. The change in energy is

∆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.

More complex argument for 1D


We now look at an arbitrary number of domain walls, not only at one island with two
of them. We introduce the number of such domain walls M and write the free energy
in the domain wall picture:
 
N
F = 2J M − k B T ln
M

In the thermodynamic limit and with the Stirling formula we get

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.

Simple argument for 2D


We now want to make the simple argument for 2D rather than for 1D. We immediately
encounter the problem that now there are two processes we have to account for: where
to place the domain walls, and which shape to assign to them. With some intuition, we
anticipate that shape fluctuations are now more important than where the islands are
located. Thus we consider one island of down spins in a sea of up spins. The change in
energy is
∆E = L · 2J

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.

More complex argument for 2D


Another way to identify a phase transition is to investigate the effects of boundaries.
We consider a quadratic field of spins and fix all the ones at the boundary to point up.
We then consider the spin in the middle and ask if it keeps the up-preference of the
boundary in the TD-limit (p+ > 1/2 ?). One can show that for sufficiently low but
finite T, indeed this happens. This means that correlations do not decay completely
and that spontaneous magnetisation can emerge, indicating a phase transition.
We consider the quantity m = p+ − p− = 2p+ − 1, which will be finite if spontaneous
magnetization exists and vanish otherwise. We can write

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.

6.4 The 2D Ising model


Several methods of solution have been reported since Onsager’s original proof based on
the mapping to a quantum mechanical problem using transfer matrices (by graphical
solution summing over graphs, Grassmann algebra etc.). The exact solution for the free
energy is:
 Z 2π Z 2π 
1  2

F = −k B TN ln 2 + dq1 dq2 ln (1 − sinh 2K ) + sinh 2K (2 − cos q1 − cos q2 )
8π 2 0 0

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

Figure 6.12: c B as a function of temperature with a divergence at T = Tc .

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

Figure 6.13: M as a function of temperature.

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

1 H(λ = 0) = H0 reference case

2 H(λ = 1) = H case of interest


3 H1 = H − H0
 n o
⇒ − βF (λ) = ln ∑ e− βEj (λ) = ln tr e− βH(λ)
j

where ∑ j is the sum over all eigenvalues.

F (0) = F0
F (1) = F result of interest
n o
tr H e − β(H0 +λH1 )
dF 1
⇒ = = hH1 i (λ)
tr e− β(H0 +λH1 )


 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).

Application to the Ising model


We consider no external field:
H = −J ∑ Si S j
hi,ji

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

The Bogoliubov inequality then states

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 Kz < 1 ⇒ only intersection at x = 0

2 Kz > 1 ⇒ also two interactions at finite x 6= 0

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 ).

Computer simulations of the Ising model


Exact enumerations are not feasible for large systems. Consider that for a 50 x 50 2D
Ising system we would have already 22500 = 10753 states. Therefore it is better to use ‘im-
portance sampling’. Here only the relevant configurations are sampled. For the canonical
ensemble this means that we want to sample according to the Boltzmann distribution
1 − βEi
pi = e
Z
We now look for a procedure which effectively generates this distribution.
Usually such a procedure is implemented using random numbers (‘Monte Carlo simula-
tions’). A standard tool is the ‘Metropolis algorithm’. This algorithm generates a series of
configurations (‘Markov chain’) such that:
1. Any configuration can be reached in principle.

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:

1 Pick a spin i by random.

2 Calculate the energy change ∆E upon flipping the spin.

3 If ∆E < 0, accept the spin flip.

4 If ∆E > 0, accept the flip with probability e− β∆E .

By going downhill, the algorithm samples regions of high importance. By hopping


over barriers (compare 6.17), it allows to escape from metastable states. One can come
up with more complicated MC-moves than simply flipping spins. In MC-simulations
of the Ising model, it is very common to flip whole clusters. With this, one can get
very exact values, e.g. Tc = 0.221654 for the critical temperature of the 3D Ising
model (strangely
√ enough, this is close to the suggestion by Rosengreen 1986, Tc =
tanh−1 (( 5 − 2) cos(π/8)) = 0.221658).
For more details on Monte Carlo simulations, check the book by Kurt Binder and Dieter
Heermann (Springer 2010).

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’ :

H = −J ∑ ~Ji~Jj − µ~B ∑ ~Ji


hi,ji i

2 binary alloys and lattice gases:


Each lattice site is occupied either by an atom A or B. Nearest neighbour interac-
tions are e AA , eBB and e AB . We identify A with Si = 1 and B with Si = −1. The
Hamiltonian then is:
H = − ∑ Jij Si S j
hi,ji

with Jij = e AB − 21 (e AA + eBB ). Thus the Ising model describes order-disorder


transitions in regard to composition.

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).

5 neural networks representing the brain:


Si = 1 a synapse is firing, Si = −1 it is resting. The Hopfield model for neu-
ral networks is a dynamic version of the Ising model and Boltzmann machines
recognise handwriting by using the Ising model.

6 spread of opinions or diseases:


Spread of opinions, rumours or diseases in a society; these kinds of models are
used in socioeconomic physics. If nearest neighbour coupling is sufficiently strong,
the system gets ‘infected’.

140
6.5 Real gases
We consider the Hamiltonian of an ideal gas with N particles:
N p2
H= ∑ 2mi
i =1

The canonical partition sum in this case reads:


 N  N
1 
1
Z Z
2
d~p e− βp /(2m) 
 
Zid =  d~q
N! |h
 3   
{z } | {z }
=1/λ3 =V
1 VN
=
N! λ3N
h
with the thermal wavelength λ = 1 . Based on the partition sum we can calcu-
(2πmk B T ) 2
late the free energy:
⇒ Fid = −k B T ln Zid

∂Fid Nk B T
⇒ pid = − = = ρk B T
∂V T,N
V
The ideal gas is an appropriate description for diluted gases (small density ρ). At higher
densities, interactions become important.

We consider a pairwise additive and isotropic interaction potential U as the simplest


case:   
N
 p2
H= ∑  2m + ∑ U ||~ri {z
−~r j |
 
i =1 j <i }
:=rij

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).

An example for an isotropic potential is the ‘Lennard-Jones potential’ introduced by John


Lennard-Jones in 1924:   
σ 12  σ 6
U (r ) = 4e −
r r
The potential consists of two elements:

141
higher
density

isotropic orientational order

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

2 stability is provided by short-ranged ‘Born repulsion’ (∝ 1/r12 ).

σ 1.12σ r

Figure 6.19: The Lennard-Jones potential U as a function of inter-particle distance 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

In the high temperature limit z  1 we expand this expression in 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.

Obviously the virial expansion is similar to the expansion in fugacity z:

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


We note that the first order


hNi
z≈ = ρλ3
Z1
in our approximation reproduces our well known ideal gas result (p = ρk B T).

We need the next higher order to calculate B2 :

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).

Z d+δ   δd, βe1


⇒ B2 ( T ) = −2π r2 dr e βe − 1 ≈ −2πd2 δβe < 0 for e  k B T
d

145
U

d r

Figure 6.21: Potential for the hard spheres with excluded region r < d.

r

δ

Figure 6.22: Square well potential with range δ and depth e.

B2 ( T ) vanishes at high temperature. The attraction of the particles decreases the


pressure. This effect increases with increasing interaction energy e and range δ.

3 Hard core repulsion and attraction


As a course approximation to particle interactions with repulsive core and attrac-
tive well, we combine a hard hard sphere with an attractive well (Figure 6.23):

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´

Figure 6.24: B2 as a function of temperature. It vanishes at the ‘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.

the Gibbs free energy as we control temperature T and pressure p:

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.

We therefore have arrived at the following picture:


If several phases can exist at constant (T,p), the one with the lowest chemical poten-
tial µ is stable (lowest G/N). At the transition point, the chemical potentials are equal.
However, their slopes can have jumps (compare Figure 6.27).

μ(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.

G = µN, dG = −SdT + Vdp + µdN



∂µ 1 ∂G 1
⇒ = =− S
∂T p,N
N ∂T p,N
N

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).

If p, v and T are expressed in terms of their critical values:


p v e= T
pe = , ve = , T
pc vc Tc
the van der Waals equation becomes
 
3
pe + 2 v − 1) = 8 T
(3e e
v
e

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

The free volume vanishes at close packing.

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.

We first define the distribution function for absolute position:

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

r3 ...dr~N e− βW (~x1 ,x~2 ,...,r~N )


R
d~
= N ( N − 1)
r1 ...dr~N e− βW
R
d~

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 (~

We define the ‘radial distribution function’ g by:

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):

Interaction between probe at ~r and particles at {~ri }:


N
∑ U (~r − ~ri )
i =1

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.

incoming wave fluid outgoing wave

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 |.

Fermi’s golden rule for the transition rate I states:


D N E
I (~q) ∝ | ~k0 ∑ U (~r − ~ri ) ~k |2

i =1
∝ f (~q) · S (~q)
| {z } | {z }
‘form factor’ ‘structure factor’

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

You might also like