PH4202: Advanced Statistical Mechanics
Exact solution of Two Dimensional Ising
Model
Author:
Vivek Pareek
10MS46
Supervisor:
Dr. Subhasis Sinha
May 2, 2014
Introduction
Two dimensional Ising model was first exactly solved by Lars Onsager in 1944. But the solution given
is very difficult to understand involve quite a bit of mathematics. Hence, subsequently new techniques
were developed to tackle the 2D Ising problem which could ensure less mathematical rigour and could
include some interesting physics. One such technique was developed by Schultz et.al in the year 1964
which involved transforming the 2D Ising problem to an interacting fermion problem. In this report
we have tried to understand and work out the physics behind such transformation. We would first
try to formulate the transfer matrix formalism for 1D Ising model and the would extend it to the 2D
Ising model. Thereafter we would introduce the interacting fermions picture and finally would calculate
exactly the critical temperature for phase transition for spontaneous magnetization.
One Dimensional Ising Model
As discussed in the Introduction we would formulate the transfer matrix formalism for 1D Ising model
in this section. For this let us consider the following Hamiltonian,
H = −J
N
X
i=1
σi σi+1 − h
N
X
σi
(1)
i=1
where σi vcan take value either +1 or −1. The partition function for the system is given by,
X
Z=
eβhσ1 eKσ1 σ2 eβhσ2 eKσ2 σ3 · · · eβhσN eKσN σ1
(2)
(σ)
where we have taken periodic boundary condition, σN +n = σn and K = βJ. We now try to simplify the
partion fuction by writing in the matrix form as follows,
X
Z=
hµ1 | V1 |µ2 i hµ2 | V1 |µ3 i · · · hµN | V1 |µ1 i
(3)
(µ=±1)
where |µi = |+1i or |−1i. The matrices V1 and V2 are as follows,
βh
e
0
= exp {βhσz }
V1 =
0 e−βh
K
e
e−K
= eK I + e−K σx
V2 =
e−K eK
we can also write V2 as
√
∗
V2 = A(K) exp {K ∗ σx }
(4)
(5)
(6)
where A(K) = 2 sinh 2K and tanh K = exp {−2K}. So now the partition function boils down to the
1/2
1/2
trace of V1 V2 or rather trace of V2 V1 V2 i.e,
X
Z =
hµ1 | V1 |µ2 i hµ2 | V2 |µ3 i · · · hµN | V2 |µ1 i
(µ=±1)
=
=
1/2
T r(V1 V2 )N = T r(V2
N
λN
1 + λ2
1/2 N
V1 V2
)
(7)
1/2
1/2
where λ1 , λ2 are eigenvalues of transfer matrix, V = V2 V1 V2 . We only consider the largest eigenvalue as it will have the maximum contribution to the free energy. Then we calculate the average of
magnetization for 1D Ising model and find no phase transition possible for it. But the formalism introduced above will be crucial to understand the 2D Ising model. Another important point to note here is
that the 1D problem now has been transformed effectively into a single site problem (zero dimension).
Similar dimension reduction would be seen as we proceed with the 2D Ising model.
1
Two dimensional Ising Model
We straight away jump into preparing the 2D Ising model for the introducing the Interacting fermions
formalism. The transfer matrix approach discussed in the previous section will help us in this direction.
Formalism
Consider a M × M square lattice with periodic boundary condition, The Hamiltonian of the system is
given by,
X
X
H = −J
σr,c σr+1,c − J
σr,c σr,c+1
(8)
r,c
r,c
where σr+M,c = σr,c+M = σr,c and r and c represents the rows and the columns. A pictorial representation is given below in figure 1.
σM,M
σ1,1
σ1,M
σM,1
Figure 1: M × M square lattice with periodic boundary condition
Now in analogy with the 1D case, we introduce the following 2M basis states
|µi = |µ1 i |µ2 i · · · |µM i
(9)
th
considering the index µi as the orientation of the
we can see that the matrix
o
n i Pspin in a column,
elements of the matrix V1 are given by V1 = exp K j σjZ σj+1,z . Similarly the matrix elements of
V2 are given by
h{µ}| V2 |{µ′ }i
=
hµM , µM −1 · · · , µ1 |
=
exp {(M − 2n)K}
M
Y
j=1
eK I + e−K σjx |µ′1 , µ′2 , · · · µ′M i
where n of the indices {µ′ } differ from the corresponding entries {µ}.
2
(10)
In the absence of external magnetic field, partition function can be written as ,
X
Z =
hµ1 | V1 |µ2 i hµ2 | V2 |µ3 i · · · hµN | V2 |µ1 i
{µ1 },{µ2 }···{µM }
1/2
1/2
T r(V1 V2 )N = T r(V2 V1 V2 )M
(11)
n
o
PM
M/2
where V2 = (2 sinh 2K)
exp K ∗ j=1 σjx . We now have to calculate the largest eigenvalue of the
following transfer matrix V ,
=
V
=
=
1/2
V2
1/2
V1 V2
(2 sinh 2K)
M/2
exp
K ∗ /2
M
X
j=1
σjx
exp
X
K
σjz σj+1,z
j
exp
K ∗ /2
M
X
j=1
σjx
(12)
But now the problem is that the terms above do not commute and moreover the matrix dimension
becomes infinity in the thermodynamic limit. Hence we map the proble to that of an interacting fermions
problem which will be shown in the next section.
Interacting Fermion Problem
For ease of calculation let us perform the following rotation of the spin operators,
σjz
σjx
→ −σjx
→ σjz ,
∀j
(13)
These rotations leave the eigenvalues invariant. Using σjz = 2σj+ σj− − 1 and σjx = σj+ + σj− , we arrive
at,
M
X
+
−
(14)
+ σj+1
V1 = exp K
σj+ + σj− σj+1
j=1
M
X
M/2
V2 = (2 sinh 2K)
exp 2K ∗
(15)
σj+ σj− − 1/2
j=1
It would be useful if we note down the commutation and anti-commutation relations of the pauli matrices.
They are as follows,
+ −
± ±
σ m , σm = σ z ,
σm , σn = 0,
m 6= n
(16)
+
−
σm
, σm
=
+
σm
1,
2
−
= σm
2
=0
(17)
We can see that the operators here are part fermionic and part bosonic. Hence we need to choose a
transformation which would ensure that either the operators become fermionic or bosonic. We choose
Jordon-Weigner transformation which converts pauli operators to fermionic operators. The reason for
this choice is that bosons do not obey pauli exclusion principle and hence the calculation becomes more
difficult. Hence the transformation is given by,
)
( j−1
X
+
†
(18)
Cm Cm Cj†
σj = exp πi
m=1
σj−
=
(
Cj exp −πi
3
j−1
X
m=1
†
Cm
Cm
)
(19)
where C dagger and C follow
i
†
Cj† , Cm
+
†
C j , Cm +
h
=
[Cj , Cm ]+ = 0
(20)
=
δjm
(21)
We can now transform V2 using the above relations and the transformation to get,
M
X
M/2
Cj† Cj − 1/2
V2 = (2 sinh 2K)
exp 2K ∗
(22)
j=1
We would have to work a little bit for finding out V1 . We can see,
†
+
†
−
= Cj − Cj Cj+1
+ σj+1
+ Cj+1 ,
σj+ + σj− σj+1
j 6= m
(23)
For the case j = m, we have to be a bit careful. This difficulty is due to periodic boundary condition.
Finally we have for j = m, we have,
+
−
σm
+ σm
†
+
−
†
= −(−1)n Cm
− Cm Cm+1
+ Cm+1 ,
σm+1
+ σm+1
where n =
M
X
Cj† Cj
(24)
j=1
n is the total fermion number operator. It is important to note that (−1)n commutes with different
terms. In order to stitch these terms, we next divide the space into two subspace- even and odd. V1 can
be written in simplified form as
†
+ Cj+1
(25)
V1 = Cj† − Cj Cj+1
with boundary terms as
†
†
1) Anti-cyclic: CM +1 ≡ −C1 ; CM
+1 ≡ −C1 , f or n even
†
†
2) Cyclic: CM +1 ≡ C1 ; CM
+1 ≡ C1 , f or n odd
by doing this we have now achieved transnational in variance. Now we carry out the following canonical
transformation,
aq
=
a†q
=
M
1 X
√
Cj e−iqj
M j=1
M
1 X † iqj
√
Cj e
M j=1
(26)
(27)
hence the inverse transformation is given by
Cj
=
Cj†
=
1 X
√
aq eiqj
M q
M
1 X † −iqj
√
aq e
M j=1
(28)
(29)
Now being canonical transformation, they preserve the commutation rules of fermions. In order to
reproduce boundary conditions discussed above, we take q = jπ/M with
j
j
=
=
±1, ±3, · · · ± (M − 1);
0, ±2, ±4, · · · ± (M );
4
f or n even
f or n odd
(30)
(31)
Here we can choose M to be even without any loss of generality. We need to stop here a while and think
why have we done all this? Answer to this question is, by doing this we have ensured that the transfer
matrix in the momentum mode is now just n multiplications of decoupled mode. Hence, by working on
only one mode we can extract the entire behaviour of the transfer matrix. Thus, we see
)
(
M
X
M/2
†
†
∗
(32)
aq aq + a−q a−q − 1
V2 = (2 sinh 2K)
exp 2K
q>0
=
(2 sinh 2K)
M/2
Y
V2q
(33)
q>0
along with,
V1
=
=
(
exp 2K
Y
Xh
q>0
cos q a†q aq + a†−q a−q − i sin q a†q a†−q + aq a−q
V1q
i
)
(34)
(35)
q>0
The above simplification can be achieved exploiting the fact that bilinear operators with different wave
vectors commute. Because of pauli’s exclusion principle we can see that the eigenvalues of the transfer
matrix can now be written as the product of eigenvalues of atmost 4 × 4 matrices.For the case of odd n
we also need operators V1q and V2q for q = 0 and q = π. These are,
n
o
n
o
V10 = exp 2Ka†0 a0 ; V20 = exp 2K ∗ a†0 a0 − 1/2
(36)
V1π = exp −2Ka†π aπ ; V2π = exp 2K ∗ a†π aπ − 1/2
(37)
We can see they are already in diagonal form hence we do not have to work extra for them. But we still
need to tackle the case of finding the eigenvalue for the case q 6= 0, ±π. This we do in the next section.
Calculation of eigenvalues
We have found the form of V1 and V2 in the previous section which then, could be written as the product
of eigenvalues of 4 × 4 matrices. In this section we are going to find the eigenvalue of these matrices
but only for the momentum modes q 6= 0, ±π. Hence we are looking for the eigenvalue of the following
operator(for a particular momentum mode q),
1/2
1/2
Vq = V2q V1q V2q
(38)
The complete transfer matrix for the system would then be just VqM . Since, we are dealing with fermions
we have only four possible states, namely, |0i, a†q |0i, a†−q |0i and a†q a†−q |0i. For convenience we would
be using the following terminology |0i, |1i, |−1i and |2i respectively for the above described states. The
state |0i is defined by aq |0i = a−q |0i = 0, which is the zero particle state, |1i and |−1i are 1 particle
states with momentum q and −q respectively and |2i is 2 particle state, one with momentum q and
another with momentum −q.
Important observation here is that the above states are already an eigenstate for the matrix V2 . Hence
V2 in the basis of these four states is given as,
−K ∗
e
0 0
0
0
1 0
0
V2 =
(39)
0
0 1
0
K∗
0
0 0 e
We can see that the |1i and |−1i are eigenstates of V1q with eigenvalue exp 2K cos q but we also have
off diagonal terms in V1 . If we look closely then we can see that these off diagonal terms appear only
between the states |0i and |2i, i.e, states differing by two fermion number. Now the problem boils down
to finding the eigenvalues of V1q in the basis |0i and |2i.
5
We will only concentrate on V1 as we have seen earlier that V2 is diagonal. For computing the matrix
V1q in the basis of |0i and |2i, we introduce the pair annihilation and creation operators, b−
q = aq a−q
†
†
†
and bq = a−q aq . These operators are analogous to the Cooper pairs creation and annihilation operators
we see in the BCS theory of superconductivity. These operators follow the simple rules of spin lowering
and raising operators and are represented by pauli spin matrices. Consider the following relations,
a†q aq + a†−q a−q
=
z
2b†q b−
q = bq + 1
(40)
a†−q a†q
=
bxq
(41)
aq a−q +
Now we can write V1q as,
V1q
exp[2Kcos q(bzq + 1) + sin q(bxq )]
exp[2K cos q] exp[2Kbzq′ ],
=
=
where bzq′ = bzq cos q + bxq sin q and (bzq′ )2 = 1
exp[2K cos q][cosh 2K + bzq′ sinh 2K]
cosh 2K + sinh 2K cos q
sinh 2K sin q
exp[2K cos q]
sinh 2K sin q
cosh 2K − sinh 2K cos q
=
=
(42)
So finally we get the following form for Vq in |0i and |2i basis,
Vq
1/2
1/2
=
V2q V1q V2q
=
e−K
exp[2K cos q]
0
∗
0
eK
∗
cosh 2K + sinh 2K cos q
sinh 2K sin q
sinh 2K sin q
cosh 2K − sinh 2K cos q
e−K
0
∗
0
eK
∗
We have to find the eigenvalue of the Vq , which in our case, is done by the usual method1 . In order to
have some comparative idea about the eigenvalues(as are only interested in the largest eigenvalue), we
write the eigenvalues in the following form,
λ±
q = exp[2K cos q ± ǫ(q)]
(43)
cosh ǫ(q) = cosh 2K cosh 2K ∗ + cos q sinh 2K sinh 2K ∗
(44)
where ǫ(q) was found out to be
by convention we choose ǫ(q) to be positive. We note that the minimum of ǫ(q) occurs for q → π, and
maximum occurs for q → 0,
ǫmin
=
ǫmax
=
lim ǫ(q) = 2|K − K ∗ |
(45)
lim ǫ(q) = 2(K + K ∗ )
(46)
q→π
q→0
Until now we have considered only the case of q 6= 0, ±π, and we can see that λ+
q is the largest eigenvalue
among all. Hence, for the case where we have to consider the case q = 0, ±π , we have to fall back to
the Vq=0,π described in the earlier section.
We now have to consider two different cases
1: Subspace containing even number of fermions.
2: Subspace containing odd number of fermions.
1 One
can also introduce Bogoulibov Transformations and try to find the normal modes.
6
Case 1: Even Fermions: As discussed earlier, the allowed values of q for this case is q = ±1π/M, ±3π/M...±
(M − 1)π/M . We can see that this case does not include the case q = 0, ±π. The largest eigenvalue for
each momentum mode is λ+
q . Hence, largest eigenvalue in the even subspace, given by Λe , is
Λe
=
(2 sinh 2K)M/2
Y
λ+
q
q>0
=
=
(2 sinh 2K)
(2 sinh 2K)
M/2
M/2
exp
"
"
X
[2K cos q + ǫ(q)]
q>0
1X
exp
ǫ(q)
2 q
#
#
(47)
(48)
It is important to note that in equation (11), the summation runs over all the momentum mode, i.e, we
have also included the contribution from q = 0 & ± π as their contribution in this case is zero. Moreover
we can see that ǫ(q) is an even function, hence a factor of 1/2 appears in equation (11) and we know
P
q cos q = 0.
Case 2: Odd Fermions: The allowed values of q for this case is q = 0, ±2π/M, ±4π/M... ± π. For
momentum modes q 6= 0, ±π, λ+
q is still the largest eigenvalue. We have to incorporate the case of
q = 0, ±π separately. As discussed in the previous section we have to deal with the following form of
Vq′ s for q = 0, ±π,
V10 = exp(2Ka†0 a0 )
V20 = exp[2K ∗ (a†0 a0 − 1/2)]
(49)
V1π = exp(−2Ka†π aπ )
V2π = exp[2K ∗ (a†π aπ − 1/2)]
(50)
It is important to note that the case for q 6= 0, ±π includes even number of fermions. Now we have
additional momentum mode q = 0 and q = π available. We only need to fill one momentum mode to
get odd fermions. Now the important question is, which one should we choose? The answer to this is
that, we have already seen from equation (8) and (9) that the maximum contribution of ǫ(q) to λ+
q is
for the case q → 0 and minimum contribution is for the case q → π. Hence, we can afford to neglect the
contribution from q = π. We can see from equation (12) and (13) the contribution is 2K. Hence the
largest eigenvalue for the cse of odd fermions is given by,
X
1
Λo = (2 sinh 2K)M/2 exp 2K +
(51)
ǫ(q)
2
q6=0,π
For thermodynamic limit we can see that the discrete sum over momentum modes would be replaced
by an integration over all momentum modes, which would be accompanied by a volume term, i.e, M in
this case. Hence in the thermodynamic limit Λe and Λo must converge to a single eigenvalue. Consider
the following,
1
1
lim ǫ(q) + lim ǫ(q)
2 q→0
2 q→π
=
|K − K ∗ | + (K + K∗)
=
=
2K, K > K ∗
2K ∗ , K < K ∗
(52)
(53)
Now, we see that for K > K ∗ , we can have Λo = Λe . This is because the extra 2K in equation (14)
can be incorporated within the summation as it is the contribution from q = 0, π, which would make
the summation in equation (14) over all the momentum modes and in the thermodynamic limit, we will
have degeneracy. This degeneracy, at a finite temperature marks the spontaneous magnetization in the
system. Hence we can see that Phase transition indeed occurs for 2D Ising Model! This degeneracy is
7
the signature for the case T < Tc . For K < K ∗ , this degeneracy is lifted and we have Λe > Λo . Hence
there is no spontaneous magnetization in the system. This is the signature for the case T > Tc . Hence
we can see that critical behavior occurs for K = K ∗ . Hence, from the definition of K ∗ we have,
tanh K ∗
=
exp(−2K)
K∗
=
tanh−1 (exp[−2K])
8
(54)
K = K∗
K∗
K ∗ = tanh−1 (exp[−2K])
K = 0.441
K
Figure 2: Plotting K ∗ = tanh−1 (exp[−2K]) and K = K ∗
From figure 1, we get
⇒
⇒
K = 0.4414
J
= 0.4414
KB T c
Tc = 2.26
J
KB
From Mean field theory calculations we had found that Tc =
solution we have obtained here.
(55)
2J
which is pretty close to the exact
KB
Acknowledgment
I would like to thank Dr. Subhasis Sinha for his guidance in this project as well as I would like to thank
Deep Chatterjee , Mohit Pandey and all my friends with whom I have had long and fruitful discussions,
which has led to a successful understanding of the model presented.
References
• Equilibrium Statistical Mechanics by Michael Plischke and Birger Bergersen
• Review of Modern Physics 1964 Two dimensional Ising model as a soluble problem of Many fermoins by T.D.Schulz, D.C.Mattis and E.H.Lieb.
9