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