1 Introduction

Since Onsager^[1] gave the exact solution of 2-dimensional Ising model, 3-dimensional Ising model has attracted the attention of many researchers^[2-6], because the Onsager solution can explain the phase-transfer phenomenon. There are various techniques to investigate 3-dimensional Ising model: the renormalization group technique^[7], cluster variational method^[8], linear chain approximation^[9], molecular field approximation^[10], etc. Liu et al. ^[11] have discussed nonlinear magnetization dynamics of a classical ferromagnet with anisotropies in an external magnetic field. Jiang et al. ^[12] have examined the phase diagrams of a single spin-1 Ising system and given the general expressions of the second-order phase transition. Recently, El-Showk et al. ^[13] have tried to use the constraints of crossing symmetry and unitarity in general 3D conformal field theories to solve 3-dimension Ising model. They found that 3-dimensional Ising model lies at a corner point on the boundary of the allowed parameter space which makes it difficult to be solved. Jin et al. ^[14] studied the phase transition of Ising model by using Monte Carlo simulations. Wu et al. ^[15] used the bond-propagation algorithm to study the finite-size Ising model. They found that boundary conditions have remarkable effects on the free energy and specific heat. Deviren et al. ^[16] have studied Ising model in a time-dependent external magnetic field by using effective-field theory. They identified the phase transition point as the minimum-correlation point and studied the critical condition of the system. Heyl et al. ^[17] have shown that 3D Ising model generates periodic nonanalytic behavior during phase transition. Campo et al. ^[18] used the method of controlling asymptotic parameter to study phase transition of Ising model. Ates et al. ^[19] identified a first-order dynamical phase transition of 3D Ising model and studied the behavior of thermal parameters during phase transition. Zhang et al.^[20] used Ising model to explain the frustration in magnetic systems. Barnett et al. ^[21] compared the long-range interaction and shor-range interaction of Ising model near critical point. Schoelz et al. ^[22] used Ising model to investigate the graphene system and studied the changing of thermal parameters during phase transition. Dublenych^[23] found a complete and exact solution of the ground-state problem for the Ising model on the Shastry-Sutherland lattice in an applied magnetic field. Lau et al. ^[24] used numerical calculation to study the thermal parameters for Ising model. They used a bond propagation algorithm and calculate the change of entropy during phase transition. Langer^[25] used numerical simulations to study the correlation lengths of Ising model at critical points. Fytas et al. ^[26] solved the partition function of Ising model by performing statistics simulation. They found that specific heat rise exponentially near the critical point. Wu et al. ^[27] used the bond-propagation algorithm to calculate a rectangle Ising model. They have determined the exact expansions of the critical free energy and specific heat. Corberi et al. ^[28] studied the nonconserved phase-ordering dynamics of random-field Ising model. They confirmed the crossover to an asymptotic logarithmic behavior of thermal parameters during phase transition. Powalski et al. ^[29] investigated the transverse field Ising model on a triangular lattice and discovered a second-order quantum phase transition in the high-field limit. Damski^[30] derived an exact closed-form expression for fidelity susceptibility of the quantum Ising model. Bastidas et al. ^[31] established a set of nonequilibrium phase transitions in the Ising model. They showed that the system exhibits and anisotropic transition. Henry et al. ^[32] studied the ground-state properties of Ising model by using spin-wave theory. They found that the ground-state system becomes stable in zero field. Jin et al. ^[33] investigated the thermal phase transition of a square-lattice Ising model. They discovered that the system is stabilized below a critical temperature.

Many theoretical works have been studied by the effective-field theory with correlations, but the core of Ising model is to obtain the largest eigenvalue of the transfer matrix and there is lack of understanding of the mathematical structure of the transfer matrix.

In this paper, we study the eigenvalues of the 3-dimensional Ising model of N₁×N₂×N₃. In Section 2, we build a 2^N₁N₂×2^N₁N₂ transfer matrix V₁, and then we build a much smaller matrix V₂ which contains some eigenvalues including the largest eigenvalue λ_max of V₁. If N₃ is infinite, we only need to obtain the largest eigenvalue of V₁.But in this paper, we set N₁ and N₂ at finite values.N₁ and N₂ are the finite width of our model. Then we use λ_max to calculate the specific heat per atom C_pa and the magnetic momentum per atom M_pa in Section 3.In Section 4, the important conclusions are drawn.

2 Basic Model

The 3-dimensional Ising model is a cubic lattice of N₁×N₂×N₃ sites. These sites are occupied by atoms with positive or negative spin. We use σ_{j, k, l}=±1 to represent the spin sign of atoms. The total energy of the lattice is:

$\begin{gathered} E = \sum\limits_{j = 0}^{{N_1} - 1} {\sum\limits_{h = 0}^{{N_2} - 1} {\sum\limits_{l = 0}^{{N_3} - 1} {\left( { - {J_1}{\sigma _{j,h,l}}{\sigma _{j + 1,h,l}} - } \right.} } } \\ \left. {{J_2}{\sigma _{j,h,l}}{\sigma _{j,h + 1,l}} - {J_3}{\sigma _{j,h,l}}{\sigma _{j,h,l + 1}} - \mu B{\sigma _{j,h,l}}} \right) \\ \end{gathered} $

(1)

where J₁, J₂, J₃ are the interaction energy between neighbors; μ is the spin magnetic moment; B is the magnetic field. The boundary condition is σ_{0, h, l}=σ_{N1, h, l}, σ_{j, 0, l}=σ_{j, N2, l}, σ_{j, h, 0}=σ_{j, h, N3}.

The partition function is:

$Z = \sum\limits_{{\sigma _{j,h,l}} = \pm 1} {{e^{ - \frac{E}{{kT}}}}} $

(2)

The free energy F can be derived from Z. The specific heat C, the magnetic moment M can be derived from F:

$F = - kT\ln Z$

(3)

$C = - T\frac{{{\partial ^2}F}}{{\partial {T^2}}}$

(4)

$M = - \frac{{\partial F}}{{\partial B}}$

(5)

To calculate Z, we divide the lattice into N₃ layers. Each layer contains N₁ rows and N₂ sites per row. We use a number ν_l to represent the state of the lth layer, 1≤ν_l≤2^N₁N₂. The total energy of the lattice can be written as:

$E = \sum\limits_{l = 0}^{{N_3} - 1} {{E_1}\left( {{v_1}} \right) + \sum\limits_{l = 0}^{{N_3} - 1} {{E_2}\left( {{v_l},{v_{l + 1}}} \right)} } $

(6)

$\begin{gathered} {E_1}\left( {{v_l}} \right) = \sum\limits_{j = 0}^{{N_1} - 1} {\sum\limits_{h = 0}^{{N_2} - 1} {\left( { - {J_1}{\sigma _{j,h,l}}{\sigma _{j + 1,h,l}} - } \right.} } \hfill \\ {J_2}{\sigma _{j,h,l}}{\sigma _{j,h + 1,l}} - \mu B{\sigma _{j,h,l}} \hfill \\ \end{gathered} $

(7)

${E_2}\left( {{v_l},{v_{l + 1}}} \right) = \sum\limits_{j = 0}^{{N_1} - 1} {\sum\limits_{h = 0}^{{N_2} - 1} {{J_3}{\sigma _{j,h,l}}{\sigma _{j,h,l + 1}}} } $

(8)

Then the partition function turns into:

$Z = \sum\limits_{{v_l}} {\prod\limits_{l = 0}^{{N_3} - 1} {{e^{ - {E_1}\left( {{v_l}} \right) - {E_2}\left( {{v_l},{v_{l + 1}}} \right)}}} } $

(9)

where ${e^{\frac{{ - {E_1}\left( {{\nu _l}} \right) - {E_2}\left( {{\nu _l},{\nu _{l + 1}}} \right)}}{{kT}}}}$ can be seen as the elements of matrix V₁; ν_l is the row number; ν_l+1 is the column number. Then Z turns into:

$Z = Tr\left( {V_1^{{N_3}}} \right) = \lambda _1^{{N_3}} + \lambda _2^{{N_3}} + \lambda _3^{{N_3}} + \ldots $

(10)

where λ₁, λ₂, λ₃, ... are the different eigenvalues of V₁. Assuming λ₁ is the largest eigenvalue λ_max, if N₃ is infinite, then:

$\mathop {\lim }\limits_{{N_3} \to \infty } \left( {\lambda _1^{{N_3}} + \lambda _2^{{N_3}} + \lambda _3^{{N_3}} + \ldots } \right) = \lambda _1^{{N_3}}$

(11)

Now we only need to obtain the largest eigenvalue λ_max of matrix V₁ and Z=λ_max^N₃.

We will use a special case to explain how to obtain λ_max. When N₁=2, N₂=1, V₁ can be expressed as:

$\begin{array}{l} {V_1} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B + 2{J_3}}}{{kT}}}}{D_1} + {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}{D_2} + \\ {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B - 2{J_3}}}{{kT}}}}{D_3} + {e^{\frac{{ - 2{J_1} + 2{J_2}}}{{kT}}}}\left( {{D_4} + {D_6}} \right) + \\ {e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}}\frac{1}{2}\left( {{D_5} + {D_{10}}} \right) + {e^{\frac{{ - 2{J_1} + 2{J_2} + 2\mu B}}{{kT}}}}{D_8} + {{\rm{e}}^{\frac{{{\rm{2}}{{\rm{J}}_{\rm{1}}}{\rm{ + 2}}{{\rm{J}}_{\rm{2}}}{\rm{ + 2}}\mu {\rm{B + 2}}{{\rm{J}}_{\rm{3}}}}}{{{\rm{kT}}}}}}{D_9}12 \end{array}$

(12)

${\boldsymbol{D}_1} = \left[{\begin{array}{*{20}{c}} 1&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \end{array}} \right]$

(13)

${\boldsymbol{D}_2} = \left[{\begin{array}{*{20}{c}} 0&1&1&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \end{array}} \right]$

(14)

${\boldsymbol{D}_3} = \left[{\begin{array}{*{20}{c}} 0&0&0&1 \\ 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \end{array}} \right]$

(15)

${\boldsymbol{D}_4} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 1&0&0&0 \\ 1&0&0&0 \\ 0&0&0&0 \end{array}} \right]$

(16)

${\boldsymbol{D}_5} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&1&1&0 \\ 0&1&1&0 \\ 0&0&0&0 \end{array}} \right]$

(17)

${\boldsymbol{D}_6} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&0&0&1 \\ 0&0&0&1 \\ 0&0&0&0 \end{array}} \right]$

(18)

${\boldsymbol{D}_7} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ 1&0&0&0 \end{array}} \right]$

(19)

${\boldsymbol{D}_8} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ 0&1&1&0 \end{array}} \right]$

(20)

${\boldsymbol{D}_9} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&1 \end{array}} \right]$

(21)

${{\mathbf{\boldsymbol{D}}}_{10}} = \left[{\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&1&-1&0 \\ 0&-1&0&0 \\ 0&0&0&0 \end{array}} \right]$

(22)

If D_j multiples V₁, the results are:

${\boldsymbol{D}_1}{\boldsymbol{V}_1} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B + 2{J_3}}}{{kT}}}}{\boldsymbol{D}_1} + {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}{\boldsymbol{D}_2} + {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B - 2{J_3}}}{{kT}}}}{\boldsymbol{D}_3}$

(23)

$\begin{gathered} {\boldsymbol{D}_2}{\boldsymbol{V}_1} = 2{e^{\frac{{ - 2{J_1} + 2{J_2}}}{{kT}}}}\left( {{\boldsymbol{D}_1} + {\boldsymbol{D}_3}} \right) + \\ \left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}}{e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right){\boldsymbol{D}_2} \\ \end{gathered} $

(24)

${\boldsymbol{D}_3}{\boldsymbol{V}_1} = {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B - 2{J_3}}}{{kT}}}}{\boldsymbol{D}_1} + {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B}}{{kT}}}}{\boldsymbol{D}_2} + {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B + 2{J_3}}}{{kT}}}}{\boldsymbol{D}_3}$

(25)

${\boldsymbol{D}_4}{\boldsymbol{V}_1} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B + 2{J_3}}}{{kT}}}}{\boldsymbol{D}_4} + {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}{\boldsymbol{D}_5} + {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B - 2{J_3}}}{{kT}}}}{\boldsymbol{D}_6}$

(26)

$\begin{gathered} {\boldsymbol{D}_5}{\boldsymbol{V}_1} = 2{e^{\frac{{ - 2{J_1} + 2{J_2}}}{{kT}}}}\left( {{\boldsymbol{D}_4} + {\boldsymbol{D}_6}} \right) + \hfill \\ \left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} + {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right){\boldsymbol{D}_5} \hfill \\ \end{gathered} $

(27)

$\begin{gathered} {\boldsymbol{D}_6}{\boldsymbol{V}_1} = {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B - 2{J_3}}}{{kT}}}}{\boldsymbol{D}_4} + {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B}}{{kT}}}}{\boldsymbol{D}_5} + \\ {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B + 2{J_3}}}{{kT}}}}{\boldsymbol{D}_6} \\ \end{gathered} $

(28)

$\begin{gathered} {\boldsymbol{D}_7}{\boldsymbol{V}_1} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B + 2{J_3}}}{{kT}}}}{\boldsymbol{D}_7} + {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}{\boldsymbol{D}_8} + \\ {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B - 2{J_3}}}{{kT}}}}{\boldsymbol{D}_9} \\ \end{gathered} $

(29)

$\begin{gathered} {\boldsymbol{D}_8}{\boldsymbol{V}_1} = {e^{\frac{{ - 2{J_1} + 2{J_2}}}{{kT}}}}2\left( {{\boldsymbol{D}_7} + {\boldsymbol{D}_9}} \right) + \hfill \\ \left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} + {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right){\boldsymbol{D}_8} \hfill \\ \end{gathered} $

(30)

$\begin{gathered} {D_9}{V_1} = {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B - 2{J_3}}}{{kT}}}}{D_7} + \hfill \\ {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B}}{{kT}}}}{D_8} + {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B + 2{J_3}}}{{kT}}}}{D_9} \hfill \\ \end{gathered} $

(31)

${\boldsymbol{D}_{10}}{\boldsymbol{V}_1} = \left( {{e^{\frac{{2{J_1} + 2{J_2} + 2\mu B}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right){\boldsymbol{D}_{10}}$

(32)

From above, it can be seen that D₁V₁-D₃V₁ can be expressed as the linear combinations of D₁-D₃. D₄V₁-D₆V₁ can be expressed as the linear combinations of D₄-D₆. D₇V₁-D₉V₁ can be expressed as the linear combinations of D₇-D₉. D₁₀V₁ can be expressed with D₁₀. Now D₁-D₁₀ are divided into 4 groups: D₁-D₃, D₄-D₆, D₇-D₉ and D₁₀. D_jV₁ can be expressed as the linear combinations of D_k from the same group. This means we can calculate V₁^N₃ based on separate groups of D₁-D₁₀. We can establish a transfer matrix V₂ to describe the relationship between D_j from the same group:

${\boldsymbol{V}_2} = {e^{\frac{{2{J_2}}}{{kT}}}}\left[{\begin{array}{*{20}{c}} {{e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_1} - 2\mu B}}{{kT}}}}}&{{e^{\frac{{ - 2{J_3} + 2{J_1} - 2\mu B}}{{kT}}}}} \\ {2{e^{\frac{{ - 2{J_1}}}{{kT}}}}}&{{e^{\frac{{2{J_3} - 2{J_1}}}{{kT}}}} + {e^{\frac{{ - 2{J_3} - 2{J_1}}}{{kT}}}}}&{{e^{\frac{{ - 2{J_1}}}{{kT}}}}} \\ {{e^{\frac{{ - 2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_1} + 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}} \end{array}} \right]$

(33)

$\left( {{a_1}{\boldsymbol{D}_1} + {a_2}{\boldsymbol{D}_2} + {a_3}{\boldsymbol{D}_3}} \right)\boldsymbol{V}_1^n = {a_{10}}{\boldsymbol{D}_1} + {a_{11}}{\boldsymbol{D}_2} + {a_{12}}{\boldsymbol{D}_3}$

(34)

$\left[{\begin{array}{*{20}{c}} {{a_{10}}}&{{a_{11}}}&{{a_{12}}} \end{array}} \right]{\text{ = }}\left[{\begin{array}{*{20}{c}} {{a_1}}&{{a_2}}&{{a_3}} \end{array}} \right]\boldsymbol{V}_2^n$

(35)

$\left[{{a_4}{\boldsymbol{D}_4} + {a_5}{\boldsymbol{D}_5} + {a_6}{\boldsymbol{D}_6}} \right]\boldsymbol{V}_1^n = {a_{13}}{\boldsymbol{D}_4} + {a_{14}}{\boldsymbol{D}_5} + {a_{15}}{\boldsymbol{D}_6}$

(36)

$\left[{\begin{array}{*{20}{c}} {{a_{13}}}&{{a_{14}}}&{{a_{15}}} \end{array}} \right]{\text{ = }}\left[{\begin{array}{*{20}{c}} {{a_4}}&{{a_5}}&{{a_6}} \end{array}} \right]\boldsymbol{V}_2^n$

(37)

$\left[{{a_7}{\boldsymbol{D}_7} + {a_8}{\boldsymbol{D}_8} + {a_9}{\boldsymbol{D}_9}} \right]\boldsymbol{V}_1^n = {a_{16}}{\boldsymbol{D}_7} + {a_{17}}{\boldsymbol{D}_8} + {a_{18}}{\boldsymbol{D}_9}$

(38)

$\left[{\begin{array}{*{20}{c}} {{a_{16}}}&{{a_{17}}}&{{a_{18}}} \end{array}} \right]{\text{ = }}\left[{\begin{array}{*{20}{c}} {{a_7}}&{{a_8}}&{{a_9}} \end{array}} \right]\boldsymbol{V}_2^n$

(39)

${\boldsymbol{D}_{10}}\boldsymbol{V}_1^n = {\left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right)^n}{\boldsymbol{D}_{10}}$

(40)

Now we can use V₂ to calculate V₁^N₃:

$\begin{gathered} \boldsymbol{V}_1^{{N_3}} = {\boldsymbol{V}_1}\boldsymbol{V}_{^1}^{{N_3} - 3} = \hfill \\ \left( {{a_1}{\boldsymbol{D}_1} + {a_2}{\boldsymbol{D}_2} + {a_3}{\boldsymbol{D}_3} + {a_4}{\boldsymbol{D}_4} + {a_5}{\boldsymbol{D}_5} + } \right. \hfill \\ \left. {{a_6}{\boldsymbol{D}_6} + {a_7}{\boldsymbol{D}_7} + {a_8}{\boldsymbol{D}_8} + {a_9}{\boldsymbol{D}_9} + {a_{19}}{\boldsymbol{D}_{19}}} \right)\boldsymbol{V}_1^{{N_3} - 1} = \hfill \\ {a_{10}}{\boldsymbol{D}_1} + {a_{11}}{\boldsymbol{D}_2} + {a_{12}}{\boldsymbol{D}_3} + {a_{13}}{\boldsymbol{D}_4} + {a_{14}}{\boldsymbol{D}_5} + \hfill \\ {a_{15}}{\boldsymbol{D}_6} + {a_{16}}{\boldsymbol{D}_7} + {a_{17}}{\boldsymbol{D}_8} + {a_{18}}{\boldsymbol{D}_9} + {a_{19}} \hfill \\ {\left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right)^{{N_3} - 1}}{\boldsymbol{D}_{10}} \hfill \\ \end{gathered} $

(41)

where ${a_1} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B + 2{J_3}}}{{kT}}}}$, ${a_2} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}$, ${a_3} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B + 2{J_3}}}{{kT}}}}$, ${a_4} = {e^{\frac{{ - 2{J_1} + 2{J_2}}}{{kT}}}}$, ${a_5} = \frac{1}{2}{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} + \frac{1}{2}{e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}$, ${a_6} = {e^{\frac{{ - 2{J_1} + 2{J_2}}}{{kT}}}}$, ${a_7} = {e^{\frac{{2{J_1} + 2{J_2} - 2\mu B - 2{J_3}}}{{kT}}}}$, ${a_8} = {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B}}{{kT}}}}$, ${a_9} = {e^{\frac{{2{J_1} + 2{J_2} + 2\mu B + 2{J_3}}}{{kT}}}}$, ${a_{19}}\frac{1}{2}\left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right)$. By using Eqs. (35), (37), (39), we can get a₁₀-a₁₈:

$\left[{\begin{array}{*{20}{c}} {{a_{10}}}&{{a_{11}}}&{{a_{12}}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {{a_1}}&{{a_2}}&{{a_3}} \end{array}} \right]\boldsymbol{V}_2^{{N_3} - 1}$

(42)

$\left[{\begin{array}{*{20}{c}} {{a_{13}}}&{{a_{14}}}&{{a_{15}}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {{a_4}}&{{a_5}}&{{a_6}} \end{array}} \right]\boldsymbol{V}_2^{{N_3} - 1}$

(43)

$\left[{\begin{array}{*{20}{c}} {{a_{17}}}&{{a_{18}}}&{{a_{19}}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {{a_7}}&{{a_8}}&{{a_9}} \end{array}} \right]\boldsymbol{V}_2^{{N_3} - 1}$

(44)

The trace of V₁^N₃ is:

$\begin{gathered} Tr\left( {V_2^{{N_3} - 1}} \right) = {\left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right)^{{N_3} - 1}} = \hfill \\ \left[{\begin{array}{*{20}{c}} {{a_1}}&{{a_2}}&{{a_3}} \end{array}} \right]V_2^{{N_3} - 1}\left[{\begin{array}{*{20}{c}} 1 \\ 0 \\ 0 \end{array}} \right] + 2\left[{\begin{array}{*{20}{c}} {{a_4}}&{{a_5}}&{{a_6}} \end{array}} \right]V_2^{{N_3} - 1}\begin{array}{*{20}{c}} 0 \\ {\left[1 \right]} \\ 0 \end{array} + \hfill \\ \left[{\begin{array}{*{20}{c}} {{a_7}}&{{a_8}}&{{a_9}} \end{array}} \right]V_2^{{N_3} - 1}\left[{\begin{array}{*{20}{c}} 0 \\ 0 \\ 1 \end{array}} \right] + \hfill \\ 2{a_{19}}{\left( {{e^{\frac{{ - 2{J_1} + 2{J_2} + 2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right)^{{N_3} - 1}} \hfill \\ \end{gathered} $

(45)

Here V₂^N₃-1 can be written as:

$\boldsymbol{V}_2^{{N_3} - 1} = G\left[{\begin{array}{*{20}{c}} {\alpha _1^{{N_3} - 1}}&{}&{} \\ {}&{\alpha _2^{{N_3} - 1}}&{} \\ {}&{}&{\alpha _3^{{N_3} - 1}} \end{array}} \right]{\boldsymbol{G}^{ - 1}}$

(46)

where α₁, α₂, α₃ are the eigenvalues of V₂; G is a 3×3 matrix consisted of 3 eigenvectors corresponding to α₁, α₂, α₃.

$\begin{gathered} {\alpha _1} = {e^{\frac{{2{J_2}}}{{kT}}}}\left[{{a_{21}} - \frac{{{a_{20}}}}{{3{a_{21}}}} + \frac{1}{3}\left( {{e^{\frac{{2{J_3} + 2{J_1} - 2\mu B}}{{kT}}}} + } \right.} \right.{e^{\frac{{2{J_3} - 2{J_1}}}{{kT}}}} + \hfill \\ \left. {\left. {{e^{\frac{{ - 2{J_3} - 2{J_1}}}{{kT}}}} + {e^{\frac{{2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}} \right)} \right] \hfill \\ \end{gathered} $

(47)

$\begin{gathered} {\alpha _2} = {e^{\frac{{2{J_2}}}{{kT}}}}\left[{{a_{21}}{e^{i\frac{{ - 2\pi }}{3}}} - \frac{{{a_{20}}}}{{3{a_{21}}{e^{i\frac{{ - 2\pi }}{3}}}}} + \frac{1}{3}\left( {{e^{\frac{{2{J_3} + 2{J_1} - 2\mu B}}{{kT}}}} + } \right.} \right.{e^{\frac{{2{J_3} - 2{J_1}}}{{kT}}}} + \hfill \\ \left. {\left. {{e^{\frac{{ - 2{J_3} - 2{J_1}}}{{kT}}}} + {e^{\frac{{2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}} \right)} \right] \hfill \\ \end{gathered} $

(48)

$\begin{gathered} {\alpha _3} = {e^{\frac{{2{J_2}}}{{kT}}}}\left[{{a_{21}}{e^{i\frac{{ - 2\pi }}{3}}} - \frac{{{a_{20}}}}{{3{a_{21}}{e^{i\frac{{ - 2\pi }}{3}}}}}} \right. + \frac{1}{3}\left( {{e^{\frac{{2{J_3} + 2{J_1} - 2\mu B}}{{kT}}}} + } \right. \hfill \\ \left. {\left. {{e^{\frac{{2{J_3} - 2{J_1}}}{{kT}}}} + {e^{\frac{{ - 2{J_3} - 2{J_1}}}{{kT}}}} + {e^{\frac{{2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}} \right)} \right] \hfill \\ \end{gathered} $

(49)

$\begin{gathered} {a_{20}} = {e^{\frac{{4{J_1}}}{{kT}}}}\left[{\left( {{e^{\frac{{4{J_3}}}{{kT}}}} - {e^{\frac{{ - 4{J_3}}}{{kT}}}}} \right) - \frac{1}{3}{e^{\frac{{4{J_3}}}{{kT}}}}{{\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}}} \right)}^2}} \right] + \hfill \\ \frac{1}{3}{e^{\frac{{2{J_3}}}{{kT}}}}\left( {{e^{\frac{{2\mu B}}{{kT}}}} + {e^{\frac{{ - 2\mu B}}{{kT}}}}} \right)\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - 5{e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right) - \hfill \\ \frac{1}{3}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} + {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}{e^{\frac{{4{J_1}}}{{kT}}}} \hfill \\ \end{gathered} $

(50)

${a_{21}} = \sqrt[3]{{ - \frac{1}{2}{a_{22}} \pm \sqrt {{a_{23}}} }}$

(51)

$\begin{gathered} {a_{22}} = {e^{\frac{{6{J_1}}}{{kT}}}}\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}}} \right)\left[{\frac{1}{3}} \right.{e^{\frac{{2{J_3}}}{{kT}}}}\left( {{e^{\frac{{4{J_3}}}{{kT}}}} - {e^{\frac{{ - 4{J_3}}}{{kT}}}}} \right) - \hfill \\ \left. {\frac{2}{{27}}{e^{\frac{{6{J_3}}}{{kT}}}}\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}}} \right)} \right] + \hfill \\ \frac{1}{9}{e^{\frac{{2{J_1}}}{{kT}}}}{e^{\frac{{4{J_3}}}{{kT}}}}\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - 5{e^{{e^{\frac{{ - 2{J_3}}}{{kT}}}}}}} \right){\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}}} \right)^2} + \hfill \\ \frac{1}{3}{e^{\frac{{2{J_1}}}{{kT}}}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} + {e^{{e^{\frac{{ - 2{J_3}}}{{kT}}}}}}} \right)^2}\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{{e^{\frac{{ - 2{J_3}}}{{kT}}}}}}} \right) + \hfill \\ {\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{{e^{\frac{{ - 2{J_3}}}{{kT}}}}}}} \right)^3} + \frac{1}{9}{e^{\frac{{2{J_1}}}{{kT}}}}{e^{\frac{{2{J_3}}}{{kT}}}}\left( {{e^{\frac{{2\mu B}}{{kT}}}} + } \right. \hfill \\ \left. {{e^{\frac{{ - 2\mu B}}{{kT}}}}} \right)\left( {{e^{\frac{{2{J_3}}}{{kT}}}} + {e^{{e^{\frac{{ - 2{J_3}}}{{kT}}}}}}} \right)\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - 5{e^{{e^{\frac{{ - 2{J_3}}}{{kT}}}}}}} \right) - \hfill \\ {e^{\frac{{ - 6{J_1}}}{{kT}}}}\frac{2}{{27}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^3} \hfill \\ \end{gathered} $

(52)

$\begin{array}{l} {a_{23}}{ = ^1}108{e^{\frac{{12{J_1}}}{{kT}}}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}\left[ { - 4{e^{\frac{{ - 4{J_3}}}{{kT}}}} - } \right.\\ \left. {{e^{\frac{{4{J_3}}}{{kT}}}}{{\left( {{e^{\frac{{2\mu B}}{{kT}}}} + {e^{\frac{{ - 2\mu B}}{{kT}}}}} \right)}^2}} \right]{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2} + \\ \frac{1}{{54}}{e^{\frac{{8{J_1}}}{{kT}}}}{e^{\frac{{6{J_3}}}{{kT}}}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - } \right.\\ \left. {3{e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}}} \right) \times \left[ {{{\left( {{e^{\frac{{2\mu B}}{{kT}}}} + {e^{\frac{{ - 2\mu B}}{{kT}}}}} \right)}^2} - } \right.\\ \left. {4{e^{\frac{{ - 4{J_3}}}{{kT}}}}\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - 3{e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)\left( {{e^{\frac{{2{J_3}}}{{kT}}}} + {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)} \right] - \\ \frac{1}{{108}}{e^{\frac{{4{J_1}}}{{kT}}}}{e^{\frac{{8{J_3}}}{{kT}}}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}{\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}}} \right)^4} + \\ \frac{1}{{54}}{e^{\frac{{4{J_1}}}{{kT}}}}{e^{\frac{{4{J_3}}}{{kT}}}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}\left( {{e^{\frac{{4{J_3}}}{{kT}}}} + 14 - } \right.\\ \left. {24 - 14{e^{\frac{{ - 4{J_3}}}{{kT}}}} + {e^{\frac{{ - 8{J_3}}}{{kT}}}}} \right) + \\ \frac{1}{{54}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}{\left( {{e^{\frac{{2\mu B}}{{kT}}}} + {e^{\frac{{ - 2\mu B}}{{kT}}}}} \right)^2} \times \\ \left[ {{e^{\frac{{6{J_3}}}{{kT}}}}\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - 3{e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)} \right.{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - 2{e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2} + \\ \frac{1}{{108}}{e^{\frac{{ - 4{J_1}}}{{kT}}}}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}{\left( {{e^{\frac{{2{J_3}}}{{kT}}}} + {e^{\frac{{ - 2{J_3}}}{{kT}}}}} \right)^2}\left[ {4\left( {{e^{\frac{{4{J_3}}}{{kT}}}}} \right.} \right. - \\ \left. {\left. {{e^{\frac{{ - 4{J_3}}}{{kT}}}}} \right) - {e^{\frac{{4{J_3}}}{{kT}}}}{{\left( {{e^{\frac{{2\mu B}}{{kT}}}} + {e^{\frac{{ - 2\mu B}}{{kT}}}}} \right)}^2}} \right] \end{array}$

(53)

If we insert Eqs. (46)-(53) into Eq. (45), then we can get the expression of Tr(V₁^N₃):

$\begin{gathered} Tr\left( {\boldsymbol{V}_1^{{N_3}}} \right) = \alpha _1^{{N_3}} + \alpha _2^{{N_3}} + \alpha _3^{{N_3}} + \hfill \\ {\left( {{e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}} - {e^{\frac{{ - 2{J_1} + 2{J_2} - 2{J_3}}}{{kT}}}}} \right)^{{N_3}}} \hfill \\ \end{gathered} $

(54)

The eigenvalues of V₁ are λ₁, λ₂, λ₃, λ₄, the trace of V₁^N₃ is:

$Tr\left( {\boldsymbol{V}_1^{{N_3}}} \right) = \lambda _1^{{N_3}} + \lambda _2^{{N_3}} + \lambda _3^{{N_3}} + \lambda _4^{{N_3}}$

(55)

We can see the similarity between two forms of Tr(V₁^N₃), and α₁, α₂, α₃ are also the eigenvalues of V₁. Then we need to prove that V₂ contains the largest eigenvalue of V₁. To prove this, we set J₁=J₂=J₃=B=0, and then V₂ in Eq. (33) and V₁ becomes:

${\boldsymbol{V}_2}\left| {_{{J_1} = {J_2} = {J_3} = B = 0} = \left[{\begin{array}{*{20}{c}} 1&1&1 \\ 2&2&2 \\ 1&1&1 \end{array}} \right]} \right.$

(56)

${\boldsymbol{V}_1}\left| {_{{J_1} = {J_2} = {J_3} = B = 0} = \left[{\begin{array}{*{20}{c}} 1&1&1&1 \\ 1&1&1&1 \\ 1&1&1&1 \\ 1&1&1&1 \end{array}} \right]} \right.$

(57)

Obviously, the eigenvalues of V₂ in Eq. (56) are 4 and 0. The eigenvalues of V₁ in Eq.(57) are also 4 and 0.For other combinations of N₁, N₂, the first row elements of V₂|_J1=J2=J3=B=0 are 1, and the second row elements are f₁>0, and the third row elements are f₂>0, the jth row elements are f_j>0.The summation of f₁f₂f₃… is: f₁+f₂+f₃+…=2^N₁N₂-1. This means the eigenvalues of V₂|_J1=J2=J3=B=0 are 2^N₁N₂ and 0. The eigenvalues of V₁|_J1=J2=J3=B=0 are also 2^N₁N₂ and 0. Now we know V₂|_J1=J2=J3=B=0 contains the largest eigenvalue of V₁|_J1=J2=J3=B=0. When J₁, J₂, J₃, B change to non-zero values, the eigenvalues of 2^N₁N₂ and 0 will change too, the eigenvalue that originates from 2^N₁N₂ will keep being larger than other eigenvalues that originates from zero. For example, when N₁=N₂=1, the two eigenvalues are:

$\begin{gathered} {\lambda _1} = \frac{1}{2}{e^{\frac{{2{J_1} + 2{J_2}}}{{kT}}}}\left[{{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}} + } \right. \hfill \\ \left. {\sqrt {{e^{\frac{{2{J_3}}}{{kT}}}}\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}} - 2} \right) + 4{e^{\frac{{ - 2{J_3}}}{{kT}}}}} } \right] \hfill \\ \end{gathered} $

(58)

$\begin{gathered} {\lambda _2} = \frac{1}{2}{e^{\frac{{2{J_1} + 2{J_2}}}{{kT}}}}\left[{{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}} - } \right. \hfill \\ \left. {\sqrt {{e^{\frac{{2{J_3}}}{{kT}}}}\left( {{e^{\frac{{ - 2\mu B}}{{kT}}}} + {e^{\frac{{2\mu B}}{{kT}}}} - 2} \right) + 4{e^{\frac{{ - 2{J_3}}}{{kT}}}}} } \right] \hfill \\ \end{gathered} $

(59)

where λ₁ remains the largest eigenvalue when J₁, J₂, J₃, B change to non-zero values. Therefore V₂ still contains the largest eigenvalue of V₁. Obviously the order of V₂ is smaller than the order of V₁.

From above, we can summarize the steps to establish a smaller transfer matrix V₂. Firstly, we change the form of V₁ into:

${\boldsymbol{V}_1} = \sum\limits_n {{e^{\frac{{{p_n}{J_1} + {q_n}{J_2} + {r_n}\mu B + {s_n}{J_3}}}{{kT}}}}} {\boldsymbol{D_n}}$

(60)

where (p₁, q₁, r₁, s₁), (p₂, q₂, r₂, s₂) … are different combinations of integers; D₁, D₂, D₃… are a series of 2^N₁N₂×2^N₁N₂ matrixes. The elements of D_n are either 0 or 1. Then we try to establish a group of matrixes F₁, F₂, F₃…. The product of F_kV₁ can be expressed as the linear combinations of F₁, F₂, F₃… from the same group:

${\boldsymbol{F_k}}{\boldsymbol{V}_1} = \sum\limits_1 {{\boldsymbol{F}_l}{b_{k,l}}} $

(61)

where b_{k, l} are obviously the elements of the transfer matrix V₂.

To start off, we use F₁ as the first member of this group and then we will add new members to the group later:

${\boldsymbol{F}_1} = \left[{\begin{array}{*{20}{c}} 1&0& \cdots \\ 0&0& \cdots \\ \vdots & \vdots & \vdots \end{array}} \right]$

(62)

All elements of F₁ are zero except the leftmost element in the first row.

Assuming there are w matrixes of the group at some points: F₁, F₂, …F_w, if for a certain pair of F_k and D_l, F_kD_l cannot be expressed as the linear combination of F₁, F₂, …F_w, and then we construct a new matrix F_w+1 to add to the group as a new member:

${\boldsymbol{F}_{w + 1}} = {\boldsymbol{F}_k}{\boldsymbol{D}_l} - \sum\limits_{m = 1}^w {\frac{{{\boldsymbol{F}_m}\left[{\left( {{\boldsymbol{F}_k}{D_l}} \right) \cdot {F_m}} \right]}}{{{\boldsymbol{F}_m} \cdot {\boldsymbol{F}_m}}}} $

(63)

where the sign · means the point product between two matrixes ${\boldsymbol{F}_m} \cdot {\boldsymbol{F}_n} = \sum\limits_{p,q} {{\boldsymbol{F}_{m.p,q}}{\boldsymbol{F}_{n.p,q}}} $. We should keep in mind that F_k is from the group F₁, F₂, …F_w, and D_l can be any matrix from the series D₁, D₂, D₃….

Then we try to find another F_kD_l that cannot be expressed as the linear combination of F₁, F₂, …F_w+1. Repeat this process over and over. Eventually we will come to a point that for any pair of F_k and D_l, F_kD_l can be expressed as the linear combinations of F₁, F₂, …F_Nc, this means the group of matrixes is complete. Here we have used N_c to denote the total number of the complete group. N_c is also the order of the transfer matrix V₂.

All elements of F₁, F₂, …F_Nc are zero except the first row. For simplicity, we can use 1×2^N₁N₂ vectors F₁^′, F₂^′, …F_Nc^′ instead of the 2^N₁N₂×2^N₁N₂ matrixes. For the convenience of calculation, we make the elements of F₁^′, F₂^′, …F_Nc^′ either 0 or 1. For example, assuming there is a pair of F_j^′ and F_k^′ like this:

${{F'}_j} = \left[{\begin{array}{*{20}{c}} 1&9&1&1 \end{array}} \right]$

(64)

${{F'}_k} = \left[{\begin{array}{*{20}{c}} 3&{ - 1}&3&3 \end{array}} \right]$

(65)

For these two vectors, the elements of the first, third and fourth spots have the same value, but the second spot has a different value. So for F_j^′, we set the value of the second spot at zero. But for F_k^′, we set the value of the second spot at 1, we turn the elements of other spots into zero:

${{\boldsymbol{F}'}_j} = \left[{\begin{array}{*{20}{c}} 1&0&1&1 \end{array}} \right]$

(66)

${{\boldsymbol{F}'}_k} = \left[{\begin{array}{*{20}{c}} 0&1&0&0 \end{array}} \right]$

(67)

The principle is: the 2^N₁N₂ spots of the vectors are divided into N_c areas. Each vector has a different corresponding area. For a certain vector, all elements of that corresponding area are 1, other spots are zero. For a certain spot, there is only one F_k^′ of which the element of that spot is 1.

After this modification, all elements of F₁^′, F₂^′, …F_Nc^′ are either 0 or 1.

The elements of V₂ are:

${b_{k,m}} = \left[{{{\left( {{{\boldsymbol{F}'}_{\rm{k}}}{\boldsymbol{V}_1}} \right)}^{\rm{T}}}{{\boldsymbol{F}'}_m}} \right]/\left( {{{\boldsymbol{F}'}_m}^T{{\boldsymbol{F}'}_m}} \right)$

(68)

In the case of N₁=N₂=1, V₂ is identical to V₁. If N₁=2, N₂=1, V₂ is:

$\begin{array}{l} {V_2}\left| {_{{N_1} = 2,{N_2} = 1} = } \right.\\ {e^{\frac{{2{J_2}}}{{kT}}}}\left[{\begin{array}{*{20}{c}} {{e^{\frac{{2{J_1} + 2{J_2} - 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_1} - 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_3} + 2{J_1} - 2\mu B}}{{kT}}}}}\\ {2{e^{\frac{{ - 2{J_1}}}{{kT}}}}}&{{e^{\frac{{2{J_3} - 2{J_1}}}{{kT}}}} + {e^{\frac{{ - 2{J_3} - 2{J_1}}}{{kT}}}}}&{2{e^{\frac{{ - 2{J_1}}}{{kT}}}}}\\ {{e^{\frac{{ - 2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_1} + 2\mu B}}{{kT}}}}}&{{e^{\frac{{2{J_3} + 2{J_1} + 2\mu B}}{{kT}}}}} \end{array}} \right]69\\ \end{array}$

(69)

If N₁=3, N₂=1, V₂ is:

${V_2}\left| {_{{N_1} = 2,{N_2} = 1} = } \right.\\{e^{\frac{{3{J_2}}}{{kT}}}}\left[{\begin{array}{*{20}{c}} {{e^{\frac{{2{J_1} + 2{J_2} + 3{J_3}}}{{kT}}}}}&{{e^{\frac{{3{J_1} - 3\mu B + {J_3}}}{{kT}}}}}&{{e^{\frac{{3{J_1} - 3\mu B - {J_3}}}{{kT}}}}}&{{e^{\frac{{3{J_1} - 3\mu B - 3{J_3}}}{{kT}}}}}\\ {{e^{\frac{{ - {J_1} - \mu B + {J_3}}}{{kT}}}}}&{2{e^{\frac{{ - {J_1} - \mu B - {J_3}}}{{kT}}}} + {e^{\frac{{ - {J_1} - \mu B + 3{J_3}}}{{kT}}}}}&{2{e^{\frac{{ - {J_1} - \mu B + {J_3}}}{{kT}}}} + {e^{\frac{{ - {J_1} - \mu B - 3{J_3}}}{{kT}}}}}&{3{e^{\frac{{ - {J_1} - \mu B - {J_3}}}{{kT}}}}}\\ {3{e^{\frac{{ - {J_1} + \mu B + {J_3}}}{{kT}}}}}&{2{e^{\frac{{ - {J_1} + \mu B + {J_3}}}{{kT}}}} + {e^{\frac{{ - {J_1} + \mu B - 3{J_3}}}{{kT}}}}}&{2{e^{\frac{{ - {J_1} + \mu B - {J_3}}}{{kT}}}} + {e^{\frac{{ - J + \mu B + 3{J_3}}}{{kT}}}}}&{3{e^{\frac{{ - {J_1} + \mu B + {J_3}}}{{kT}}}}}\\ {{e^{\frac{{3{J_1} + 3\mu B - 3{J_3}}}{{kT}}}}}&{{e^{\frac{{3{J_1} + 3\mu B - {J_3}}}{{kT}}}}}&{{e^{\frac{{3{J_1} + 3\mu B + {J_3}}}{{kT}}}}}&{{e^{\frac{{3{J_1} + 3\mu B + 3{J_3}}}{{kT}}}}} \end{array}} \right]$

(70)

Here are the orders of V₂ in some cases:

N₁=3, N₂=1: N_C=4

N₁=4, N₂=1: N_C=6

N₁=2, N₂=2: N_C=7

N₁=5, N₂=1: N_C=8

N₁=6, N₂=1: N_C=13

N₁=3, N₂=2: N_C=13

N₁=7, N₂=1: N_C=18

N₁=8, N₂=1: N_C=30

N₁=4, N₂=2: N_C=34

N₁=9, N₂=1: N_C=46

N₁=3, N₂=3: N_C=36

N₁=10, N₂=1: N_C=78

N₁=5, N₂=2: N_C=78

N₁=11, N₂=1: N_C=126

N₁=12, N₂=1: N_C=224

N₁=6, N₂=2: N_C=237

N₁=4, N₂=3: N_C=158

N₁=13, N₂=1: N_C=380

N₁=14, N₂=1: N_C=687

N₁=7, N₂=2: N_C=687

N₁=15, N₂=1: N_C=1 224

N₁=5, N₂=3: N_C=707

N₁=16, N₂=1: N_C=2 250

N₁=8, N₂=2: N_C=2 299

N₁=4, N₂=4: N_C=1 459

N₁=17, N₂=1: N_C=4 112

N₁=18, N₂=1: N_C=7 685

N₁=9, N₂=2: N_C=7 685

N₁=6, N₂=3: N_C=4 236

N₁=19, N₂=1: N_C=14 310

From above, we can see that N_C is much smaller than 2^N₁N₂.

3 Calculation

The thermal-physical parameters we use in this paper are: J₁=J₂=J₃=1 kJ/mol, μ=9.274 009 49×10^-24 J/T, k=1.380 648 8×10^-23 J/K, here μ is the magnetic momentum of an atom, and k is the Boltzmann's constant.

Fig. 1 shows how the specific heat per atom C_pa=C/N₁N₂N₃ changes with temperature T and magnetic field intensity B. When B=0, C_pa gets to its largest value 1.94×10^-23 J/K at T=47.966 K. But the overall maximum value for C_pa is 2.142 907×10^-23J/K at T=49.178 K, B=±4.1715 T. We use T_C to represent the critical temperature where C_pa gets to its maximum value. T_C becomes larger when the value of magnetic field intensity B increases. If B=0, T_C has its smallest value.

Since we only calculate 3-dimensional Ising model of finite width, the value of C_pa is finite too. If N₁ or N₂ is infinite, the value of C_pa at T_C will become infinite, but T_C will approach to a certain value which is finite. T_C is the temperature where phase transition happens. We use T_X and B_X to represent the temperature and magnetic field intensity where C_pa get to its maximum value. T₀ is the temperature where C_pa get to its maximum value when B=0. Here's T_X, B_X and T₀ in some cases:

N₁=17, N₂=1: T_X=27.795 K, B_X=±1.372T, T₀=27.249 K

N₁=16, N₂=1: T_X=27.598 K, B_X=±1.135T, T₀=27.245 K

N₁=8, N₂=2: T_X=44.031 K, B_X=±3.508T, T₀=43.177 K

N₁=4, N₂=4: T_X=49.080 K, B_X=±4.426T, T₀=48.063 K

N₁=15, N₂=1: T_X=27.835 K, B_X=±1.603T, T₀=27.232 K

N₁=5, N₂=3: T_X=48.466 K, B_X=±5.222T, T₀=47.192 K

N₁=7, N₂=2: T_X=44.048 K, B_X=±4.168T, T₀=43.156 K

If N₂=1, the model turns into 2-dimensional model. T_X≈27.8 K, B_X≈1.3 K, T₀≈27.25 K. If the value of N₁N₂ does not change and N₂ increases, T_X, B_X and T₀ will become larger. We can see that in 3 dimensions, it needs higher temperature and magnetic field intensity to create phase transition. This is because the 3-dimensional interaction between atoms is stronger than that of 2-dimensional.

Fig. 2 shows the relationship between the magnetic momentum per atom M_pa=M/N₁N₂N₃ and the magnetic field intensity B. We can see that M_pa becomes smaller when temperature increases. If B is large enough, M_pa will increase slower and approach a certain value. This is called magnetic saturation.

Fig. 3 shows the maximum value of C_pa in different cases when B=0. In this picture, we use T^′=T-T₀ instead of normal temperature T so that the maximum point of C_pa in different cases can meet together to show the difference, and T₀ is the temperature where C_pa gets to its maximum value. We can see that when N₁N₂ increases, the maximum value of C_pa|_B=0 becomes larger. If N₁N₂ is infinite, C_pa will becomes infinite at critical temperature T_C. The value of T_C changes with magnetic field intensity.

4 Conclusions

In this paper, we use the 2^N₁N₂×2^N₁N₂ matrix V₁ to build a much smaller matrix V₂ which contains the largest eigenvalue of V₁. Then we obtain the largest eigenvalue of V₂ and calculate the specific heat per atom C_pa and the magnetic momentum per atom M_pa. There exists a pair of T and B to make C_pa get to its maximum value. When N₁N₂ increases, the maximum value of C_pa|_B=0 increases too.