Exact Expressions of Spin-Spin Correlation Functions of the Two-Dimensional Rectangular Ising Model on a Finite Lattice

Abstract

We employ the spinor analysis method to evaluate exact expressions of spin-spin correlation functions of the two-dimensional rectangular Ising model on a finite lattice, special process enables us to actually carry out the calculation process. We first present some exact expressions of correlation functions of the model with periodic-periodic boundary conditions on a finite lattice. The corresponding forms in the thermodynamic limit are presented, which show the short-range order. Then, we present the exact expression of the correlation function of the two farthest pair of spins in a column of the model with periodic-free boundary conditions on a finite lattice. Again, the corresponding form in the thermodynamic limit is discussed, from which the long-range order clearly emerges as the temperature decreases.

1. Introduction

Since the exact solution of the partition function in the absence of a magnetic field of the two-dimensional rectangular Ising model with periodic-periodic boundary conditions is obtained in the thermodynamic limit [1] and in finite-size systems [2], many authors have contributed to the knowledge of various aspects of this model, such as different boundary conditions, the arrangement modes of the spin lattice, surfaces, or mathematical methods, etc. [3,4,5,6,7].

Besides the partition function of the model, the calculations of spin-spin correlation functions are an important subject in the research of the two-dimensional Ising model. Some expressions of correlation functions in the thermodynamic limit have been obtained [3,4,5,8,9], and the case in a finite lattice has been studied [10,11,12].

The determination of exact expressions of the partition function and spin-spin correlation functions of the model on a finite lattice is not only a theoretical subject; the results obtained can also be used in the research of finite-size scaling, finite-size corrections, and boundary effects [7,13,14,15,16].

In this paper, we present some exact expressions of spin-spin correlation functions of the two-dimensional rectangular Ising model on a finite lattice by employing the spinor analysis method [2].

In Section 2, for the model with $L$ rows and $N$ columns and periodic-periodic boundary conditions (Onsager’s lattice), we calculate some exact expressions of the correlation function $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 1+Q}\rangle$ and compare the corresponding forms in the thermodynamic limit obtained here with known results presented in Reference [9]. The investigation in Section 2 shows the main steps, key points, and problems of the approach used in this paper. Since the whole process are complex, here we outline the approach.

(1) Although any spin-spin correlation functions $\langle {\sigma }_{l , n}{\sigma }_{l^{\prime } , n^{\prime }}\rangle$ can be expressed by matrices, and the matrices belong to spin representatives [8], here we only consider $\langle {\sigma }_{l , 1}{\sigma }_{l , 1+Q}\rangle$, i.e., the correlation functions of pairwise spins in one column, of which exact expressions can be obtained by the spinor analysis method.

(2) When we employ the spinor analysis method to evaluate $\langle {\sigma }_{l , 1}{\sigma }_{l , 1+Q}\rangle$, it is very difficult to find the exact eigenvalues of the corresponding rotation matrix; however, the operators of the first derivative and limit in the expressions of $\langle {\sigma }_{l , 1}{\sigma }_{l , 1+Q}\rangle$ (see (8)) allow us to obtain exact expressions of $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 1+Q}\rangle$ by only finding approximate eigenvalues of the rotation matrix.

(3) We employ Rayleigh-Schrodinger Perturbation Theory (RSPT) in quantum mechanics and change “finding eigenvalue up to $Q$-th order” to “finding eigenvalue through $Q$ times first-order approximation” (see the discussions in Section 2.5) to find approximate eigenvalues of the rotation matrix. On the one hand, this approximate method enables us to actually carry out the calculation process. On the other hand, since RSPT is irregular, the approximate method is in fact incapable when Q is very larger. Hence, by this approach we can only obtain exact expressions of correlation functions when Q is a small number, for example, that of $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 2}\rangle ,$ $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 3}\rangle ,$ $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 4}\rangle , \cdots$, etc., which belong to the short-range order.

What is more interesting is the long-range order, as it is closely related to properties of the phase transformation of the system. To obtain the correlation function that can reveal the long-range order of the model on a finite lattice, we turn to the model with $L$ rows and $N$ columns and periodic boundary condition in the horizontal direction and free boundary condition in the vertical direction. For this model, because of the free boundary condition in a column, and ${\sigma }_{l , N}$ is the farthest spin of ${\sigma }_{l , 1}$ in the column, we forecast that the correlation functions $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$, $\langle {\sigma }_{l , 1}{\sigma }_{l , N-1}\rangle$, $\langle {\sigma }_{l , 1}{\sigma }_{l , N-2}\rangle , \cdots$ will display the long-range order.

For the model with periodic-free boundary conditions, if we write $\langle {\sigma }_{l , 1}{\sigma }_{l , 1+Q}\rangle$ in matrix forms and use the method presented in Section 2, then we still can only obtain exact expressions of $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 2}\rangle ,$ $\langle {\sigma }_{1 , 1}{\sigma }_{1 , 3}\rangle , \cdots$ but cannot obtain that of $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$, $\langle {\sigma }_{l , 1}{\sigma }_{l , N-1}\rangle , \cdots$. On the other hand, if we write $\langle {\sigma }_{1 , 1}{\sigma }_{1 , N-n}\rangle$ in matrix forms directly (see Formulas (52) and (53) in this paper), then the method presented in Section 2 is feasible to deal with the matrix forms of $\langle {\sigma }_{1 , 1}{\sigma }_{1 , N-n}\rangle$, and we therefore can obtain some exact expressions of $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$, $\langle {\sigma }_{l , 1}{\sigma }_{l , N-1}\rangle$, $\langle {\sigma }_{l , 1}{\sigma }_{l , N-2}\rangle , \cdots$.

However, to save space, in this paper we no longer discuss $\langle {\sigma }_{l , 1}{\sigma }_{l , N-1}\rangle$, $\langle {\sigma }_{l , 1}{\sigma }_{l , N-2}\rangle , \cdots$ but only evaluate $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$, for which all of the matrix forms, the corresponding rotation matrices, and the eigenvalue equations have been given by Reference [12]. Therefore, we only need to derive the exact expression of $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$ by employing the method presented in Section 2. (The reason why the determination of the exact expression of $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$ fails in Reference [12] is explained in Section 3.3).

After obtaining the exact expression of $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$, we discuss the properties of the expression of $\langle {\sigma }_{l , 1}{\sigma }_{l , N}\rangle$ in the thermodynamic limit, from which the long-range order emerges as the temperature decreases, as shown clearly.

2. Short-Range Order in Onsager’s Lattice

2.1. The Definition of the Spin-Spin Correlation Functions and Their Basic Properties

According to the definition of the spin-spin correlation functions, $\langle {\sigma }_{l, n}{\sigma }_{l+P, n+Q}\rangle$ of pairwise spins ${\sigma }_{l , n}$ and ${\sigma }_{l+P , n+Q}$ of Onsager’s lattice read:

$$\langle {\sigma }_{l, n}{\sigma }_{l+P, n+Q}\rangle =\frac{1}{Z}{\displaystyle \sum _{\left\{{\sigma }_{h , v}=\pm 1\right\}}{\sigma }_{l, n}{\sigma }_{l+P, n+Q}\exp \left(\frac{J^{\prime }}{kT}{\displaystyle \sum _{l^{\prime }=1}^L{\displaystyle \sum _{n^{\prime }=1}^N{\sigma }_{l^{\prime } , n^{\prime }}{\sigma }_{l^{\prime }+1 , n^{\prime }}}}\right)}\exp \left(\frac{J}{kT}{\displaystyle \sum _{l^{\prime }=1}^L{\displaystyle \sum _{n^{\prime }=1}^N{\sigma }_{l^{\prime } , n^{\prime }}{\sigma }_{l^{\prime } , n^{\prime }+1}}}\right) ,$$

(1)

where ${\sigma }_{L+1 , n}={\sigma }_{1 , n}$, ${\sigma }_{l, N+1}={\sigma }_{l , 1}$; $J^{\prime }(>0)$, and $J(>0)$ are the interaction constants for the horizontal and vertical directions, respectively;

$$Z=\sum _{\left\{{\sigma }_{h , v}=\pm 1\right\}}\exp \left(\frac{J^{\prime }}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\sigma }_{l , n}{\sigma }_{l+1 , n}}}\right)\exp \left(\frac{J}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\sigma }_{l , n}{\sigma }_{l , n+1}}}\right)$$

(2)

is the partition function in absence of a magnetic field of the model.

According to the periodic-periodic boundary conditions of Onsager’s lattice, it is easy to prove $\langle {\sigma }_{l,n}{\sigma }_{l+P,n+Q}\rangle =\langle {\sigma }_{1,1}{\sigma }_{1+P,1+Q}\rangle$; in this paper, we calculate $\langle {\sigma }_{l,n}{\sigma }_{l,n+Q}\rangle =\langle {\sigma }_{1,1}{\sigma }_{1,1+Q}\rangle$, i.e., we only calculate correlation functions of pairwise spins in one column.

Further, because of periodic-periodic boundary conditions, both ${\sigma }_{1 , 2}$ and ${\sigma }_{1 , N}$ are the closest spins of ${\sigma }_{1 , 1}$, and we therefore have $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle =\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$. Generally speaking, it is easy to prove $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle =\langle {\sigma }_{1, 1}{\sigma }_{1, 1+N-Q}\rangle$ in terms of periodic-periodic boundary conditions. Hence, ${\sigma }_{1 , 1+[N/2]}$ is the farthest spin of ${\sigma }_{1 , 1}$, where $\left[x\right]$ denotes the greatest integer not exceeding $x$. Thus, we only need to calculate $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ for $Q=1, 2, \cdots , \left[\frac{N}{2}\right]$.

2.2. Some Results Concerning the Partition Function $Z$

From (1), we see that to obtain $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$, we need knowledge about the partition function $Z$ given by (2), from which we summarize some results presented in Reference [2] as follows.

$$Z={\left(2\sinh \frac{2J^{\prime }}{kT}\right)}^{\frac{LN}{2}}\tilde{Z}, \tilde{Z}=\mathrm{Tr}\left(\frac{1+{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{U}}^{(\uparrow )}+\frac{1-{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{U}}^{(\downarrow )}\right).$$

(3)

$\tilde{Z}$ can thus be obtained by finding the trace of a matrix, where ${\mathbf{\Gamma }}_{2N+1}$ is the matrix U defined by (15.68) in Reference [17], and the matrices $\frac{1+{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{U}}^{(\uparrow )}$ and $\frac{1-{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{U}}^{(\downarrow )}$ can be diagonalized at the same time.

On the other hand, both matrices ${\mathbf{U}}^{(\uparrow )}$ and ${\mathbf{U}}^{(\downarrow )}$ are spin representatives. By ${\mathsf{\omega }}^{(\uparrow )}$ and ${\mathsf{\omega }}^{(\downarrow )}$ we denote the corresponding rotation matrices of ${\mathbf{U}}^{(\uparrow )}$ and ${\mathbf{U}}^{(\downarrow )}$, respectively; both ${\mathsf{\omega }}^{(\uparrow )}$ and ${\mathsf{\omega }}^{(\downarrow )}$ can be diagonalized:

$$\begin{array}{l}{({\mathsf{\zeta }}^{(\uparrow \downarrow )})}^{+}{\mathsf{\omega }}^{(\uparrow \downarrow )}{\mathsf{\zeta }}^{(\uparrow \downarrow )}={\mathsf{\lambda }}^{(\uparrow \downarrow )};\\ {\mathsf{\lambda }}^{(\uparrow \downarrow )}=\mathrm{diag}\left[\begin{array}{cccc}{\mathsf{\lambda }}_1^{(\uparrow \downarrow )}& {\mathsf{\lambda }}_2^{(\uparrow \downarrow )}& \cdots & {\mathsf{\lambda }}_N^{(\uparrow \downarrow )}\end{array}\right] , {\mathsf{\lambda }}_n^{(\uparrow \downarrow )}=\mathrm{diag}\left[\begin{array}{cc}{\mathrm{e}}^{{\gamma }_n^{(\uparrow \downarrow )}}& {\mathrm{e}}^{-{\gamma }_n^{(\uparrow \downarrow )}}\end{array}\right](n=1, 2, \cdots , N) ,\end{array}$$

(4)

where ${\mathbf{A}}^{+}$ indicates a complex conjugate to a matrix $\mathbf{A}$,

$$\begin{array}{l}\cosh {\gamma }_n=\cosh 2K^{\prime }\cosh 2K-\cos {\phi }_n\sinh 2K^{\prime }\sinh 2K , {\mathrm{e}}^{-2K^{\prime }}=\tanh \frac{J^{\prime }}{kT} , K=\frac{J}{kT};\\ {\phi }_n^{(\uparrow )}=\frac{(2n-1)\mathsf{\pi }}{N}, {\phi }_n^{(\downarrow )}=\frac{2n\mathsf{\pi }}{N} (n=1, 2, \cdots , N=2M) .\end{array}$$

(5)

In the above formulas, ${\gamma }_n$ and ${\phi }_n$ are as the abbreviations for ${\gamma }_n^{(\uparrow \downarrow )}$ and ${\phi }_n^{(\uparrow \downarrow )}$, respectively. To save space, we use these abbreviation as far as possible in this paper.

Since properties of the partition function $Z$ vary between even and odd numbers $N$, we must calculate $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ separately in terms of whether $N$ is an even or odd number. In this paper we only consider the case that $N=2M$ as an even number (the case of $N=2M+1$ can be dealt with by the same approach). From Reference [2], we have:

$$\begin{array}{l}\tilde{Z}=\frac{1}{2}({\left({\displaystyle \prod _{m=1}^{M}2\cosh \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2+{\left({\displaystyle \prod _{m=1}^{M}2\sinh \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2\\ +\left(2\cosh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\right) \left(2\cosh \frac{L{\gamma }_M^{(\downarrow )}}{2}\right) {\left({\displaystyle \prod _{m=1}^{M-1}2\cosh \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)}^2\\ +\left(2\sinh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\right) \left(2\sinh \frac{L{\gamma }_M^{(\downarrow )}}{2}\right) {\left({\displaystyle \prod _{m=1}^{M-1}2\sinh \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)}^2).\end{array}$$

(6)

When $N=2M$, ${\gamma }_n$ and ${\phi }_n$ have properties:

$$\begin{array}{l}{\gamma }_m={\gamma }_{m^{\prime }}, {\phi }_m=2\mathsf{\pi }-{\phi }_{m^{\prime }}, m^{\prime }=\{\begin{array}{ll}2M-(m-1) ,& \mathrm{for} (\uparrow ) ;\\ 2M-m ,& \mathrm{for} (\downarrow ) ,\end{array}m=\{\begin{array}{ll}1, 2, \cdots , M ,& \mathrm{for} (\uparrow ) ;\\ 1, 2, \cdots , M-1 , & \mathrm{for} (\downarrow ) ,\end{array}\\ {\gamma }_M^{(\downarrow )}=2(K+K^{\prime }) , {\gamma }_{2M}^{(\downarrow )}=2(K-K^{\prime }) , {\phi }_M^{(\downarrow )}=\mathsf{\pi } , {\phi }_{2M}^{(\downarrow )}=2\mathsf{\pi } ,\\ 0<{\gamma }_1^{(\uparrow )}<{\gamma }_2^{(\uparrow )}<\cdots <{\gamma }_M^{(\uparrow )} , 0<\left|{\gamma }_{2M}^{(\downarrow )}\right|<{\gamma }_1^{(\downarrow )}<{\gamma }_2^{(\downarrow )}<\cdots <{\gamma }_{M-1}^{(\downarrow )}<{\gamma }_M^{(\downarrow )} .\end{array}$$

(7)

2.3. Writing $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ in Matrix Form

Since ${\sigma }_{h , v}^2=1$, we have ${\sigma }_{1 , 1}{\sigma }_{1 , 1+Q}={\sigma }_{1 , 1}{\sigma }_{1 , 2}\cdot {\sigma }_{1 , 2}{\sigma }_{1 , 3}\cdot \cdots \cdot {\sigma }_{1 , Q-1}{\sigma }_{1 , Q}\cdot {\sigma }_{1 , Q}{\sigma }_{1 , 1+Q}$, and, further, considering ${\sigma }_{1 , q}{\sigma }_{1 , 1+q}=\underset{{\varphi }_q\to 0}{\lim }\frac{\partial {\mathrm{e}}^{{\varphi }_q{\sigma }_{1 , q}{\sigma }_{1 , 1+q}}}{\partial {\varphi }_q}$, according to (1), $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ can be written in the form:

$$\begin{array}{l}\langle {\sigma }_{1 , 1}{\sigma }_{1 , 1+Q}\rangle =\frac{1}{Z}\left({\displaystyle \prod _{q=1}^Q\underset{{\varphi }_q\to 0}{\lim }\frac{\partial }{\partial {\varphi }_q}}\right)Y_Q ,\\ Y_Q={\displaystyle \sum _{\left\{{\sigma }_{h , v}=\pm 1\right\}}\exp \left({\displaystyle \sum _{q=1}^Q{\varphi }_q{\sigma }_{1 , q}{\sigma }_{1 , 1+q}}\right)}\exp \left(\frac{J^{\prime }}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\sigma }_{l , n}{\sigma }_{l+1 , n}}}\right)\exp \left(\frac{J}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\sigma }_{l , n}{\sigma }_{l , n+1}}}\right) ,\end{array}$$

(8)

Also, by a standard method [1,17] we can further write $Y_Q$ in the form:

$$Y_Q={\left(2\sinh \frac{2J^{\prime }}{kT}\right)}^{\frac{LN}{2}}{\tilde{Y}}_Q, {\tilde{Y}}_Q=\mathrm{Tr}\left(\frac{1+{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{V}}^{(\uparrow )}+\frac{1-{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{V}}^{(\downarrow )}\right).$$

(9)

${\tilde{Y}}_Q$ can thus be obtained by finding the trace of a matrix. Here we are not going to write out the explicit expressions of ${\mathbf{V}}^{(\uparrow )}$ and ${\mathbf{V}}^{(\downarrow )}$, but only point out that:

• The matrices $\frac{1+{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{V}}^{(\uparrow )}$ and $\frac{1-{\mathbf{\Gamma }}_{2N+1}}{2}{\mathbf{V}}^{(\downarrow )}$ can be diagonalized at the same time.
• Both ${\mathbf{V}}^{(\uparrow )}$ and ${\mathbf{V}}^{(\downarrow )}$ are spin representatives, whose corresponding rotation matrices are ${\mathsf{\zeta }}^{(\uparrow )}{\mathbf{H}}^{(\uparrow )}{({\mathsf{\zeta }}^{(\uparrow )})}^{+}$ and ${\mathsf{\zeta }}^{(\downarrow )}{\mathbf{H}}^{(\downarrow )}{({\mathsf{\zeta }}^{(\downarrow )})}^{+}$, respectively, where $\mathsf{\zeta }$ are introduced by (4),

$$\mathbf{H}={\mathbf{H}}^{(0)}+{\displaystyle \sum _{q=1}^Q\sinh 2{\varphi }_q{\mathbf{H}}^{(q)}}+2{\displaystyle \sum _{q=1}^Q{\sinh }^2{\varphi }_q{{\mathbf{H}}^{\mathbf{\prime }}}^{(q)}},$$

(10)

the forms of the matrices ${\mathbf{H}}^{(0)}$, ${\mathbf{H}}^{(q)}$, and ${{\mathbf{H}}^{\mathbf{\prime }}}^{(q)}$ can be expressed in terms of $2\times 2$ blocks ${\mathbf{H}}_{lm}^{(0)}$, ${\mathbf{H}}_{lm}^{(q)}$ and ${{\mathbf{H}}^{\mathbf{\prime }}}_{lm}^{(q)}$, $1\le l, m\le N(=2M)$, given by:

$${\mathbf{H}}_{lm}^{(0)}=\left[\begin{array}{cc}{\mathrm{e}}^{L{\gamma }_l}& 0\\ 0& {\mathrm{e}}^{-L{\gamma }_l}\end{array}\right]{\delta }_{lm} ,$$

(11)

$${\mathbf{H}}_{lm}^{(q)}={\mathrm{e}}^{-\mathrm{i}q\frac{(l-m)\mathsf{\pi }}{M}}{\mathbf{H}}_{lm} , {\mathbf{H}}_{lm}=\frac{1}{N}{\mathrm{e}}^{-\mathrm{i}\frac{{\theta }_l-{\theta }_m}{2}}\left[\begin{array}{cc}{\mathrm{e}}^{\frac{L({\gamma }_l+{\gamma }_m)}{2}}\cos \frac{{\theta }_l+{\theta }_m}{2}& -{\mathrm{ie}}^{\frac{L({\gamma }_l-{\gamma }_m)}{2}}\sin \frac{{\theta }_l+{\theta }_m}{2}\\ {\mathrm{ie}}^{-\frac{L({\gamma }_l-{\gamma }_m)}{2}}\sin \frac{{\theta }_l+{\theta }_m}{2}& -{\mathrm{e}}^{-\frac{L({\gamma }_l+{\gamma }_m)}{2}}\cos \frac{{\theta }_l+{\theta }_m}{2}\end{array}\right] ,$$

(12)

$${{\mathbf{H}}^{\mathbf{\prime }}}_{lm}^{(q)}={\mathrm{e}}^{-\mathrm{i}q\frac{(l-m)\mathsf{\pi }}{M}}{{\mathbf{H}}^{\mathbf{\prime }}}_{lm} , {{\mathbf{H}}^{\mathbf{\prime }}}_{lm}=\frac{1}{N}{\mathrm{e}}^{-\mathrm{i}\frac{{\theta }_l-{\theta }_m}{2}}\left[\begin{array}{cc}{\mathrm{e}}^{\frac{L({\gamma }_l+{\gamma }_m)}{2}}\cos \frac{{\theta }_l-{\theta }_m}{2}& {\mathrm{ie}}^{\frac{L({\gamma }_l-{\gamma }_m)}{2}}\sin \frac{{\theta }_l-{\theta }_m}{2}\\ {\mathrm{ie}}^{-\frac{L({\gamma }_l-{\gamma }_m)}{2}}\sin \frac{{\theta }_l-{\theta }_m}{2}& {\mathrm{e}}^{-\frac{L({\gamma }_l+{\gamma }_m)}{2}}\cos \frac{{\theta }_l-{\theta }_m}{2}\end{array}\right].$$

(13)

In (12) and (13), the quantities ${\theta }_n$ are defined by:

$$\begin{array}{l}\sinh {\gamma }_n\cos {\theta }_n=\cosh 2K^{\prime }\sinh 2K-\cos {\phi }_n\sinh 2K^{\prime }\cosh 2K ,\\ \sinh {\gamma }_n\sin {\theta }_n=\sin {\phi }_n\sinh 2K^{\prime }\hspace{1em}\hspace{1em}(n=1, 2, \cdots , N=2M) .\end{array}$$

(14)

When $N=2M$, ${\theta }_n$ have properties:

$${\theta }_m=2\mathsf{\pi }-{\theta }_{m^{\prime }}, {\theta }_M^{(\downarrow )}={\theta }_{2M}^{(\downarrow )}=0 .$$

(15)

In (15), the values of $m^{\prime }$ and $m$ are exactly the same as those in (7).

We can first evaluate the eigenvalues of the rotation matrices $\mathsf{\zeta }\mathbf{H}{\mathsf{\zeta }}^{+}$, and then obtain ${\tilde{Y}}_Q$ in terms of the spinor analysis method. Finally, we obtain $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ according to (8) and (9).

From (10) to (13), we see that $\mathbf{H}$ is a Hermitian conjugate matrix. Hence, the eigenvalues of both matrices $\mathsf{\zeta }\mathbf{H}{\mathsf{\zeta }}^{+}$ and $\mathbf{H}$ are the same, and we therefore can only evaluate the eigenvalues of the rotation matrix $\mathbf{H}$.

2.4. Basic Properties of the Eigenvalues and Eigenvectors of the Matrix $\mathbf{H}$

Any rotation matrix A has the following property: If $\tau$ is an eigenvalue of A, then ${\tau }^{-1}$ is also an eigenvalue of A [2]. Since $\mathsf{\zeta }\mathbf{H}{\mathsf{\zeta }}^{+}$ is a rotation matrix and the eigenvalues of $\mathsf{\zeta }\mathbf{H}{\mathsf{\zeta }}^{+}$ and $\mathbf{H}$ are the same, $\mathbf{H}$ should have the same property. In this sub-section we prove this conclusion.

The eigen equation of $\mathbf{H}$ reads:

$$\mathbf{H}\mathbf{\Psi }=\tau \mathbf{\Psi }.$$

(16)

By $\mathbf{\Psi }={\left[\begin{array}{ccccc}{\mathsf{\psi }}_1& \cdots & {\mathsf{\psi }}_n& \cdots & {\mathsf{\psi }}_N\end{array}\right]}^{\mathrm{T}},{\mathsf{\psi }}_n=\left[\begin{array}{c}{\psi }_n^{\Delta }\\ {\psi }_n^{\nabla }\end{array}\right]$, where ${\mathbf{A}}^{\mathrm{T}}$ means the transpose of a matrix A, we denote the eigenvector of $\mathbf{H}$, and, introducing $C_q^{\Delta }$, $C_q^{\nabla }$ $(q=1, 2, \cdots , Q)$ in terms of:

$$\begin{array}{l}\left[\begin{array}{c}{\psi }_n^{\Delta }\\ {\psi }_n^{\nabla }\end{array}\right]={\mathrm{e}}^{-\mathrm{i}\frac{{\theta }_n}{2}}\left[\begin{array}{cc}{\mathrm{e}}^{-\frac{L{\gamma }_n}{2}}& 0\\ 0& {\mathrm{e}}^{\frac{L{\gamma }_n}{2}}\end{array}\right] \left[\begin{array}{cc}\cos \frac{{\theta }_n}{2}& \mathrm{isin}\frac{{\theta }_n}{2}\\ \mathrm{isin}\frac{{\theta }_n}{2}& \cos \frac{{\theta }_n}{2}\end{array}\right] \\ \times \left[\begin{array}{cc}\frac{\tau \cos \mathrm{h}(L{\gamma }_n)-1+\tau \cos {\theta }_n\sin \mathrm{h}(L{\gamma }_n)}{{\tau }^2-2\tau \cosh (L{\gamma }_n)+1}& \mathrm{i}\frac{\tau \sin {\theta }_n\sin \mathrm{h}(L{\gamma }_n)}{{\tau }^2-2\tau \cosh (L{\gamma }_n)+1}\\ -\mathrm{i}\frac{\tau \sin {\theta }_n\sin \mathrm{h}(L{\gamma }_n)}{{\tau }^2-2\tau \cosh (L{\gamma }_n)+1}& \frac{\tau \cos \mathrm{h}(L{\gamma }_n)-1+\tau \cos {\theta }_n\sin \mathrm{h}(L{\gamma }_n)}{{\tau }^2-2\tau \cosh (L{\gamma }_n)+1}\end{array}\right] \left[\begin{array}{c}{\displaystyle \sum _{q=1}^Q{\mathrm{e}}^{-\mathrm{i}q{\phi }_n}{C}_q^{\Delta }}\\ {\displaystyle \sum _{q=1}^Q{\mathrm{e}}^{-\mathrm{i}q{\phi }_n}{C}_q^{\nabla }}\end{array}\right] ,\end{array}$$

we can prove that the eigen Equation (16) is equivalent to:

$$\coth {\varphi }_q\left[\begin{array}{c}{C}_q^{\Delta }\\ C_q^{\nabla }\end{array}\right]={\displaystyle \sum _{q^{\prime }=1}^Q\left[\begin{array}{cc}{A}_{+}(q-q^{\prime })& -B(q-q^{\prime })\\ -B(q-q^{\prime })& -A_{-}(q-q^{\prime })\end{array}\right]\left[\begin{array}{c}{C}_{q^{\prime }}^{\Delta }\\ C_{q^{\prime }}^{\nabla }\end{array}\right]}\hspace{1em}\hspace{1em}(q=1, 2, \cdots , Q) ;$$

(17)

where

$$\begin{array}{l}{A}_{\pm }(k)=\frac{1}{N}{\displaystyle \sum _{n=1}^N{\mathrm{e}}^{\mathrm{i}k\frac{2n\mathsf{\pi }}{N}}}\frac{{\tau }^2\pm 2\tau \cos {\theta }_n\sin \mathrm{h}(L{\gamma }_n)-1}{{\tau }^2-2\tau \cosh (L{\gamma }_n)+1}\\ =\{\begin{array}{ll}\frac{1}{M}{\displaystyle \sum _{m=1}^M\cos \frac{km\mathsf{\pi }}{M}}\frac{{\tau }^2\pm 2\tau \cos {\theta }_m^{(\uparrow )}\sin \mathrm{h}(L{\gamma }_m^{(\uparrow )})-1}{({\mathrm{e}}^{L{\gamma }_m^{(\uparrow )}}\tau -1) ({\mathrm{e}}^{-L{\gamma }_m^{(\uparrow )}}\tau -1)} ,& \mathrm{for} (\uparrow ) ;\\ \frac{1}{M}(\frac{{(-1)}^k}{2}\frac{{\mathrm{e}}^{\mp L{\gamma }_M^{(\downarrow )}}\tau +1}{{\mathrm{e}}^{\mp L{\gamma }_M^{(\downarrow )}}\tau -1}+\frac{1}{2}\frac{{\mathrm{e}}^{\mp L{\gamma }_{2M}^{(\downarrow )}}\tau +1}{{\mathrm{e}}^{\mp L{\gamma }_{2M}^{(\downarrow )}}\tau -1}+{\displaystyle \sum _{m=1}^{M-1}\cos \frac{km\mathsf{\pi }}{M}\frac{{\tau }^2\pm 2\tau \cos {\theta }_m^{(\downarrow )}\sin \mathrm{h}(L{\gamma }_m^{(\downarrow )})-1}{({\mathrm{e}}^{L{\gamma }_m^{(\downarrow )}}\tau -1) ({\mathrm{e}}^{-L{\gamma }_m^{(\downarrow )}}\tau -1)}}) ,& \mathrm{for} (\downarrow ),\end{array}\end{array}$$

$$B(k)=-\mathrm{i}\frac{1}{N}{\displaystyle \sum _{n=1}^N{\mathrm{e}}^{\mathrm{i}k\frac{2n\mathsf{\pi }}{N}}}\frac{2\tau \sin {\theta }_n\sin \mathrm{h}(L{\gamma }_n)}{{\tau }^2-2\tau \cosh (L{\gamma }_n)+1}=\{\begin{array}{ll}\frac{1}{M}{\displaystyle \sum _{m=1}^M\sin \frac{km\mathsf{\pi }}{M}}\frac{2\tau \sin {\theta }_m^{(\uparrow )}\sin \mathrm{h}(L{\gamma }_m^{(\uparrow )})}{{\tau }^2-2\tau \cosh (L{\gamma }_m^{(\uparrow )})+1},& \mathrm{for} (\uparrow ) ;\\ \frac{1}{M}{\displaystyle \sum _{m=1}^{M-1}\sin \frac{km\mathsf{\pi }}{M}}\frac{2\tau \sin {\theta }_m^{(\downarrow )}\sin \mathrm{h}(L{\gamma }_m^{(\downarrow )})}{{\tau }^2-2\tau \cosh (L{\gamma }_m^{(\downarrow )})+1},& \mathrm{for} (\downarrow ),\end{array}$$

where we have used (7) and (15). We see that all $A_{\pm }(k)$ and $B(k)$ are real functions and satisfy:

$$A_{\pm }(-k)=A_{\pm }(k), B(-k)=-B(k).$$

According to the above expressions and the properties of $A_{\pm }(k)$ and $B(k)$, we can conclude that if $\tau$ and $C_q^{\Delta }(\tau ) , C_q^{\nabla }(\tau ) (q=1, 2, \cdots , Q)$ satisfy (17), then ${\tau }^{\prime }={\tau }^{-1}$ and

$$C_q^{\Delta }({\tau }^{\prime })=C_0(\tau )C_q^{\nabla }(\tau ) , C_q^{\nabla }({\tau }^{\prime })=C_0(\tau )C_q^{\Delta }(\tau )\hspace{1em}(q=1, 2, \cdots , Q)$$

(18)

also satisfy (17), where $C_0(\tau )$ is an arbitrary function. From this discussion we not only prove the conclusion “If $\tau$ is an eigenvalue of the matrix H, then ${\tau }^{-1}$ is also an eigenvalue”, but also obtain the relation (18) between $\left\{ C_q^{\Delta }(\tau ), C_q^{\nabla }(\tau )\right\}$ and $\left\{ C_q^{\Delta }({\tau }^{-1}) , C_q^{\nabla }({\tau }^{-1})\right\}$. The conclusion and the relation (18) are useful to determine the forms of approximate eigenvalues and the expressions of the normalized eigenvectors of $\mathbf{H}$, as well as to calculate the determinant of the matrix consisting of the eigenvectors in the actual calculation process.

2.5. Approximate Method for Solving the Eigen Equation (16)

It is very difficult to find the exact eigenvalues of $\mathbf{H}$ by solving the eigen Equation (16). On the other hand, the operator ${\displaystyle \prod _{q=1}^Q\underset{{\varphi }_q\to 0}{\lim }\frac{\partial }{\partial {\varphi }_q}}$ in (8) allows us to ignore all terms whose orders are higher than ${\varphi }_q^1(={\varphi }_q)$ $(q=1, 2, \cdots , Q)$ in all eigenvalues of $\mathbf{H}$. According to this key property, we can obtain the exact expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ by only finding approximate eigenvalues of $\mathbf{H}$.

Concretely, as the first step, the term $2{\displaystyle \sum _{q=1}^Q{\sinh }^2{\varphi }_q{{\mathbf{H}}^{\mathbf{\prime }}}^{(q)}}$ with the factors ${\sinh }^2{\varphi }_q$ in (10) can be ignored, since ${\sinh }^2{\varphi }_q\approx {\varphi }_q^2$ have ${\varphi }_q^2$ order. Then, from (11) we see that ${\mathbf{H}}^{(0)}$ in (10) is a diagonal matrix, whose eigenvalues and eigenvectors are summarized in the following formulas:

$$\begin{array}{l}{\mathbf{H}}^{(0)}{\mathbf{\Psi }}_{n, \pm }^{(0)}={\mathrm{e}}^{\pm L{\gamma }_n}{\mathbf{\Psi }}_{n, \pm }^{(0)}\hspace{1em}(n=1, 2, \cdots , N=2M);\\ {\mathbf{\Psi }}_{n, \pm }^{(0)}={\left[\begin{array}{ccccc}{\mathsf{\psi }}_{n , 1 , \pm }^{(0)}& \cdots & {\mathsf{\psi }}_{n , m , \pm }^{(0)}& \cdots & {\mathsf{\psi }}_{n , N , \pm }^{(0)}\end{array}\right]}^{ \mathrm{T}}, {\mathsf{\psi }}_{n, m, \pm }^{(0)}=\left[\begin{array}{c}0\\ 0\end{array}\right] (m\ne n) ,\\ {\mathsf{\psi }}_{n, n, +}^{(0)}=\{\begin{array}{l}\left[\begin{array}{c}1\\ 0\end{array}\right] , \mathrm{for} \mathrm{the} \mathrm{eigenvalue} {\mathrm{e}}^{L{\gamma }_n};\\ \left[\begin{array}{c}0\\ 0\end{array}\right] , \mathrm{for} \mathrm{the} \mathrm{eigenvalue} {\mathrm{e}}^{-L{\gamma }_n},\end{array}{\mathsf{\psi }}_{n, n, -}^{(0)}=\{\begin{array}{l}\left[\begin{array}{c}0\\ 0\end{array}\right] , \mathrm{for} \mathrm{the} \mathrm{eigenvalue} {\mathrm{e}}^{L{\gamma }_n};\\ \left[\begin{array}{c}0\\ 1\end{array}\right] , \mathrm{for} \mathrm{the} \mathrm{eigenvalue} {\mathrm{e}}^{-L{\gamma }_n}.\end{array}\end{array}$$

(19)

which are as the zeroth order approximation of the eigen Equation (16).

The term ${\displaystyle \sum _{q=1}^Q\sinh 2{\varphi }_q{\mathbf{H}}^{(q)}}$ with the factors $\sinh 2{\varphi }_q\approx 2{\varphi }_q$ $(q=1, 2, \cdots , Q)$ in (10) can be regarded as a perturbation term. Then, by using RSPT, we can obtain the approximate eigenvalues of $\mathbf{H}$.

However, although what eigenvalues we need are only corrected to the ${\varphi }_q^1(={\varphi }_q)$ order $(q=1, 2, \cdots , Q)$, we must calculate the perturbation terms up to the $Q$-th order, not only for the first-order approximation, because all of the terms with the factor ${\displaystyle \prod _{q=1}^Q\sinh 2{\varphi }_q}$ appear in the $Q$-th order eigenvalues and are needed, which only include the ${\varphi }_q^1$ order for every ${\varphi }_q$.

However, if we calculate the eigenvalues up to the $Q$-th order by using RSPT, then not only the actual calculation process is very complex, but there are also many unwanted terms with factors ${\varphi }_q^k (k\ge 2)$, for example, the term with the factor ${\sinh }^{3}2{\varphi }_1{\displaystyle \prod _{q=4}^Q\sinh 2{\varphi }_q}$, in the $Q$-th order eigenvalues.

To take out those terms with factors ${\varphi }_q^k (k\ge 2)$, we change “finding the eigenvalue up to the $Q$-th order” to “finding the eigenvalue through $Q$ times first-order approximation”.

Concretely, since now $\mathbf{H}={\mathbf{H}}^{(0)}+{\displaystyle \sum _{q=1}^Q\sinh 2{\varphi }_q{\mathbf{H}}^{(q)}}$, we first consider the matrix ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$, in which the eigenvalues and eigenvectors ${\left\{ {\tau }^{(0)}, {\mathbf{\Psi }}^{(0)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}$ are given by (19) and $\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$ is as perturbation term. By only calculating first-order approximation we obtain all eigenvalues and eigenvectors ${\left\{ {\tau }^{(1)}, {\mathbf{\Psi }}^{(1)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$; therefore, all terms in ${\left\{ {\tau }^{(1)}, {\mathbf{\Psi }}^{(1)}\right\}}^{}$ only correct to the ${\varphi }_1^1$ order.

Then, we consider the matrix ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}+\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}$. Since now all eigenvalues and eigenvectors ${\left\{ {\tau }^{(1)}, {\mathbf{\Psi }}^{(1)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$ are known, we regard $\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}$ as a perturbation term, and, by only calculating the first-order approximation, obtain all eigenvalues and eigenvectors ${\left\{ {\tau }^{(2)}, {\mathbf{\Psi }}^{(2)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}+\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}$, in which all terms are only of the ${\varphi }_1^1$ and ${\varphi }_2^1$ orders. In particular, all of the terms with the factor $\sinh 2{\varphi }_1\sinh 2{\varphi }_2\approx {\varphi }_1{\varphi }_2$ remain.

Then, we consider the matrix ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}+\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}+\sinh 2{\varphi }_3{\mathbf{H}}^{(3)}$. Since now all eigenvalues and eigenvectors ${\left\{ {\tau }^{(2)}, {\mathbf{\Psi }}^{(2)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}+\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}$ are known, we regard $\sinh 2{\varphi }_3{\mathbf{H}}^{(3)}$ as a perturbation term, and, by only calculating the first-order approximation, obtain all eigenvalues and eigenvectors ${\left\{ {\tau }^{(3)}, {\mathbf{\Psi }}^{(3)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}+\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}+\sinh 2{\varphi }_3{\mathbf{H}}^{(3)}$, in which all terms are only of the ${\varphi }_1^1$, ${\varphi }_2^1$ and ${\varphi }_3^1$ orders. In particular, all of the terms with the factor $\sinh 2{\varphi }_1\sinh 2{\varphi }_2\sinh 2{\varphi }_3\approx {\varphi }_1{\varphi }_2{\varphi }_3$ remain, and, many unwanted terms with factors ${\sinh }^{2}2{\varphi }_1\sinh 2{\varphi }_2$, $\sinh 2{\varphi }_2{\sinh }^{2}2{\varphi }_3$, etc., do not appear in the eigenvalues of ${\tau }^{(3)}$.

We follow this approach up to $\sinh 2{\varphi }_Q{\mathbf{H}}^{(Q)}$ and every time we only calculate thr first-order approximation, which leads to the eigenvalues and eigenvectors ${\left\{ {\tau }^{(Q)}, {\mathbf{\Psi }}^{(Q)}\right\}}^{}$ of all terms being only of the ${\varphi }_q^1(={\varphi }_q)$ order. All of the terms with ${\displaystyle \prod _{q=1}^Q\sinh 2{\varphi }_q}$ remain, and at the same time those unwanted terms with ${\varphi }_q^k (k\ge 2)$ do not appear.

On the one hand, the above approximate method allows us to actually carry out the calculation process to find the eigenvalues and eigenvectors of $\mathbf{H}$. In particular, once we obtain ${\left\{ {\tau }^{(1)}, {\mathbf{\Psi }}^{(1)}\right\}}^{}$, we can obtain ${\tilde{Y}}_1$ by the spinor analysis method, as well as obtain $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle$ in terms of (8) and (9). Once we obtain ${\left\{ {\tau }^{(2)}, {\mathbf{\Psi }}^{(2)}\right\}}^{}$, we can obtain ${\tilde{Y}}_2$ by the spinor analysis method, and, further, obtain $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ in terms of (8) and (9),$\cdots$. Generally speaking, once we obtain ${\left\{ {\tau }^{(q)}, {\mathbf{\Psi }}^{(q)}\right\}}^{}$, we can obtain ${\tilde{Y}}_q$, and, further, obtain $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+q}\rangle$.

On the other hand, since RSPT is irregular, when Q is very large, e.g., $Q\approx \left[\frac{N}{2}\right]$, the above approach no longer functions. Hence, by this approach we can only obtain the exact expressions of correlation functions when Q is a small number, for example, $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle ,$ $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle ,$ $\langle {\sigma }_{1, 1}{\sigma }_{1, 4}\rangle , \cdots$, etc., which belong to the short-range order, but we cannot obtain the exact expressions of correlation functions when Q is larger, for example, $\langle {\sigma }_{1, 1}{\sigma }_{1, [N/2]+1}\rangle$, $\langle {\sigma }_{1, 1}{\sigma }_{1, [N/2]}\rangle$, $\langle {\sigma }_{1, 1}{\sigma }_{1, [N/2]-1}\rangle , \cdots$, etc., which belong to the long-range order.

2.6. Recurrence Formulas of the Eigenvalues and Eigenvectors $\left\{{\tau }^{(Q)}, {\mathbf{\Psi }}^{(Q)}\right\}$

According to the discussions in the above sub-section, we first regard $\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$ as a perturbation term, and, by using RSPT, evaluate eigenvalues and eigenvectors ${\left\{ {\tau }^{(1)}, {\mathbf{\Psi }}^{(1)}\right\}}^{}$ of the matrix ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$ up to the first-order approximation. However, according to (7), all eigenvalues ${\mathrm{e}}^{\pm L{\gamma }_n^{(\uparrow )}}$ are doubly-degenerate; and, except ${\mathrm{e}}^{\pm L{\gamma }_{2M}^{(\downarrow )}}$ and ${\mathrm{e}}^{\pm L{\gamma }_M^{(\downarrow )}}$, all the remaining eigenvalues ${\mathrm{e}}^{\pm L{\gamma }_n^{(\downarrow )}}$ are also doubly-degenerate. Hence, for doubly-degenerate eigenvalues of ${\mathbf{H}}^{(0)}$, we must use the degenerate perturbation theory; the results obtained up to ${\varphi }_1^1$ order are as follows.

$${\alpha }_{m , \uparrow \downarrow }^{(1)}=\frac{\sinh 2{\varphi }_1}{M}{\cos }^2\frac{{\theta }_m^{(\uparrow \downarrow )}}{2} , {\beta }_{m , \uparrow \downarrow }^{(1)}=\frac{\sinh 2{\varphi }_1}{M}{\sin }^2\frac{{\theta }_m^{(\uparrow \downarrow )}}{2} , {\delta }_M^{(1)}={\delta }_{2M}^{(1)}=\frac{\sinh 2{\varphi }_1}{2M},$$

(20)

$$\begin{array}{l}\Delta {\mathbf{\Psi }}_{m,\pm ,\uparrow \downarrow ,\mathrm{I}}^{(0)}=\frac{\cos \frac{{\theta }_m^{(\uparrow \downarrow )}}{2}}{\sqrt{2}}\left({\displaystyle \sum _{\begin{array}{l}l=1\\ l\ne m,l\ne m^{\prime }\end{array}}^{2M}\frac{\cos \frac{{\theta }_l^{(\uparrow \downarrow )}}{2}}{\sin \mathrm{h}\frac{L\left({\gamma }_l^{(\uparrow \downarrow )}-{\gamma }_m^{(\uparrow \downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\pm }^{(0)}}-{\displaystyle \sum _{l=1}^{2M}\frac{\mathrm{i}\sin \frac{{\theta }_l^{(\uparrow \downarrow )}}{2}}{\sinh \frac{L\left({\gamma }_l^{(\uparrow \downarrow )}+{\gamma }_m^{(\uparrow \downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\mp }^{(0)}}\right),\\ \Delta {\mathbf{\Psi }}_{m,\pm ,\uparrow \downarrow ,\mathrm{II}}^{(0)}=\frac{\mathrm{i}\sin \frac{{\theta }_m^{(\uparrow \downarrow )}}{2}}{\sqrt{2}}\left({\displaystyle \sum _{\begin{array}{l}l=1\\ l\ne m,l\ne m^{\prime }\end{array}}^{2M}\frac{\mathrm{i}\sin \frac{{\theta }_l^{(\uparrow \downarrow )}}{2}}{\sin \mathrm{h}\frac{L\left({\gamma }_l^{(\uparrow \downarrow )}-{\gamma }_m^{(\uparrow \downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\pm }^{(0)}}-{\displaystyle \sum _{l=1}^{2M}\frac{\cos \frac{{\theta }_l^{(\uparrow \downarrow )}}{2}}{\sinh \frac{L\left({\gamma }_l^{(\uparrow \downarrow )}+{\gamma }_m^{(\uparrow \downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\mp }^{(0)}}\right),\\ \Delta {\mathbf{\Psi }}_{M,\pm ,\downarrow }^{(0)}={\displaystyle \sum _{\begin{array}{l}l=1\\ l\ne M\end{array}}^{2M}\frac{\cos \frac{{\theta }_l^{(\downarrow )}}{2}}{\sin \mathrm{h}\frac{L\left({\gamma }_l^{(\downarrow )}-{\gamma }_M^{(\downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\pm }^{(0)}}-{\displaystyle \sum _{l=1}^{2M}\frac{\mathrm{i}\sin \frac{{\theta }_l^{(\downarrow )}}{2}}{\sinh \frac{L\left({\gamma }_l^{(\downarrow )}+{\gamma }_M^{(\downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\mp }^{(0)}},\\ \Delta {\mathbf{\Psi }}_{2M,\pm ,\downarrow }^{(0)}={\displaystyle \sum _{\begin{array}{l}l=1\\ l\ne 2M\end{array}}^{2M}\frac{\cos \frac{{\theta }_l^{(\downarrow )}}{2}}{\sin \mathrm{h}\frac{L\left({\gamma }_l^{(\downarrow )}-{\gamma }_{2M}^{(\downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\pm }^{(0)}}-{\displaystyle \sum _{l=1}^{2M}\frac{\mathrm{i}\sin \frac{{\theta }_l^{(\downarrow )}}{2}}{\sinh \frac{L\left({\gamma }_l^{(\downarrow )}+{\gamma }_M^{(\downarrow )}\right)}{2}}{\mathbf{\Psi }}_{l,\mp }^{(0)}}.\end{array}$$

(21)

In

Table 1. The eigenvalues and eigenvectors of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$ corrected to the ${\varphi }_1^1(={\varphi }_1)$ order.

Eigenvalue ${\left\{ {\mathit{\tau }}^{(1)}\right\}}^{}$	Eigenvector ${\left\{ {\mathbf{\Psi }}^{(1)}\right\}}^{}$
${\mathrm{e}}^{\pm (L{\gamma }_m^{(\uparrow \downarrow )}+{\alpha }_{m,\uparrow \downarrow }^{(1)})}$	$\frac{1}{\sqrt{2}}({\mathbf{\Psi }}_{m,\pm }^{(0)}-{\mathbf{\Psi }}_{m^{\prime },\pm }^{(0)})-\frac{\sinh 2{\varphi }_1}{2M}\Delta {\mathbf{\Psi }}_{m,\pm ,\uparrow \downarrow ,\mathrm{I}}^{(0)}$
${\mathrm{e}}^{\pm (L{\gamma }_m^{(\uparrow \downarrow )}-{\beta }_{m,\uparrow \downarrow }^{(1)})}$	$\frac{1}{\sqrt{2}}({\mathbf{\Psi }}_{m,\pm }^{(0)}+{\mathbf{\Psi }}_{m^{\prime },\pm }^{(0)})-\frac{\sinh 2{\varphi }_1}{2M}\Delta {\mathbf{\Psi }}_{m,\pm ,\uparrow \downarrow ,\mathrm{II}}^{(0)}$
${\mathrm{e}}^{\pm (L{\gamma }_M^{(\downarrow )}+{\delta }_M^{(1)})}$	${\mathbf{\Psi }}_{M,\pm }^{(0)}-\frac{\sinh 2{\varphi }_1}{4M}\Delta {\mathbf{\Psi }}_{M,\pm ,\downarrow }^{(0)}$
${\mathrm{e}}^{\pm (L{\gamma }_{2M}^{(\downarrow )}+{\delta }_{2M}^{(1)})}$	${\mathbf{\Psi }}_{2M,\pm }^{(0)}-\frac{\sinh 2{\varphi }_1}{4M}\Delta {\mathbf{\Psi }}_{2M,\pm ,\downarrow }^{(0)}$

, (20) and (21), the values of $m^{\prime }$ and $m$ are exactly the same as those in (7).

From Table 1, we see that all eigenvalues ${\left\{ {\tau }^{(1)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$ are nondegenerate. Hence, all degenerate eigenvalues of ${\mathbf{H}}^{(0)}$ are relieved by $\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}$. Thus, when we calculate the eigenvalues and eigenvectors ${\left\{ {\tau }^{(2)}, {\mathbf{\Psi }}^{(2)}\right\}}^{}$ of ${\mathbf{H}}^{(0)}+\sinh 2{\varphi }_1{\mathbf{H}}^{(1)}+\sinh 2{\varphi }_2{\mathbf{H}}^{(2)}$, we only need use nondegenerate perturbation theory; this is applicable up to $\sinh 2{\varphi }_Q{\mathbf{H}}^{(Q)}$. Further, since from ${\left\{ {\tau }^{(q)}, {\mathbf{\Psi }}^{(q)}\right\}}^{}$ to ${\left\{ {\tau }^{(q+1)}, {\mathbf{\Psi }}^{(q+1)}\right\}}^{}$ $(q=1, 2, \cdots , Q-1)$, we need only to calculate the first-order approximation in terms of $\sinh 2{\varphi }_{q+1}\approx 2{\varphi }_{q+1}$, and the corresponding recurrence formulas are:

$$\begin{array}{l}{\tau }_m^{(q+1)}={\tau }_m^{(q)}+\frac{\sinh 2{\varphi }_{q+1}}{2M}{({\mathbf{\Psi }}_m^{(q)})}^{ +}{\mathbf{H}}^{(q+1)}{\mathbf{\Psi }}_m^{(q)},\\ {\mathbf{\Psi }}_m^{(q+1)}={\mathbf{\Psi }}_m^{(q)}-\frac{\sinh 2{\varphi }_{q+1}}{2M}{\displaystyle \sum _{\begin{array}{l} l=1\\ l\ne m\end{array}}^{4M}\frac{{({\mathbf{\Psi }}_l^{(q)})}^{+}{\mathbf{H}}^{(q+1)}{\mathbf{\Psi }}_m^{(q)}}{{\tau }_l^{(q)}-{\tau }_m^{(q)}}{\mathbf{\Psi }}_l^{(q)}}\\ (m=1, 2, \cdots , 4M ; q=1, 2, \cdots , Q-1) ,\end{array},$$

(22)

In the calculation, all terms including the ${\varphi }_q^k(k\ge 2)$ order are ignored.

In principle, by following the above approach we obtain the eigenvalues ${\left\{ {\tau }^{(Q)}\right\}}^{}$. Furthermore, considering that up to the first-order approximation for ${\varphi }_q$, we have $1+C\sinh 2{\varphi }_q\approx {\mathrm{e}}^{C\sinh 2{\varphi }_q}$, the eigenvalues ${\left\{ {\tau }^{(Q)}\right\}}^{}$ of $\mathbf{H}$ can be denoted by the forms:

$${\mathrm{e}}^{\pm (L{\gamma }_m^{(\uparrow \downarrow )}+{\alpha }_{m , \uparrow \downarrow }^{(Q)})}{, \mathrm{e}}^{\pm (L{\gamma }_m^{(\uparrow \downarrow )}-{\beta }_{m , \uparrow \downarrow }^{(Q)})}{ , \mathrm{e}}^{\pm (L{\gamma }_M^{(\downarrow )}+{\delta }_M^{(Q)})}{, \mathrm{e}}^{\pm (L{\gamma }_{2M}^{(\downarrow )}+{\delta }_{2M}^{(Q)})},$$

where the value of $m$ is exactly the same as that in (7).

Based on the above forms of the eigenvalues ${\left\{ {\tau }^{(Q)}\right\}}^{}$ and using the spinor analysis method, we obtain:

$$\begin{array}{ll}{\tilde{Y}}_Q& =\left({\displaystyle \prod _{l=1}^{M}2\cosh \frac{L{\gamma }_l^{(\uparrow )}+{\alpha }_{l , \uparrow }^{\left(Q\right)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{M}2\cosh \frac{L{\gamma }_m^{(\uparrow )}-{\beta }_{m , \uparrow }^{\left(Q\right)}}{2}}\right)\\ & +\left({\displaystyle \prod _{l=1}^{M}2\sinh \frac{L{\gamma }_l^{(\uparrow )}+{\alpha }_{l , \uparrow }^{\left(Q\right)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{M}2\sinh \frac{L{\gamma }_m^{(\uparrow )}-{\beta }_{m , \uparrow }^{\left(Q\right)}}{2}}\right)\\ & +\left(2\cosh \frac{L{\gamma }_{2M}^{(\downarrow )}+{\delta }_{2M}^{\left(Q\right)}}{2}\right)\left(2\cosh \frac{L{\gamma }_M^{(\downarrow )}+{\delta }_M^{\left(Q\right)}}{2}\right)\left({\displaystyle \prod _{l=1}^{M-1}2\cosh \frac{L{\gamma }_l^{(\downarrow )}+{\alpha }_{l , \downarrow }^{\left(Q\right)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{M-1}2\cosh \frac{L{\gamma }_m^{(\downarrow )}-{\beta }_{m , \downarrow }^{\left(Q\right)}}{2}}\right)\\ & +\left(2\sinh \frac{L{\gamma }_{2M}^{(\downarrow )}+{\delta }_{2M}^{\left(Q\right)}}{2}\right)\left(2\sinh \frac{L{\gamma }_M^{(\downarrow )}+{\delta }_M^{\left(Q\right)}}{2}\right)\left({\displaystyle \prod _{l=1}^{M-1}2\sinh \frac{L{\gamma }_l^{(\downarrow )}+{\alpha }_{l , \downarrow }^{\left(Q\right)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{M-1}2\sinh \frac{L{\gamma }_m^{(\downarrow )}-{\beta }_{m , \downarrow }^{\left(Q\right)}}{2}}\right).\end{array}$$

(23)

Finally, according to (8) and (9), we obtain:

$$\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle =\frac{1}{\tilde{Z}}({\displaystyle \prod _{q=1}^Q\underset{{\varphi }_q\to 0}{\lim }\frac{\partial }{\partial {\varphi }_q}}){\tilde{Y}}_Q,$$

(24)

where $\tilde{Z}$ and ${\tilde{Y}}_Q$ are given by (6) and (23), respectively.

2.7. The Exact Expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle$ and $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ on a Finite Lattice

Although in Section 2.5 we presented a simplified approximate method, the actual calculation process of $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ is still complex; here, we only present the expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle$ and $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ directly.

When $Q=1$, substituting ${\alpha }_{m , \uparrow \downarrow }^{(1)} , {\beta }_{m , \uparrow \downarrow }^{(1)} , {\delta }_M^{(1)}$, and ${\delta }_{2M}^{(1)}$ given by (20) into (23), we obtain ${\tilde{Y}}_1$, and, further, we have:

$$\begin{array}{ll}\underset{{\varphi }_1\to 0}{\lim }\frac{\partial {\tilde{Y}}_1}{\partial {\varphi }_1}& =\left({\displaystyle \sum _{m=1}^M\frac{\cos {\theta }_m^{(\uparrow )}}{M}\tanh \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right) {\left({\displaystyle \prod _{l=1}^{M}2\cosh \frac{L{\gamma }_l^{(\uparrow )}}{2}}\right)}^2\\ & +\left({\displaystyle \sum _{m=1}^M\frac{\cos {\theta }_m^{(\uparrow )}}{M}\coth \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right) {\left({\displaystyle \prod _{l=1}^{M}2\sinh \frac{L{\gamma }_l^{(\uparrow )}}{2}}\right)}^2\\ & +\left(\frac{1}{2M}\tanh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}+\frac{1}{2M}\tanh \frac{L{\gamma }_M^{(\downarrow )}}{2}+{\displaystyle \sum _{m=1}^{M-1}\frac{\cos {\theta }_m^{(\downarrow )}}{M}\tanh \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)\\ & \times \left(2\cosh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\right) \left(2\cosh \frac{L{\gamma }_M^{(\downarrow )}}{2}\right) {\left({\displaystyle \prod _{l=1}^{M-1}2\cosh \frac{L{\gamma }_l^{(\downarrow )}}{2}}\right)}^2\\ & +\left(\frac{1}{2M}\coth \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}+\frac{1}{2M}\coth \frac{L{\gamma }_M^{(\downarrow )}}{2}+{\displaystyle \sum _{m=1}^{M-1}\frac{\cos {\theta }_m^{(\downarrow )}}{M}\coth \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)\\ & \times \left(2\sinh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\right) \left(2\sinh \frac{L{\gamma }_M^{(\downarrow )}}{2}\right) {\left({\displaystyle \prod _{l=1}^{M-1}2\sinh \frac{L{\gamma }_l^{(\downarrow )}}{2}}\right)}^2 .\end{array}$$

(25)

Substituting $\tilde{Z}$ given by (6) and the above expression into (24), we obtain the expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle$ of the model on a finite lattice:

$$\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle =\frac{1}{\tilde{Z}}\underset{{\varphi }_1\to 0}{\lim }\frac{\partial {\tilde{Y}}_1}{\partial {\varphi }_1}.$$

(26)

Then, using ${\left\{ {\tau }^{(1)}, {\mathbf{\Psi }}^{(1)}\right\}}^{}$ presented in Table 1, (20) and (21), and according to (22), we obtain ${\left\{ {\tau }^{(2)}\right\}}^{}$, which can be denoted by the forms:

$${\mathrm{e}}^{\pm (L{\gamma }_m^{(\uparrow \downarrow )}+{\alpha }_{m , \uparrow \downarrow }^{(2)})}{, \mathrm{e}}^{\pm (L{\gamma }_m^{(\uparrow \downarrow )}-{\beta }_{m , \uparrow \downarrow }^{(2)})}{ , \mathrm{e}}^{\pm (L{\gamma }_M^{(\downarrow )}+{\delta }_M^{(2)})}{, \mathrm{e}}^{\pm (L{\gamma }_{2M}^{(\downarrow )}+{\delta }_{2M}^{(2)})},$$

where

$$\begin{array}{ll}{\alpha }_{m, \uparrow \downarrow }^{(2)}& =\frac{\sinh 2{\varphi }_1+\sinh 2{\varphi }_2}{M}{\cos }^2\frac{{\theta }_m}{2}-\frac{\sinh 2{\varphi }_2}{M}{\sin }^2{\phi }_m\\ & +\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{4M^2}{\displaystyle \sum _{\begin{array}{c}k=1\\ k\ne m , k\ne m^{\prime }\end{array}}^{2M}\frac{{4\cos }^2\frac{{\theta }_m}{2}{\cos }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k-\sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{{\mathrm{e}}^{L({\gamma }_m-{\gamma }_k)}-1}}\\ & +\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{4M^2}{\displaystyle \sum _{k=1}^{2M}\frac{{4\cos }^2\frac{{\theta }_m}{2}{\sin }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k+\sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{{\mathrm{e}}^{L({\gamma }_m+{\gamma }_k)}-1}} ,\end{array}$$

$$\begin{array}{ll}{\beta }_{m, \uparrow \downarrow }^{(2)}& =\frac{\sinh 2{\varphi }_1+\sinh 2{\varphi }_2}{M}{\sin }^2\frac{{\theta }_m}{2}-\frac{\sinh 2{\varphi }_2}{M}{\sin }^2{\phi }_m\\ & +\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{4M^2}{\displaystyle \sum _{\begin{array}{c}k=1\\ k\ne m , k\ne m^{\prime }\end{array}}^{2M}\frac{- 4{\sin }^2\frac{{\theta }_m}{2}{\sin }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k+ \sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{1-{\mathrm{e}}^{-L({\gamma }_m-{\gamma }_k)}}}\\ & +\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{4M^2}{\displaystyle \sum _{k=1}^{2M}\frac{-4{\sin }^2\frac{{\theta }_m}{2}{\cos }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k- \sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{1-{\mathrm{e}}^{-L({\gamma }_m+{\gamma }_k)}}} ,\end{array}$$

$$\begin{array}{ll}{\delta }_M^{(2)}& =\frac{\sinh 2{\varphi }_1+\sinh 2{\varphi }_2}{2M}-\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{2M^2}{\displaystyle \sum _{\begin{array}{c} k=1\\ k\ne M\end{array}}^{2M}\frac{1}{{\mathrm{e}}^{L({\gamma }_M^{(\downarrow )}-{\gamma }_k^{(\downarrow )})}-1}{\cos }^2\frac{{\theta }_k^{(\downarrow )}}{2}\cos {\phi }_k^{(\downarrow )}}\\ & -\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{2M^2}{\displaystyle \sum _{k=1}^{2M}\frac{1}{{\mathrm{e}}^{L({\gamma }_M^{(\downarrow )}+{\gamma }_k^{(\downarrow )})}-1}{\sin }^2\frac{{\theta }_k^{(\downarrow )}}{2}\cos {\phi }_k^{(\downarrow )}},\end{array}$$

$$\begin{array}{ll}{\delta }_{2M}^{(2)}& =\frac{\sinh 2{\varphi }_1+\sinh 2{\varphi }_2}{2M}+\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{2M^2}{\displaystyle \sum _{k=1}^{2M-1}\frac{1}{{\mathrm{e}}^{L({\gamma }_{2M}^{(\downarrow )}-{\gamma }_k^{(\downarrow )})}-1}{\cos }^2\frac{{\theta }_k^{(\downarrow )}}{2}\cos {\phi }_k^{(\downarrow )}}\\ & +\frac{\sinh 2{\varphi }_1\sinh 2{\varphi }_2}{2M^2}{\displaystyle \sum _{k=1}^{2M}\frac{1}{{\mathrm{e}}^{L({\gamma }_{2M}^{(\downarrow )}+{\gamma }_k^{(\downarrow )})}-1}{\sin }^2\frac{{\theta }_k^{(\downarrow )}}{2}\cos {\phi }_k^{(\downarrow )}}.\end{array}$$

In the above expressions, the values of $m^{\prime }$ and $m$ are exactly the same as those in (7). We see that all terms with the factor $\sinh 2{\varphi }_1\sinh 2{\varphi }_2$ remain in the above expressions.

Substituting the above expressions of ${\alpha }_{m , \uparrow \downarrow }^{(2)}$, ${\beta }_{m , \uparrow \downarrow }^{(2)}$, ${\delta }_M^{(2)}$, and ${\delta }_{2M}^{(2)}$ into (23), we obtain ${\tilde{Y}}_2$, and, further, we have:

$$\begin{array}{l}\underset{{\varphi }_2\to 0}{\lim }\frac{\partial }{\partial {\varphi }_2}\underset{{\varphi }_1\to 0}{\lim }\frac{\partial {\tilde{Y}}_2}{\partial {\varphi }_1}\\ =\left({\displaystyle \sum _{m=1}^M\frac{W_m^{(\uparrow )}}{M}\tanh \frac{L{\gamma }_m^{(\uparrow )}}{2}}+\frac{1}{2M}{\displaystyle \sum _{m=1}^M\frac{1}{M}\frac{2{\cos }^2{\phi }_m^{(\uparrow )}-{\sin }^2{\theta }_m^{(\uparrow )}}{{\cosh }^2\frac{L{\gamma }_m^{(+)}}{2}}}+{\left({\displaystyle \sum _{m=1}^M\frac{1}{M}\cos {\theta }_m^{(\uparrow )}\tanh \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2\right)\\ \times {\left({\displaystyle \prod _{l=1}^{M}2\cosh \frac{L{\gamma }_l^{(\uparrow )}}{2}}\right)}^2\\ +\left({\displaystyle \sum _{m=1}^M\frac{W_m^{(\uparrow )}}{M}\coth \frac{L{\gamma }_m^{(\uparrow )}}{2}}-\frac{1}{2M}{\displaystyle \sum _{m=1}^M\frac{1}{M}\frac{2{\cos }^2{\phi }_m^{(\uparrow )}-{\sin }^2{\theta }_m^{(\uparrow )}}{{\sinh }^2\frac{L{\gamma }_m^{(+)}}{2}}}+{\left({\displaystyle \sum _{m=1}^M\frac{1}{M}\cos {\theta }_m^{(\uparrow )}\coth \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2\right)\\ \times {\left({\displaystyle \prod _{l=1}^{M}2\sinh \frac{L{\gamma }_l^{(\uparrow )}}{2}}\right)}^2\\ +(\left(\frac{W_{2M}^{(\downarrow )}}{2M}\tanh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}+\frac{W_M^{(\downarrow )}}{2M}\tanh \frac{L{\gamma }_M^{(\downarrow )}}{2}+{\displaystyle \sum _{m=1}^{M-1}\frac{W_m^{(\downarrow )}}{M}\tanh \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)\\ +\frac{1}{2M}\left(\frac{1}{2M}\frac{1}{{\cosh }^2\frac{L{\gamma }_{2M}^{(\downarrow )}}{2}}+\frac{1}{2M}\frac{1}{{\cosh }^2\frac{L{\gamma }_M^{(\downarrow )}}{2}}+{\displaystyle \sum _{m=1}^{M-1}\frac{1}{M}\frac{2{\cos }^2{\phi }_m^{(\downarrow )}-{\sin }^2{\theta }_m^{(\downarrow )}}{{\cosh }^2\frac{L{\gamma }_m^{(\downarrow )}}{2}}}\right)\\ +{\left(\frac{1}{2M}\tanh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}+\frac{1}{2M}\tanh \frac{L{\gamma }_M^{(\downarrow )}}{2}+{\displaystyle \sum _{m=1}^{M-1}\frac{1}{M}\cos {\theta }_m^{(\downarrow )}\tanh \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)}^2) \\ \times \left(2\cosh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\right)\left(2\cosh \frac{L{\gamma }_M^{(\downarrow )}}{2}\right){\left({\displaystyle \prod _{l=1}^{M-1}2\cosh \frac{L{\gamma }_l^{(\downarrow )}}{2}}\right)}^2\\ +(\left(\frac{W_{2M}^{(\downarrow )}}{2M}\coth \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}+\frac{W_M^{(\downarrow )}}{2M}\coth \frac{L{\gamma }_M^{(\downarrow )}}{2}+{\displaystyle \sum _{m=1}^{M-1}\frac{W_m^{(\downarrow )}}{M}\coth \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)\\ -\frac{1}{2M}\left(\frac{1}{2M}\frac{1}{{\sinh }^2\frac{L{\gamma }_{2M}^{(\downarrow )}}{2}}+\frac{1}{2M}\frac{1}{{\sinh }^2\frac{L{\gamma }_M^{(\downarrow )}}{2}}+{\displaystyle \sum _{m=1}^{M-1}\frac{1}{M}\frac{2{\cos }^2{\phi }_m^{(\downarrow )}-{\sin }^2{\theta }_m^{(\downarrow )}}{{\sinh }^2\frac{L{\gamma }_m^{(-)}}{2}}}\right)\\ +{\left(\frac{1}{2M}\coth \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}+\frac{1}{2M}\coth \frac{L{\gamma }_M^{(\downarrow )}}{2}+{\displaystyle \sum _{m=1}^{M-1}\frac{1}{M}\cos {\theta }_m^{(\downarrow )}\coth \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)}^2)\\ \times \left(2\sinh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\right)\left(2\sinh \frac{L{\gamma }_M^{(\downarrow )}}{2}\right){\left({\displaystyle \prod _{l=1}^{M-1}2\sinh \frac{L{\gamma }_l^{(\downarrow )}}{2}}\right)}^2 ,\end{array}$$

(27)

where $W_m^{(\uparrow \downarrow )}$ is introduced by:

$$\begin{array}{l}{W}_m=\frac{1}{2M}({\displaystyle \sum _{\begin{array}{c}k=1\\ k\ne m , k\ne m^{\prime }\end{array}}^{2M}\frac{{4\cos }^2\frac{{\theta }_m}{2}{\cos }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k-\sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{{\mathrm{e}}^{L({\gamma }_m-{\gamma }_k)}-1}}\\ +{\displaystyle \sum _{k=1}^{2M}\frac{{4\cos }^2\frac{{\theta }_m}{2}{\sin }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k+\sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{{\mathrm{e}}^{L({\gamma }_m+{\gamma }_k)}-1}}\\ +{\displaystyle \sum _{\begin{array}{c}k=1\\ k\ne m , k\ne m^{\prime }\end{array}}^{2M}\frac{ 4{\sin }^2\frac{{\theta }_m}{2}{\sin }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k- \sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{1-{\mathrm{e}}^{-L({\gamma }_m-{\gamma }_k)}}}\\ +{\displaystyle \sum _{k=1}^{2M}\frac{4{\sin }^2\frac{{\theta }_m}{2}{\cos }^2\frac{{\theta }_k}{2}\cos {\phi }_m\cos {\phi }_k+ \sin {\theta }_m\sin {\theta }_k\sin {\phi }_m\sin {\phi }_k}{1-{\mathrm{e}}^{-L({\gamma }_m+{\gamma }_k)}}}) ,\end{array}$$

(28)

where the value of $m^{\prime }$ is exactly the same as that in (7).

Substituting $\tilde{Z}$ given by (6) and (27) into (24), we obtain the expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ of the model on a finite lattice:

$$\langle {\sigma }_{1,1}{\sigma }_{1,3}\rangle =\frac{1}{\tilde{Z}}\underset{{\varphi }_2\to 0}{\lim }\frac{\partial }{\partial {\varphi }_2}\underset{{\varphi }_1\to 0}{\lim }\frac{\partial {\tilde{Y}}_2}{\partial {\varphi }_1}$$

(29)

2.8. The Expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle$ and $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ in the Thermodynamic Limit

We now consider the thermodynamic limit. First, if $L$ is very large, then according to (7) we have:

$$2\cosh \frac{L{\gamma }_m^{(\uparrow \downarrow )}}{2}\approx 2\sinh \frac{L{\gamma }_m^{(\uparrow \downarrow )}}{2}\approx \exp \frac{L{\gamma }_m^{(\uparrow \downarrow )}}{2},\tanh \frac{L{\gamma }_m^{(\uparrow \downarrow )}}{2}\approx \coth \frac{L{\gamma }_m^{(\uparrow \downarrow )}}{2}\approx 1\hspace{1em}(m=1,2\cdots M);$$

However, when the system crosses its critical temperature, ${\gamma }_{2M}^{(\downarrow )}=2(K-K^{\prime })$ changes sign, following which we therefore have:

$$2\cosh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\approx \exp \frac{L\left|{\gamma }_{2M}^{(\downarrow )}\right|}{2}, 2\sinh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\approx \frac{{\gamma }_{2M}^{(\downarrow )}}{\left|{\gamma }_{2M}^{(\downarrow )}\right|}\exp \frac{L\left|{\gamma }_{2M}^{(\downarrow )}\right|}{2},\tanh \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\approx \coth \frac{L{\gamma }_{2M}^{(\downarrow )}}{2}\approx \frac{{\gamma }_{2M}^{(\downarrow )}}{\left|{\gamma }_{2M}^{(-)}\right|}.$$

Hence, for $\tilde{Z}$ and $\underset{{\varphi }_1\to 0}{\lim }\frac{\partial {\tilde{Y}}_1}{\partial {\varphi }_1}$ given by (6) and (25), respectively, when $L$ is very large, we obtain:

$$\begin{array}{l}\tilde{Z}\approx \frac{1}{2}({\left({\displaystyle \prod _{m=1}^M\exp \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2+{\left({\displaystyle \prod _{m=1}^M\exp \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2\\ +\exp \frac{L\left|{\gamma }_{2M}^{(\downarrow )}\right|}{2}\exp \frac{L{\gamma }_M^{(\downarrow )}}{2}{\left({\displaystyle \prod _{m=1}^{M-1}\exp \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)}^2+\frac{{\gamma }_{2M}^{(\downarrow )}}{\left|{\gamma }_{2M}^{(\downarrow )}\right|}\exp \frac{L\left|{\gamma }_{2M}^{(\downarrow )}\right|}{2}\exp \frac{L{\gamma }_M^{(\downarrow )}}{2}{\left({\displaystyle \prod _{m=1}^{M-1}\exp \frac{L{\gamma }_m^{(\downarrow )}}{2}}\right)}^2)\\ \approx {\left({\displaystyle \prod _{m=1}^M\exp \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2\times \{\begin{array}{l}1 , K<K^{\prime } ;\\ 2 , K>K^{\prime } ,\end{array}\end{array}$$

(30)

$$\begin{array}{l}\underset{{\varphi }_1\to 0}{\lim }\frac{\partial {\tilde{Y}}_1}{\partial {\varphi }_1}\approx \left({\displaystyle \sum _{n=1}^M\frac{\cos {\theta }_n^{(\uparrow )}}{M}}\right) {\left({\displaystyle \prod _{l=1}^M\exp \frac{L{\gamma }_l^{(\uparrow )}}{2}}\right)}^2\\ +\frac{1}{2}\left(\frac{1}{2M}\frac{{\gamma }_{2M}^{(\downarrow )}}{\left|{\gamma }_{2M}^{(-)}\right|}+\frac{1}{2M}+{\displaystyle \sum _{n=1}^{M-1}\frac{\cos {\theta }_n^{(\downarrow )}}{M}}\right) \left(1+\frac{{\gamma }_{2M}^{(\downarrow )}}{\left|{\gamma }_{2M}^{(-)}\right|}\right) \left(\exp \frac{L\left|{\gamma }_{2M}^{(\downarrow )}\right|}{2}\right) \left(\exp \frac{L{\gamma }_M^{(\downarrow )}}{2}\right) {\left({\displaystyle \prod _{l=1}^{M-1}\exp \frac{L{\gamma }_l^{(\downarrow )}}{2}}\right)}^2\\ \approx \left({\displaystyle \sum _{n=1}^M\frac{\cos {\theta }_n^{(\uparrow )}}{M}}\right) {\left({\displaystyle \prod _{m=1}^M\exp \frac{L{\gamma }_m^{(\uparrow )}}{2}}\right)}^2\times \{\begin{array}{l}1 , K<K^{\prime } ;\\ 2 , K>K^{\prime } .\end{array}\end{array}$$

Substituting the above two expressions into (26), we obtain:

$$\begin{array}{l}\underset{\begin{array}{l}L\to \infty \\ N\to \infty \end{array}}{\lim }\langle {\sigma }_{1 , 1}{\sigma }_{1 , 2}\rangle =\underset{M\to \infty }{\lim }{\displaystyle \sum _{m=1}^M\frac{\cos {\theta }_m^{(\uparrow )}}{M}}=\underset{M\to \infty }{\lim }\left(\frac{1}{2}{\displaystyle \sum _{m=1}^M\frac{\cos {\theta }_m^{(\uparrow )}}{M}}+\frac{1}{2}{\displaystyle \sum _{m=M+1}^{2M}\frac{\cos {\theta }_m^{(\uparrow )}}{M}}\right)\\ =\underset{N\to \infty }{\lim }{\displaystyle \sum _{n=1}^N\frac{\cos {\theta }_n^{(\uparrow )}}{N}}={\displaystyle {\int }_0^1\mathrm{d}x}\cos \theta (\mathsf{\pi }x) ,\end{array}$$

(31)

where the function $\theta (\mathsf{\pi }x)$ in terms of (14) is defined by:

$$\begin{array}{l}\cos \theta (\mathsf{\pi }x)=\frac{\cosh 2K^{\prime }\sinh 2K-\cos (\mathsf{\pi }x)\sinh 2K^{\prime }\cosh 2K}{\sqrt{{(\cosh 2K^{\prime }\cosh 2K-\cos (\mathsf{\pi }x)\sinh 2K^{\prime }\sinh 2K)}^2-1}} ,\\ \sin \theta (\mathsf{\pi }x)=\frac{\sin (\mathsf{\pi }x)\sinh 2K^{\prime }}{\sqrt{{(\cosh 2K^{\prime }\cosh 2K-\cos (\mathsf{\pi }x)\sinh 2K^{\prime }\sinh 2K)}^2-1}} .\end{array}$$

(32)

The result (31) is in accordance with that in Reference [8].

According to the similar discussions, for $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ we first have:

$$\underset{\begin{array}{l}L\to \infty \\ N\to \infty \end{array}}{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle \approx \underset{\begin{array}{l}L\to \infty \\ M\to \infty \end{array}}{\lim }{\displaystyle \sum _{m=1}^M\frac{W_m^{(\uparrow )}}{M}}+{(\underset{M\to \infty }{\lim }{\displaystyle \sum _{m=1}^M\frac{\cos {\theta }_m^{(\uparrow )}}{M}})}^2.$$

(33)

We discuss the first term in (28) as an example to show how to calculate $\underset{L\to \infty }{\lim }{W}_m^{(+)}$. First, using (7) and (15), the first term in (28) can be written in the form:

$$\begin{array}{l}\frac{1}{2M}{\displaystyle \sum _{\begin{array}{c}k=1\\ k\ne m , k\ne m^{\prime }\end{array}}^{2M}\frac{{4\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}{\cos }^2\frac{{\theta }_k^{(\uparrow )}}{2}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}-\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}{{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}-1}}\\ =\frac{1}{M}{\displaystyle \sum _{\begin{array}{c} k=1\\ k\ne m\end{array}}^M\frac{{4\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}{\cos }^2\frac{{\theta }_k^{(\uparrow )}}{2}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}-\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}{{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}-1}}\\ =\frac{1}{M}{\displaystyle \sum _{k=1}^{m-1}\frac{{4\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}{\cos }^2\frac{{\theta }_k^{(\uparrow )}}{2}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}-\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}{{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}-1}}\\ +\frac{1}{M}{\displaystyle \sum _{k=m+1}^M\frac{{4\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}{\cos }^2\frac{{\theta }_k^{(\uparrow )}}{2}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}-\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}{{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}-1}} .\end{array}$$

According to (7), when $k<m$, ${\gamma }_k^{(\uparrow )}<{\gamma }_m^{(\uparrow )}$, $\underset{L\to \infty }{\lim }{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}=\infty$, and, thus, the first term in the above expression vanishes; when $k>m$, ${\gamma }_k^{(\uparrow )}>{\gamma }_m^{(\uparrow )}$, $\underset{L\to \infty }{\lim }{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}=0$, we therefore obtain:

$$\begin{array}{l}\underset{L\to \infty }{\lim }\frac{1}{2M}{\displaystyle \sum _{\begin{array}{c}k=1\\ k\ne m , k\ne m^{\prime }\end{array}}^{2M}\frac{{4\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}{\cos }^2\frac{{\theta }_k^{(\uparrow )}}{2}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}-\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}{{\mathrm{e}}^{L({\gamma }_m^{(\uparrow )}-{\gamma }_k^{(\uparrow )})}-1}}\\ =\frac{1}{M}{\displaystyle \sum _{k=m+1}^M\frac{{4\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}{\cos }^2\frac{{\theta }_k^{(\uparrow )}}{2}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}-\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}{-1}} .\end{array}$$

Using this method to deal with the remaining terms in $W_m^{(\uparrow )}$, we finally obtain:

$$\begin{array}{l}\underset{L\to \infty }{\lim }{W}_m^{(\uparrow )}=\frac{2}{M}({\displaystyle \sum _{k=m+1}^M\sin {\theta }_m^{(\uparrow )}\sin {\theta }_k^{(\uparrow )}\sin {\phi }_m^{(\uparrow )}\sin {\phi }_k^{(\uparrow )}}-{\displaystyle \sum _{k=m+1}^M\cos {\theta }_m^{(\uparrow )}\cos {\theta }_k^{(\uparrow )}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}}\\ +{\displaystyle \sum _{k=1}^{m+1}\cos {\phi }_m^{(\uparrow )}\cos {\phi }_k^{(\uparrow )}})+\frac{2}{M}\left(2{\sin }^2{\frac{{\theta }_m^{(\uparrow )}}{2}}^2\left({\cos }^2\frac{{\theta }_m^{(\uparrow )}}{2}-{\cos }^2{\phi }_m^{(\uparrow )}\right)-\cos {\phi }_m^{(\uparrow )}\cos {\phi }_{m+1}^{(\uparrow )}\right).\end{array}$$

As $M\to \infty$, the second term in the above expression vanishes, and, according to the definition of the Riemann integral, we have:

$$\begin{array}{ll}\underset{\begin{array}{l}L\to \infty \\ M\to \infty \end{array}}{\lim }{W}_m^{(\uparrow )}& =2{\displaystyle {\int }_x^1\mathrm{d}y}(\sin \theta (\mathsf{\pi }x)\sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }x)\sin (\mathsf{\pi }y)-\cos \theta (\mathsf{\pi }x)\cos \theta (\mathsf{\pi }y)\cos (\mathsf{\pi }x)\cos (\mathsf{\pi }y))\\ & +2{\displaystyle {\int }_0^x\mathrm{d}y} \cos (\mathsf{\pi }x)\cos (\mathsf{\pi }y) ,\end{array}$$

where the function $\theta (\mathsf{\pi }x)$ is introduced by (32), $x=\frac{m+1}{M}$. Further,

$$\begin{array}{ll}\underset{\begin{array}{l}L\to \infty \\ M\to \infty \end{array}}{\lim }{\displaystyle \sum _{m=1}^M\frac{W_m^{(\uparrow )}}{M}}& =2{\displaystyle {\int }_0^1\mathrm{d}x}{\displaystyle {\int }_x^1\mathrm{d}y}\sin \theta (\mathsf{\pi }x)\sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }x)\sin (\mathsf{\pi }y)\\ & -2{\displaystyle {\int }_0^1\mathrm{d}x}{\displaystyle {\int }_x^1\mathrm{d}y}\cos \theta (\mathsf{\pi }x)\cos \theta (\mathsf{\pi }y)\cos (\mathsf{\pi }x)\cos (\mathsf{\pi }y)+2{\displaystyle {\int }_0^1\mathrm{d}x}{\displaystyle {\int }_0^x\mathrm{d}y} \cos (\mathsf{\pi }x)\cos (\mathsf{\pi }y) .\end{array}$$

(34)

Generally speaking, for the function $f(u, v)$ and the domain $D$ of the integration shown in Figure 1, we have:

$${\displaystyle \underset{D}{\iint }\mathrm{d}u\mathrm{d}v}f(u, v)={\displaystyle {\int }_a^b \mathrm{d}u}{\displaystyle {\int }_a^{ u} \mathrm{d}v}f(u, v)={\displaystyle {\int }_a^b\mathrm{d}v}{\displaystyle {\int }_v^b \mathrm{d}u}f(u, v) .$$

(35)

Using (35), for the first term in (34) we obtain:

$$\begin{array}{l}2{\displaystyle {\int }_0^1\mathrm{d}x}{\displaystyle {\int }_x^1\mathrm{d}y} \sin \theta (\mathsf{\pi }x)\sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }x)\sin (\mathsf{\pi }y)=2{\displaystyle {\int }_0^1\mathrm{d}x}\sin \theta (\mathsf{\pi }x)\sin (\mathsf{\pi }x){\displaystyle {\int }_x^1\mathrm{d}y} \sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }y)\\ ={\displaystyle {\int }_0^1\mathrm{d}x} \sin \theta (\mathsf{\pi }x)\sin (\mathsf{\pi }x){\displaystyle {\int }_x^1\mathrm{d}y} \sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }y)+{\displaystyle {\int }_0^1\mathrm{d}y} \sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }y){\displaystyle {\int }_0^y\mathrm{d}x} \sin \theta (\mathsf{\pi }x)\sin (\mathsf{\pi }x)\\ ={\displaystyle {\int }_0^1\mathrm{d}x} \sin \theta (\mathsf{\pi }x)\sin (\mathsf{\pi }x){\displaystyle {\int }_0^1\mathrm{d}y} \sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }y) .\end{array}$$

Likewise, the second term in (34) becomes:

$$-2{\displaystyle {\int }_0^1\mathrm{d}x}{\displaystyle {\int }_x^1\mathrm{d}y} \cos \theta (\mathsf{\pi }x)\cos \theta (\mathsf{\pi }y)\cos (\mathsf{\pi }x)\cos (\mathsf{\pi }y)=-{\displaystyle {\int }_0^1\mathrm{d}x} \cos \theta (\mathsf{\pi }x)\cos (\mathsf{\pi }x){\displaystyle {\int }_0^1\mathrm{d}y} \cos \theta (\mathsf{\pi }y)\cos (\mathsf{\pi }y) .$$

Further, the third term in (34) vanishes due to:

$${\displaystyle {\int }_0^1\mathrm{d}x}{\displaystyle {\int }_0^x\mathrm{d}y} \cos (\mathsf{\pi }x)\cos (\mathsf{\pi }y)=\frac{1}{\mathsf{\pi }}{\displaystyle {\int }_0^1\mathrm{d}x}\cos (\mathsf{\pi }x)\sin (\mathsf{\pi }x)=0$$

Therefore, (34) becomes:

$$\begin{array}{l}\underset{\begin{array}{l}L\to \infty \\ M\to \infty \end{array}}{\lim }{\displaystyle \sum _{m=1}^M\frac{W_m^{(\uparrow )}}{M}}={\displaystyle {\int }_0^1\mathrm{d}x} \sin \theta (\mathsf{\pi }x)\sin (\mathsf{\pi }x){\displaystyle {\int }_0^1\mathrm{d}y} \sin \theta (\mathsf{\pi }y)\sin (\mathsf{\pi }y)\\ -{\displaystyle {\int }_0^1\mathrm{d}x} \cos \theta (\mathsf{\pi }x)\cos (\mathsf{\pi }x){\displaystyle {\int }_0^1\mathrm{d}y} \cos \theta (\mathsf{\pi }y)\cos (\mathsf{\pi }y)\\ =-{\displaystyle {\int }_0^1\mathrm{d}x} \cos (\theta (\mathsf{\pi }x)-\mathsf{\pi }x){\displaystyle {\int }_0^1\mathrm{d}y} \cos (\theta (\mathsf{\pi }y)+\mathsf{\pi }y) .\end{array}$$

(36)

Substituting (31) and (36) into (33), we obtain the form of $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle$ in the thermodynamic limit:

$$\underset{\begin{array}{l}L\to \infty \\ N\to \infty \end{array}}{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle =-{\displaystyle {\int }_0^1\mathrm{d}x} \cos (\theta (\mathsf{\pi }x)-\mathsf{\pi }x){\displaystyle {\int }_0^1\mathrm{d}y} \cos (\theta (\mathsf{\pi }y)+\mathsf{\pi }y)+{\left({\displaystyle {\int }_0^1\mathrm{d}x}\cos \theta (\mathsf{\pi }x)\right)}^2 .$$

(37)

On the other hand, the expressions of $\langle {\sigma }_{1, 1}{\sigma }_{1, 1+Q}\rangle$ in the thermodynamic limit have been obtained [3,5,9]. Thus, we here cite the formulas (B6) and (B7) in Reference [9] for comparison. According to those two formulas:

$$\underset{\begin{array}{l}L\to \infty \\ N\to \infty \end{array}}{\lim }\langle {\sigma }_{1 , 1}{\sigma }_{1 , 2}\rangle =a_0,\underset{\begin{array}{l}L\to \infty \\ N\to \infty \end{array}}{\lim }\langle {\sigma }_{1 , 1}{\sigma }_{1 , 3}\rangle =\left|\begin{array}{cc}{a}_0& a_1\\ a_{-1}& a_0\end{array}\right|=a_0^2-a_1{a}_{-1};a_r=\frac{1}{\mathsf{\pi }}{\displaystyle {\int }_0^{\mathsf{\pi }} \mathrm{d}\omega }\cos (\theta (\omega )-r\omega ).$$

(38)

where $\theta (\omega )$ is the function $\delta \ast (\omega )$ in Reference [9]. We see that (31) and (37) obtained here are exactly the same as (38).

3. Long Range-Order in the Model with Periodic-Free Boundary Conditions

For the model with $L$ rows and $N$ columns and periodic boundary condition in the horizontal direction and free boundary condition in the vertical direction, we consider $\langle {\sigma }_{l , n}{\sigma }_{l , n^{\prime }}\rangle$, i.e., correlation functions of pairwise spins in one column, periodic boundary condition in the horizontal direction leads to:

$$\langle {\sigma }_{l , n}{\sigma }_{l , n^{\prime }}\rangle =\langle {\sigma }_{1 , n}{\sigma }_{1 , n^{\prime }}\rangle =\frac{1}{Z_0}{\displaystyle \sum _{\left\{{\sigma }_{h , v}=\pm 1\right\}}{\sigma }_{1 , 1}{\sigma }_{1 , n^{\prime }}\exp \left(\frac{J^{\prime }}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^N{\sigma }_{l , m}{\sigma }_{l+1 , m}}}\right)\exp \left(\frac{J}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^{N-1}{\sigma }_{l , m}{\sigma }_{l , m+1}}}\right)} ,$$

(39)

where

$$Z_0={\displaystyle \sum _{\left\{{\sigma }_{h , v}=\pm 1\right\}}\exp \left(\frac{J^{\prime }}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^N{\sigma }_{l , m}{\sigma }_{l+1 , m}}}\right)\exp \left(\frac{J}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^{N-1}{\sigma }_{l , m}{\sigma }_{l , m+1}}}\right)}$$

(40)

is the partition function of the system in absence of a magnetic field, where ${\sigma }_{L+1 , m}={\sigma }_{1 , m}$, $J^{\prime }(>0)$ and $J(>0)$ are the interaction constants for the horizontal and vertical directions, respectively.

3.1. Some Results Concerning the Partition Function $Z_0$

We summarize some results concerning $Z_0$ given by (40), some of which are obtained in Reference [12]. However, the approximate values of some quantities presented here show improvement over those given by Reference [12].

By using the spinor analysis method, $Z_0$ is obtained in Reference [12]:

$$Z_0={\left(2\sinh \frac{2J^{\prime }}{kT}\right)}^{\frac{LN}{2}}{\displaystyle \prod _{n=1}^N\left(2\cosh \frac{L{\gamma }_{n-1}}{2}\right)} ,$$

(41)

where ${\gamma }_{n-1}$ $(n=1, 2, \cdots , N)$ are determined by:

$$\cosh {\gamma }_{n-1}=\cosh 2K^{\prime }\cosh 2K-x_{n-1}\sinh 2K^{\prime }\sinh 2K,{\mathrm{e}}^{-2K^{\prime }}=\tanh \frac{J^{\prime }}{kT} , K=\frac{J}{kT},$$

(42)

where $(n=2, 3, \cdots , N)$ $(n=1, 2, \cdots , N)$ are N roots of the N-th order algebraic equation in x:

$$g_N(x)-2g_{N-1}(x)\coth 2K^{\prime }\tanh K+g_{N-2}(x){\tanh }^{2}K=0,$$

(43)

where

$$g_n(x)={\displaystyle \sum _{k=0}^{\left[n/2\right]}\frac{(n+1)!}{(2k+1)!(n-2k)!}{x}^{n-2k}{(x^2-1)}^{ k}}$$

(44)

is an n-th degree polynomial in $x$. If by $x\equiv \frac{d+d^{-1}}{2}$ we introduce the quantity $d$, then $g_n(x)$ can be written in the form:

$$g_n(x)=\frac{d^{n+1}-d^{-(n+1)}}{d-d^{-1}}.$$

(45)

The expression in (45) is not only simple but also convenient for investigating the properties of $g_n(x)$, especially if we assume $x=\cos \phi$, then $g_n(x)=\frac{\sin (n+1)\phi }{\sin \phi }$. Substituting these forms of $g_n(x)$ into (43), for the $N-1$ roots of the N-th order algebraic Equation (43) we obtain:

$$x_{n-1}=\cos {\phi }_{n-1},{\phi }_{n-1}=\frac{(n-1)\mathsf{\pi }+{\theta }_{n-1}}{N},0<{\theta }_{n-1}<\mathsf{\pi }(n=2, 3, \cdots , N),$$

(46)

Further, ${\theta }_{n-1}$ can be determined by solving a transcendental equation about $\theta$; if N is finite, then the evaluation of ${\theta }_{n-1}$ is complex because of the so-called “finite size effect”; for the limit case $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$, we can assume ${\theta }_{n-1}={\displaystyle \sum _{k=0}^{\infty }\frac{{\theta }_{n-1}^{(k)}}{N^k}}$ and obtain ${\theta }_{n-1}$ by the iterative method. Further, we obtain ${\gamma }_1 , {\gamma }_2 , \cdots , {\gamma }_{N-1}$ in terms of (42). Concretely, correcting to $\frac{1}{N}$ order, we have:

$$\begin{array}{c}{x}_{n-1}=\cos \frac{(n-1)\mathsf{\pi }+{\theta }_{n-1}^{(0)}}{N} , {\gamma }_{n-1}\approx {\gamma }_{n-1}^{(0)}+2\sin \frac{(n-1)\mathsf{\pi }+{\theta }_{n-1}^{(0)}}{N}\sin \frac{{\theta }_{n-1}^{(0)}}{2N}\frac{\sinh 2K^{\prime }\sinh 2K}{\sinh {\gamma }_{n-1}^{(0)}} ,\\ (n=2, 3, \cdots , N)\end{array}$$

(47)

where ${\gamma }_{n-1}^{(0)}$ and ${\theta }_{n-1}^{(0)}$ $(n=2, 3, \cdots , N)$ are introduced by:

$$\begin{array}{l}\cosh {\gamma }_{n-1}^{(0)}=\cosh 2K^{\prime }\cosh 2K-\cos \frac{(n-1)\mathsf{\pi }}{N}\sinh 2K^{\prime }\sinh 2K , {\gamma }_{n-1}^{(0)}>0;\\ \cosh 2K^{\prime }\sinh 2K-\cos \frac{(n-1)\mathsf{\pi }}{N}\sinh 2K^{\prime }\cosh 2K=\sinh {\gamma }_{n-1}^{(0)}\cos {\theta }_{n-1}^{(0)} ,\\ \sin \frac{(n-1)\mathsf{\pi }}{N}\sinh 2K^{\prime }=\sinh {\gamma }_{n-1}^{(0)}\sin {\theta }_k^{(0)} , 0<{\theta }_{n-1}^{(0)}<\mathsf{\pi } .\end{array}$$

(48)

More important are the values of $x_0$ and ${\gamma }_0$; to present $x_0$ and ${\gamma }_0$, we first introduce a temperature $T_c$ in terms of:

$$N=\frac{\sinh 2{K^{\prime }}_{\mathrm{c}}}{\sinh 2(K_{\mathrm{c}}-{K^{\prime }}_{\mathrm{c}})}, \tanh {K^{\prime }}_{\mathrm{c}}={\mathrm{e}}^{-\frac{2J^{\prime }}{k{T}_{\mathrm{c}}}}, K_{\mathrm{c}}=\frac{J}{k{T}_{\mathrm{c}}}.$$

(49)

When $T\ge T_{\mathrm{c}}$, $0<x_0\le 1$; however, when $T<T_{\mathrm{c}}$, $K^{\prime }<K$ and $1<x_0<\frac{\tanh 2K}{\tanh 2K^{\prime }}$. For the limit case $N\to \infty$, we can obtain the approximate values of $x_0$ and ${\gamma }_0$, whose low-order approximations are:

$$x_0\approx \{\begin{array}{ll}\cos \frac{\mathsf{\pi }}{N},& T\ge T_{\mathrm{c}}, K^{\prime }>K ;\\ \cos \frac{\mathsf{\pi }}{2N},& T\ge T_{\mathrm{c}}, K^{\prime }=K ;\\ \frac{1}{2}\left(\frac{\tanh K}{\tanh K^{\prime }}+\frac{\tanh K^{\prime }}{\tanh K}\right)-2{\left(\frac{\tanh K^{\prime }}{\tanh K}\right)}^{2N}\frac{{(\cosh 2K-\cosh 2K^{\prime })}^2}{\sin \mathrm{h}2K^{\prime }{\sinh }^{3}2K},& T<T_{\mathrm{c}},\end{array}$$

(50)

$${\gamma }_0\approx \{\begin{array}{ll}2(K^{\prime }-K)+2{\sin }^2\frac{\mathsf{\pi }}{2N}\frac{\sinh 2K^{\prime }\sinh 2K}{\sinh 2(K^{\prime }-K)},& T\ge T_{\mathrm{c}}, K^{\prime }>K ;\\ 2\sin \frac{\mathsf{\pi }}{4N}\sinh 2K,& T\ge T_{\mathrm{c}}, K^{\prime }=K ;\\ 2{\left(\frac{\tanh K^{\prime }}{\tanh K}\right)}^N\frac{\cosh 2K-\cosh 2K^{\prime }}{\sinh 2K},& T<T_{\mathrm{c}} .\end{array}$$

(51)

We can thus make a comparison between Onsager’s lattice and the model with periodic-free boundary conditions. For Onsager’s lattice, when the system crosses its critical temperature, ${\gamma }_N^{(-)}={\gamma }_{2M}^{(-)}=2(K-K^{\prime })$ given by (7) changes sign; however, from (51) we see that, for the model with periodic-free boundary conditions, when $T\ge T_{\mathrm{c}}$, ${\gamma }_0\approx 2(K^{\prime }-K)$. Once the system crosses its critical temperature $T_{\mathrm{c}}$, ${\gamma }_0$ becomes exponentially smaller and then vanishes rapidly as $N\to \infty$. This property of ${\gamma }_0$ plays a key role for the correlation function $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$.

3.2. The Matrix Forms of $\langle {\sigma }_{1, 1}{\sigma }_{1, N-Q}\rangle$ and Some Results Concerning $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ Obtained in Reference [12]

If we write (39) in forms similar to (8), then by employing the method presented in Section 2, the exact expressions we can obtain are still $\langle {\sigma }_{1, 1}{\sigma }_{1, 2}\rangle ,$ $\langle {\sigma }_{1, 1}{\sigma }_{1, 3}\rangle , \cdots$, which belong to the short-range order. We still cannot obtain the exact expressions of $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$, $\langle {\sigma }_{l, 1}{\sigma }_{l, N-1}\rangle , \cdots$.

To obtain the exact expressions of $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$, $\langle {\sigma }_{l, 1}{\sigma }_{l, N-1}\rangle , \cdots$, we consider the forms:

$$\langle {\sigma }_{1, 1}{\sigma }_{1, N-n}\rangle =\frac{1}{Z_0}{\displaystyle \sum _{\left\{{\sigma }_{h , v}=\pm 1\right\}}{\sigma }_{1, 1}{\sigma }_{1, N-n}\exp \left(\frac{J^{\prime }}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^N{\sigma }_{l , m}{\sigma }_{l+1 , m}}}\right)\exp \left(\frac{J}{kT}{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^{N-1}{\sigma }_{l , m}{\sigma }_{l , m+1}}}\right)}.$$

(52)

Taking advantage of ${\sigma }_{l , m}^2=1$ and ${\sigma }_{1 , l}{\sigma }_{1 , m}=\underset{\alpha \to 0}{\lim }\frac{\partial {\mathrm{e}}^{\alpha {\sigma }_{1 , l}{\sigma }_{1 , m}}}{\partial \alpha }$, we have:

$$\begin{array}{l}{\sigma }_{1, 1}{\sigma }_{1, N-n}={\sigma }_{1, 1}{\sigma }_{1, N}{\sigma }_{1, N}{\sigma }_{1, N-1}\cdots {\sigma }_{1, N-(n-2)}{\sigma }_{1, N-(n-1)}{\sigma }_{1, N-(n-1)}{\sigma }_{1, N-n}\\ =\left({\displaystyle \prod _{k=1}^n\underset{{\beta }_k\to 0}{\lim }\frac{\partial }{\partial {\beta }_k}}\right)\underset{{\beta }_N\to 0}{\lim }\frac{\partial }{\partial {\beta }_N}\left({\mathrm{e}}^{{\beta }_N{\sigma }_{1 , 1}{\sigma }_{1 , N}}{\displaystyle \prod _{k=0}^{n-1}{\mathrm{e}}^{{\beta }_{k+1}{\sigma }_{1 , N-k}{\sigma }_{1 , N-(k+1)}}}\right) ,\end{array}$$

Further, (52) can be written in the form:

$$\langle {\sigma }_{1, 1}{\sigma }_{1, N-n}\rangle =\frac{1}{Z_0}{\left(2\sinh \frac{2J^{\prime }}{kT}\right)}^{\frac{LN}{2}}\left({\displaystyle \prod _{k=1}^n\underset{{\beta }_k\to 0}{\lim }\frac{\partial }{\partial {\beta }_k}}\right)\underset{{\beta }_N\to 0}{\lim }\frac{\partial }{\partial {\beta }_N}\mathrm{Tr}(\mathbf{W}),$$

(53)

where the matrix W belongs to the spin representative, and, by employing the method presented in Section 2, we can obtain the exact expressions of $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$, $\langle {\sigma }_{l, 1}{\sigma }_{l, N-1}\rangle$, $\langle {\sigma }_{l, 1}{\sigma }_{l, N-2}\rangle , \cdots$.

However, to save space, here we no longer discuss $\langle {\sigma }_{l, 1}{\sigma }_{l, N-1}\rangle$, $\langle {\sigma }_{l, 1}{\sigma }_{l, N-2}\rangle , \cdots$, but only consider $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$, for which a closed formula was given by Reference [12]:

$$\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle =\frac{1}{Z_0}{\left(2\sinh \frac{2J^{\prime }}{kT}\right)}^{\frac{LN}{2}}\underset{\varphi \to 0}{\lim }\frac{\partial Y}{\partial \varphi },$$

(54)

$$\begin{array}{ll}Y=& \frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\cosh \frac{{\chi }_{2(l-1)}^{(+)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\cosh \frac{-{\chi }_{2m-1}^{(+)}}{2}}\right)+\frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\sinh \frac{{\chi }_{2(l-1)}^{(+)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\sinh \frac{-{\chi }_{2m-1}^{(+)}}{2}}\right)\\ & +\frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\cosh \frac{{\chi }_{2(l-1)}^{(-)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\cosh \frac{-{\chi }_{2m-1}^{(-)}}{2}}\right)-\frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\sinh \frac{{\chi }_{2(l-1)}^{(-)}}{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\sinh \frac{-{\chi }_{2m-1}^{(-)}}{2}}\right) ,\end{array}$$

(55)

where ${\chi }_{2(l-1)}^{(\pm )}$ and ${\chi }_{2m-1}^{(\pm )}$ are determined by:

$${\mathrm{e}}^{{\chi }_{2(l-1)}^{(\pm )}}={\tau }_{2(l-1)}^{(\pm )}\left(l=1, 2, \cdots , \left[\frac{N+1}{2}\right]\right){; \mathrm{e}}^{{\chi }_{2m-1}^{(\pm )}}={\tau }_{2m-1}^{(\pm )}\left(m=1, 2, \cdots , \left[\frac{N}{2}\right]\right) ,$$

(56)

${\tau }_n^{(\pm )}$ $(n=1, 2, \cdots , N)$ are N roots of the N-th order algebraic equation $F_{\pm }(\tau )=0$, where

$$\begin{array}{ll}{F}_{\pm }\left(\tau \right)=& \left({\displaystyle \prod _{l=1}^{[(N+1)/2]}\left(\tau { \mathrm{e}}^{-L{\gamma }_{2(l-1)}}-1\right)}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}\left(\tau { \mathrm{e}}^{L{\gamma }_{2m-1}}-1\right)}\right)\\ & \mp \tanh \varphi \left({\displaystyle \prod _{l=1}^{[(N+1)/2]}\left(\tau { \mathrm{e}}^{-L{\gamma }_{2(l-1)}}-1\right)}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}\left(\tau { \mathrm{e}}^{L{\gamma }_{2m-1}}-1\right)}\right)\\ & \mp 4\tanh \varphi {\displaystyle \sum _{n=1}^{[(N+1)/2]}\left({\Omega }_{2(n-1)}^2\left({\displaystyle \prod _{\begin{array}{l}l=1\\ l\ne n\end{array}}^{[(N+1)/2]}\left(\tau { \mathrm{e}}^{-L{\gamma }_{2(l-1)}}-1\right)}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}\left(\tau { \mathrm{e}}^{L{\gamma }_{2m-1}}-1\right)}\right)\right)}\\ & \mp 4\tanh \varphi {\displaystyle \sum _{n=1}^{[N/2]}\left({\Omega }_{2n-1}^2\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}\left(\tau { \mathrm{e}}^{-L{\gamma }_{2(l-1)}}-1\right)}\right)\left({\displaystyle \prod _{\begin{array}{l}m=1\\ m\ne n\end{array}}^{[N/2]}\left(\tau { \mathrm{e}}^{L{\gamma }_{2m-1}}-1\right)}\right)\right)} ,\end{array}$$

(57)

where

$${\Omega }_{n-1}=\sinh 2K^{\prime }\cosh K\sqrt{\frac{1-x_{n-1}^2}{N{\sinh }^2{\gamma }_{n-1}+\cosh {\gamma }_{r-1}\cosh 2K^{\prime }-\cosh 2K}}\hspace{1em}(n=1, 2, \cdots , N)$$

(58)

are the normalization constants of the eigenvectors of a rotation matrix [12], as $N\to \infty$. Thus, according to (47)~(51), we obtain:

$$\begin{array}{l}\underset{N\to \infty }{\lim }{\Omega }_{n-1}\approx \frac{1}{\sqrt{N}}\sin \frac{(n-1)\mathsf{\pi }}{N}\frac{\sinh 2K^{\prime }\cosh K}{\sinh {\gamma }_{n-1}}\hspace{1em}(n=2, 3, \cdots , N) ,\\ \underset{N\to \infty }{\lim }{\Omega }_0^2\approx \{\begin{array}{l}\frac{1}{N}{\sin }^2\frac{\mathsf{\pi }}{N}\frac{{\sinh }^{2}2K^{\prime }{\cosh }^{2}K}{{\sinh }^{2}2(K^{\prime }-K)} ,T\ge T_{\mathrm{c}} , K^{\prime }>K ;\\ \frac{1}{N}\frac{{\sinh }^{2}2K^{\prime }{\cosh }^{2}K}{{\sinh }^{2}2K} ,T\ge T_{\mathrm{c}} , K^{\prime }=K ;\\ \frac{\cosh 2K-\cosh 2K^{\prime }}{4{\sinh }^{2}K} .T<T_{\mathrm{c}} .\end{array}\end{array}$$

(59)

3.3. An Exact Expression of $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ on a Finite Lattice

In Reference [12], all roots of the equation $F_{\pm }(\tau )=0$ are obtained by correcting to the ${\mathrm{e}}^{-L{C}_0}$ order of magnitude ($C_0$ is a positive constant). These approximate roots can lead to the exact expression of $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$ in the thermodynamic limit, since $\underset{L\to \infty }{\lim }{\mathrm{e}}^{-L{C}_0}=0$, but cannot lead to the exact expression of $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$ on a finite lattice. Hence, the expression of $\langle {\sigma }_{l, 1}{\sigma }_{l, N}\rangle$ presented in Reference [12] is only an approximate result.

On the other hand, similar to the analysis in Section 2.5, the operator $\underset{\varphi \to 0}{\lim }\frac{\partial }{\partial \varphi }$ in (54) allows us to ignore all terms whose order is higher than ${\varphi }^1(=\varphi )$ in all roots of the equation $F_{\pm }(\tau )=0$. Hence, to obtain the exact expression of $\underset{\varphi \to 0}{\lim }\frac{\partial Y}{\partial \varphi }$, we need only to find all roots of the equation $F_{\pm }(\tau )=0$ corrected to the ${\varphi }^1(=\varphi )$ order, The corresponding calculations are in fact simpler than those required of find the roots corrected to the ${\mathrm{e}}^{-L{C}_0}$ order of magnitude in Reference [12]; concretely, we obtain:

$${\chi }_{2(l-1)}^{(\pm )}\approx L{\gamma }_{2(l-1)}\pm 4{\Omega }_{2(l-1)}^2\tanh \varphi , {\chi }_{2m-1}^{(\pm )}\approx -L{\gamma }_{2m-1}\pm 4{\Omega }_{2m-1}^2\tanh \varphi .$$

Substituting the above results into (55), we obtain $Y$ correcting to $\tanh \varphi (\approx {\varphi }^1)$ order:

$$\begin{array}{ll}Y\approx & \frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\cosh \frac{L{\gamma }_{2(l-1)}+4{\Omega }_{2(l-1)}^2\tanh \varphi }{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\cosh \frac{L{\gamma }_{2m-1}-4{\Omega }_{2m-1}^2\tanh \varphi }{2}}\right)\\ & +\frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\sinh \frac{L{\gamma }_{2(l-1)}+4{\Omega }_{2(l-1)}^2\tanh \varphi }{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\sinh \frac{L{\gamma }_{2m-1}-4{\Omega }_{2m-1}^2\tanh \varphi }{2}}\right)\\ & +\frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\cosh \frac{L{\gamma }_{2(l-1)}-4{\Omega }_{2(l-1)}^2\tanh \varphi }{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\cosh \frac{L{\gamma }_{2m-1}+4{\Omega }_{2m-1}^2\tanh \varphi }{2}}\right)\\ & -\frac{1}{2}\left({\displaystyle \prod _{l=1}^{[(N+1)/2]}2\sinh \frac{L{\gamma }_{2(l-1)}-4{\Omega }_{2(l-1)}^2\tanh \varphi }{2}}\right)\left({\displaystyle \prod _{m=1}^{[N/2]}2\sinh \frac{L{\gamma }_{2m-1}+4{\Omega }_{2m-1}^2\tanh \varphi }{2}}\right).\end{array}$$

Substituting the above result and (41) into (54), we obtain the exact expression of $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ of the model on a finite lattice:

$$\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle =2\left({\displaystyle \sum _{l=1}^{[(N+1)/2]}{\Omega }_{2(l-1)}^2\coth \frac{L{\gamma }_{2(l-1)}}{2}}-{\displaystyle \sum _{m=1}^{[N/2]}{\Omega }_{2m-1}^2\coth \frac{L{\gamma }_{2m-1}}{2}}\right)\left({\displaystyle \prod _{n=1}^N\tanh \frac{L{\gamma }_{n-1}}{2}}\right)$$

(60)

Although the whole calculation process is complex, the final result (60) is simple.

3.4. The Expression of $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ in the Thermodynamic Limit

To derive the expression of $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ in the thermodynamic limit, we first discuss some properties of ${\gamma }_{n-1}$ $(n=1, 2, \cdots , N)$.

For ${\gamma }_{n-1}$ $(n=2, 3, \cdots , N)$ given by (47), we have:

$$\underset{L\to \infty }{\lim }L{\gamma }_{n-1}=\infty , \underset{L\to \infty }{\lim }\tanh \frac{L{\gamma }_{n-1}}{2}=\underset{L\to \infty }{\lim }\coth \frac{L{\gamma }_{n-1}}{2}=1\hspace{1em}(n=2, 3, \cdots , N).$$

(61)

As for ${\gamma }_0$, when $T\ge T_{\mathrm{c}}$, from (51) we see that (61) still holds for ${\gamma }_0$; however, when $T<T_{\mathrm{c}}$, from the last expression in (51) we see that maybe $\underset{L\to \infty }{\lim }L{\gamma }_0=\infty$ does not hold. For example, if $L=N^a$, where a is a positive integer, then:

$$\underset{L\to \infty }{\lim }L{\gamma }_0=\underset{N\to \infty }{\lim }{N}^a\cdot 2{\left(\frac{\tanh K^{\prime }}{\tanh K}\right)}^N\csc \mathrm{h}2K(\cosh 2K-\cosh 2K^{\prime })=0,$$

since now $0<\frac{\tanh K^{\prime }}{\tanh K}<1$, for $0<b<1$, $\underset{N\to \infty }{\lim }{N}^a{b}^N=0$. Hence, for ${\gamma }_0$ we have:

$$\underset{L\to \infty , N\to \infty }{\lim } L{\gamma }_0=\{\begin{array}{l}\infty , T\ge T_{\mathrm{c}} ;\\ 0 , T<T_{\mathrm{c}} ,\end{array}\hspace{1em}\underset{L\to \infty , N\to \infty }{\lim } \tanh \frac{L{\gamma }_0}{2}=\{\begin{array}{l}1 , T\ge T_{\mathrm{c}} ;\\ 0 , T<T_{\mathrm{c}} ,\end{array}\hspace{1em}(L=N^a)$$

(62)

According to the above discussions, to obtain the expression of $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ in the thermodynamic limit, we first write (60) in the form:

$$\begin{array}{ll}\langle {\sigma }_{1 , 1}{\sigma }_{1 , N}\rangle =& 2{\Omega }_0^2\left({\displaystyle \prod _{n=2}^N\tanh \frac{L{\gamma }_{n-1}}{2}}\right)\\ & +2\tanh \frac{L{\gamma }_0}{2}\left({\displaystyle \sum _{l=2}^{[(N+1)/2]}{\Omega }_{2(l-1)}^2\coth \frac{L{\gamma }_{2(l-1)}}{2}}-{\displaystyle \sum _{m=1}^{[N/2]}{\Omega }_{2m-1}^2\coth \frac{L{\gamma }_{2m-1}}{2}}\right)\left({\displaystyle \prod _{n=2}^N\tanh \frac{L{\gamma }_{n-1}}{2}}\right) ,\end{array}$$

as $L\to \infty$, according to (61), the above expression becomes:

$$\underset{L\to \infty }{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle \approx 2{\Omega }_0^2+2\tanh \frac{L{\gamma }_0}{2}\left({\displaystyle \sum _{l=2}^{[(N+1)/2]}{\Omega }_{2(l-1)}^2}-{\displaystyle \sum _{m=1}^{[N/2]}{\Omega }_{2m-1}^2}\right) .$$

(63)

When $T\ge T_{\mathrm{c}}$, according to (59) and (62), Equation (63) becomes:

$$\underset{L\to \infty }{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle \approx 2\left({\displaystyle \sum _{l=1}^{[(N+1)/2]}{\Omega }_{2(l-1)}^2}-{\displaystyle \sum _{m=1}^{[N/2]}{\Omega }_{2m-1}^2}\right) ,$$

Further, as $N\to \infty$,

$${\displaystyle \sum _{l=1}^{[(N+1)/2]}{\Omega }_{2(l-1)}^2}\approx {\displaystyle \sum _{m=1}^{[N/2]}{\Omega }_{2m-1}^2}\approx \frac{1}{2}{\displaystyle \sum _{n=1}^N{\Omega }_{n-1}^2},$$

(64)

hence, for this case $\underset{L\to \infty , N\to \infty }{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle =0$.

When $T<T_{\mathrm{c}}$, according to (62) and (64), the second term in (63) vanishes, and (63) thus becomes $\underset{L\to \infty }{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle \approx 2{\Omega }_0^2$, where ${\Omega }_0^2$ is given by the last expression in (59) as $N\to \infty$.

Summarizing the above results, in the thermodynamic limit, if $L=N^a$, then (60) becomes:

$$\underset{L\to \infty , N\to \infty }{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle = \{\begin{array}{ll}0,& T\ge T_{\mathrm{c}} ;\\ \frac{\cosh 2K-\cosh 2K^{\prime }}{2{\sinh }^{2}K},& T<T_{\mathrm{c}} .\end{array}$$

The above result was obtained in Reference [12] in terms of an approximate result of $\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$. Some further discussions about the above result can be found in Reference [12].

From the above discussions, it is revealed how the changes of the values of ${\gamma }_{n-1}$ $(n=1, 2, \cdots , N)$, especially the change of the value of ${\gamma }_0$, lead to the change of $\underset{L\to \infty , N\to \infty }{\lim }\langle {\sigma }_{1, 1}{\sigma }_{1, N}\rangle$ when the system crosses its critical temperature $T_{\mathrm{c}}$, as well as how the long-range order emerges as the temperature decreases.