Exact Solution for Three-Dimensional Ising Model

Abstract

The three-dimensional Ising model in a 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.

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 $T=T_c$ given by the condition: $sinh\frac{2J}{k_B{T}_c}sinh\frac{2J^{\prime }}{k_B{T}_c}=1$, where $\left(J^{\prime }J\right)$ 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 $H=-{\sum }_{i,j,k=1}^{m,n,l}(J_1{\sigma }_{ijk}^z{\sigma }_{i+1kj}^z+J_2{\sigma }_{ijk}^z{\sigma }_{ij+1k}^z+J{\sigma }_{ijk}^z{s}_{ijk+1}^z)$, where ${\sigma }_{ijk}^z=\pm 1$ is the spin on the site $\left[ijk\right]$. Assume that ${\nu }_k$ labels the spin configurations in the kth layer, we have $1\le {\nu }_k\le 2^{mn}$. As a result, the energy of a spin configuration of the crystal $E_{sc}={\sum }_{k=1}^l{E}_1\left({\nu }_k\right)+{\sum }_{k=1}^l{E}_2\left({\nu }_k\right)+{\sum }_{k=1}^{l}E({\nu }_k,{\nu }_{k+1})$, where $E_1\left({\nu }_k\right)$ and $E_2\left({\nu }_k\right)$ are the energies along two perpendicular directions in the kth layer, respectively, and $E({\nu }_k,{\nu }_{k+1})$ is the energy between two adjacent layers. Now we define ${\left(V_1{V}_2\right)}_{{\nu }_k{\nu }_k}={\sum }_{{\nu }_k^{\prime }}{\left(V_1\right)}_{{\nu }_k{\nu }_k^{\prime }}{\left(V_2\right)}_{{\nu }_k^{\prime }{\nu }_k}={\left(V_1\right)}_{{\nu }_k{\nu }_k}{\left(V_2\right)}_{{\nu }_k{\nu }_k}\equiv exp[-E_1\left({\nu }_k\right)/\left(k_{B}T\right)]\times exp[-E_2\left({\nu }_k\right)/\left(k_{B}T\right)]$ and ${\left(V_3\right)}_{{\nu }_k{\nu }_{k+1}}\equiv exp[-E({\nu }_k,{\nu }_{k+1})/\left(k_{B}T\right)]$. Here we use the periodic boundary conditions along both $\left(010\right)$ and $\left(001\right)$ directions and the screw boundary condition along the $\left(100\right)$ direction for simplicity [3] (see Figure 1). Therefore, the spin configurations along the $\mathbf{X}$ direction in a layer can be described by the spin variables ${\sigma }_1^z,{\sigma }_2^z,\cdots ,{\sigma }_{mn}^z$. Because the probability of a spin configuration is proportional to $exp[-E_{sc}/\left(k_{B}T\right)]={\left(V_1{V}_2\right)}_{{\nu }_1{\nu }_1}{\left(V_3\right)}_{{\nu }_1{\nu }_2}{\left(V_1{V}_2\right)}_{{\nu }_2{\nu }_2}{\left(V_3\right)}_{{\nu }_2{\nu }_3}\cdots {\left(V_1{V}_2\right)}_{{\nu }_l{\nu }_l}{\left(V_3\right)}_{{\nu }_l{\nu }_1}$, the partition function of the 3D Ising model is

$$\begin{array}{ccc}Z\hfill & =\hfill & {\sum }_{{\nu }_1,{\nu }_2,\cdots ,{\nu }_l}{\left(V_1{V}_2\right)}_{{\nu }_1{\nu }_1}{\left(V_3\right)}_{{\nu }_1{\nu }_2}\cdots {\left(V_1{V}_2\right)}_{{\nu }_l{\nu }_l}{\left(V_3\right)}_{{\nu }_l{\nu }_1}\hfill \\ & \equiv \hfill & \mathrm{tr}{\left(V_1{V}_2{V}_3\right)}^l.\hfill \end{array}$$

(1)

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

$$\begin{array}{ccc}{V}_1\hfill & =\hfill & exp\left(H_1\sum _{\tau =1}^{mn}{\sigma }_{\tau }^z{\sigma }_{\tau +1}^z\right)\equiv exp\left(H_1{H}_x\right),\hfill \\ V_2\hfill & =\hfill & exp\left(H_2\sum _{\tau =1}^{mn}{\sigma }_{\tau }^z{\sigma }_{\tau +m}^z\right)\equiv exp\left(H_2{H}_y\right),\hfill \\ V_3\hfill & =\hfill & {[2sinh\left(2H\right)]}^{mn/2}exp\left(H^{*}\sum _{\tau =1}^{mn}{\sigma }_{\tau }^x\right)\hfill \\ & \equiv \hfill & {[2sinh\left(2H\right)]}^{mn/2}exp\left(H^{*}{H}_z\right),\hfill \end{array}$$

(2)

where $H_1=J_1/\left(k_{B}T\right),H_2=J_2/\left(k_{B}T\right),H=J/\left(k_{B}T\right)$, and $H^{*}=\frac{1}{2}lncothH={tanh}^{-1}$$\left(e^{-2H}\right)$.

To diagonalize the transfer matrix $V\equiv V_1{V}_2{V}_3$, following Onsager’s famous work in two dimensions, we first introduce the operators

$$L_{a,a}=-{\sigma }_a^x,L_{a,b}={\sigma }_a^z{\sigma }_{a+1}^x{\sigma }_{a+2}^x\cdots {\sigma }_{b-1}^x{\sigma }_b^z$$

(3)

in spin space $\Gamma$ along the X direction under the boundary conditions mentioned above. Here $a,b=1,2,\cdots ,2mn$, ${\sigma }_a^x$, ${\sigma }_a^y$ and ${\sigma }_a^z$ are the Pauli matrices at site a, respectively. Then we have $L_{a,b}^2=1$ and

$$L_{a,b+mn}=L_{a+mn,b}=-Q{L}_{a,b}=-L_{a,b}Q$$

(4)

with $Q\equiv {\prod }_{a=1}^{nm}{\sigma }_a^x=\pm 1$. It is obvious that the period of $L_{a,b}$ is 2$mn$. 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

$$\begin{array}{c}{H}_x=\sum _{a=1}^{mn}{L}_{a,a+1},\hfill \\ H_z=\sum _{a=1}^{mn}{\sigma }_a^x=-\sum _{a=1}^{mn}{L}_{a,a}.\hfill \end{array}$$

(5)

Following Onsager’s idea [4], we introduce the operators

$$\begin{array}{ccc}{\alpha }_r\hfill & =\hfill & -\frac{1}{4mn}\sum _{a,b=1}^{2mn}{L}_{a,b}cos\frac{(a-b)r\pi }{mn},\hfill \\ {\beta }_r\hfill & =\hfill & -\frac{1}{4mn}\sum _{a,b=1}^{2mn}{L}_{a,b}sin\frac{(a-b)r\pi }{mn},\hfill \\ {\gamma }_r\hfill & =\hfill & \frac{i}{8mn}\sum _{a,b=1}^{2mn}(L_{a,x}{L}_{b,x}-L_{x,a}{L}_{x,b})sin\frac{(a-b)r\pi }{mn}\hfill \end{array}$$

(6)

where x is an arbitrary index. Obviously, we have ${\alpha }_{-r}={\alpha }_r$, ${\beta }_{-r}=-{\beta }_r$, ${\beta }_0={\beta }_{mn}=0$, ${\gamma }_{-r}=-{\gamma }_r$, and ${\gamma }_0={\gamma }_{mn}=0$. Equation (6) can be rewritten as

$$\begin{array}{ccc}{\alpha }_r\hfill & =\hfill & -\frac{1}{2mn}\sum _{s=1}^{2mn}{A}_{s}cos\frac{rs\pi }{mn},\hfill \\ {\beta }_r\hfill & =\hfill & \frac{1}{2mn}\sum _{s=1}^{2mn}{A}_{s}sin\frac{rs\pi }{mn},\hfill \\ {\gamma }_r\hfill & =\hfill & -\frac{i}{2mn}\sum _{s=1}^{2mn}{G}_{s}sin\frac{rs\pi }{mn},\hfill \end{array}$$

(7)

where $A_s={\sum }_{a=1}^{mn}{L}_{a,a+s}$ and $G_s=\frac{1}{2}{\sum }_{a=1}^{mn}(L_{a,x}{L}_{a+s,x}-L_{x,a}{L}_{x,a+s})$. According to the orthogonal properties of the coefficients, we obtain

$$\begin{array}{ccc}{A}_s\hfill & =\hfill & \sum _{r=1}^{2mn}[-{\alpha }_{r}cos\frac{rs\pi }{mn}+{\beta }_{r}sin\frac{rs\pi }{mn}],\hfill \\ G_s\hfill & =\hfill & i\sum _{r=1}^{2mn}{\gamma }_{r}sin\frac{rs\pi }{mn}.\hfill \end{array}$$

(8)

From Equations (5)–(8), $H_x$ and $H_z$ have the expansions

$$\begin{array}{c}{H}_x=A_1=-2\sum _{r=1}^{mn-1}({\alpha }_{r}cos\frac{r\pi }{mn}-{\beta }_{r}sin\frac{r\pi }{mn})\hfill \\ -{\alpha }_0+{\alpha }_{mn},\hfill \\ H_z=-A_0={\alpha }_0+2\sum _{r=1}^{mn-1}{\alpha }_r+{\alpha }_{mn}.\hfill \end{array}$$

(9)

Because $A_{mn+s}=-Q{A}_s=-A_{s}Q$ and $G_{mn+s}=-Q{G}_s=-G_{s}Q$, and combining with Equation (8), we have

$$[1+{(-1)}^{r}Q]{\alpha }_r=[1+{(-1)}^{r}Q]{\beta }_r=[1+{(-1)}^{r}Q]{\gamma }_r=0.$$

(10)

When $Q=1$, ${\alpha }_{2r}={\beta }_{2r}={\gamma }_{2r}=0$ while ${\alpha }_{2r+1}={\beta }_{2r+1}={\gamma }_{2r+1}=0$ if $Q=-1$. Therefore, we can investigate the algebra (8) with $Q=1$ 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 ${\alpha }_r$, ${\beta }_r$ and ${\gamma }_r$. Similar to those calculations in Ref. [4], we obtain

$$\begin{array}{c}[A_i,A_j]=4G_{i-j},[G_i,G_j]=0,\hfill \\ [G_i,A_j]=2(A_{j+i}-A_{j-i}).\hfill \end{array}$$

(11)

Substituting Equation (8) into Equation (11), we arrive at

$$[{\alpha }_r,{\beta }_r]=2i{\gamma }_r,[{\beta }_r,{\gamma }_r]=2i{\alpha }_r,[{\gamma }_r,{\alpha }_r]=2i{\beta }_r,$$

(12)

where $r=1,2,\cdots ,mn-1$, and all the other commutators vanish. Obviously, the algebra in (12) is associated with the site r, and hence is local. Because ${\alpha }_r$, ${\beta }_r$, and ${\gamma }_r$ obey the same commutation relations with $-X_r$, $-Y_r$, and $-Z_r$ in Ref. [4], we have the further relations

$$\begin{array}{c}{\alpha }_0^2=\frac{1}{2}(1-Q)=R_0,\hfill \\ {\alpha }_{mn}^2=\frac{1}{2}[1-{(-1)}^{mn}Q]=R_{mn},\hfill \\ {\alpha }_r{\beta }_r=i{\gamma }_r,{\beta }_r{\gamma }_r=i{\alpha }_r,{\gamma }_r{\alpha }_r=i{\beta }_r,\hfill \\ {\alpha }_r^2={\beta }_r^2={\gamma }_r^2=R_r^2=R_r,{\alpha }_r=R_r{\alpha }_r={\alpha }_r{R}_r,\hfill \\ {\beta }_r=R_r{\beta }_r={\beta }_r{R}_r,{\gamma }_r=R_r{\gamma }_r={\gamma }_r{R}_r.\hfill \end{array}$$

(13)

We note that $A_{sm}={\sum }_{p=1}^m{A}_{p,s}={\sum }_{p=1}^m{\sum }_{a=1}^n{L}_{a,a+s}^p={\sum }_{p=1}^m{\sum }_{a=1}^n$$L_{p+(a-1)m,p+(a-1+s)m}$ and $G_{sm}={\sum }_{p=1}^m{G}_{p,s}=\frac{1}{2}{\sum }_{p=1}^m{\sum }_{a=1}^n[L_{p+(a-1)m,x}\times L_{p+(a-1+s)m,x}-L_{x,p+(a-1)m}{L}_{x,p+(a-1+s)m}]$, where $s=1,2,\cdots ,2n$, and

$$\begin{array}{ccc}{A}_{p,s}\hfill & =\hfill & \sum _{q=1}^{2n}\{-{\alpha }_{p+(q-1)m}cos\frac{[p+(q-1\left)m\right]s\pi }{n}\hfill \\ & & +{\beta }_{p+(q-1)m}sin\frac{[p+(q-1\left)m\right]s\pi }{n}\},\hfill \\ G_{p,s}\hfill & =\hfill & i\sum _{q=1}^{2n}{\gamma }_{p+(q-1)m}sin\frac{[p+(q-1\left)m\right]s\pi }{n}.\hfill \end{array}$$

(14)

When $m=p=1$, Equation (14) recovers the results in two dimensions [4]. It is obvious that $A_{p,i}$ and $G_{p,j}$ also satisfy the commutation relations (11). When $p\ne p^{\prime }$, $[A_{p,i},A_{p^{\prime },i^{\prime }}]=[G_{p,j},G_{p^{\prime },j^{\prime }}]=[A_{p,i},G_{p^{\prime },i^{\prime }}]=0$.

We have obtained the expressions of $H_x$ and $H_z$ in terms of the operators ${\alpha }_r$, ${\beta }_r$ and ${\gamma }_r$ in the space $\Gamma$. To obtain the Hamiltonian in the third dimension, we project the operator algebra in the space $\Gamma$ into the $\mathbf{Y}$ direction. Then we have m subspaces ${\Gamma }_p(p=1,2,\cdots ,m)$, in which the operator algebra with period $2n$ is same with that in $\Gamma$. In ${\Gamma }_p$, we define

$$\begin{array}{c}{\mathcal{L}}_{a,a}^p=-{\sigma }_{p+(a-1)m}^x,\hfill \\ {\mathcal{L}}_{a,b}^p={\sigma }_{p+(a-1)m}^z{\sigma }_{p+am}^x\cdots {\sigma }_{p+(b-2)m}^x{\sigma }_{p+(b-1)m}^z\hfill \end{array}$$

(15)

along the $\mathbf{Y}$ direction. Then we have ${\mathcal{A}}_{p,s}={\sum }_{a=1}^n{\mathcal{L}}_{a,a+s}^p$ and ${\mathcal{G}}_{p,s}=\frac{1}{2}{\sum }_{a=1}^n[{\mathcal{L}}_{p+(a-1)m,x}{\mathcal{L}}_{p+(a-1+s)m,x}-{\mathcal{L}}_{x,p+(a-1)m}\times {\mathcal{L}}_{x,p+(a-1+s)m}]$, which also obey the same commutation relations (11) and (12), similar to $A_{p,s}$ and $G_{p,s}$. Then the Hamiltonian $H_y={\sum }_{p=1}^m{\mathcal{A}}_{p,1}$.

Because $[{\mathcal{L}}_{a,a+s}^p,L_{b,b+s}^p]=0$ (see Figure 2), we have $[{\mathcal{A}}_{p,s},A_{p,s}]=0$, which leads to ${\mathcal{A}}_{p,s}\equiv A_{p,s}$ due to their common local algebra (12). This is a renormalization of operator, which means that ${\mathcal{A}}_{p,s}$ and $A_{p,s}$ have same eigenfunctions and eigenvalues in ${\Gamma }_p$ or Γ space. We note that $V_2$ is the transfer matrix along $\mathbf{Y}$ direction, which must be calculated in $\Gamma$ rather than ${\Gamma }_p$ space by mapping ${\mathcal{A}}_{p,1}\equiv A_{p,1}$ to diagonalize total transfer matrix V. Therefore, we have

$$\begin{array}{ccc}{H}_y\hfill & =\hfill & \sum _{p=1}^m{\mathcal{A}}_{p,1}=\sum _{p=1}^m{A}_{p,1}\equiv A_m\hfill \\ & =\hfill & -{\alpha }_0-2\sum _{r=1}^{mn-1}({\alpha }_{r}cos\frac{r\pi }{n}-{\beta }_{r}sin\frac{r\pi }{n})\hfill \\ & & -{(-1)}^m{\alpha }_{mn}.\hfill \end{array}$$

(16)

Here, we would like to mention that $H_z=-{\sum }_{p=1}^m{\mathcal{A}}_{p,0}\equiv -A_0$, which is same with that in (9). This means that when $J_1=0$, the Hamiltonian of the 2D Ising model is recovered immediately.

Because $[Q,H_x]=[Q,H_y]=[Q,H_z]=[Q,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 $\pm 1$ of Q.

The transfer matrix V with Equations (9) and (16) becomes

$$\begin{array}{ccc}V\hfill & =\hfill & {[2sinh\left(2H\right)]}^{\frac{mn}{2}}{e}^{H_1{A}_1}{e}^{H_2{A}_m}{e}^{-H^{*}{A}_0}\hfill \\ & =\hfill & {[2sinh\left(2H\right)]}^{\frac{mn}{2}}{e}^{(H^{*}-H_1-H_2){\alpha }_0}\hfill \\ & & \times \prod _{r=1}^{mn-1}{U}_r{e}^{[H^{*}+H_1-{(-1)}^m{H}_2]{\alpha }_{mn}},\hfill \end{array}$$

(17)

where

$$\begin{array}{ccc}{U}_r\hfill & =\hfill & e^{-2H_1({\alpha }_{r}cos\frac{r\pi }{mn}-{\beta }_{r}sin\frac{r\pi }{mn})}\hfill \\ & & \times e^{-2H_2({\alpha }_{r}cos\frac{r\pi }{n}-{\beta }_{r}sin\frac{r\pi }{n})}{e}^{2H^{*}{\alpha }_r}.\hfill \end{array}$$

To obtain the eigenvalues of the transfer matrix V, we first diagonalize $U_r$ by employing the general unitary transformation:

$$\begin{array}{c}{e}^{\frac{i}{2}{\eta }_r{\gamma }_r}{e}^{a_r({\alpha }_{r}cos{\theta }_r+{\beta }_{r}sin{\theta }_r)}{U}_r\hfill \\ \times e^{-a_r({\alpha }_{r}cos{\theta }_r+{\beta }_{r}sin{\theta }_r)}{e}^{-\frac{i}{2}{\eta }_r{\gamma }_r}=e^{{\xi }_r{\alpha }_r}.\hfill \end{array}$$

(18)

Here ${\theta }_r$ is an arbitrary constant and can be taken to be zero without loss of generality, and

$$\begin{array}{c}cosh{\xi }_r={\mathcal{D}}_r,\hfill \\ sinh{\xi }_{r}cos{\eta }_r={\mathcal{A}}_r,tanh\left(2a_r\right)=\frac{{\mathcal{C}}_r}{{\mathcal{B}}_r},\hfill \\ sinh{\xi }_{r}sin{\eta }_r={\mathcal{B}}_{r}cosh\left(2a_r\right)-{\mathcal{C}}_{r}sinh\left(2a_r\right),\hfill \end{array}$$

(19)

where

$$\begin{array}{ccc}{\mathcal{A}}_r\hfill & =\hfill & cosh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)\hfill \\ & & -sinh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)cos\frac{r\pi }{mn}\hfill \\ & & -cosh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)cos\frac{r\pi }{n}\hfill \\ & & +sinh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{(m-1)r\pi }{mn},\hfill \\ {\mathcal{B}}_r\hfill & =\hfill & sinh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)sin\frac{r\pi }{mn}\hfill \\ & & +cosh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)sin\frac{r\pi }{n}\hfill \\ & & +sinh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{(m-1)r\pi }{mn},\hfill \\ {\mathcal{C}}_r\hfill & =\hfill & sinh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{r\pi }{mn}\hfill \\ & & +cosh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{r\pi }{n}\hfill \\ & & +sinh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)sin\frac{(m-1)r\pi }{mn},\hfill \\ {\mathcal{D}}_r\hfill & =\hfill & cosh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)\hfill \\ & & -sinh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{r\pi }{mn}\hfill \\ & & -cosh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{r\pi }{n}\hfill \\ & & +sinh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)cos\frac{(m-1)r\pi }{mn}.\hfill \end{array}$$

We note that ${\mathcal{D}}_r^2+{\mathcal{C}}_r^2-{\mathcal{A}}_r^2-{\mathcal{B}}_r^2\equiv 1$, which ensures that the 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^{*}$. Therefore, Equation (19) recovers Onsager’s results in the 2D Ising model [4].

Then the transfer matrix V has a diagonal form

$$\begin{array}{c}{e}^{{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r{\gamma }_r}{e}^{{\sum }_{r=1}^{mn-1}{a}_r{\alpha }_r}V{e}^{-{\sum }_{r=1}^{mn-1}{a}_r{\alpha }_r}\hfill \\ \times e^{-{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r{\gamma }_r}={[2sinh\left(2H^{*}\right)]}^{\frac{mn}{2}}\hfill \\ \times e^{(H^{*}-H_1-H_2){\alpha }_0+{\sum }_{r=1}^{mn-1}{\xi }_r{\alpha }_r+[H^{*}+H_1-{(-1)}^m{H}_2]{\alpha }_{mn}}.\hfill \end{array}$$

(20)

3. Transformations

3.1. Transformation 1

To explore the symmetries in the 3D Ising model, we take the transformation

$$\begin{array}{ccc}{\alpha }_r^{*}\hfill & =\hfill & -{\alpha }_{r}cos\frac{r\pi }{mn}+{\beta }_{r}sin\frac{r\pi }{mn},\hfill \\ {\beta }_r^{*}\hfill & =\hfill & {\alpha }_{r}sin\frac{r\pi }{mn}+{\beta }_{r}cos\frac{r\pi }{mn},{\gamma }_r^{*}=-{\gamma }_r.\hfill \end{array}$$

(21)

It is easy to prove that ${\alpha }_r^{*}$, ${\beta }_r^{*}$ and ${\gamma }_r^{*}$ satisfy the same commutation relations with ${\alpha }_r$, ${\beta }_r$ and ${\gamma }_r$. Then we have

$$\begin{array}{c}{H}_x={\alpha }_0^{*}+2\sum _{r=1}^{mn-1}{\alpha }_r^{*}+{\alpha }_{mn}^{*},\hfill \\ H_y={\alpha }_0^{*}+2\sum _{r=1}^{mn-1}[{\alpha }_r^{*}cos\frac{(m-1)r\pi }{mn}+{\beta }_r^{*}sin\frac{(m-1)r\pi }{mn}]\hfill \\ -{(-1)}^m{\alpha }_{mn}^{*},\hfill \\ H_z=-{\alpha }_0^{*}-2\sum _{r=1}^{mn-1}[{\alpha }_r^{*}cos\frac{r\pi }{mn}-{\beta }_r^{*}sin\frac{r\pi }{mn}]+{\alpha }_{mn}^{*}.\hfill \end{array}$$

(22)

Obviously, compared to Equations (9) and (16), such a transformation (21) exchanges the interaction forms in $(1,0,0)$ and $(0,0,1)$ directions (i.e., $H_x$ and $H_z$), but changes the interaction form in $(0,1,0)$ direction (i.e., $H_y$). Therefore, the 3D Ising model has no dual transformation, and the critical point cannot be fixed by Kramers and Wannier’s approach [3].

The transfer matrix can be expressed as

$$\begin{array}{ccc}V\hfill & =\hfill & {[2sinh\left(2H\right)]}^{\frac{mn}{2}}{e}^{H_1{A}_1}{e}^{H_2{A}_m}{e}^{-H^{*}{A}_0}\hfill \\ & =\hfill & {[2sinh\left(2H\right)]}^{\frac{mn}{2}}{e}^{(H_1+H_2-H^{*}){\alpha }_0^{*}}\hfill \\ & & \times \prod _{r=1}^{mn-1}{U}_r^{*}{e}^{[H_1-{(-1)}^m{H}_2+H^{*}]{\alpha }_{mn}^{*}},\hfill \end{array}$$

(23)

where

$$\begin{array}{ccc}{U}_r^{*}\hfill & =\hfill & e^{2H_1{\alpha }_r^{*}}{e}^{2H_2[{\alpha }_r^{*}cos\frac{(m-1)r\pi }{mn}+{\beta }_r^{*}sin\frac{(m-1)r\pi }{mn}]}\hfill \\ & & \times e^{-2H^{*}({\alpha }_r^{*}cos\frac{r\pi }{mn}-{\beta }_r^{*}sin\frac{r\pi }{mn})}.\hfill \end{array}$$

Following the procedure above, we can diagonalize the transfer matrix V, i.e.,

$$\begin{array}{c}{e}^{{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r^{*}{\gamma }_r^{*}}{e}^{{\sum }_{r=1}^{mn-1}{a}_r^{*}{\alpha }_r^{*}}V{e}^{-{\sum }_{r=1}^{mn-1}{a}_r^{*}{\alpha }_r^{*}}\hfill \\ \times e^{-{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r^{*}{\gamma }_r^{*}}={[2sinh\left(2H^{*}\right)]}^{\frac{mn}{2}}\hfill \\ \times e^{(H_1+H_2-H^{*}){\alpha }_0^{*}+{\sum }_{r=1}^{mn-1}{\xi }_r{\alpha }_r^{*}+[H_1-{(-1)}^m{H}_2+H^{*}]{\alpha }_{mn}^{*}},\hfill \end{array}$$

(24)

where

$$\begin{array}{c}sinh{\xi }_{r}cos{\eta }_r^{*}={\mathcal{A}}_r^{*},tanh\left(2a_r^{*}\right)=-\frac{{\mathcal{C}}_r}{{\mathcal{B}}_r^{*}},\hfill \\ sinh{\xi }_{r}sin{\eta }_r^{*}={\mathcal{B}}_r^{*}cosh\left(2a_r^{*}\right)+{\mathcal{C}}_{r}sinh\left(2a_r^{*}\right),\hfill \end{array}$$

(25)

and

$$\begin{array}{ccc}{\mathcal{A}}_r^{*}\hfill & =\hfill & sinh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)\hfill \\ & & -cosh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{r\pi }{mn}\hfill \\ & & -sinh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{r\pi }{n}\hfill \\ & & +cosh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)cos\frac{(m-1)r\pi }{mn},\hfill \\ {\mathcal{B}}_r^{*}\hfill & =\hfill & cosh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{r\pi }{mn}\hfill \\ & & +sinh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{r\pi }{n}\hfill \\ & & +cosh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)sin\frac{(m-1)r\pi }{mn}.\hfill \end{array}$$

We also have ${\mathcal{D}}_r^2+{\mathcal{C}}_r^2-{\mathcal{A}}_r^{*2}-{\mathcal{B}}_r^{*2}\equiv 1.$

3.2. Transformation 2

Let

$$\begin{array}{ccc}{\alpha }_r^{\prime }\hfill & =\hfill & -{\alpha }_{r}cos\frac{r\pi }{n}+{\beta }_{r}sin\frac{r\pi }{n},\hfill \\ {\beta }_r^{\prime }\hfill & =\hfill & {\alpha }_{r}sin\frac{r\pi }{n}+{\beta }_{r}cos\frac{r\pi }{n},{\gamma }_r^{\prime }=-{\gamma }_r,\hfill \end{array}$$

(26)

then we have

$$\begin{array}{c}{H}_x={\alpha }_0^{\prime }+2\sum _{r=1}^{mn-1}[{\alpha }_r^{\prime }cos\frac{(m-1)r\pi }{mn}-{\beta }_r^{\prime }sin\frac{(m-1)r\pi }{mn}]\hfill \\ -{(-1)}^m{\alpha }_{mn}^{\prime },\hfill \\ H_y={\alpha }_0^{\prime }+2\sum _{r=1}^{mn-1}{\alpha }_r^{\prime }+{\alpha }_{mn}^{\prime },\hfill \\ H_z=-{\alpha }_0^{\prime }-2\sum _{r=1}^{mn-1}[{\alpha }_r^{\prime }cos\frac{r\pi }{n}-{\beta }_r^{\prime }sin\frac{r\pi }{n}]\hfill \\ -{(-1)}^m{\alpha }_{mn}^{\prime }.\hfill \end{array}$$

(27)

By also comparing with Equations (9) and (16), the transformation (26) exchanges the interaction forms in $(0,1,0)$ and $(0,0,1)$ directions (i.e., $H_y$ and $H_z$), but changes the interaction form in $(1,0,0)$ direction (i.e., $H_x$). Therefore, such the transformation is not a dual transformation yet, which cannot be used to fix the critical point [3].

The transfer matrix reads

$$\begin{array}{ccc}V\hfill & =\hfill & {[2sinh\left(2H\right)]}^{\frac{mn}{2}}{e}^{H_1{A}_1}{e}^{H_2{A}_m}{e}^{-H^{*}{A}_0}\hfill \\ & =\hfill & {[2sinh\left(2H\right)]}^{\frac{mn}{2}}{e}^{(H_1+H_2-H^{*}){\alpha }_0^{\prime }}\hfill \\ & & \times \prod _{r=1}^{mn-1}{U}_r^{\prime }{e}^{[-{(-1)}^m{H}_1+H_2-{(-1)}^m{H}^{*}]{\alpha }_{mn}^{\prime }},\hfill \end{array}$$

(28)

where

$$\begin{array}{ccc}{U}_r^{\prime }\hfill & =\hfill & e^{2H_1[{\alpha }_r^{\prime }cos\frac{(m-1)r\pi }{mn}-{\beta }_r^{\prime }sin\frac{(m-1)r\pi }{mn}]}{e}^{2H_2{\alpha }_r^{\prime }}\hfill \\ & & \times e^{-2H^{*}({\alpha }_r^{\prime }cos\frac{r\pi }{n}-{\beta }_r^{\prime }sin\frac{r\pi }{n})}.\hfill \end{array}$$

Similarly, we have

$$\begin{array}{c}{e}^{{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r^{\prime }{\gamma }_r^{\prime }}{e}^{{\sum }_{r=1}^{mn-1}{a}_r^{\prime }{\alpha }_r^{\prime }}V{e}^{-{\sum }_{r=1}^{mn-1}{a}_r^{\prime }{\alpha }_r^{\prime }}\hfill \\ \times e^{-{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r^{\prime }{\gamma }_r^{\prime }}={[2sinh\left(2H^{*}\right)]}^{\frac{mn}{2}}\hfill \\ \times e^{(H_1+H_2-H^{*}){\alpha }_0^{\prime }+{\sum }_{r=1}^{mn-1}{\xi }_r{\alpha }_r^{\prime }+[H_2-{(-1)}^m(H_1+H^{*})]{\alpha }_{mn}^{\prime }}.\hfill \end{array}$$

(29)

Here,

$$\begin{array}{c}sinh{\xi }_{r}cos{\eta }_r^{\prime }={\mathcal{A}}_r^{\prime },tanh\left(2a_r^{\prime }\right)=-\frac{{\mathcal{C}}_r}{{\mathcal{B}}_r^{\prime }},\hfill \\ sinh{\xi }_{r}sin{\eta }_r^{\prime }={\mathcal{B}}_r^{\prime }cosh\left(2a_r^{\prime }\right)+{\mathcal{C}}_{r}sinh\left(2a_r^{\prime }\right),\hfill \end{array}$$

(30)

and

$$\begin{array}{ccc}{A}_r^{\prime }\hfill & =\hfill & cosh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)\hfill \\ & & -cosh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{r\pi }{n}\hfill \\ & & -sinh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{(2m-1)r\pi }{mn}\hfill \\ & & +sinh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)cos\frac{(m-1)r\pi }{mn},\hfill \\ B_r^{\prime }\hfill & =\hfill & cosh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{r\pi }{n}\hfill \\ & & -sinh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)sin\frac{(m-1)r\pi }{mn}\hfill \\ & & +sinh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)sin\frac{(2m-1)r\pi }{mn}.\hfill \end{array}$$

The identity ${\mathcal{D}}_r^2+{\mathcal{C}}_r^2-{\mathcal{A}}_r^{\prime 2}-{\mathcal{B}}_r^{\prime 2}\equiv 1$ also holds.

4. Results

Because ${\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_{mn}$ have the common eigenvectors ${\chi }_0$ with the corresponding eigenvalues ${\Delta }_0,{\Delta }_1,\cdots ,{\Delta }_{mn}$, from Equation (20), we have $V\psi =\lambda \psi$, where

$$\begin{array}{cc}\hfill \psi =& e^{-{\sum }_{r=1}^{mn-1}{a}_r{\alpha }_r}{e}^{-{\sum }_{r=1}^{mn-1}\frac{i}{2}{\eta }_r{\gamma }_r}{\chi }_0,\hfill \\ \hfill ln\lambda =& \frac{1}{2}mnln[2sinh\left(2H\right)]+(H^{*}-H_1-H_2){\Delta }_0\hfill \\ & +\sum _{r=1}^{mn-1}{\xi }_r{\Delta }_r+[H_1-{(-1)}^m{H}_2+H^{*}]{\Delta }_{mn}.\hfill \end{array}$$

(31)

At the critical point, we have ${\xi }_0=H^{*}-H_1-H_2=0$ [4]. This leads to a critical temperature $T=T_c$ given by the condition

$$sinh\left(2H\right)sinh(2H_1+2H_2)=1.$$

(32)

If $H_2=0$ or $H_1=0$, 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., $\eta =(H_1+H_2)/H\to 0$, the critical temperature determined by Equation (32) also agrees perfectly with the asymptotically exact value $H=2{[\ln {\eta }^{-1}-\mathrm{lnln}{\eta }^{-1}+0\left(1\right)]}^{-1}$ shown in Refs. [18,19].

When $H_1=H_2=H$, the critical value $H_c=J/\left(k_B{T}_c\right)=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 ${\lambda }_{max}$ of the transfer matrix V. Following Ref. [4], we have

$$\begin{array}{c}ln{\lambda }_{max}-\frac{1}{2}mnln[2sinh\left(2H\right)]\hfill \\ =\{\begin{array}{c}{\xi }_1+{\xi }_3+\cdots +{\xi }_{2L-1}\mathrm{for}mn=2L;\hfill \\ {\xi }_1+{\xi }_3+\cdots +{\xi }_{2L-1}\hfill \\ +H_1-{(-1)}^m{H}_2+H^{*}\mathrm{for}mn=2L+1.\hfill \end{array}\hfill \end{array}$$

(33)

Here ${\Delta }_1={\Delta }_3=\cdots ={\Delta }_{mn-1}=1$, which are same as the eigenvalues of the operators $X_r$ in Ref. [4]. We note that these two results above can be combined due to ${\xi }_{-r}={\xi }_r$ and ${\xi }_{mn}=2[H_1-{(-1)}^m{H}_2+H^{*}]$. Therefore, Equation (33) has the compact form

$$\begin{array}{c}ln{\lambda }_{max}-\frac{1}{2}mnln[2sinh\left(2H\right)]=\frac{1}{2}\sum _{r=1}^{mn}{\xi }_{2r-1}\hfill \\ =\frac{1}{2}\sum _{r=1}^{mn}{cosh}^{-1}[cosh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)\hfill \\ -sinh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{(2r-1)\pi }{2mn}\hfill \\ -cosh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)cos\frac{(2r-1)\pi }{2n}\hfill \\ +sinh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)cos\frac{(m-1)(2r-1)\pi }{2mn}].\hfill \end{array}$$

(34)

To calculate the partition function per atom ${\lambda }_{\infty }={lim}_{m,n\to \infty }{\left({\lambda }_{\max }\right)}^{\frac{1}{mn}}$ for the infinite crystal, we replace the sum in Equation (34) by the integral

$$ln{\lambda }_{\infty }=\frac{1}{2}ln[2sinh\left(2H\right)]+\frac{1}{2\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }{\xi }_m\left(\omega \right)d\omega ,$$

(35)

where

$$\begin{array}{c}cosh{\xi }_m\left(\omega \right)=\mathcal{D}\left(\omega \right)=cosh\left(2H_1\right)cosh\left(2H_2\right)cosh\left(2H^{*}\right)\hfill \\ -sinh\left(2H_1\right)cosh\left(2H_2\right)sinh\left(2H^{*}\right)cos\omega \hfill \\ -cosh\left(2H_1\right)sinh\left(2H_2\right)sinh\left(2H^{*}\right)cos\left(m\omega \right)\hfill \\ +sinh\left(2H_1\right)sinh\left(2H_2\right)cosh\left(2H^{*}\right)cos\left[(m-1)\omega \right].\hfill \end{array}$$

(36)

Similarly, the continuous $\mathcal{A}\left(\omega \right)$, ${\mathcal{A}}^{*}\left(\omega \right)$, ${\mathcal{A}}^{\prime }\left(\omega \right)$, $\mathcal{B}\left(\omega \right)$, ${\mathcal{B}}^{*}\left(\omega \right)$, ${\mathcal{B}}^{\prime }\left(\omega \right)$, $\mathcal{C}\left(\omega \right)$, ${\xi }_m\left(\omega \right)$, $\eta \left(\omega \right)$, ${\eta }^{*}\left(\omega \right)$, and ${\eta }^{\prime }\left(\omega \right)$ replace the discrete ${\mathcal{A}}_r$, ${\mathcal{A}}_r^{*}$, ${\mathcal{A}}_r^{\prime }$, ${\mathcal{B}}_r$, ${\mathcal{B}}_r^{*}$, ${\mathcal{B}}_r^{\prime }$, ${\mathcal{C}}_r$, ${\xi }_r$, ${\eta }_r$, ${\eta }_r^{*}$, and ${\eta }_r^{\prime }$, respectively, by letting $\omega =\frac{r\pi }{mn}$. Here we emphasize that when $H_2=0$, or $H_1=0$, 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 $N=mnl$, the free energy

$$F=U-TS=-N{k}_{B}Tln{\lambda }_{\infty },$$

(37)

the internal energy

$$\begin{array}{ccc}U\hfill & =\hfill & F-T\frac{dF}{dT}=N{k}_B{T}^2\frac{ln{\lambda }_{\infty }}{dT}\hfill \\ & =\hfill & -N{k}_{B}T[H_1\frac{\partial ln{\lambda }_{\infty }}{\partial H_1}+H_2\frac{\partial ln{\lambda }_{\infty }}{\partial H_2}+H\frac{\partial ln{\lambda }_{\infty }}{\partial H}],\hfill \end{array}$$

(38)

and the specific heat

$$\begin{array}{cc}C\hfill & =\frac{dU}{dT}=N{k}_B[H_1^2\frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_1^2}+H_2^2\frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_2^2}+H^2\frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H^2}\hfill \\ & +2H_1{H}_2\frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_1\partial H_2}+2H_{1}H\frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_1\partial H}+2H_{2}H\frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_2\partial H}].\hfill \end{array}$$

(39)

Here,

$$\begin{array}{c}\frac{\partial ln{\lambda }_{\infty }}{\partial H_1}=\frac{1}{\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }cos{\eta }^{*}d\omega ,\hfill \\ \frac{\partial ln{\lambda }_{\infty }}{\partial H_2}=\frac{1}{\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }\frac{\frac{\partial \mathcal{D}}{\partial H_2}}{sinh{\xi }_m}d\omega ,\hfill \\ \frac{\partial ln{\lambda }_{\infty }}{\partial H}=cosh\left(2H^{*}\right)-\frac{1}{\pi }sinh\left(2H^{*}\right)\underset{m\to \infty }{lim}{\int }_0^{\pi }cos\eta d\omega ,\hfill \\ \frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_1^2}=\frac{2}{\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }{sin}^2{\eta }^{*}coth{\xi }_{m}d\omega ,\hfill \\ \frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_2^2}=\frac{1}{2\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }[4-\frac{1}{{sinh}^2{\xi }_m}{\left(\frac{\partial \mathcal{D}}{H_2}\right)}^2]coth{\xi }_{m}d\omega ,\hfill \\ \frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H^2}=2{sinh}^2\left(2H^{*}\right)[\frac{1}{\pi }coth\left(2H^{*}\right)\underset{m\to \infty }{lim}{\int }_0^{\pi }cos\eta d\omega \hfill \\ +\frac{1}{\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }{sin}^2\eta coth{\xi }_{m}d\omega -1],\hfill \\ \frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_1\partial H_2}=\frac{1}{\pi }\underset{m\to \infty }{lim}{\int }_0^{\pi }\frac{d\omega }{sinh{\xi }_m}(\frac{\partial {\mathcal{A}}^{*}}{\partial H_2}-\frac{\partial \mathcal{D}}{\partial H_2}cos{\eta }^{*}coth{\xi }_m),\hfill \\ \frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_1\partial H}=-\frac{1}{\pi }sinh\left(2H^{*}\right)\hfill \\ \times \underset{m\to \infty }{lim}{\int }_0^{\pi }\frac{d\omega }{sinh{\xi }_m}(\frac{\partial {\mathcal{A}}^{*}}{\partial H^{*}}-2cosh{\xi }_{m}cos\eta cos{\eta }^{*}),\hfill \\ \frac{{\partial }^{2}ln{\lambda }_{\infty }}{\partial H_2\partial H}=-\frac{1}{\pi }sinh\left(2H^{*}\right)\hfill \\ \times \underset{m\to \infty }{lim}{\int }_0^{\pi }\frac{d\omega }{sinh{\xi }_m}(\frac{\partial A}{\partial H_2}-\frac{\partial \mathcal{D}}{\partial H_2}cos\eta coth{\xi }_m).\hfill \end{array}$$

We note that at the critical point, ${lim}_{\omega \to 0}{\xi }_m\to 0$. However, ${lim}_{\omega \to 0}\frac{\partial \mathcal{D}}{\partial H_2}/sinh{\xi }_m\to -cos\eta \left(0\right)$. 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 $J_1=J_2$, where the calculation of the thermodynamic functions can be simplified considerably. After integrating, Equation (36) can be rewritten as

$$\begin{array}{c}cosh{\xi }_{\infty }\left(\omega \right)=cosh\left(2H_1\right)cosh\left(2H_{2D}^{*}\right)\hfill \\ -sinh\left(2H_1\right)sinh\left(2H_{2D}^{*}\right)cos\omega ,\hfill \end{array}$$

(40)

where

$$H_{2D}^{*}=H^{*}-H_1.$$

(41)

It is surprising that Equation (40) is nothing but that in 2D Ising model with the interaction energies $(J_1,J_{2D})$ and $H_{2D}=\frac{J_{2D}}{k_{B}T}$. Therefore, $ln{\lambda }_{\infty }-\frac{1}{2}ln[2sinh\left(2H\right)]$ in three dimensions can be obtained from $ln{\lambda }_{\infty }^{2D}-\frac{1}{2}ln[2sinh\left(2H_{2D}\right)]$ 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., $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 the 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 ${\alpha }_{2D}=0$. Because the 3D Ising model with the simple cubic symmetry (i.e., $H_1=H_2=H$) 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 ${\alpha }_{3D}$ of the 3D Ising model is identical to ${\alpha }_{2D}$, i.e. ${\alpha }_{3D}=0$. According to the scaling laws $d\nu =2-\alpha$ and $\mu +\nu =2-\alpha$ [6], we have ${\nu }_{3D}=\frac{2}{3}$ and ${\mu }_{3D}=\frac{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\left(y\right)$ 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 $x\left(y\right)$ 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.

5. High-Temperature Expansions

Now we calculate the high-temperature expansions of the partition function per atom when $J_1=J_2=J$. According to the identity

$${\int }_0^{2\pi }\ln (2\cosh x-2\cos {\omega }^{\prime })d{\omega }^{\prime }=2\pi x,$$

(42)

from Equation (35) and (36), we obtain

$$\begin{array}{ccc}ln\frac{{\lambda }_{\infty }}{2}\hfill & =\hfill & \frac{1}{2{\pi }^2}{\int }_0^{\pi }{\int }_0^{\pi }ln\{{cosh}^3\left(2H\right)\hfill \\ & & -sinh\left(2H\right)cosh\left(2H\right)[cos\omega +cos(m\omega \left)\right]\hfill \\ & & +{sinh}^2\left(2H\right)cosh\left(2H\right)cos\left[(m-1)\omega \right]\hfill \\ & & -sinh\left(2H\right)cos{\omega }^{\prime }\}d\omega d{\omega }^{\prime }\hfill \\ & =\hfill & 3*ln(cosh\left(H\right))+\frac{3}{2}ln(1+k^2)\hfill \\ & & +\frac{1}{2{\pi }^2}{\int }_0^{\pi }{\int }_0^{\pi }ln\{1-\frac{2k(1-k^2)}{{(1+k^2)}^2}[cos\omega +cos\left(m\omega \right)]\hfill \\ & & +\frac{4k^2}{{(1+k^2)}^2}cos\left[(m-1)\omega \right]\hfill \\ & & -\frac{2k{(1-k^2)}^2}{{(1+k^2)}^3}cos{\omega }^{\prime }\}d\omega d{\omega }^{\prime }\hfill \\ & =\hfill & 3*ln(cosh\left(H\right))-3k^4-62k^6-\frac{2081}{2}{k}^8\hfill \\ & & -21024k^{10}-\cdots ,\hfill \end{array}$$

(43)

where $k=tanhH$. Therefore, the partition function per atom at high temperatures is

$${\lambda }_{\infty }=2{cosh}^{3}H(1-3k^4-62k^6-1036k^8-20838k^{10}-\cdots ).$$

(44)

We note that for PBCs, the high-temperature partition function per atom reads [20]

$${\lambda }_{\infty }^p=2{cosh}^{3}H(1+3k^4+22k^6+192k^8+2046k^{10}+\cdots ).$$

(45)

Obviously, the difference between ${\lambda }_{\infty }$ and ${\lambda }_{\infty }^p$ comes from the screw boundary condition along the $\mathbf{X}$ direction (see Figure 1). We note that the $k^2$ term in Equations (44) and (45) vanishes, which can be seen as a feature of the 3D Ising model.

6. Conclusions

We have exactly solved the 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 $\left(100\right)$ direction and the periodic boundary conditions along both $\left(010\right)$ and $\left(001\right)$ 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 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 $H_1=H_2$ 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.