Chapter 3

The Ising Model

3.1 Introduction
3.1.1 Definition of the Ising model
Given some regular lattice in d dimensions (d=1,2,3,...)
(e.g., in 3 dimensions, cubic, fcc, bcc, ...)

Ingredients:

(i) Each lattice site carries a ”Spin” Si “ ˘1, i.e., a variable that can
take one out of two values (not a quantum mechanical spin)
(ii) Cooperativity: The value of one spin influences the neighbor spins
(iii) Possibly an external ”field” H that favors a certain value of Si
No kinetic degrees of freedom (momentum etc.)
ÿ ÿ
ñ Energy function H rtSi us “ ´J Si Sj ´ H Si
xijy i
(”Hamiltonian”) nearest neighbors
Cooperativity (ii) Field (iii)

The most common choice is J ą 0 (”ferromagnetic”),

but J ă 0 is also possible (”antiferromagnetic”).

Based on the energy function, one calculates

– The partition function Z “ tSi “˘1u e´βH rtSi us
ř

– The free energy F “ ´kB T ln Z

etc.

Extensions and generalizations are possible and have been studied,

e.g., additional interactions between spins that are further apart,
anisotropic or spatially varying interactions Jij or fields Hi ,
Ising models on irregular lattices or other graphs, etc.

19
20 CHAPTER 3. THE ISING MODEL

3.1.2 Motivation
‚ Original motivation: Designed as simple model for magnetism, therefore
”magnetic” language (spins, field, etc.)
However, the Ising model is actually not a good model for a magnetic
system, since magnetic moments (”Spins”) are quantum objects and three
~
dimensional vectors (S)

‚ Can be a good model for certain binary alloys in the context of order-disorder
transitions (a ”spin” S “ ˘1 then indicates the occupation of a lattice
site with an atom of type A or B).

‚ The Ising model at H “ 0 is one of the simplest model systems that exhibits
a phase transition
; used to study fundamental properties of phase transitions

‚ Universality: Many practically important phase transitions are in the so-

called ”Ising universality class”, e.g., the liquid-gas transition, binary
mixtures, ... (Ñ see sections ?? and ?? )

(‚ Historically: Has been studied extensively. Many techniques and arguments

that were developed in this context are simply fun!
First case, where a non-trivial phase transition has been calculated ex-
actly from first principles, i.e., starting from the partition function of a
microscopic model)

3.1.3 History
1925: Introduction by Ising (in his PhD thesis, supervised by Lenz)
Model for magnetism
Exact solution in one dimension, unfortunately no phase transition :-(
(Ising believed/argued, that there would be no phase transition in higher
dimensions either)

1934, 35: Bragg, Williams, Bethe

Took interest in Ising model as model for binary mixtures
Ñ Developed approximate solution methods (”mean-field approximations”)
that gave a phase transition

1936: Bethe
Argument, why there should be a phase transition in 2 dimensions
(not entirely accurate, later completed in 1964 by Griffiths)

1941: Kramers, Wannier

Symmetry considerations Ñ exact expression for Tc in two dimensions

1944: Onsager
Exact solution in two dimensions for the case H “ 0
(Solutions for H ‰ 0 also became available later.)
3.1. INTRODUCTION 21

3.1.4 Remark: Identifying phase transitions in the Ising model

ř
Problem: For symmetry reasons, one always has M “ x i Si y “ 0 at H “ 0.
Ñ How can one identify a phase transition then?

Possible approaches:

(i) Calculate free energy F and search for singularities

(ii) Introduce a symmetry breaking infinitesimal field H Ñ 0˘
Symmetry breaking phase transition: lim M pT, Hq ‰ lim M pT, Hq
HÑ0` HÑ0´
(iii) Symmetry breaking boundary conditions
++++++ −−−−−−
+ + − −
+ + Ñ M0 pT q Ñ ´M0 pT q
+ ++
−
− − −
−
+ + − −
++++++ −−−−−−

(iv) Analyze histogram of M

P(M) P(M)

(disordered) (ordered)
M 0 =0 −M0 M0
22 CHAPTER 3. THE ISING MODEL

3.2 One-Dimensional Ising Model

We first consider the one dimensional Ising chain:
Nř´1
+ - - + - + + - - H “ ´J Si Si`1
1 N i“1

3.2.1 Why there cannot be a phase transition

‚ ”Ground state” (State with lowest energy)
All Spins have the same sign, e.g., + + + + + + +
Energy: E “ ´JpN ´ 1q “: E0

‚ Lowest excitation:
One ”kink”: + + + +|- - -
Energy: E “ ´JpN ´ 2q ` J “ E0 ` 2J
; Energy costs compared to ground state: ∆E “ 2J
; Boltzmann probability of such a kink
Pkink 9 e´β∆E “ e´2βJ : Finite number (0 ă Pkink ă 1).
But: Every kink destroys the order
Probability, that no kink is present (i.e., order persists):
N Ñ8
Pno kink 9 p1 ´ Pkink qN ´1 ÝÑ 0
Ó
possible positions of kinks

ñ Ising chain is always disordered !

NB: Argument does not work in two dimensions, +++++ +++++

+++++ +++++
since lowest excitation (one flipped spin) +++++ +−+++
+++++ +++++
does not yet destroy global order +++++ +++++
Global order is destroyed by an excitation of the form
+++++ ++−−−
+++++ ++−−−
+++++ ++−−− L (boundary line between domains)
+++++ +++−−
+++++ +++−−

However, the energy costs of this excitation are: ∆E ě 2JL

ñ P 9 e´2βLJ Ñ 0 for L Ñ 8
3.2. ONE-DIMENSIONAL ISING MODEL 23

3.2.2 Solution of the one-dimensional Ising model

N
Consider one dimensional Ising chain H “ ´J
ř ř
Si Si`1 ´ H i Si
with periodic boundary conditions: SN `1 :“ S1 i“1

‹ Free energy: Exact calculation via transfer matrix method

Starting point: Partition function Z “ tSi u e´βH rtSi us
ř

Notation (motivated by quantum

` ˘ mechanical` bras
˘ and kets):
S fl |Sy with S “ `1 fl 10 , S “ ´1 fl 01
H H
eβJS1 S2 `β 2 pS1 `S2 q looooooooooomooooooooooon
ñ e´βH “ looooooooooomooooooooooon eβJS2 S3 `β 2 pS2 `S3 q ¨ ¨ ¨
fl ˆ xS1 |V |S2 y ˙ xS2 |V |S3 y ¨ ¨ ¨
eβJ`βH e´βJ
with V “ (check by inserting!)
e´βJ eβJ´βH
ñZ “ xS1 |V |S2 y¨¨ xSN |V |S1 y “ xS1 |V N |S1 y “ TrpV N q
ř ř
S1¨¨ SN S1
N N N N N Ñ8
N
“ pλ 1 ` λ2 q “ λ1 p1 ` pλ2 {λ1 q q ÝÑ λ1
looooomooooon
Eigenvalues, λ1 ąλ2
ñ F “ ´kB T ln Z “ ´kB T N ln λ1
Specifically: Eigenvalues of a
V are given by
λ1,2 “ eβJ coshpβHq ˘ e2βJ sinh2 pβHq ` e´2βJ
(at H “ 0: λ1,2 “ eβJ ˘ e´βJ )
” b ı
ñ F “ ´kB T N ln eβJ coshpβHq ` e2βJ sinh2 pβHq ` e´2βJ
“ ‰
In particular, at H “ 0: F “ ´kB T N ln 2 coshpβJq
Analytical function for all temperatures ; No phase transition!
ř BF
‹ ”Magnetization” : M “ x i Si y “ ´ BH
ř ř
(Since: F “ ´kB T ln e´βH “ ´kMB T lnp eβJ Si Sj `βH Si
ř ř
BF
ř ´βH ř ř ´βH ř
ñ BH “ ´k BT
///// e / Si
β e “ ´x Si y X q

Mb
ñ M “ ¨ ¨ ¨ “ N sinhpβHq cosh2 pβHq ´ 1 ` e´2βJ
M/N

βJ Ñ 0 : M {N « tanhpβHq H

βJ Ñ 8 : M {N « ? 2sinhpβHq M/N
sinh pβHq“signpβJq
H
step function
; ”Phase transition” at T “ 0?

‹ Correlations: Gij “ xSi Sj y ´ xSi yxSj y

Specifically H “ 0 ; xSi y “ xSj y “ 0

Longer calculation (homework) ñ Gij 9 e´|i´j|{ξ

1
with ξpH “ 0, T q “ ´ lnptanhpβJqq „ e2βJ Ñ 8 for β Ñ 8
; Correlation length ξ diverges at T Ñ 0
; Characteristic feature of critical behavior
24 CHAPTER 3. THE ISING MODEL

3.3 Two-Dimensional Ising Model: Exact Results

Now we discuss the two dimension Ising model at H “ 0,
H “ ´J
ř
Si Sj (with xijy: nearest neighbors)
xijy
In this section we will restrict ourselves to the square lattice. (Generalizations of
the arguments below to other lattices are sometimes possible, but not always.)

3.3.1 Peierl’s argument for the existence of a phase transition

(Peierls 1936, Griffiths 1964, Dobrushin 1965)

Idea of the argument:

+++++
– Consider Ising system Ω
p +++++ fixed
L ++−−+
with fixed boundary condition: ++−++
+++++
Surrounded by spins S “ `1.
– Show, that for sufficiently low temperatures, there exists an α ą 0
1 ř
such that, at H “ 0, the magnetization per spin MΩp {N “ N x i Si y
is always larger than α, MΩp {N ě α, independent of system size
(i.e., this system has positive magnetization).
– At these temperatures, limN Ñ8,HÑ0` M {N ě α holds generally, inde-
pendent of the boundary condition. (H Ñ 0` is introduced to break
the S “ ˘1 symmetry in the absence of the boundary condition.)
(Reason: Free energy F pT, Hq ´ ¯
Magnetization at H Ñ 0` : M “ ´ BF
BH
HÑ0`
But: F is extensive (additive) at L Ñ 8 ?
; Boundary only contributes with surface term to F , Fsurf 9 L 9 N
; Fgeneral “ FΩp ` OpLq ñ lim MΩp “ lim Mgeneral ` OpLq
HÑ0` HÑ0`

M
Goal therefore: Search for lower bound for NΩp “ 1 ´ 2 xNN´ y in the system Ω
p
where N´ : Number of sites with Si “ ´1

‹ Consider some (arbitrary) configuration C

+++++++
+−−−−++ Draw domain wall lines between spins of different sign
+−+−−++
+−−−−−+ ; Lines form closed polygons
+++−−++ ; Every spin S “ ´1 lies inside at least one polygon
+−+++−+
+++++++ Label all
" possible polygons (Ñ label p)
1 : polygon p appears in configuration C
Define variables Xp “
0 : otherwise
lp :“ contour length of polygon p
Polygon p contains at most plp {4q2 spins
ř ` l ˘2
ñ Estimate: N´ ď p Xp 4p
ř ` l ˘2
‹ Also holds in the thermal average ñ xN´ y ď p xXp y 4p
ř1 ´βH
e configurations C that contain polygon p
with xXp y “ ř ´βH
e all configurations
3.3. TWO-DIMENSIONAL ISING MODEL: EXACT RESULTS 25

‹ Upper bound for xXp y and xN´ y

+++++++ +++++++
Trick: For each configuration C (in Σ1 ) +−−−−++ +++++++
+−+−−++ ++−++++
that contains polygon p, construct a +−−−−−+ +++++++
+++−−++ +++++++
partner configuration C ˚ by invert- +−+++−+ +−+++−+
+++++++ +++++++
ing all spins inside the polygon.
; Energy difference: H pC ˚ q “ H pC q ´ 2lp J
“ e2lp βJ e´βH “ e2lp βJ 1 e´βH
ř ´βH ř ´βH ř ř
ñ Estimate : e ě e
all confs C˚ C
ř1 ´βH
ñ xXp y “ ř e´βH ď ´2l
e p βJ
e number of polygons of length l
hkkikkj
βJ lp 2
ř ` ˘ ř ` l ˘2
ñ xN´ y ď p e´2lp
4 “ l 4 e´2lβJ
nplq
‹ Upper bound for nplq (number of polygons of length l)
Construction of a polygon: *
First line: 2N possibilities (+ border)
ñ nplq ď 2N 3l´1
Second line: At most 3 possibilities
‹ Combine everything:
` ˘2 ` ´2βJ ˘l 2
xN´ y ď l e´2lβJ 4l 2N 3l´1 “ 24 N ř8 N
ř
l“0 3e l “ 24 f p3e´2βJ q
xp1`xq ř l2 d 2 ř ˇ
α lˇ xp1`xq
with f pxq “ p1´xq3 (since: l x l “ dα 2 l pxe q α“0 “ ¨ ¨ ¨ “ p1´xq3 )

f pxq becomes arbitrarily small at x Ñ 0:

For example, choose β large enough f p3e´2βJ q ď 21 ñ xN´ y ď 48 N
MΩ 2xN´ y 1
ñ N
p
“1´ N ě1´ 24 independent of system size X

3.3.2 Kramers-Wannier method to determine Tc

Consider Ising model, square lattice, H “ 0 : +++++
+++++
Number of sites: N ++−−+
+−+++
Number of edges (”bonds”): K +++++
This subsection: Exact method to determine Tc (under certain assumptions)
based on a comparison of different series expansions of the partition function
; Smart approach (fun), will also teach us about series expansions
(a) Series Expansions of the partition function: Two approaches

(i) Low temperature expansion

+++++ ř
‹ T “ 0: + + + + + ZN “ tSi u eβJ xijy Si Sj “ eβJK
ř
+++++
+++++
+++++
+++++
‹T finite: + + + + + ; Polygons P with total contour length l
++−−+
+−+++ Costs: ∆H “ 2βJl
ř+ + + + + ř
βJ
ZN “ xijy i Sj “ eβJK
S ř ´2βJ lpP q
tSi u e e |
P summed over all polygons
all polygon
configurations

(NB: Uniqueness: Polygons do not cross each other: =

)
‹ At low temperatures, the dominant contributions correspond to
polygons with short contour lengths.
ñ ”Low temperature expansion”: Expansion in powers of e´2βJ
26 CHAPTER 3. THE ISING MODEL

‹ Construction of the coefficients (lowest order terms)

l “ 4: + − + Ñ N possibilities
l “ 6: − − Ñ 12 N ¨ 4 possibilities
l “ 8: − − 1
,
2 N pN ´ 5q .
“ 1 N pN ` 9q possibilities
−−− , −−
− Ñ `6N
−−
−− - 2
`N
; ZN “ e βJK p1 ` N e´2βJ¨4 ` 2N e´2βJ¨6 ` 12 N pN ` 9qe´2βJ¨8 ` ¨ ¨ ¨ q
ÿÿ
ñ General result: ZN “ eβJK e´2βJl
l tP ul
|
Possible polygon configurations made of boundary lines
with total contour length l (all polygons)
NB: Expansion is finite (l ď 2N ) for finite systems
Infinite series at N Ñ 8, possibly with convergence radius

(ii) High temperature expansion

‹ T Ñ 8 or β “ 0 ñ ZN “ tSi u 1 “ 2N “: ZN,0
ř
ř
ř βJ xijy Si Sj
ř βJ ř Si Sj e
‹ T finite: ZN “ e xijy “ ZN,0 tSi u
ř
1 tSi u
tSi u
ř
βJ xijy Si Sj
“ ZN,0 e
@ D
0–average at β “ 0
βJSi Sj
e “ coshpβJq ` Si Sj sinhpβJq

“ ZN,0 pcoshpβJqqK
@ś D
p1 ` v Si Sj q 0
xijy
with v “ tanhpβJq
‹ At high temperatures, the dominant contributions correspond to
terms of low order in v “ tanhpβJq .
ñ ”High temperature expansion”: Expansion in powers of v
‹ Construction of the coefficients: Graphical approach
Every term ”vSi Sj ” corresponds to a line i —j along the edge xijy.
; First Order: — ; Second order: —¨— + — etc.
˛ For unequal indices i, j, one has xSin1 Sjn2 y0 “ xSin1 y0 xSjn2 y0
For odd powers n, one has xSin y0 “ 0
; Graphs with free line ends or graphs, where an odd
number of lines meet at one vertex, do not contribute
(e.g., first order terms: - fl xSi Sj y0 “ xSi y0 xSj y0 “ 0).
; Only graphs consisting of closed polygons contribute.
˛ Lowest nonvanishing order: v 4 ( )
Next orders: v 6 ( )
v8 ( + + + )
ÿÿ
ñ General result: ZN “ 2N pcoshpβJqqK ptanhpβJqql
l tP ul
|
Possible polygon configurations made of edges
with total length l (all polygons)
3.3. TWO-DIMENSIONAL ISING MODEL: EXACT RESULTS 27

(b) Duality
Compare (i) and (ii): Very similar expressions
Z fl in both cases a sum over polygon configurations,
however, different alignment with the underlying lattice
(i) Low temperature expansion
Polygon lines perpendicular to bonds between lattice sites
(ii) High temperature expansion
Polygon lines lie on to bonds between lattice sites

Formal relation: Dual lattice

Original lattice: N0 lattice sites (vertices)
N1 edges (bonds) or
N2 faces (plaquettes)
Dual lattice:
Plaquette centers Ñ vertices pN0˚ “ N2 q
Bonds Ñ transverse bonds pN1˚ “ N1 q ,
Vertices Ñ plaquettes pN2˚ “ N0 q
Euler relation: N0 ´ N1 ` N2 “ N0˚ ´ N1˚ ` N2˚ “ 2
(NB: The outside of the graph counts as one face.)
Consequence for partition function
High temperature
L expansion: Polygon expansion on original lattice:
ZN0 ,HT pβq 2N0 pcoshpβJqqN1 “ ptanhpβJqql
ř ř
l tP ul
Low temperature expansion:
ř ř Polygon expansion on dual lattice:
˚
ZN2˚ ,LT pβq e
L 2βJN ´2βJ l
1 “ pe q
l tP ul
; Consider now an Ising model on the dual lattice with partition
function ZN˚˚ . Then the expansions of Z , Z ˚ can be mapped
0
onto each other: ZN˚˚ ,HT pβ ˚ q 9 ZN2˚ ,LT pβq
0

with tanhpβ ˚ Jq “ e´2βJ ô sinhp2βJq sinhp2β ˚ Jq “ 1

; Duality relation: High temperatures map onto low temperatures
Specifically, the square lattice is self dual: Z ˚ pβ ˚ q “ Z pβ ˚ q
Consequence for phase transition
Phase transition Ø Singularity of F “ ´kB T ln Z at N Ñ 8
` ř ř ´2βJ l ˘
Here: F pβq “ ´kB T ln pe q `N ¨ analytic terms
l tP ul
` ˘
ptanhpβ ˚ Jql `N ¨ analytic terms
ř ř
F ˚ pβ ˚ q “ ´kB T ln
l tP ul
ñ If F pβq is nonalytic at βc , then F ˚ pβ ˚ q is nonalytic at βc˚
Self dual lattice: F ˚ pβ ˚ q “ F pβ ˚ q
ñ F pβq is also nonanalytic at βc˚ , singularities come in pairs
Assume only one singularity ñ βc “ βc˚ ñ psinhp2βc Jqq2 “ 1
1 1 ?
ñ βc J “ arsinhp1q “ lnp1 ` 2q “ 0.4407
2 2
28 CHAPTER 3. THE ISING MODEL

Summary and Conclusions

‚ We have used symmetry considerations to determine Tc without actu-

ally ”solving” the Ising model (i.e., calculating the partition function)
‚ On this occasion, we also introduced the method of series expansions
for β Ñ 0 and β Ñ 8. They generally play an important role
independent of this argument here.

Remarks

‚ Duality trick does not work in three dimensions, since the lattice and
the dual lattice are too different
‚ The trick can also be used for the triangular lattice / honeycomb lattice
; ”star-triangle transformation” ?
"
´1 2`? 3 : triangular lattice
Result: ptanhpβc Jqq “
3 : honeycomb lattice

3.3.3 Exact solution of the Ising model on the square lattice

In his original solution of the 2D Ising model, Onsager (1944) used the transfer
matrix method (see Sec. 3.2.2).

Here, we present an alternative approach due to Samuel (1980), which is based

on the high temperature expansion (Sec. 3.3.2) and the mathematical
framework of ”Grassmann variables”.
For a more detailed discussion see also the book by F. Wegner:
”Supermathematics and its Applications to Statistical Physics”

We consider an Ising model with N “ Lx ˆ Ly spins on a square lattice

(1) General properties of Grassmann variables

‚ Symbols ξ with ξi ξj “ ´ξj ξi pñ ξi2 “ 0q

Application example: exppξi q “ 1 ` ξi
‚ Generate so-called ”Grassmann algebra”
ř ř
with elements A “ a ` i ai ξi ` i,j aij ξi ξj ` ¨ ¨ ¨
ş
‚ Formally define ”Integration” via dξN dξN ´1 ¨ ¨ ¨ dξ1 ξ1 ξ1 ¨ ¨ ¨ ξN “ 1,
where integral becomes zero, if one of the ξi is missing.
¯ exppξaξq ¯ “ dξdξp1
¯ ¯ “a
ş ş
ñ dξdξ ` ξaξq
¯ expp ij ξi Aij ξ¯j q “ detpAq
ş ř
drξsdrξs
¯ :“ dξ¯N dξN dξ¯N ´1 dξN ´1 ¨ ¨ ¨ dξ¯1 dξ1
ş ş
with drξsdrξs
ş
‚ Change of variables: Given a Grassmann integral drξsf rξs
ř
Linear substitution ξi Ñ θj “ k Jjk ξk
1
ş ş
Then, one has (no proof): drξsf rξs “ drθsf rξrθss detpJq
NB: For nonlinear substitution rules, the transformation rule may
be more complicated, but we will not need that here!
3.3. TWO-DIMENSIONAL ISING MODEL: EXACT RESULTS 29

(2) Strategy for solving the two dimensional Ising model (Samuel, 1980)

‹ Starting point: High temperature expansion (Sec. 3.3.2)

ZN “ 2N pcoshpβJqq2N Zx N with ZN “
lP
ř
CP v , v “ tanhpβJq
x
ř
with CP : Sum over all configurations of closed polygons
lP : Total length of polygons
Idea: Generate sum over all polygon configurations
via Grassmann integral
‹ Procedure
1
piq x
‚ Assign four Grassmann variables ξnm 4x x 3
x
(i “ 1¨¨ 4) to each lattice site pn, mq 2
‚ Define linker ˆl: Pairs of Grassmann variables
pijq piq pjq x x
Local: ênm “ ξnm ξnm for i ă j fl x x , x x
x x x
pxq p3q p4q
Bond: b̂nm “ ξnm ξpn`1qm x x
x x
pyq p1q p2q x
b̂nm “ ξnm ξnpm`1q , x x x x
x x x
x x
x
‚ Define ”action” (quadratic in ξ)
pxq pyq pijq ‰
Srξs “ nm v pax b̂nm ` ay b̂nm q ` iăj aij ênm ” cα ˆlα
ř “ ř ř
with ax , ay , aij “ ˘1 (sign to be determined below).
ş
‚ Consider Grassmann integral I “ drξs exppSrξsq
ş ş p4q p1q p4q p1q p4q p1q
(Order: drξs “ dξLx ,Ly¨¨ dξLx ,Ly dξLx ´1,Ly¨¨ dξ1,Ly dξLx ,Ly ´1¨¨ dξ1,1 )
; I has additive contributions from all configurations, in which
2Lx Ly linkers ˆlα are distributed such that every position
pn, m, iq is occupied by precisely one Grassmann variable.
Example: Configurations
ş ś contribute
ˆ l
(additively)
x
x
x x
x
xwith drξs cα lα 9 ˘ v
x
x
x x
x
x
P
x x x x
α
x
x

x
x x
x

x
x x
x

x
x x
x

x
x Sign: Depends on the power of ax , ay , aij
x x x x
in I and the number of transpositions
piq
needed to sort α ˆlα by ascending ξnm
x x x x x x x x ś
x x x x

ñ - I sums automatically over all configurations with closed

polygons. Free ends are not possible :-)
- Contribution of each polygon configuration is 9 v lp :-)
- But: Prefactor could be negative :-(
; Goal: Choose ax , ay , aij “ ˘1 such that every polygon configuration
has the weight v lP . Then we have ZˆN “ I, and I can be calculated
following the integration rules described in (1).

(3) Determination of coefficients ax , ay , aij

First consider configurations, in which polygons don’t touch, i.e., don’t
share corners (For touching polygons see step (iii) below). Calculate their
weight by rearranging, reorienting and reassigning linkers.
30 CHAPTER 3. THE ISING MODEL

Notation: kx , ky : Number of bonds in x,y direction (lP “ kx ` ky )

Np : Number of polygons
NB: Linkers commute, since they are pairs of Grassmann variables.
kx and ky are even, since polygons are closed.
Steps:
(i) Calculate contribution of isolated lattice sites
Three possibilities, x x x
x x , x x , x x
x x x x x x x
x x x x x x x x

Factor: pa12 a34 ` a23 a14 ´ a13 a24 q per site. x x x x

!
; Postulate pa12 a34 ` a23 a14 ´ a13 a24 q “ 1
x x x x
x x x x x x x x
x x x x

(ii) Now consider polygons. Assemble all bonds x

x

x
x x
x

x
x x
x

x
x x
x

x
x

belonging to same polygons together in chains.

x x x x

(iia) Reorient bonds in polygons such that they

x x x x x x x x
x x x x

run counterclockwise. x
x

x
x x
x

x
x x
x

x
x x
x

x
x

Costs: Factor p´1qpkx `ky q{2 x x x x

x x x x x x x x

(half of all bonds must be reoriented) x x x x

(iib) Shift linkers in polygon chains by one ξ x

x
x x
x
x x
x
x x
x
x

; New linkers are all local.

x x x x

x x x x
pijq x x x x x x x x

nij new linkers of the form ênm . x x x x

Costs: p´1qNp (one cyclic permutation x

x

x
x x
x

x
x x
x

x
x x
x

x
x

of ξ variables per polygon)

x x x x

(iic) Reorient new linkers pijq such that i ă j

ř ê
x x x x x x x x
x x x x

Costs: Factor p´1q iąj nij . x

x

x
x x
x

x
x x
x

x
x x
x

x
x

(iid) Determine weight of polygon elements: x

x
x x
x
x x
x
x x
x
x

k
- Bonds: v lp akxx ayy “ v lp (kx , ky even)
x x x x

- Joints between bonds:

x x
x
x
x
x x
x
x x
x
x
x x
x x
x x (n̄ij “ nij ` nji )
x x
x x x x

an̄1234 an̄3412 an̄2314 an̄1423 p´a13 qn̄24 p´a24 qn̄13

ş
(Sign: dξ p4q dξ p3q dξ p2q dξ p1q êpijq êpklq for site with linkers êpijq êpklq )
(iie) Summarize: Polygon configuration has weight W vřlP with
1
W “ an̄1234 an̄3412 an̄2314 an̄1423 an̄1324 an̄2413 p´1qNp ` 2 pkx `ky q` iąj nij `n̄24 `n̄13
Exploit relations between kx,y and nij :
– Every straight polygon line has corners at both ends
ñ Lines up: ky {2 ´ n21 “ n24 ` n23 “ n41 ` n31
Lines down: ky {2 ´ n12 “ n42 ` n32 “ n14 ` n13
Lines right: kx {2 ´ n43 “ n41 ` n42 “ n13 ` n23
Lines left: kx {2 ´ n34 “ n14 ` n24 “ n31 ` n32
– If one runs through a polygon in a counterclockwise way, one has
four more left corners than right corners.
ñ n24 ` n32 ` n13 ` n41 “ 4Np ` n42 ` n14 ` n31 ` n23
– Collect all this:
n24 “ Np ` n31 , n13 “ Np ` n42 , n32 “ Np ` n14 . n41 “ Np ` n23
n̄14 “ n̄23 , n̄13 “ n̄24 , ky ´ n̄12 “ kx ´ n̄34 “ n̄13 ` n̄14
ñ W “ pa13 a24 qn̄13 p´a14 a23 qn̄14 pa12 a34 q´n̄13 ´n̄14
3.3. TWO-DIMENSIONAL ISING MODEL: EXACT RESULTS 31

(iii) Possible choices for ax , ay , aij

‹ Recall |ax | “ |ay | “ |aij | “ 1
‹ Further conditions from (i) and (ii)
!
(i) a12 a34 ` a23 a14 ´ a13 a24 “ 1
!
(ii) W “ pa13 a24 qn̄13 p´a14 a23 qn̄14 pa12 a34 q´n̄13 ´n̄14 “ 1 @ n̄12 , n̄14
ñ a12 a34 “ ´1, a13 a24 “ ´1, a14 a23 “ 1, ax.y “ ˘1

Consistency check: What happens, if polygons touch each other?

x
Corresponds to bond constellation x
x
x
(”cross”)
x x
Steps in (ii) turn this into x
x
x
or x
x
x
(equivalent)
; In W (iii), two ”corners” are replaced by one ”cross”
x x
x
x
x
x
x x fl p´a13 qp´a24 q “ ´1 Ñ x
x
x
“ ´1 X
x
x
x x
fl pa14 qpa23 q “ 1 Ñ “1 X
x x
x x x x
x
x x

Conclusion: Conditions (iii) make sure that I corresponds to a sum

over all polygon configurations with weight v lP .
ñ I “ Zx“ CP v lP , hence I can be used to calculate Zx.
ř
In practice, we still have some freedom and choose
ax “ ´1, ay “ 1, a12 “ a24 “ ´1, a13 “ a14 “ a23 “ a34 “ 1

(4) Calculation of the partition function

‹ Remaining task: Calculate Zx

ş
N “ drξs exp Srξs
ř “ pxq pyq ř pijq ‰
with Srξs “ nm vpax b̂nm `ay b̂nm q` iăj aij ênm pv “ tanhpβJqq
and coefficients ax,y , aij from (3)
pijq piq pjq pxq p3q p4q pyq p1q p2q
(ênm “ ξnm ξnm , b̂nm “ ξnm ξpn`1qm , b̂nm “ ξnm ξnpm`1q )
‹ Fourier transform: Assume periodic boundary conditions (as in 3.2.2).
Then, configurations exist where domain interfaces span the whole
system, i.e., they do not form closed polygons. However, the sta-
tistical weight of such configurations decreases exponentially with
increasing system size, and we will thus neglect them.
2πk Lx Lx L L
Define pk “ ql “ 2πl
Lx ,
y y
Ly with k P r´ 2 , 2 s, l P r´ 2 , 2 s
piq ´ippk n`ql mq ξ piq piq piq˚
“?1
ř
and ζkl mn e nm (ñ ξ´k,´l “ ξkl )
Lx Ly
piq a ´1 ř ippk n`ql mq piq
ñ ξnm
“ Lx Ly kl e ζ
ř pijq ř piq pjq 1
ř kl` piq pjq pjq piq ˘
nm ênm “ ζ ζ
kl kl ´k,´l “ 2 kl ζkl ζ´k,´l ´ ζkl ζ´k,´l
pxq ř p1q p3q ` p1q p3q p3q p1q ˘
´iql
“ 12 kl ζkl ζ´k,´l e´iql ´ζkl ζ´k,´l eiql
ř ř
nm b̂nm “ kl ζkl ζ´k,´l e
pyq ř p2q p4q ř ` p2q p4q p4q p2q ˘
´ipk
“ 21 kl ζkl ζ´k,´l e´ipk ´ζkl ζ´k,´l eipk
ř
nm b̂nm “ kl ζkl ζ´k,´l e

1 ř ř piq pklq pjq ř1 ´ ř piq pklq pjq˚

¯
ñ S“ 2 kl ij ζkl Aij ζ´k,´l “ kl ij ζkl Aij ζkl “: Srζ, ζ ˚ s,
ř1
where klsums only over half of the pklq,
such that pklq and p´k, ´lq are both fully covered
32 CHAPTER 3. THE ISING MODEL
¨ ˛
0 ´1 ´ ve´iql 1 1
˚1 ` veiql 0 1 ´1
and Apklq “
‹
˚ ‹
´ipk ‚
˝ ´1 ´1 0 1 ` ve
´1 1 ´1 ´ veipk 0
` pklq p´k,´lq pklq˚ ˘
NB: Aij “ ´Aji “ ´Aji
ñ Partition function (using equations from (1))
Zx
ş
N “ drξs exp Srξs
(Jacobi determinant is one)
´ř `ř ¯
ş ş 1 piq pklq pjq˚ ˘
“ drζ ˚ sdrζs exp Srζ, ζ ˚ s “ drζ ˚ sdrζs exp kl ζ A
ij kl ij ζ kl
ś1 ś ?
“ kl det A pklq “ p kl q det A pklq
all
ś a
“ kl p1 ` v 2 q2 ´ 2vp1 ´ v 2 qpcos pk ` cos qk qq

(5) Conclusion: Exact solution of the Ising model

Free energy: (From ZN “ 2N pcoshpβJqq2N Zx N)

F “ ´kB T ln!Z
‰)
“ ´N kB T lnp2 cosh2 βJq ` 2N
1 ř
“ 2 q2 ´ 2vp1 ´ v 2 qpcos p ` cos q q
kl ln p1 ` v k l

v “ tanh βJ
” ı
1 ř 2
“ ´N kB T 2N kl ln 4 cosh p2βJq ´ 4 sinhp2βJqpcos p k ` cos q l q
ř 1
ťπ
Thermodynamic limit: kl Ñ p2πq2 ´π dpx dpy
ijπ
1 ” ı
ñ F “ ´N kB T 2 dpx dpy ln 4 cosh2 p2βJq ´ 4 sinhp2βJqpcos px ` cos py q
8π
´π
Corresponds to the result of Onsager!
Analysis:
A phase transition is expected, if the argument of lnr¨ ¨ ¨ s is zero.
; cosh2 p2βJq “ sinhp2βJqpcos px ` cos py q for one ppx , py q
; Possible for pcos px ` cos py q “ 2, i.e., ppx , py q “ p0, 0q
!
Then, one has: cosh2 2βc J “ 1 ` sinh2 2βc J “ 2 sinh 2βc J
ñ p1 ´ sinh2 2βc Jq “ 0 ?
ñ βc J “ 12 arsinh1 “ 12 lnp1 ` 2q
; Same result as in Sec. 3.3.2!
But: from the exact solution, one can also calculate other quantities,
such as, e.g., the specific heat Ñ Exercise
(One obtains c „ lnpT ´ Tc q: Logarithmic divergence)
3.4. SERIES EXPANSIONS: GENERAL REMARKS 33

3.4 Series Expansions: General Remarks

Last subsection (3.3.2): Introduction of the concept of series expansions – an
important technique when studying phase transitions analytically

In particular, the high temperature expansion turns out to be a powerful and

highly versatile tool in statistical physics.

Basic idea (in quantum mechanics notation)

‚ For arbitrary statistical averages (canonical ensemble), we have
TrpAe´βH q (Classically, ”Tr” refers to the suitable phase
xAyβ “ space
ř integral or sum, e.g. in the Ising model,
Trpe´βH q ” tSi u ¨ ¨ ¨ ”, and H fl H to the Hamiltonian.)
‚ Define Z0 “ Trp1q and rewrite:
Numerator of xAyβ : TrpAe´βH q “ Z10 xAe´βH y0
Denominator of xAyβ : Trpe´βH q “ Z10 xe´βH y0 ,
where x¨ ¨ ¨ y0 : Statistical average at β “ 0 (T Ñ 8)
xAe´βH y0
‚ Then expand xAyβ “ xe´βH y0
in powers of β
´βH 1 2 2
(e “ 1 ´ βH ` 2
β H ` ¨¨¨)
Leading terms:
xAy0 ´βxAHy0 ` 12 β 2 xAH 2 y0 `¨¨¨
xAyβ “ 1´βxHy0 ` 12 β 2 xH 2 y0 `¨¨¨
` ˘
= xAy0 ´´β xAHy0 ´ xAy0 xHy0 ¯
` 21 β 2 xAH 2 y0 ´ 2xAHy0 xHy0 ´ xAy0 xH 2 y0 ` 2xAy0 xHy20
` Opβ 3 q
Free energy (βF “ ´ lnpTrpe´βH qq) can be expanded analogeously.
Very general approach

Low temperature expansions are also possible (see, e.g., Sec. 3.3.2), but the
design principles are less generic (Setting up such an expansion requires
the knowledge of the elementary excitations in the system).

General remarks on series expansions

‚ In general, graphical methods are useful for the construction, see, e.g.,
the graphical expansions in Sec. 3.3.2 in polygon configurations
(The ”diagrams” of the expansion are the polygon configurations.)
‚ Simplifications can often be identified beforehand based on general
considerations (e.g., symmetry considerations).
‚ Important example: Linked Cluster Theorem: Only connected dia-
grams (configurations with connected polygons) contribute to F 9 ln Z ,
diagrams with unconnected components cancel out.
(Heuristic ”proof”: Every unconnected component comes with a
combinatorial factor 9 N . However, F is extensive, therefore, they
must all cancel each other!)
34 CHAPTER 3. THE ISING MODEL

Analysis of series expansions

Starting point: Series f pzq “ n an z n

ř
Only a finite number of coefficients are known.
Question: Assume that f pzq has a singularity, f pzq „ pz ´ zc q´γ
What can we learn from the series about the singularity?
Example: Consider simple function f pzq “ p1 ´ z{zc q´γ (with γ ą 1)
ř ` ˘ ´1 n n
ñ Expansion f pzq “ n ´γ n p zc q z
ñ an “ p z1c qn p´1qn n “ p z1c qn γpγ`1q¨¨¨pγ`n´1q
` ´γ
˘
n!
ñ an´1an
“ z1c γ`n´1
n “ 1
zc p1 ` γ´1
n q
ñ Possible strategies for determining zc :
‚ Simply plot rn “ an´1 an
versus n1 ; Axis intercept gives 1{zc !
(Generally, limnÑ8 an´1 an gives the radius of convergence of
the series. Therefore, this method works, if the convergence
radius is determined by the singularity at zc )
‚ More efficient method: Eliminate term γ´1 n in our example
by choosing rn “ n an´1 an
´ pn ´ 1q aan´2
n´1

(Gives rn ” z1c in our example. In general, corrections ap-

ply.)
‚ There exist numerous other, much more sophisticated ap-
proaches, e.g., Padé approximants. Analyzing series expan-
sions is an art in itself
If zc is known, similar techniques can be applied to determine γ.
an
For example, a simple estimator is Sn “ 1 ` np an´1 zc ´ 1q
3.5. MEAN-FIELD APPROXIMATION 35

3.5 Mean-Field Approximation

Often less involved than series expansions, more general approach, not restricted
to regions without singularities. ”Sufficient” for many purposes.
But: Uncontrolled approximation

3.5.1 Simplest approach: Spins in mean fields

3.5.1.1 Approach via effective field (intuitive approach)

H “ ´J
ř ř
+++ Si Sj ´ H i Si
xijy
−++ nearest neighbors
+−− i “ ´ BH “ H ` J
Interactions Ø Effective field Heff BSi
ř
j Sj

Mean-field approximation: Replace Sj by xSy “ m.

Heff “ H ` Jqm, with q: coordination number (cubic lattice: 2 ¨ D)
piq
Consider single spin in the external field Heff : Heff “ ´Heff S
; Partition function Z “ eβHeff ` e´βHeff
βH ´e´βHeff
Magnetization: m “ eeβHeffeff `e´βHeff
“ tanhpβHeff q
` ˘
ñ Implicit, self consistent equation for m: m “ tanh βpJqm ` Hq

At H “ 0, the equation can be solved graphically:

β large β small β = βc
m m
m
tanh(βJqm)
tanh(βJqm)
tanh(βJqm)

Three crossings Only one crossing Tangent crossing

At the critical point β “ βc

ñ Only one crossingˇ point, butˇ slopes are equal
d d ˇ
tanhpβc Jqmqˇ m“1
ˇ
ñ dm “ dm ˇ
m“0 m“0
ñ βc Jq “ 1

Close to the critical point, tanhp¨ ¨ ¨ q can be expanded in powers of m.

ñ m “ tanhpβqJmq « βqJm ´ 13 pβqJq3 m 3 ` ¨¨¨
a ? b
ñ mpT q « 3pβqJ ´ 1q{pβqJq3 “ 3 TTc 1 ´ TTc 9 pTc ´ T q1{2
ñ Critical exponent β: β “ 1{2

Within this approach, one can calculate the spontaneous magnetization and
the susceptibility, but not the entropy or the free energy.
36 CHAPTER 3. THE ISING MODEL

3.5.1.2 Approach via free energy (Bragg-Williams approximation)

Starting point: Free energy F “ U ´ T S

”Mean-field” approximation: Spins are not correlated:

; Joint probability function factorizes: P pS1 , S2 ,¨¨ , SN q « N
ś p1q
j“1 p pSj q
with pp1q pSq: probability distribution for single spin
A ř E
(a) Energy: N1 U “ N1 H “ N1 ´ J
@ D ř
Si Sj ´ H Si
xijy i
´ ř @ D@ D ř @ D¯
1
« N ´J Si Sj ´ H Si “ ´ 12 Jq m2 ´ H m
xijy i

(b) Entropy: S “ ´kB tS1 ,¨¨ ,SN u P pS1 ,¨¨ , SN q lnpP pS1 ,¨¨ , SN qq
ř
ř
« ´kB N S pp1q pSq ln pp1q pSq

Given magnetization m, uncorrelated spins Sj

; Construct probability function pp1q pSq such that xSy “ m
p1q p1q
Notation: p` :“ pp1q p`1q, p´ :“ pp1q p´1q
p1q p1q p1q p1q
ñ xSy “ p` ´ p´ and p` ` p´ “ 1
p1q
ñ xSy “ 2p` ´ 1
p1q p1q
ñ p` “ p1 ` mq{2, p´ “ p1 ´ mq{2
` p1q p1q p1q p1q ˘
ñ N1 S “ ´k“B p` ln p` ` p´ ln p´ ‰
“ kB ´ 1`m 2 ln 1`m
2 ´ 1´m
2 ln 1´m
2

F “ ‰
ñ Free energy: “ ´ 12 Jqm2 ´ Hm ` kB T 1`m 1`m
2 ln 2 `
1´m
2 ln 1´m
2
N
` ˘ !
BF
Minimization: Bm “ ´qJm ´ H ` kB T 21 ln 1`m
1´m “ 0
`1 ` m˘
ñ ln “ 2βpqJm ` Hq ñ m “ tanh βpqJm ` Hq
1´m
; Approximation equivalent to the approximation of 3.5.1.1

3.5.1.3 Problems with this mean-field aproximation

– Geometry enters only via coordination number q
; no dependence on dimension, local structure etc.

– Predicts a phase transition also for one dimensional Ising model

(clearly wrong)

– Wrong critical point Tc , wrong critical exponents

Example: Two dimensional Ising model on the square lattice:
Mean-field Ñ βc J “ 1q , Critical exponent β “ 1{2
?
Exact Ñ βc J “ 1q 2 lnp1 ` 2q “ 1q ¨ 1.76
Critical exponent β “ 1{8
3.5. MEAN-FIELD APPROXIMATION 37

3.5.2 Improved theory: Clusters in mean fields

3.5.2.1 Approach via effective field (Bethe approximation)

Cluster with central particle

ř ”0”
S0
ñ HCluster “ ´J S0 qj“1 Sj ´ Heff qj“1 Sj
ř
(Effective field Heff acts on outer spins only
Central spin is treated exactly)
˘q
ñ ZCluster “ eβpJ`Heff q ` e´βpJ`Heff q
`
: Contribution of S0 “ `1
` βp´J`H q ´βp´J`H q
˘q
` e eff `e eff : Contribution of S0 “ ´1
”´ ¯q ´ ¯q ı
1 βpJ`Heff q ´βpJ`Heff q βp´J`Heff q ´βp´J`Heff q
xS0 y “ Z
e ` e ´ e ` e
Cluster
„´ ¯q´1 ´ ¯
xSj y “ Z 1 eβpJ`Heff q ` e´βpJ`Heff q ¨ eβpJ`Heff q ´ e´βpJ`Heff q
Cluster
´ ¯q´1 ´ ¯
βp´J`Heff q ´βp´J`Heff q βp´J`Heff q ´βp´J`Heff q
` e `e ¨ e ´e

!
; Condition for Heff : xS0 y “ xSj y p“ mq
cosh βpJ ` Heff q
ñ “ e2βHeff {pq´1q
cosh βp´J ` Heff q

Solutions:

(i) Heff “ 0: Disordered state

(ii) Heff ‰ 0 (if β is not too small): Ordered state

Transition point: Expand about small Heff

” ı ” ı
cosh βpJ`Heff q !
; cosh βp´J`Heff q « 1`2βHeff tanh βJ “ e2βHeff {pq´1q « 1`2β Heff
pq´1q
` q ˘
ñ cothpβc Jq “ q ´ 1 ñ 2βc J “ ln
q´2

Remarks:

– In one dimensions, one has q “ 2 ; no phase transition X

– Two dimensional Ising model on the square lattice:
Exact: βc J “ 0.44
Bragg-Williams: βc J “ 1{4 “ 0.25
Bethe: βc J “ lnp2q “ 0.35: Significant improvement!
– Higher coordination numbers:
Bragg-Williams: 2βc J “ 2{q
q
Bethe: 2βc J “ ln q´2 “ ´ lnp1 ´ 2{qq “ 2{q ` ¨ ¨ ¨
; Bragg-Williams and Bethe approximation agree at lowest order
of 1{q. Results are never identical!
– Critical exponents in Bethe and Bragg-Williams approximation are
the same: No improvement in this respect.
38 CHAPTER 3. THE ISING MODEL

3.5.2.2 Approach via free energy (Guggenheim approximation)

Main approximation in Bragg-Williams theory: Independent spins

Guggenheim approximation: Independent clusters, neglect cluster correlations

; Improved treatment of pairs of neighbor spins

; Improved treatment of pairs of neighbor spins

• Probability for one single cluster:
pcluster pS0 ; tSj uq “ pp1q pS0 q qj“1 P pS0 Sj |S0 q
ś
conditional probability for pS0 Sj q given S0
S0
P pS0 Sj |S0 q “ pp2q pS0 Sj q{pp1q pS0 q
with pp2q pS0 Sj q: Pair probability
śq pp2q pS0 Sj q
“ pp1q pS0 q j“1 pp1q pS0 q
• Two neighbor clusters: Must correct for double counting of bonds
; pcluster S0 ; tSj u Pcluster S01 ; tSj u|pS0 S01M
` ˘ ` ˘
q
` ˘ ` ˘ S0 S‘0
“ pcluster S0 ; tSj u pcluster S01 ; tSj u pp2q pS0 S01 q

ź Mź
; Whole system: P pS1 ,¨¨ , SN q « pp2q pSi Sj q pp1q pSi qq´1
xijy i

Construct probability functions such that xSy “ m

p1q p2q p2q p2q
Notation: p˘ :“ pp1q p˘1q, p˘˘ :“ pp2q p˘1, ˘1q, p`´ “ p´` “: a
p1q p2q p2q p1q p2q p2q
ñ p` “ p`` ` p`´ “ 1`m 2 , p´ “ p`´ ` p´´ “ 1´m 2
p2q 1`m p2q 1´m
ñ p`` “ 2 ´ a, p´´ “ 2 ´a

ñ Entropy: S “
ř
tS1 ,¨¨ ,SN u P pS1 ,¨¨ , SN q ln P pS1 ,¨¨ , SN q
´ ř ¯
q
NS
1 ř
ñ « ´kB 2 pp2q pSS 1 q ln pp2q pSS 1 q ´ pq ´ 1q pp1q pSq ln pp1q pSq
SS 1 S
! “
q
“ ´kB 2p 1`m 1`m 1´m
2 ´ aq lnp 2 ´ aq ` p 2 ´ aq lnp 2 ´ aq
1´m

‰ “ ‰)
` 2a ln a ´ pq ´ 1q p 1`m
2 q lnp 1`m
2 q ` 1´m
2 q lnp 1´m
2 q
“ p2q p2q p2q p2q ‰
Energy: 1
NU “ 1
N xH y “ ´J 2q pp`` ` p´´ q ´ pp`´ ` p´` q “ J 2q p4a ´ 1q

Free energy: F “ U ´ β1 S

BF ! BF !
Minimize free energy: Ba “ 0, Bm “ 0
” ˘ı !
1q
1 BF
`
• N Ba “ 2qJ ` β2 2 lnp2aq ´ ln p1 ` m ´ 2aqp1 ´ m ´ 2aq “ 0
´ p1 ` m ´ 2aqp1 ` m ` 2aq ¯
ñ 4βJ “ ln (i)
p2aq2
 3.5. MEAN-FIELD APPROXIMATION 39
” ı
1 BF 1 q q´1 !
• N Bm “ βln 1`m´2a
4
1´m´2a ´ 2 ln 1`m
1´m “ 0

q ´ 1 ` m ´ 2a ¯ ´1 ` m¯
ñ ln “ pq ´ 1q ln (ii)
2 1 ´ m ´ 2a 1´m

Critical point: m Ñ 0
” ı ” ı
!
(ii) ñ 2q lnp1 ` 1´2a
2m
q« q 2m
2 1´2a “ pq ´ 1q lnp1 ` 2mq « pq ´ 1q2m
q´2 q
ñ 2a “ 2pq´1q , 1 ´ 2a “ 2pq´1q

q
` q ˘
(i): ñ 4βc J “ 2 ln 1´2a
2a “ 2 ln q´2 ñ 2βc J “ ln
q´2

Remarks:

– Same result as in Bethe approximation X

– Systematic generalization to larger clusters is possible
; Cluster variation method
Popular method in the context of order/disorder phase transitions
For large clusters: Very good phase diagrams
But still: wrong critical behavior (see next chapter)

3.5.3 Critical behavior in mean-field theory

As already mentioned earlier, one often observes critical behavior at con-
tinuous transitions: Many properties exhibit singularities when plotted
against intensive variables such as temperature and magnetic field, which
are often characterized by power laws.

At a qualitative level, the same behavior can already be seen in mean-field

approximation. This will be shown in the present section.

Preview (defining t “ pT ´ Tc q{Tc and d: Spatial dimension)

Quantity Expo- Power law Value
nent Mean-field Ising exact
2D 3D
β m „ p´tqβ β “ 1{2 1{8 0.33
Magnetization m δ mδ „ H δ“3 15 4.8
at t “ 0
Bm
Susceptibility χ “ BH γ χ „ |t|´γ γ“1 7{4 1.24
Specific heat
BS
cH “ T p BT qH “ p BE BT qH α cH „ |t|´α α“0 0plogq 0.1
Correlations
Gij “ xSi Sj y ´ xSi y xSj y η Gp~rq „ r2´d´η η“0 1{4 0.04
“: Gp~ri ´ ~rj q at t “ 0 1{4 0.04
Correlation length
Gp~rq „ e´r{ξ ν ξ „ |t|´ν ν “ 1{2 1 0.63
40 CHAPTER 3. THE ISING MODEL

Calculation: In Bragg-Williams approximation

Bethe-Guggenheim calculation shall not be shown here, but the results
are the same. In the next section we will see, why.
Define h “ βH and t “ pT ´ Tc q{Tc “ βc {β ´ 1
with βc “ 1{qJ (q: Coordination number)

3.5.3.1 Magnetization

Starting point: (see Sec. 3.5.1.2): m “ tanh βpqJm ` Hq “ tanhp ββc m ` hq

β
At |t|, |h| ! 1, one has |m| ! 1 and hence m « βc m ` h ´ 13 p ββc m ` hq3

Consider limit t Ñ 0´ , h “ 0
β
ñ m« βc m ´ 31 p ββc mq3
β
?
ñ βc m “ 3 p´tq1{2 ñ m „ p´tq1{2 ñ β “ 1{2

Consider limit t “ 0, h Ñ 0 (i.e., β{βc “ 1)

ñ m « m ` h ´ 31 pm ` hq3
ñ h « 13 m3 ` Opm2 h, mh2 q ñ h „ m3 ñ δ“3

3.5.3.2 Susceptibility

Starting point: Same as before in Section 3.5.3.1

Consider limit t Ñ 0, h Ñ 0
Define gpm, hq :“ tanhp ββc m ` hq ´ m « mp ββc ´ 1q ` h ´ 13 p ββc mq3
ˇ M
ñ gpm, hq ” 0 ñ Bm ˇ “ ´ Bg Bg Bh t Bh Bm
with “ 1, Bg
Bh “ Bg
Bm p ββc ´ 1q ´ β β
βc p βc mq
2 “ ´ ββc pt ` p ββc mq2 q
βc
#
Bm ˇ
ˇ
β βc t´1 : t ą 0 pm “ 0q
ñ ˇ “ “
Bh t β 2
t ` p βc mq q β p´2tq´1 : t ă 0 pp ββc mq2 “ ´3tq
ñ γ“1

3.5.3.3 Specific heat

Starting point: U “ xH y “ ´J xijy xSi y xSj y “ ´N J 2q m2 “ ´N 2β1 c m2

ř
#
0 : T ą Tc
with m2 “ βc pT ´Tc q T 2
p´tq3 β “ Tc p Tc q : T ă Tc
"
BU ˇˇ 0 : T ą Tc
ñ cH “ ˇ “ N
BT H ´ 2 kB : T À Tc with (T Ñ Tc´ )

Consider limit t Ñ 0: Finite jump ñ α“0

3.5. MEAN-FIELD APPROXIMATION 41

3.5.3.4 Correlation functions

Less straightforward, since correlations are ignored in mean-field theory.

Starting point: Consider Ising model at H “ 0, regular lattice with simple unit
cell and lattice vectors ~ri , d dimensions, periodic boundary conditions.
The interaction range is characterized by a set of neighbor vectors t~τ u,
i.e., spins Si , Sj interact if ~rij :“ p~rj ´ ~ri q P t~τ u). NB: If τ is a neighbor
vector, then (- τ ) is a neighbor vector as well.

Goal: Calculate Gp~rij q “ xSi Sj y ´ xSi y xSj y

Trick: Use general relation between fluctuations and susceptibilities

Consider generally an energy function of the form B: H “ H0 ´ HB B

1 BxAy
Then we have the general relation xABy ´ xAy xBy “ β BHB
ř ´βH `βH B
1 e 0 B A
e´βH A “
ř
(Proof: xAy “ Z
ř ´βH `βH B
e 0 B
Configurations
ř ´βH `βH B
Z B βBA´p e´βH0 `βHB B βBqp e´βH0 `βHB B Aq
ř ř
BxAy e 0
ñ BHB
“ Z2
“ βxABy ´ βxAy xBy X)

Of course also valid in the case HB “ 0

Here: Consider H “ ´J xijy Si Sj ´ i Hi Si
ř ř
; Inhomogeneous system with xS`i y “ řmi different for all i ˘
Solution as before: mi “ tanh βpJ neighbors j of i mj ` Hˇ i q
Can be evaluated ñ G~rij “ xSi Sj y ´ xSi y xSj y “ β1 BH
Bmi ˇ
j
ˇ
hk “0@k

ż
1 1 ~
Solution for T ě Tc : Gp~rq “ dd k eik¨~r
p2πqd ~ τ qq
ř
1 ´ βJ τ cospk ¨ ~
~
1st Brillouin
zone
(Calculation:
` ř ˘ ř
First linearize: mi “ tanh βpJ mj ` Hi q « βpJ mj ` Hj q
neighbors neighbors
"
J : ~rij P t~τ u
ñ j Bij mj “ Hj with Bij “ β1 δij ´
ř
0 : otherwise
Bmi
ñ mi “ j pB ´1 qij Hj ñ Gp~rij q “ β1 BH “ β1 pB ´1 qij
ř
j

Then diagonalize and invert B by Fourier transform. Define Bij “: Bp~rij q

ř ´i~k¨~r
ñ B̃p~kq “ Bp~rq “ β1 ´ J cosp~k ¨ ~τ q
ř
e
lattice ~
τ
vectors ~ r
~
d k pβ B̃p~kqq
ş d
ñ 1
β
pB ´1
q ij “ 1
p2πqd
´1
eip~ri ´~rj q¨k Xq

ż
1 βc {β ~
Simplification for t Ñ 0`: Gp~rq “ dd k eik¨~r
p2πqd 8 t ` k 2 vp~ e ~k q
~k 1 ř
with ~e~k “ k and vp~eq “ 2q ~τ p~e ¨ ~τ q2

(Calculation:
Rewrite Bp~kq “ ` cosp~k ¨ ~τ q
1
` 1
ř ˘
βc
t´ q τ p1
~
t Ñ 0: Main“ contribution to integral‰ stems ~
“ from small‰ k!
1
ř ~ 2
ñ Expand t ` q ~τ p1 ´ cospk ¨ ~τ qq « t ` k vp~e~k q Xq
42 CHAPTER 3. THE ISING MODEL

Consider t “ 0 (~e~r “ ~r{r)

βc i~ r k̃“~
kr βc
ş d k¨~
ş d
ñ Gp~rq “ 1
p2πqd β 8
d k 1
k2 vp~e~ q
e “ 1
p2πqd β 8
d k̃ 1
k̃2 vp~
r
e~ q
eik̃¨~e~r 2´d
k loooooooooooooooooomoooooooooooooooooon k

independent of r

2´d
ñ Gp~rq „ r ; ñ η“0

Now consider case t Ñ 0

βc
ş ~ k̃“~
kr βc
ş
ñ Gp~rq “ 1
p2πqd β 8
dd k 1
t`k2 vp~
e~ q
eik¨~r “ 1
p2πqd β 8
dd k̃ 1
r 2 t`k̃2 vp~
e~ q
eik̃¨~e~r r2´d
k k

‚ r2 t ! 1 ñ r2 t « 0 ñ Gp~rq „ r2´d as before

‚ r2 t " 1 : Choose x direction in direction of ~r. Other directions: ~kK
şπ β
ş
ñ Gp~rq “ 1
2π ´π
dkx eirkx 1
βc p2πqd´1
dd´1 kk t`pk2 `k12 qvp~e q
x K ~
k
loooooooooooooooooooooomoooooooooooooooooooooon
“:gpkx q
Use theorem of residues, search for poles pj of gpkx q in the upper complex plane.
1 8
ş rÑ8
dkx eirkx gpkx q “ i j lim rpz ´ k̄j qf pzq eirz s ÝÑ const eirp̄
ř
ñ 2π 8 zÑpj
where p̄ is the pole that is closest to the real axis.
Main contribution to the integral stems from kK « 0
a
ñ Pole at t “ ´k̄x2 vp~ex q ñ k̄x “ i t{vp~ex q “: i{ξp~ex q
?
ñ Gp~rq „ e´r{ξp~e~r q with ξ „ 1{ t ñ ν “ 1{2

3.5.4 Validity of mean-field theory, Ginzburg criterion

3.5.4.1 Compare two methods for determining specific heat

Consider specifically the case t ą 0, t Ñ 0`

Recall Sec. 3.5.3.3: cH “ constř“ 0

Calculated from U “ ´J xijy xSi y xSj y
ř ř
Now: Alternative calculation from U “ ´J xijy xSi Sj y “ ´J xijy Gij
using the results from Sec. 3.5.3.4
1
ş 1 J ~
dd k eik¨~τ
ř
U “ ´N p2πq d ~ τ qq 2
ř
1´βJ τ cospk¨~
~
1st Brillouin ~
τ
ˇ zone
ˇ
ˇ Expansion about k “ 0
ş d 1´k2 vp~e~k q
« ´ N k2B T p2πq
1
d 8 d k t`k 2 vp~
e~k q
ˇ ?
ˇ
ˇ ~
k̃ “ k{ t, t Ñ 0 `

N kB T d{2´1 1 1
ş d
«´ 2 t p2πqd 8
d k̃ 1`vp~ e 9 T td{2´1
k̃ q

1 BU
ñ cH “ N BT 9 td{2´1 (at t Ñ 0` )
"
const for d ą 4 : consistent with Sec. 3.5.3.3
ñ cH „
divergent for d ą 4 : not consistent with Sec. 3.5.3.3

ñ Apparently, mean-field approximation breaks down for dimensions d ă 4.

; ”Upper critical dimension”

3.5. MEAN-FIELD APPROXIMATION 43

3.5.4.2 Alternative argument: Ginzburg criterion

Mean-field theory neglects correlations.
; should be oK, if the fluctuations within the correlation length ξ are small
compared to the magnetization!
ÿ ÿ ” ı ÿ ÿ
ñ xSi Sj y ´ xSi y xSj y ! xSi y xSj y
V pξq V pξq V pξq V pξq
loooooooooooooooooomoooooooooooooooooon looooooooomooooooooon
|| ||

ξd χ pξ d q2 m2
χ
ñ !1 with ξ „ |t|´ν , m „ |t|β , γ „ |t|´γ
ξ d m2
ñ R |t|´γ`dν´2β ! 1
R : System dependent factor (range of interactions etc.)

At the critical point t Ñ 0

ñ Condition dν ´ 2β ´ γ ą 0 with γ “ 1, β “ ν “ 1{2
ñ Fulfilled for d2 ´ 2 ą 0 ñ dą4

; Mean-field approximation captures correct critical behavior at d ą 4.

However: Fails for t Ñ 0` at d ď 4
44 CHAPTER 3. THE ISING MODEL

3.6 The Monte Carlo method

Problem: Calculate partition functions, statistical expectation values, phase
transitions in the Ising model or other ”microscopic” models

Looking back: Approaches we have discussed so far

3.3: Exact techniques

; Exact solutions for special cases (1D, 2D Ising)
3.4: Series expansions
; Also exact, but limited applicability (convergence radius)
3.5: Mean-field approximation
; More generally applicable, but uncontrolled approximation

Question: How can one obtain an ”exact” solution in the general case?

Answer: Up to now - Only numerically

”Sledgehammer approach”:
ř ´βH
Calculate xAy “ řee´βHA directly for finite systems
But: Use ”smart” sledgehammer ; Monte Carlo Simulations

Very general method

Broad applications in all areas of statistical physics and beyond. Also
heavily used in Mainz. Shall therefore be briefly illustrated here at the
example of the Ising model. (For more details see class ”Computer simu-
lations in statistical physics”).

3.6.1 Main idea of Monte Carlo integration

ř ´βH
Task: Calculate statistical expectation values xAy “ ře ´βHA
e

Solution strategies

(a) Exact enumeration : Full Sum over all configurations

Pros: Exact
Cons: very time consuming, only possible for tiny systems
Inefficient at small temperatures, since most configurations don’t
contribute much to the sum
Way out: Monte Carlo Integration
Sum only over a random sample of configurations, not all!
(b) Simple sampling Entirely random sample (every configuration
has equal probability)
pjq
(e.g.: Sample#j Ø N random numbers ri P r0 : 1s,
pjq
pjq `1 : ri ą 1{2
Si “ pjq
´1 : ri ă 1{2
pjq
pjq
ArtSi use´βH rtSi us
řn
j“1
Analysis: xAy “ limnÑ8 pjq
ArtSi use´βH rtSi us
řn
j“1
3.6. THE MONTE CARLO METHOD 45

Pro: Results can already be obtained with small samples

Can be improved systematically by increasing sample size
; One is less restricted with respect to system size
Cons: – No longer exact
(however, accuracy can be controlled via sample size)
– Still inefficient at low temperatures
(c) Importance sampling Draw sample according to the distribution
P tSi u 9 e´βH rtSi us
Analysis: xAy “ limnÑ8 nj“1 ArtSi s
ř pjq

Pros: – First results can already be obtained with small samples

Can be improved systematically by increasing sample size
(as in (b))
– Efficient: Configurations that contribute to the sum with
higher weight are drawn more often
(Cons: Not fully exact, but that’s life!)

ñ Importance sampling seems to be the method of choice, but ...

Question: How generate a sample with a prescribed probability distribution?
Solution: Generate a Markov Chain

‹ Stochastic process without memory defined by

– State space (Here: Configuration space Γ “ tSi u
– Transition probability WΓÑΓ1
W W
; generates chain of states Γ0 ùñ Γ2 ùñ Γ2 ¨ ¨ ¨
or, respectively, chain of probabilities Pn pΓq
ř 1
(
Master equation: Pn`1 pΓq “ Pn pΓq` Γ1 W Γ1 ÑΓ Pn pΓ q ´ loooooomoooooon
looooooomooooooon WΓÑΓ1 Pn pΓq
flow in flow out
‹ For Markov chains with finite state space, one has a central limit theorem
(stated without proof)
If the Markov chain is irreducible, i.e., every state can be reached
from every other states (possibly by more than one step), then there
exists a unique stationary limit distribution P̄ pΓq with limnÑ8 Pn pΓq “
P̄ pΓq, independent of the initial distribution P0 pΓq
ÿ ÿ
‹ The stationary limit distribution fulfills WΓ1 ÑΓ P̄ pΓ1 q “ WΓÑΓ1 P̄ pΓq
Γ1 Γ1

‹ ñ Trick: (Metropolis, Rosenbluth, Teller)

Construct the transition function WΓÑΓ1 such that the limit distri-
bution is just the target distribution function. This can be achieved
with the following sufficient (but not necessary) conditions:
(i) irreducible (every state can be reached from every other state)
WΓÑΓ1 P̄ pΓ1 q
(ii) detailed balance “
WΓ1 ÑΓ P̄ pΓq
46 CHAPTER 3. THE ISING MODEL

One example of a popular implementation is the Metropolis algorithm

P̄ pΓq
WΓÑΓ1 “ NΓΓ1 minp1, q with NΓΓ1 “ NΓ1 Γ
P̄ pΓ1 q
‹ In particular, to obtain the Boltzmann distribution, one requires
WΓÑΓ1 1
“ e´βpH pΓq´H pΓ qq “ e´β∆E
WΓ1 ÑΓ
which in the Metropolis algorithm results in
WΓÑΓ1 “ NΓΓ1 minp1, e´β∆E q with NΓΓ1 “ NΓ1 Γ

Remarks:

• The Metropolis algorithm is the most popular algorithm, but every

other algorithm works too, as long as it fulfills the conditions (i)
and
ř (ii) or, instead of
ř (ii), at least the condition of ”global balance”
1
Γ1 WΓ1 ÑΓ P̄ pΓ q “ Γ WΓÑΓ1 P̄ pΓq
• It is not strictly necessary to target the Boltzmann distribution in
the Markov chain (i.e., choose P̄ pΓq „ e´βH . In some cases, it may
be more convenient to choose another target distributions and then
”reweight” the data when calculating the expectation values.
(; ”Reweighting” methods such as ”multicanonical” sampling,
”Wang/Landau” sampling, ”Metadynamics” etc.)

3.6.2 Examples of Monte Carlo algorithms

(a) Simple ”Single Flip” Metropolis algorithm

Algorithm
İ
§ (0) Initial configuration Γ “ tSi u
§
§ Ó
§
§ (i) Choose randomly a spin site j
§
§ Ó
§
§ (ii) Calculate energy difference ∆E between configuration Γ
§ §
and a configuration Γ˚ where Sj Ñ ´Sj
§ §
§ §
(2D cubic: ∆E “ 0, ˘2J, ˘4J)
§ §
§ §
§ đ
§
§
§ (iii) Pick a random number r P r0 : 1s
§
§ Ó
§
§ (iv)
§ Adopt Γ1 “ Γ˚ if r ă e´β∆E , otherwise keep Γ1 “ Γ
§ Ó
§
§ (v) New configuration Γ1

Remarks:
– Similar algorithms can be designed easily also for other systems.
– Close to the critical point, spin clusters become very large
; Dynamics become very slow (critical slowing down)
; Sampling becomes inefficient!
3.6. THE MONTE CARLO METHOD 47

(b) Ising model: Cluster algorithm (Wolff algorithm)

Algorithm
İ
§ (0) Initial configuration Γ
§
§ Ó
§
§ (i) Choose randomly a spin site j
§
§ Ó
§
§ (ii) Identify neighbors k of j with same spin direction Sk “ Sj
§ §
§ §
§ § and bonds pjkq that connect them
§ đ
§
§ (iii) Assign bond variables ujk to these bonds,
§ §
§ §
§ § choosing uij “ ´1 with probability e´2βJ .
Spins connected by bonds with uij “ 1 form a ”cluster”
§ §
§ §
§ đ
§
§
§ (iv) Identify neighbors of cluster with the same spin value.
§ §
§ § Assign bond variables to unoccupied connecting bonds.
§ §
§ § Extend cluster accordingly.
§ §
§ đ
§
§ (v) Continue until cluster can no longer grow.
§
§ Ó
§
§ (vi) Change sign of all spins in the cluster (Sk Ñ ´Sk )
§
§ Ó
§
§ (vi) New configuration Γ1

Proof that this algorithm fulfills detailed balance

Γ, Γ1 : Configurations, in which light pink cluster − − − − − − −
− + + + + − −
contains spins `1 or ´1, respectively − + − + + + −
Cluster is bounded by the contour L “ l` ` l − + + + + + −
with l` : Boundary to spins +1 (dashed line) − − + + + − −
+ − + + + − −
l´ : Boundary to spins -1 (solid line) + − − − − + −
Transition Γ Ñ Γ1 : Energy difference ∆E “ 2Jpl` ´ l´ q
`
WΓÑΓ1 “ W inside
loomo on ¨ e´2βJl
loomoon
Probability of having Probability that cluster
selected pink cluster does nothkkikkj
grow further
hkkikkj
´
WΓ1 ÑΓ M
“ Winside ¨ e´2βJl
` ´
ñ WΓÑΓ1 WΓ1 ÑΓ “ e´2βJpl ´l q “ e´2β∆E X

Remark: Global dynamics, totally ’unrealistic’,

but correlations break up much faster. ; more efficient sampling!
48 CHAPTER 3. THE ISING MODEL

3.6.3 ”Problems” with the Monte Carlo method

If one could invest an infinite amount of computing time, the Monte Carlo
method would be exact, on principle. Nevertheless, one has to apply caution
when analyzing the data.

‹ In fact, the computing time is never infinite

; Statistical error

‹ The systems have finite sie

; Systematic error
Causes problems in particular in the vicinity of critical points, where the
corrrelation length diverges.
Way out: Finite size scaling (Chapter 5)

‹ When using importance sampling, the entropy and free energy cannot be
calculated directly. (NB: Similar to experiments: Only observables can
be calculated!)
; Special methods must be developed, e.g., ”thermodynamic integra-
tion”, determination of free energy differences from histograms etc.
(See textbooks on simulation methods)