1. Introduction
The exact solution for the 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 to study the order–disorder transition in ferromagnets [
1]. Subsequently the so-called Ising model has been widely applied in condensed matter physics. Unfortunately, the 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 the two-dimensional (2D) Ising model at finite temperature by employing the dual transformation [
3]. About two and a half years later Onsager solved the exactly 2D Ising model using an algebraic approach [
4] and calculated the thermodynamic properties. Contrary to continuous internal energy, specific heat becomes infinite at the transition temperature
given by the condition:
, where
are the interaction energies along two perpendicular directions in a plane, respectively. Later, the partition function of the 2D Ising model was also re-evaluated by spinor analysis [
5]. Up to now, many 2D statistical systems have been exactly solved [
6].
Since Onsager exactly solved the 2D Ising model in 1944, much attention has been paid to the investigation of the 3D Ising model. In Ref. [
7], Griffiths presented the first rigorous proof of an order–disorder phase transition in the 3D Ising model at finite temperature by extending Peierls’s argument in the 2D case [
2]. In 2000, Istrail proved that solving the 3D Ising model on the lattice is an NP-complete problem [
8]. We also note that the critical properties of the 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], principal component analysis [
15], etc. However, despite great efforts, the 3D Ising model has not been solved exactly yet due to its complexity. There is no question that an exact solution of the 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 the 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 the 3D Ising model with periodic boundary conditions (PBCs) in terms of Onsager’s operators. In addition, due to the existence of nonlocal rotation, the 3D Ising model with PBCs seems not to be also solved by the spinor analysis [
5]. Therefore, the key to solve the 3D Ising model is to find out the operator expression of the interaction along the third dimension. We note that the transfer matrix in the 3D Ising model is constructed by the spin configurations on a plane, in which the boundary conditions (BCs) play an important role to solve exactly the 3D Ising model. In this paper, we introduce a set of operators, which is similar to that in solving the 2D Ising model [
4]. Under suitable BCs, the 3D Ising model with vanishing external field can be described by the operator algebras, and thus can be solved exactly.
2. 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
, where
is the spin on the site
. Assume that
labels the spin configurations in the
kth layer, we have
. As a result, the energy of a spin configuration of the crystal
, where
and
are the energies along two perpendicular directions in the
kth layer, respectively, and
is the energy between two adjacent layers. Now we define
and
. Here we use the periodic boundary conditions along both
and
directions and the screw boundary condition along the
direction for simplicity [
3] (see
Figure 1). Therefore, the spin configurations along the
direction in a layer can be described by the spin variables
. Because the probability of a spin configuration is proportional to
, the partition function of the 3D Ising model is
We note that
,
and
are
-dimensional matrices, and both
and
are diagonal. Following Ref. [
4], we obtain
where
, and
.
To diagonalize the transfer matrix
, following 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
,
,
and
are the Pauli matrices at site
a, respectively. Then we have
and
with
. It is obvious that the period of
is 2
. We note that these operators
are identical to
in Ref. [
4] except
replaces
n.
and
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
,
,
,
, and
. Equation (
6) can be rewritten as
where
and
. According to the orthogonal properties of the coefficients, we obtain
From Equations (5)–(8),
and
have the expansions
Because
and
, and combining with Equation (
8), we have
When
,
while
if
. Therefore, we can investigate the algebra (8) with
or −1 independently. However, we keep them together for convenience. To diagonalize the transfer matrix
V, we must first determine the commutation relations among the operators
,
and
. Similar to those calculations in Ref. [
4], we obtain
Substituting Equation (
8) into Equation (
11), we arrive at
where
, and all the other commutators vanish. Obviously, the algebra in (12) is associated with the site
r, and hence is local. Because
,
, and
obey the same commutation relations with
,
, and
in Ref. [
4], we have the further relations
We note that
and
, where
, and
When
, Equation (
14) recovers the results in two dimensions [
4]. It is obvious that
and
also satisfy the commutation relations (11). When
,
.
We have obtained the expressions of
and
in terms of the operators
,
and
in the space
. To obtain the Hamiltonian in the third dimension, we project the operator algebra in the space
into the
direction. Then we have
m subspaces
, in which the operator algebra with period
is same with that in
. In
, we define
along the
direction. Then we have
and
, which also obey the same commutation relations (11) and (12), similar to
and
. Then the Hamiltonian
.
Because
(see
Figure 2), we have
, which leads to
due to their common local algebra (12).
This is a renormalization of operator, which means that and have same eigenfunctions and eigenvalues in or Γ space. We note that
is the transfer matrix along
direction, which must be calculated in
rather than
space by mapping
to diagonalize total transfer matrix
V. Therefore, we have
Here, we would like to mention that , which is same with that in (9). This means that when , the Hamiltonian of the 2D Ising model is recovered immediately.
Because , 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 of Q.
The transfer matrix
V with Equations (9) and (16) becomes
where
To obtain the eigenvalues of the transfer matrix
V, we first diagonalize
by employing the general unitary transformation:
Here
is an arbitrary constant and can be taken to be zero without loss of generality, and
where
We note that
, which ensures that the 3D Ising model can be solved exactly in the whole parameter space. When
(i.e.,
) and
, or
(i.e.,
) and
, we have
. Therefore, Equation (
19) recovers Onsager’s results in the 2D Ising model [
4].
Then the transfer matrix
V has a diagonal form
4. Results
Because
have the common eigenvectors
with the corresponding eigenvalues
, from Equation (
20), we have
, where
At the critical point, we have
[
4]. This leads to a critical temperature
given by the condition
If
or
, we obtain the critical temperature in the 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.,
, the critical temperature determined by Equation (
32) also agrees perfectly with the asymptotically exact value
shown in Refs. [
18,
19].
When
, the critical value
, 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
of the transfer matrix
V. Following Ref. [
4], we have
Here
, which are same as the eigenvalues of the operators
in Ref. [
4]. We note that these two results above can be combined due to
and
. Therefore, Equation (
33) has the compact form
To calculate the partition function per atom
for the infinite crystal, we replace the sum in Equation (
34) by the integral
where
Similarly, the continuous
,
,
,
,
,
,
,
,
,
, and
replace the discrete
,
,
,
,
,
,
,
,
,
, and
, respectively, by letting
. Here we emphasize that when
, or
, Equation (
35) is nothing but Onsager’s famous result in the 2D case [
4]. We also note that very different from the 2D case, the partition function of the 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.
For a crystal of
, the free energy
the internal energy
and the specific heat
We note that at the critical point, . However, . Therefore, we can see from Equations (37) and (38) that at the critical point, the internal energy U is continuous while the specific heat C becomes infinite, similar to the 2D case.
We consider the special case of
, where the calculation of the thermodynamic functions can be simplified considerably. After integrating, Equation (
36) can be rewritten as
where
It is surprising that Equation (
40) is nothing but that in 2D Ising model with the interaction energies
and
. Therefore,
in three dimensions can be obtained from
in two dimensions by taking the transformation (41). In other words, the thermodynamic properties of the 3D Ising model originate from those in 2D case. We can also see from Equation (
41) that both the 2D and 3D Ising systems approach simultaneously the critical point, i.e.,
and
. 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 the 2D Ising model with the quadratic symmetry (i.e.,
) 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
. Because the 3D Ising model with the simple cubic symmetry (i.e.,
) can be mapped exactly into the 2D one by Equation (
41), the expressions of
U and
C in three dimensions have similar forms with those in two dimensions. Therefore, the critical exponent
of the 3D Ising model is identical to
, i.e.
. According to the scaling laws
and
[
6], we have
and
.
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
axis is parallel to the transfer matrix direction, the corresponding partition function per site can be achieved from Equation (
35) and (36) by exchanging the interaction constants along the
and
z axes. Therefore, the total physical quantity in the 3D Ising model, such as free energy, internal energy, specific heat, etc., can be calculated by taking the average over three directions. We note that the average of a physical quantity naturally holds for the 2D Ising model.
6. Conclusions
We have exactly solved the 3D Ising model by an algebraic approach. The critical temperature
, at which an order transition occurs, is determined. The expression of
is consistent with the exact formula in Ref. [
17]. At
, the internal energy is continuous while the specific heat diverges. We note that if and only if the screw boundary condition along the
direction and the periodic boundary conditions along both
and
directions are imposed, the Onsager operators (15) along
Y direction can form a closed Lie algebra, and then the Hamiltonian
(16) is obtained rigorously. For PBCs, the Onsager operators along
X or
Y direction cannot construct a Lie algebra, and hence the 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) heavily affect the results on the 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 the 3D Ising model with
is equivalent to that of a the 2D Ising model, the thermodynamic properties in three dimensions are highly correlated with those of the 2D Ising system. When the interaction energy in the third dimension vanishes, Onsager’s exact solution of the 2D Ising model is recovered immediately. This guarantees the correctness of the exact solution of the 3D Ising model.