Solving the 2D Ising Model

Solving the 2D Ising Model Nina Kuklisova PHYS 35200 March 17, 2013 Abstract Finding a solution to the 2-dimensional Ising model has necessitated the effort of many physicists. In this expository paper, we consider the development of the analytic solution to the Ising Model and further evolution of the solution such as the applicability of Monte Carlo techniques. 1 Introduction The Ising model was one of the first attempts to simulate the structure of a physical ferromagnetic substance by approximatng the system as an N -dimensional periodic lattice. At each site of the lattice, there is a point with a spin up (+1) or down (-1). This problem was first stated by Lenz, who proposed it to his PhD student Ernst Ising. It turned out that the problem is much more complex than it would initially seem, and was introduced Ising’s thesis which appeared in 1925. In this thesis, Ising solved the problem in 1 dimension. A few years later, Peierls showed that a phase transition occurs in 2 dimensions. Still, it wasn’t until 1943 when Lars Onsager published his paper “Crystal Statistics I: A 2dimensional Model with an Order-Disorder Transition” [1], where the 2-dimensional model was fully solved with the assumption of a zero magnetic field. This solution was later simplified and became more useful for subsequent uses. This problem has applications beyond ferromagnetic materials: its solution is equivalent to that of lattice gas (a collection of atoms whose positions can take only discrete values) and binary alloy (alloy only composed of two types of atoms). Therefore, finding a solution for a ferromagnetic material would also reveal properties of these systems. 2 L Figure 1: Examples of an arrangement of spins in a 2D Ising model. As we can see, they can form clusters of the same spin. Source: ref.[8] 2 Systems studied in the Ising Model For a 2-dimensional lattice of spins in a magnetic field, each spin has 4 neighbors and is subject to the magnetic force. If we assume it only interacts with its nearest neighbors(not the diagonals), we find the Hamiltonian of this system to be � � H({si }) = −J si sj − h si , <i,j> i where J is the strength of interaction between the 2 spins, the first sum over all nearest neighbor pairs, and h is the external magnetic field, with the second sum over all spins. In this system, the partition function is � � � ZI = · ·... e−βH(si ) , s1 s2 sN which is generated by the interactions of each spin with its neighbors. Solving this problem consists in finding the transfer matrix of this distribution. This will allow predicting how will the spins in this system change when some parameters change. 3 3.1 Earlier Analytic solutions Onsager’s solution Onsager’s solution is valid in a situation where external magnetic field H = 0; for other values, no analytical solution has been found yet. 3 In his paper published in 1943, Onsager [1] was the first one to state the problem so that it could be solved through extensive use of properties of algebras. He defined an operator matrix that describes the addition of a new tier to the system by its effect on the interaction energy, which remained used in a simpler form in all subsequent solutions. Despite the complexity of these operators, he found the solution of this eigenvalue problem. First, he showed that in a 1-dimensional model, solving this problem is equivalent to finding the eigenvalue of this transition matrix. This solution revealed a certain type of operators. In the 2-dimensional problem, he showed that these operators are linearly combined to become a generating basis elements of an algebra, and this algebra can be decomposed into a direct product of quarternion algebras. From these, he determined the representation matrices of the operators. Afterwards, the quadratic equations that were used in this procedure were used to find the roots. So, these roots are determined from general properties of the matrix operators. Once the largest roots were found, they can be used for a qualitative investigation of the spectrum. The thermodynamic properties of the model are determined by finding the largest root of the characteristic equation as a function of the temperature. After finding the eigenvalue, he also found the critical temperature where the system starts behaving in a different way. At this “Curie point” of T = TC , the transition operators do not produce the same distribution. Here, the spontaneous magnetization disappears, and the specific heat becomes logarithmically infinite. At T > TC , the spontaneous magnetization is zero. At multiple places in his derivation derivation, he used results previously derived by Kramers and Wannier [3]. We omi them for simplicity here. This way, he found the in the thermodynamic limit to be A = −kBT ln Z Z = λN 1 ln λ = ln(2 cosh 2βJ) + pi � π 2 0 � 1 dw ln {1 + (1 − K 2 sin2 w)1/2 } 2 � where A is the Helmholtz free energy, Z the partition function of the system, and λ the eigenvalue. We do not show the full procedure here, because it later turned out that it can be simplified, and the main parts of the solution will be explained in the next section. Later, this solution was also used by other authors for slightly modified sysetms: R. Peierls [5] showed that for sufficiently low temperatures, the 2-dimensional result hols for 3 dimensions. Majumdar [6] discussed further analytic properties of this solution. Yang and Lee [10] were the first ones to notice that below critical temperature at vanishing external field, the free energy function does not have analytic behavior. 4 3.2 Kac and Ward’s Solution Kac and Ward’s solution [2] was more combinatorial and therefore was an important simplification. Its main idea was easier to use in later simulations. The layout of this problem is still the same as before, but its contribution was in noting some properties of the systems that their predecessors didn’t. Therefore, we will describe this solution in slightly more detail than Onsager’s. The partition function in this system in the multiplicative form with all the terms expanded is � Π<i,j> eβJsi sj , Z= si where J is the strength of interaction. Using hyperbolic identities, Kac and Ward noticed that eβJsi sj = cosh βJ(1 + t · si sj ), where t ≡ tanh βJ, and the partition function can be rewritten as � Π<i,j> (1 + t · si sj ). Z = (cosh βJ)N � N {si } In the term following (cosh βJ) {si } , the terms in the product denote the bonds in the lattice. What made this proof simpler was its treatment of the square lattice. Kac and Ward used a theorem claiming that if they count all the closed graphs on the square lattice, the closed graphs will cancel out when different ways of traversing them are considered. Further aspects of this are discussed in the next section. Subsequently, using procedures similar to those in Onsager’s solution, they found that Z = 2N (cosh2N J)Z � , where Z � = eN z , � � 1 2π 2π dζdη log det(1 − T M ), z = 2 0 (2π)2 0 det(1 − T M ) = (T 2 + 1)2 − 2T (1 − T 2 )(cos η + cos η) T = tanh J, where M is the transition matrix. This became the most frequently studied solution, which also Feynman used in his book [9], so it became the one that the majority of advanced physics students was aware of. Its “philosophy” suggests viewing the partition function as a Markovian process. 3.3 Vdovichenko’s Solution An alternative solution which did not gain so much popularity was derived by N. V. Vdovichenko [7]. The special part of his solution was that he manipulated the operations leading to the transfer matrix so that it lead to matrices in Toeplitz form, and then used their properties to get the solution. Although many authors did not use his solution, it benefited from the most straightforward derivation out of all those listed. 5 4 Numerical Solution 4.1 Markov Chain perspective Not only was Kac and Ward’s solution mathematically relatively simple, it also suggeted an application of Monte Carlo techniques. The reason for this is the structure of the main part of the solution. It can be modified into the following: Consider all the various possible spin configurations. In an Ising model with N spins, there are 2N spin configurations. The transition from one state (A) to another state (B) is characterized by a probability matrix πAB . This πAB is a 2N × 2N matrix, with most entries being zero. If the states A and B are the same or are different in only one state, then πAB > 0. If they have two or more different spins, πAB = 0. We call pacc AB the probability of acceptance of the spin change and call αAB the trial probability with the property αAB = αBA . Then, � this transition matrix can be written in a form πAB = αAB · pacc for A = � B, and π = 1 − AA AB B�=A πAB . If we assume EB > EA , then � � EB −EA acc and pacc pAB = exp − kB T BA = 1. Therefore, πAB πBA = αAB pacc AB αBA pacc BA � � A . = exp − EBkB−E T Once the equilibrium is reached at the state A, pA = const·e E � − EA − A condition, it follows that pA = Z1 · e kB T and Z = A e kB T . 4.2 EA BT −k . From the normalization Metropolis algorithm Now, the reader can see that this procedure can be simulated by the following algorithm called a Metropolis algorithm: 0. Start with a spin configuration {sk }. 1. Pick an arbitrary spin si . 2. Try to flip it: si := −si . 3. Calculate the energy change ∆E that this flip would cause. 4. If ∆E < 0, accept the flip. 5. If ∆E > 0, accept the flip with probability pacc = e−β∆E . 6. When a trial gets rejected, the spin is put back: si := −si . 7. If the highest possible number of iterations has not been reached yet, go to 1. Just as the procedure described in the previous subsection, this algorithm generates the canonical distribution. It also is a good example of a Markov Chain, a process in which each step only depends on the previous one. 6 4.3 Spin Flipping Simulation Since this algorithm is relatively simple, it inpired many authors to produce simulations of the system evolution. One of the most simulations that show the properties of the system in the highest complexity, but are still easy to understand for the viewer is http://physics.weber.edu/schroeder/software/demos/IsingModel.html. This simulation allows to see the influence of temperature on the system. Also, it meets the predictions given by the previous analytic solutions. 4.4 Magnetization and Specific Heat The change of the other parameters shows the effects of the critical temperature even better. The two most obvious cases are the Specific heat per spin (versus temperature) and Magnetization per spin (versus temperature). L Figure 2: Magnetization per spin of the 2D Ising model calculated by Monte Carlo simulation. Source: ref.[8] The most important property of this model is that at the critical temperature, the spontaneous magnetization (average spin value) decreases almost to 0. Even if the simulation may have not shown the magnetization to be exactly zero, it is more obviously seen from the Heat capacity function. Since the heat capacity is defined as � � �� ∂ AI 2 ∂ CI (H, T ) = −kT , ∂T ∂T kT where H is the magnetization and A the Helmoltz free energy of the system; it goes to infinity at the critical temperature. This simulation shows that even with the uncertainty range, at the critical temperature, the specific heat does tend to infinity. Since these properties fit well the theoretical predictions, the algorithms can also be used for further simulations of more complex systems that may not have such relatively simple analytic solutions. 7 L Figure 3: Specific heat per spin of the 2D Ising model calculated by Monte Carlo simulation (points are with error bars). Source: ref.[8] 5 Further results At this level of understanding the problem, a solution for a simplified model has been found, and it was also shown that it can be simulated through certain type of algorithms. This allows further ways of improving the solution. We will introduce each of them through showing simple results that it produced. 5.1 Probability: Dynamics of a Glauber System The 2-dimensional Ising model is also studied by probabilists. Under the assumptions of the model described in previous sections, the system qualifies as a Markov Chain. When studying Markov Chains, it is possible to find the mixing time, which is the time until the Markov Chain is “close” to its stationary distribution π (when another applying of the transition matrix would not change the distribution). In probability theory, this system is considered as a graph with n vertices and maximal degree ∆. This allowed probabilists to find boundaries of the mixing times of this system. For example D. Levin, Y. Peres and E. Wilmer [11] showed that � � n(log n + log(1/�)) tmix (�) ≤ 1 − ∆ tanh(β) whenever β < ∆−1 . For a complete proof, the reader can see [11]. These authors have also shown that on the n-cycle of an Ising Model, for any β > 0 and fixed � > 0, 1+o(1) 2(1−tanh(2β)) ≤ tmix (�) n log n ≤ 1+o(1) . 1−tanh(2β) Recent years have produced a lot of probability research in this area. Even if these estimations do not take into account all interactions in the system, they significantly improved our estimation of the result and effects of other interactions in the system. 8 5.2 Algorithms: Wolff algorithm As more complex simulations of this system are being developed, it is also necessary to develop algorithms which will compute the next step in the simulation more efficiently. When developing better algorithms, it was necessary to overcome the following big problem. In the equilibrium (depending on the other paramenters), the spins arrange themselves into a group of size η (called correlation length) and the reduced temperature (of distance C . Near the phase transition, the correlation length behaves as from TC ) is t = T −T TC η ∼ |t|−ν , where ν is called a critical exponent. Then, we can define correlation time τ , measured in steps per lattice site, as τ ∼ |t|−zν , where z is called the dynamic exponent. Consequently, τ ∼ ηz . This means divergence of the correlation length near the critical point, which is not possible in a system of finite size, such as in the Monte Carlo simulation. L Figure 4: Cutoff in the divergence of the correlation length η near the critical temperature for a system of finite size. The solid line describes the infinite system, the solid one the finite simulated one. Source: ref.[8] For this reason, other algorithms are being developed, which should allow a simulation that would be closer to real physical systems. The most fundamental step was the Wolff algorithm, developed in 1989. Its basic idea was to look for clusters of similarly oriented spins and then flip them in entirety all at once. This reduced the time it takes to simulate one correlation time. Then, the probability of adding a spin Padd becomes reduced (derived in [8]) to Padd = 1 − e−2βJ . With these modification, this algorithm follows a procedure similar to the Metropolis algorithm. 9 As we’ve seen, different algorithms are better for predicting different properties of the system, and therefore, new algorithms are being developed for producing better simulations of systems which behave similarly to the 2-dimensional Ising Model. 5.3 Physics: system modifications Still, physicists have been improving this model. One of these advances was done by V. S. Dotsenko and V. S. Dotsenko [12], who studied effects of inhomogenous impurities in the system. Their work showed that in a system containing a low amount of impurities, the singularity point of the specific heat changes from C ∼ ln(1/|τ |), where τ = T −TC , TC to C ∼ ln ln(1/|τ |), with τ << τi , where τi ∼ exp(−const.ci , and ci is the concentration of impurities. Furthermore, they showed that the spin correlation changes more dramatically. Although this problem has not been in the center of attention of physicists in the recent years, we can expect that the work with probabilists and computer scientists will soon bring results that could predict the behaviour of real-world system with a high precision. Conclusion Now, we can see that solving the 2-dimensional Ising model allowed developing of many strategies of finding a solution for a physical system. As a result, the current reserach in this area goes beyond the boundaries of physics, this problem still allows further development of applied mathematics. References [1] Crystal Statistics I. A Two-Dimensional Model with an Order-Disorder Transition - L. Onsager, Phys. Rev. 65 , 117, (1943) [2] A Combinatorial Solution of the Two-Dimensional Ising Model - M. Kac and J. C. Ward , Phys. Rev. 8, 1332, (1952) [3] Statistics of the Two-Dimensional Ferromagnet - H. A. Kramers and G. H. Wannier, Phys. Rev. 60, 252, 263 (1941) [4] On Ising’s Model of Ferromagnetism - R. Peierls, Proc. Cambridge Phil. Soc. 32, 477 (1936) [5] Peierl’s Proof of Spontaneous Magnetization in a Two-Dimensional Ising Ferromagnet R. B. Griffiths, Phys. Rev. 136, 2A (1964) [6] Analytic Properties of the Onsager Solution of the Ising Model - C. K. Majumdar, Phys. 10 Rev. 145, 158 (1966) [7] Spontaneous Magnetization of a Plane Dipole Lattice - N. V. Vdovichenko, J. Exptl. Theoret. Phys. (U.S.S.R.), 48, 526-530 (1965) [8] Monte Carlo Methods in Statistical Physics - M. E. J. Newman and G. T. Barkema, Oxford University Press (2001) [9] Statistical Mechanics: A Set of Lectures - R. P. Feynman, W.A.Benjamin Inc (1972) [10] Statistical Theory of Equations of Satte and Phase Transitions. I. Theory of Condensation - C. N. Yang and T. D. Lee, Phys, Rev. 87, 404-410 (1952) [11] Markov Chains and Mixing times - D. Levin, Y. Peres and E. Wilmer (2005) [12] Critical behavious of the phas etransition in the 2D Ising Model with Impurities - V. S. Dotsenko and V. S. Dotsenko, Advances in Physics 32, 2, 129-172 (1983) 11