Exact Solution for Three-Dimensional Ising Model

Three-dimensional Ising model in zero external field is exactly solved by operator algebras, similar to the Onsager's approach in two dimensions. The partition function of the simple cubic crystal imposed by the periodic boundary condition along two directions and the screw boundary condition along the third direction is calculated rigorously. In the thermodynamic limit an integral replaces a sum in the formula of the partition function. The critical temperatures, at which order-disorder transitions in the infinite crystal occur along three axis directions, are determined. The analytical expressions for the internal energy and the specific heat are also presented.


I. INTRODUCTION
The exact solution for three-dimensional (3D) Ising model has been one of the greatest challenges to the physics community for decades. In 1925, Ising presented the simple statistical model in order to study the order-disorder transition in ferromagnets [1]. Subsequently the so-called Ising model has been widely applied in condensed matter physics. Unfortunately, one-dimensional Ising model has no phase transition at nonzero temperature. However, such systems could have a transition at nonzero temperature in higher dimensions [2]. In 1941, Kramers and Wannier located the critical point of two-dimensional (2D) Ising model at finite temperature by employing the dual transformation [3]. About two and a half years later Onsager solved exactly 2D Ising model by using an algebraic approach [4] and calculated the thermodynamic properties. Contrary to the continuous internal energy, the specific heat becomes infinite at the transition temperature T = T c given by the condition: sinh 2J kB Tc sinh 2J ′ kB Tc = 1, where (J ′ J) are the interaction energies along two perpendicular directions in a plane, respectively. Later, the partition function of 2D Ising model was also re-evaluated by a spinor analysis [5]. Up to now many 2D statistical systems have been exactly solved [6].
Since Onsager exactly solved 2D Ising model in 1944, much attention has been paid to the investigation of 3D Ising model. In Ref. [7], Griffiths presented the first rigorous proof of an order-disorder phase transition in 3D Ising model at finite temperature by extending the Peierls's argument in 2D case [2]. In 2000, Istrail proved that solving 3D Ising model on the lattice is an NP-complete problem [8]. We also note that the critical properties of 3D Ising model were widely explored by employing conformal field theories [9,10,11], self-consistent Ornstein-Zernike approximation [12], Renormalization group theory [13], Monte Carlo Simulations [14], the principal components analysis [15], and etc.. However, despite great efforts, 3D Ising model has not been solved exactly yet due to its complexity. It is out of question that an exact solution of 3D Ising model would be a huge jump forward, since it can be used to not only describe a broad class of phase transitions ranging from binary alloys, simple liquids and lattice gases to easy-axis magnets [16], but also verify the correctness of numerical simulations and finite-size scaling theory in three dimensions.
Because there is no dual transformation, the critical point of 3D Ising model cannot be fixed by such a symmetry. We also discover that it is impossible to write out the Hamiltonian along the third dimension of 3D Ising model with periodic boundary conditions (PBCs) in terms of the Onsager's operators. In addition, due to the existence of nonlocal rotation, 3D Ising model with PBCs seems not to be also solved by the spinor analysis [5]. Therefore, the key to solve 3D Ising model is to find out the operator expression of the interaction along the third dimension. We note that the transfer matrix in 3D Ising model is constructed by the spin configurations on a plane, which the boundary conditions (BCs) play an important role to solve exactly 3D Ising model. In this paper, we introduce a set of operators, which is similar to that in solving 2D Ising model [4]. Under suitable BCs, 3D Ising model with vanishing external field can be described by the operator algebras, and thus can be solved exactly.

II. THEORY
Consider a simple cubic lattice with l layers, n rows per layer, and m sites per row. Then the Hamiltonian of 3D Ising model is H = − m,n,l i,j,k=1 (J 1 σ z ijk σ z i+1kj + J 2 σ z ijk σ z ij+1k + Jσ z ijk s z ijk+1 ), where σ z ijk = ±1 is the spin on the site [ijk]. Assume that ν k labels the spin configurations in the kth layer, we have 1 ≤ ν k ≤ 2 mn . As a result, the and E 2 (ν k ) are the energies along two perpendicular directions in the kth layer, respectively, and E(ν k , ν k+1 ) is the energy between two adjacent layers. Now we define ( Here we use the periodic boundary conditions along both (010) and (001) directions and the screw boundary condition along the (100) direction for simplicity [3] (see Fig. 1). So the spin configurations along the X direction in a layer can be described by the spin variables σ z 1 , σ z 2 , · · · , σ z mn . Because the probability of a spin configuration is proportional to We note that V 1 , V 2 and V 3 are 2 mn -dimensional matrices, and both V 1 and V 2 are diagonal. Following Ref. [4], we obtain where , and H * = 1 2 ln coth H = tanh −1 (e −2H ). In order to diagonalize the transfer matrix V ≡ V 1 V 2 V 3 , following the Onsager's famous work in two dimensions, we first introduce the operators in spin space Γ along the X direction under the boundary conditions mentioned above. Here a, b = 1, 2, · · · , 2mn, σ x a , σ y a and σ z a are the Pauli matrices at site a, respectively. Then we have L 2 a,b = 1 and with Q ≡ nm a=1 σ x a = ±1. It is obvious that the period of L a,b is 2mn. We note that these operators L a,b are identical to P ab in Ref. [4] except mn replaces n.
H x and H z in the transfer matrix V can be expressed as Following Onsager's idea [4], we introduce the operators where x is an arbitrary index. Obviously, we have α −r = α r , β −r = −β r , β 0 = β mn = 0, γ −r = −γ r , and γ 0 = γ mn = 0. Eqs. (6) can be rewritten as where A s = mn a=1 L a,a+s and G s = 1 2 mn a=1 (L a,x L a+s,x − L x,a L x,a+s ). According to the orthogonal properties of the coefficients, we obtain From Eqs. (5)- (8), H x and H z have the expansions Because A mn+s = −QA s = −A s Q and G mn+s = −QG s = −G s Q, and combining with Eqs. (8), we have So we can investigate the algebra (8) with Q = 1 or -1 independently. However, we keep them together for convenience. In order to diagonalize the transfer matrix V , we must first determine the commutation relations among the operators α r , β r and γ r . Similar to those calculations in Ref. [4], we obtain Substituting Eqs. (8) into Eqs. (11), we arrive at where r = 1, 2, · · · , mn − 1, and all the other commutators vanish. Obviously, the algbra (12) is associated with the site r, and hence is local. Because α r , β r , and γ r obey the same commutation relations with −X r , −Y r , and −Z r in Ref. [4], we have the further relations We , where s = 1, 2, · · · , 2n, and When m = p = 1, Eqs. (14) recover the results in two dimensions [4]. It is obvious that A p,i and G p,j also satisfy the commutation relations (11).
We have obtained the expressions of H x and H z in terms of the operators α r , β r and γ r in the space Γ. In order to get the Hamiltonian in the third dimension, we project the operator algebra in the space Γ into the Y direction. Then we have m subspaces Γ p (p = 1, 2, · · · , m), in which the operator algebra with period 2n is same with that in Γ.
along the Y direction. Then we have A p,s = n a=1 L p a,a+s and G p, , which also obey the same commutation relations (11) and (12), similar to A p,s and G p,s . Then the Hamiltonian H y = m p=1 A p,1 . Because [L p a,a+s , L p b,b+s ] = 0 (see Fig. 2), we have [A p,s , A p,s ] = 0, which leads to A p,s ≡ A p,s due to their common local algebra (12). This is a renormalization of operator, which means that A p,s and A p,s have same eigenfunctions and eigenvalues in Γ p or Γ space. We note that V 2 is the transfer matrix along Y direction, which must be calculated in Γ rather than Γ p space by mapping A p,1 ≡ A p,1 in order to diagonalize total transfer matrix V . Therefore, we have Here, we would like to mention that H z = − m p=1 A p,0 ≡ −A 0 , which is same with that in (9). This means that when J 1 = 0, the Hamiltonian of 2D Ising model is recovered immediately.
V ] = 0, V and Q can be simultaneously diagonalized on the same basis. In other words, the eigenvalue problem of V can be classified by the value ±1 of Q.
The transfer matrix V with Eqs. (9) and (16) becomes where In order to obtain the eigenvalues of the transfer matrix V , we first diagonalize U r by employing the general unitary transformation: e i 2 ηr γr e ar (αr cos θr+βr sin θr) U r ×e −ar(αr cos θr+βr sin θr) e − i 2 ηrγr = e ξrαr .
Here θ r is an arbitrary constant and can be taken to be zero without loss of generality, and cosh ξ r = D r , sinh ξ r cos η r = A r , tanh(2a r ) = Cr Br , sinh ξ r sin η r = B r cosh(2a r ) − C r sinh(2a r ), where We note that D 2 r + C 2 r − A 2 r − B 2 r ≡ 1, which ensures that 3D Ising model can be solved exactly in the whole parameter space. When H 2 = 0(i.e.J 2 = 0) and n = 1, or H 1 = 0(i.e.J 1 = 0) and m = 1, we have a r = H * . So Eqs. (19) recover the Onsager's results in 2D Ising model [4].
Then the transfer matrix V has a diagonal form III. TRANSFORMATIONS
Following the procedure above, we can diagonalize the transfer matrix V , i.e. e mn−1 r=1 where sinh ξ r cos η * and We also have D 2 r + C 2 r − A * 2 r − B * 2 r ≡ 1.
The transfer matrix reads where Similarly, we have e mn−1 r=1 Here, and The identity D 2 r + C 2 r − A ′2 r − B ′2 r ≡ 1 also holds.
If H 2 = 0 or H 1 = 0, we obtain the critical temperature in 2D Ising model [3,4]. We note that the exact critical line (32) between the ferromagnetic and paramagnetic phases coincides completely with the result found in the domain wall analysis [17]. In the anisotropic limit, i.e. η = (H 1 + H 2 )/H → 0, the critical temperature determined by Eq.
(32) also agrees perfectly with the asymptotically exact value H = 2[lnη −1 − lnlnη −1 + 0(1)] −1 shown in Refs. [18,19]. When H 1 = H 2 = H, the critical value H c = J/(k B T c ) = 0.30468893, which is larger than the conjectured value about 0.2216546 from the previous numerical simulations [12,14]. We shall see from the analytical expressions (35) and (36) of the partition function per atom below that this discrepancy mainly comes from the oscillatory terms with respect to the system size m along X direction, which were not taken into account in all the previous numerical simulations.
We note that the thermodynamic properties of a large crystal are determined by the largest eigenvalue λ max of the transfer matrix V . Following Ref. [4], we have Here ∆ 1 = ∆ 3 = · · · = ∆ mn−1 = 1, which are same with the eigenvalues of the operators X r in Ref. [4]. We note that these two results above can be combined due to ξ −r = ξ r and ξ mn In order to calculate the partition function per atom λ ∞ = lim m,n→∞ (λ max ) 1 mn for the infinite crystal, we replace the sum in Eq. (34) by the integral where Similarly, the continuous A(ω), A * (ω), A ′ (ω), B(ω), B * (ω), B ′ (ω), C(ω), ξ m (ω), η(ω), η * (ω), and η ′ (ω) replace the discrete A r , A * r , A ′ r , B r , B * r , B ′ r , C r , ξ r , η r , η * r , and η ′ r , respectively, by letting ω = rπ mn . Here we emphasis that when H 2 = 0, or H 1 = 0, Eq. (35) is nothing but the Onsager's famous result in the 2D case [4]. We also note that very different from the 2D case, the partition function of 3D Ising model is oscillatory with m. Therefore, the conjectured values extrapolating to the infinite system in the numerical calculations seem to be inaccurate, and the 3D finite-size scaling theory must be modified.
We consider the special case of J 1 = J 2 , where the calculation of the thermodynamic functions can be simplified considerably. After integrating, Eq. (36) can be rewritten as where It is surprising that Eq. (40) is nothing but that in 2D Ising model with the interaction energies (J 1 , J 2D ) and H 2D = J2D kB T . Therefore, ln λ ∞ − 1 2 ln[2 sinh(2H)] in three dimensions can be obtained from ln λ 2D ∞ − 1 2 ln[2 sinh(2H 2D )] in two dimensions by taking the transformation (41). In other words, the thermodynamic properties of 3D Ising model originate from those in 2D case. We can also see from Eq. (41) that both 2D and 3D Ising systems approach simultaneously the critical point, i.e. H * 2D = H 1 and H * = 2H 1 . It is expected that the scaling laws near the critical point in two dimensions also hold in three dimensions [6].
The energy U and the specific heat C of 2D Ising model with the quadratic symmetry (i.e. H 1 = H 2D ) have been calculated analytically by Onsager and can be expressed in terms of the complete elliptic integrals [4]. The critical exponent associated with the specific heat α 2D = 0. Because 3D Ising model with the simple cubic symmetry (i.e. H 1 = H 2 = H) can be mapped exactly into 2D one by Eq. (41), the expressions of U and C in three dimensions have similar forms with those in two dimensions. So the critical exponent α 3D of the 3D Ising model is identical to α 2D , i.e. α 3D = 0. According to the scaling laws dν = 2 − α and µ + ν = 2 − α [6], we have ν 3D = 2 3 and µ 3D = 4 3 . Up to now, we have obtained the partition function per site and some physical quantities when the z axis is chosen as the transfer matrix direction. However, if the x(y) axis is parallel to the transfer matrix direction, the corresponding partition function per site can be achieved from Eqs. (35) and (36) by exchanging the interaction constants along the x(y) and z axes. Therefore, the total physical quantity in 3D Ising model, such as the free energy, the internal energy, the specific heat, and etc., can be calculated by taking the average over three directions. We note that the average of a physical quantity naturally holds for 2D Ising model.
Obviously, the difference between λ ∞ and λ p ∞ comes from the screw boundary condition along the X direction (see Fig. 1). We note that the k 2 term in Eqs. (44) and (45) vanishes, which can be seen as a feature of 3D Ising model.

VI. CONCLUSIONS
We have exactly solved 3D Ising model by an algebraic approach. The critical temperature T ic (i = 1, 2, 3), at which an order transition occurs, is determined. The expression of T ic is consistent with the exact formula in Ref. [17]. At T ic , the internal energy is continuous while the specific heat diverges. We note that if and only if the screw boundary condition along the (100) direction and the periodic boundary conditions along both (010) and (001) directions are imposed, the Onsager operators (15) along Y direction can form a closed Lie algebra, and then the Hamiltonian H y (16) is obtained rigorously. For PBCs, the Onsager operators along X or Y direction cannot construct a Lie algebra, and hence 3D Ising model is not solved exactly. Therefore, the numerical simulations on 3D finite Ising model with PBCs are unreliable due to the unclosed spin configurations on the transfer matrix plane. It is known that the BCs (the surface terms) affect heavily the results on small system, which lead to the different values extrapolating to the infinite system. However, the impact of the BCs on the critical temperatures can be neglected in the thermodynamic limit. Because the partition function per atom of 3D Ising model with H 1 = H 2 is equivalent to that of a 2D Ising model, the thermodynamic properties in three dimensions are highly correlated to those of 2D Ising system. When the interaction energy in the third dimension vanishes, the Onsager's exact solution of 2D Ising model is recovered immediately. This guarantees the correctness of the exact solution of 3D Ising model.