Influence of long-range interaction on degeneracy of eigenvalues of connection matrix of d-dimensional Ising system

Let us consider a chain of the length L and set the distance between its nodes to be equal to one. For certainty, we suppose that L is an odd number: L = 2l + 1. Let J(k) be an L × L symmetric matrix that describes the interactions only between spins spaced by the distance k (k = 1, 2, ..., l). The structure of the matrix J(k) is as follows. The ones occupy the kth and (L − k)th its diagonals which are parallel to the main diagonal and the other matrix elements are equal to zero. The central row of the matrix J(k) has the form

$$\begin{equation*}\left(0\quad \dots \quad 1\quad \dots \quad 0\quad \dots \quad 1\quad \dots \quad 0\right),\end{equation*}$$

where the ones are at the distance k from the center. We obtain all other rows by consequent cyclic shifts of the (l + 1)th row: shifting it to the left we obtain the lth row, the right shift gives the (l + 2)-row, and so on. For example, when L = 7, the matrices J(1), J(2), and J(3) are

$$\begin{align}\hfill \mathbf{J}\left(1\right)=& \left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{matrix}\right),\qquad \mathbf{J}\left(2\right)=\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right),\\ \mathbf{J}\left(3\right)=& \left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right).\end{align}$$

All the matrices J(k) commute and consequently they all have the same set of the eigenvectors ${\left\{{\mathbf{f}}_{\alpha }\right\}}_{\alpha =1}^{L}$. The components of the vectors f_α are well-known [21]:

$$\begin{equation}{f}_{\alpha }^{\left(i\right)}=\frac{\mathrm{cos}\enspace {\varphi }_{\alpha i}+\mathrm{sin}\enspace {\varphi }_{\alpha i}}{\sqrt{L}},\qquad {\varphi }_{\alpha i}=\frac{2\pi \left(\alpha -1\right)\left(i-1\right)}{L},\quad i=1,2,\dots ,L.\end{equation} \tag{ 1 }$$

Each matrix J(k) has its own set of the eigenvalues ${\left\{{\lambda }_{\alpha }\left(k\right)\right\}}_{\alpha =1}^{L}$:

$$\begin{equation*}\mathbf{J}\left(k\right)\enspace {\mathbf{f}}_{\alpha }={\lambda }_{\alpha }\left(k\right){\mathbf{f}}_{\alpha },\qquad {\lambda }_{\alpha }\left(k\right)=2\enspace \mathrm{cos}\left(\frac{2\pi k\left(\alpha -1\right)}{L}\right),\quad \alpha =1,2,\dots ,L,\quad k=1,2,\dots ,l.\end{equation*}$$

Let w(k) be the constant of interaction between spins that are at the distant k from each other, where k = 1, 2, ..., l. Then for the one-dimensional lattice, the interaction matrix taking account for an arbitrary long-range interaction has the form:

$$\begin{equation*}{\mathbf{A}}_{0}=\sum _{k=1}^{l}w\enspace \left(k\right)\enspace \mathbf{J}\left(k\right).\end{equation*}$$

The equation (1) define the eigenvectors of this matrix and its eigenvalues, A₀ f_α = μ_α f_α , are

$$\begin{equation}{\mu }_{\alpha }=\sum _{k=1}^{l}w\enspace \left(k\right){\lambda }_{\alpha }\left(k\right),\quad \alpha =1,2,\dots ,L.\end{equation} \tag{ 2 }$$

For simplicity and universality, we introduce the notations,

$$\begin{equation*}\mathbf{J}\left(0\right)=\mathbf{I},\quad {\lambda }_{\alpha }\left(0\right)=1,\quad \alpha =1,2,\dots ,L,\end{equation*}$$

where I is an L × L unit matrix. In addition, we would like to recall the product rules for matrices and vectors that are the Kronecker products. Suppose we have a matrix M = M₁ ⊗ M₂ and a vector F = F₁ ⊗ F₂. Then

$$\begin{equation}{\mathbf{F}}^{+}\mathbf{M}\mathbf{F}\equiv {\mathbf{F}}_{1}^{+}\otimes {\mathbf{F}}_{2}^{+}\enspace {\mathbf{M}}_{1}\otimes {\mathbf{M}}_{2}\enspace {\mathbf{F}}_{1}\otimes {\mathbf{F}}_{2}={\mathbf{F}}_{1}^{+}{\mathbf{M}}_{1}{\mathbf{F}}_{1}\cdot {\mathbf{F}}_{2}^{+}{\mathbf{M}}_{2}{\mathbf{F}}_{2}.\end{equation} \tag{ 3 }$$

For the following calculations, the rule (3) is very useful.

In this subsection we discuss the 2D Ising model that is a system of spins in the nods of a square lattice. By w(m, k) we denote the constant of interaction between spins that are shifted from each other at a distance m along one of the lattice axis and at a distance k along the other axis. The connection matrix of such a system is an L² × L² matrix B₀. It is convenient to present this matrix as an L × L block matrix

$$\begin{equation}{\mathbf{B}}_{0}=\left(\begin{matrix}\hfill {\mathbf{A}}_{0}\hfill & \hfill {\mathbf{A}}_{1}\hfill & \hfill {\mathbf{A}}_{2}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{l\enspace \enspace \enspace }\hfill & \hfill {\mathbf{A}}_{l\enspace \enspace }\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{2}\hfill & \hfill {\mathbf{A}}_{1}\hfill \\ \hfill {\mathbf{A}}_{1}\hfill & \hfill {\mathbf{A}}_{0}\hfill & \hfill {\mathbf{A}}_{1}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{l-1}\hfill & \hfill {\mathbf{A}}_{l\enspace \enspace }\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{3}\hfill & \hfill {\mathbf{A}}_{2}\hfill \\ \hfill {\mathbf{A}}_{2}\hfill & \hfill {\mathbf{A}}_{1}\hfill & \hfill {\mathbf{A}}_{0}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{l-2}\hfill & \hfill {\mathbf{A}}_{l-3}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{4}\hfill & \hfill {\mathbf{A}}_{3}\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\mathbf{A}}_{l}\hfill & \hfill {\mathbf{A}}_{l-1}\hfill & \hfill {\mathbf{A}}_{l-2}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{0}\hfill & \hfill {\mathbf{A}}_{1}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{l-1}\hfill & \hfill {\mathbf{A}}_{l}\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\mathbf{A}}_{2}\hfill & \hfill {\mathbf{A}}_{3}\hfill & \hfill {\mathbf{A}}_{4}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{l-1}\hfill & \hfill {\mathbf{A}}_{l-2}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{0}\hfill & \hfill {\mathbf{A}}_{1}\hfill \\ \hfill {\mathbf{A}}_{1}\hfill & \hfill {\mathbf{A}}_{2}\hfill & \hfill {\mathbf{A}}_{3}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{l\enspace \enspace }\hfill & \hfill {\mathbf{A}}_{l-1}\hfill & \hfill \dots \hfill & \hfill {\mathbf{A}}_{1}\hfill & \hfill {\mathbf{A}}_{0}\hfill \end{matrix}\right),\end{equation} \tag{ 4 }$$

where L × L matrices A_m have the form:

$$\begin{equation}{\mathbf{A}}_{m}=\sum _{k=0}^{l}w\enspace \left(m,k\right)\mathbf{J}\left(k\right),\quad \mathbf{J}\left(0\right)\equiv \mathbf{I},\quad m=0,1,\dots ,l.\end{equation} \tag{ 5 }$$

In equation (5) we set w(0, 0) = 0 and, consequently, the self-interaction is equal to zero. Note, the central block row of the matrix B₀ is

$$\begin{equation*}\left(\begin{matrix}\hfill {\mathbf{A}}_{l}{\mathbf{A}}_{l-1}\dots {\mathbf{A}}_{1}{\mathbf{A}}_{0}{\mathbf{A}}_{1}\dots {\mathbf{A}}_{l-1}{\mathbf{A}}_{l}\hfill \end{matrix}\right);\end{equation*}$$

all other block rows we obtain by evident cyclic shifts.

We can treat a two-dimensional spin system as a set of interacting one-dimensional chains (for example, the horizontal ones.) Then the matrix A₀ describes the interactions between the spins of the one horizontal chain and the matrices A_m (m ≠ 0) define interactions between the spins from different chains shifted vertically by m nods.

The matrix B₀ is a block Toeplitz matrix with the matrices A₀ on the main diagonal and the matrices A_m on its mth and (L − m)th diagonals (m ≠ 0). It is easy to show that the matrix B₀ is

$$\begin{equation}{\mathbf{B}}_{0}=\sum _{m=0}^{l}\mathbf{J}\left(m\right)\otimes {\mathbf{A}}_{m}=\sum _{m=0}^{l}\sum _{k=0}^{l}w\enspace \left(m,k\right)\cdot \mathbf{J}\left(m\right)\otimes \mathbf{J}\left(k\right).\end{equation} \tag{ 6 }$$

Since the matrices A_m commute, we can write the eigenvectors of the matrix B₀ as the Kronecker products of the eigenvectors (1):

$$\begin{equation}{\mathbf{F}}_{\alpha \beta }={\mathbf{f}}_{\alpha }\otimes {\mathbf{f}}_{\beta },\quad \alpha ,\beta =1,2,\dots ,L.\end{equation} \tag{ 7 }$$

This means that we reduce the eigenvalue problem B₀ F_αβ = μ_αβ F_αβ to calculation of the value ${\mu }_{\alpha \beta }={\mathbf{F}}_{\alpha \beta }^{+}{\mathbf{B}}_{0}{\mathbf{F}}_{\alpha \beta }$. Then substituting the matrix B₀ in the form (6) and the vector F_αβ in the form (7) with the aid of the identity (3) we obtain

$$\begin{equation}{\mu }_{\alpha \beta }=\sum _{m=0}^{l}\sum _{k=0}^{l}w\left(m,k\right)\cdot {\lambda }_{\alpha }\left(m\right)\cdot {\lambda }_{\beta }\left(k\right),\end{equation} \tag{ 8 }$$

where w(0, 0) = 0 and λ_α (0) = λ_β (0) = 1. We see that the eigenvalues μ_αβ are polynomials of the second degree of the eigenvalues λ_α (k) calculated for the one-dimensional system.

As example, let us examine a special case of an isotropic interaction only with the nearest neighbors (the interaction constants are w(0, 1) = w(1, 0) = 1) and the next nearest neighbors (the interaction constant is w(1, 1)). Then the equation (8) takes the form

$$\begin{equation}{\mu }_{\alpha \beta }={\lambda }_{\alpha }\left(1\right)+{\lambda }_{\beta }\left(1\right)+w\left(1,1\right)\cdot {\lambda }_{\alpha }\left(1\right)\enspace {\lambda }_{\beta }\left(1\right).\end{equation} \tag{ 9 }$$

The equation (9) repeats the result obtained previously in [22, 23] where we discussed this special case.

Let us discuss the three-dimensional Ising system of interacting spins that are in the nodes of a cubic lattice. By w(n, m, k) we denote the constant of interaction between spins shifted relative to each other by a distance n along one axis, by a distance m along the other axis, and by a distance k along the third axis. In such a system, it is convenient to write the connection L³ × L³ matrix C₀ in the block form:

$$\begin{equation}{\mathbf{C}}_{0}=\left(\begin{matrix}\hfill {\mathbf{B}}_{0}\hfill & \hfill {\mathbf{B}}_{1}\hfill & \hfill {\mathbf{B}}_{2}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l\enspace \enspace \enspace }\hfill & \hfill {\mathbf{B}}_{l\enspace \enspace }\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{2}\hfill & \hfill {\mathbf{B}}_{1}\hfill \\ \hfill {\mathbf{B}}_{1}\hfill & \hfill {\mathbf{B}}_{0}\hfill & \hfill {\mathbf{B}}_{1}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l-1}\hfill & \hfill {\mathbf{B}}_{l\enspace \enspace }\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{3}\hfill & \hfill {\mathbf{B}}_{2}\hfill \\ \hfill {\mathbf{B}}_{2}\hfill & \hfill {\mathbf{B}}_{1}\hfill & \hfill {\mathbf{B}}_{0}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l-2}\hfill & \hfill {\mathbf{B}}_{l-3}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{4}\hfill & \hfill {\mathbf{B}}_{3}\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\mathbf{B}}_{l}\hfill & \hfill {\mathbf{B}}_{l-1}\hfill & \hfill {\mathbf{B}}_{l-2}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{0}\hfill & \hfill {\mathbf{B}}_{1}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l-1}\hfill & \hfill {\mathbf{B}}_{l}\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill {\mathbf{B}}_{3}\hfill & \hfill {\mathbf{B}}_{4}\hfill & \hfill {\mathbf{B}}_{5}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l-2}\hfill & \hfill {\mathbf{B}}_{l-3}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{1}\hfill & \hfill {\mathbf{B}}_{2}\hfill \\ \hfill {\mathbf{B}}_{2}\hfill & \hfill {\mathbf{B}}_{3}\hfill & \hfill {\mathbf{B}}_{4}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l-1}\hfill & \hfill {\mathbf{B}}_{l-2}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{0}\hfill & \hfill {\mathbf{B}}_{1}\hfill \\ \hfill {\mathbf{B}}_{1}\hfill & \hfill {\mathbf{B}}_{2}\hfill & \hfill {\mathbf{B}}_{3}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{l\enspace \enspace }\hfill & \hfill {\mathbf{B}}_{l-1}\hfill & \hfill \dots \hfill & \hfill {\mathbf{B}}_{1}\hfill & \hfill {\mathbf{B}}_{0}\hfill \end{matrix}\right),\end{equation} \tag{ 10 }$$

where B_n are L² × L² matrices (n = 0, 1, 2, ..., l). To obtain the matrices B_n we have to generate a set of L × L matrices ${\mathbf{A}}_{m}^{\left(n\right)}$,

$$\begin{equation}{\mathbf{A}}_{m}^{\left(n\right)}=\sum _{k=0}^{l}w\left(n,m,k\right)\mathbf{J}\left(k\right),\quad \mathbf{J}\left(0\right)\equiv \mathbf{I},\quad n,m=0,1,2,\dots ,l.\end{equation} \tag{ 11 }$$

Since there is no self-interaction, we set w(0, 0, 0) = 0. With the aid of the matrices (11) we generate the matrices B_n :

$$\begin{equation*}{\mathbf{B}}_{n}=\sum _{m=0}^{l}\mathbf{J}\left(m\right)\otimes {\mathbf{A}}_{m}^{\left(n\right)}=\sum _{m=0}^{l}\sum _{k=0}^{l}w\left(n,m,k\right)\cdot \mathbf{J}\left(m\right)\otimes \mathbf{J}\left(k\right).\end{equation*}$$

By analogy with the two-dimensional system, we can consider the three-dimensional lattice as a set of interacting planar lattices. Then the matrix B₀ describes the interactions of spins belonging to one (let us say, a horizontal) plane; the matrix B_n (n ≠ 0) describes the interactions between the spins from two different planes shifted with respect to each other along the vertical axis by n nodes.

As we see, the matrix C₀ has a form of a block Toeplitz matrix with the matrices B₀ at its main diagonal and the matrices B_n (n ≠ 0) at its nth and (L − n)th diagonals. Then, we can write the matrix C₀ as

$$\begin{equation*}{\mathbf{C}}_{0}=\sum _{n=0}^{l}\mathbf{J}\left(n\right)\otimes {\mathbf{B}}_{n}=\sum _{n=0}^{l}\sum _{m}^{l}\sum _{k=0}^{l}w\left(n,m,k\right)\cdot \mathbf{J}\left(n\right)\otimes \mathbf{J}\left(m\right)\otimes \mathbf{J}\left(k\right).\end{equation*}$$

Since the matrices B_n commute, we can write the eigenvectors of the matrix C₀ as the Kronecker products of the eigenvectors (1):

$$\begin{equation*}{\mathbf{F}}_{\alpha \beta \gamma }={\mathbf{f}}_{\alpha }\otimes {\mathbf{f}}_{\beta }\otimes {\mathbf{f}}_{\gamma },\quad \alpha ,\beta ,\gamma =1,2,\dots ,L.\end{equation*}$$

This means that we reduce the eigenvalues problem C₀ F_αβγ = μ_αβγ F_αβγ to calculation of the values ${\mu }_{\alpha \beta \gamma }={\mathbf{F}}_{\alpha \beta \gamma }^{+}{\mathbf{C}}_{0}{\mathbf{F}}_{\alpha \beta \gamma }$ and with account for equation (3), we obtain:

$$\begin{equation}{\mu }_{\alpha \beta \gamma }=\sum _{n=0}^{l}\sum _{m=0}^{l}\sum _{k=0}^{l}w\left(n,m,k\right)\cdot {\lambda }_{\alpha }\left(n\right)\cdot {\lambda }_{\beta }\left(m\right)\cdot {\lambda }_{\gamma }\left(k\right),\end{equation} \tag{ 12 }$$

where, as usually, λ_α (0) = λ_β (0) = λ_γ (0) = 1. We see that in the three-dimensional case the eigenvalues are the polynomials of the third degree of the eigenvalues for the one-dimensional system.

As an example, let us discuss a special case of the three-dimensional isotropic Ising system that is w(n, m, k) = w(k, n, m) = w(m, n, k). We suppose that only the interactions with the nearest neighbors, the next nearest and the third neighbors are nonzero. We set w(0, 0, 1) = w(0, 1, 0) = w(1, 0, 0) = 1, w(0, 1, 1) = w(1, 0, 1) = w(1, 1, 0) = b₁, and w(1, 1, 1) = b₂. Then from equation (12) we obtain

$$\begin{align*}\hfill {\mu }_{\alpha \beta \gamma }& ={w}_{1}\left({\lambda }_{\alpha }\left(1\right)+{\lambda }_{\beta }\left(1\right)+{\lambda }_{\gamma }\left(1\right)\right)+{w}_{2}\left({\lambda }_{\alpha }\left(1\right)\enspace {\lambda }_{\beta }\left(1\right)+{\lambda }_{\alpha }\left(1\right)\enspace {\lambda }_{\gamma }\left(1\right)\right.\hfill \\ \hfill & \quad +\left.{\lambda }_{\beta }\left(1\right)\enspace {\lambda }_{\gamma }\left(1\right)\right)+{w}_{3}\enspace {\lambda }_{\alpha }\left(1\right)\enspace {\lambda }_{\beta }\left(1\right)\enspace {\lambda }_{\gamma }\left(1\right).\hfill \end{align*}$$

The last expression coincides with result obtained previously in [22, 23], where we examined the same special case.

A generalization to the case of a d-dimensional lattice is evident. Let us introduce the constants of interaction between spins w(k₁, k₂, ..., k_d ), where k₁ is a relative distance between the spins along the axis 1, k₂ is a relative distance between the spins along the axis 2, and so on. We generate a set of L × L matrices

$$\begin{equation*}{\mathbf{A}}_{{k}_{2}}^{\left({k}_{3},\dots ,{k}_{d}\right)}=\sum _{{k}_{1}=0}^{l}w\left({k}_{1},{k}_{2},\dots ,{k}_{d}\right)\mathbf{J}\left({k}_{1}\right),\end{equation*}$$

where k_i = 0, 1, ..., l, and i = 2, 3, ‥, d. We use these matrices to generate a set of L² × L² matrices ${\mathbf{B}}_{{k}_{3}}^{\left({k}_{4},\dots ,{k}_{d}\right)}={\sum }_{{k}_{2}=0}^{l}\mathbf{J}\left({k}_{2}\right)\otimes {\mathbf{A}}_{{k}_{2}}^{\left({k}_{3},\dots ,{k}_{d}\right)}$. In the same way the matrices ${\mathbf{B}}_{{k}_{3}}^{\left({k}_{4},\dots ,{k}_{d}\right)}$ allows us to construct the matrices ${\mathbf{C}}_{{k}_{4}}^{\left({k}_{5},\dots ,{k}_{d}\right)}={\sum }_{{k}_{3}=0}^{l}\mathbf{J}\left({k}_{3}\right)\otimes \enspace {\mathbf{B}}_{{k}_{3}}^{\left({k}_{4},\dots ,{k}_{d}\right)}$ for the three-dimensional system, and so on until we obtain the matrix ${\mathbf{U}}_{0}^{\left(d\right)}$ that defines the interactions in the d-dimensional system. This is a block matrix similar to the matrices (4) and (10), where the blocks are the L^d−1 × L^d−1 matrices that describe interactions in the (d − 1)-dimensional system. Omitting the intermediate calculations, we obtain

$$\begin{equation}{\mathbf{U}}_{0}^{\left(\mathrm{d}\right)}=\sum _{{k}_{1}=0}^{l}\sum _{{k}_{2}=0}^{l}\dots \sum _{{k}_{d}=0}^{l}w\left({k}_{1},{k}_{2},\dots ,{k}_{d}\right)\cdot \mathbf{J}\left({k}_{1}\right)\otimes \mathbf{J}\left({k}_{2}\right)\otimes \dots \otimes \mathbf{J}\left({k}_{d}\right),\end{equation} \tag{ 13 }$$

where w(0, 0, ..., 0) = 0, and J(0) ≡ I.

Since in equation (13) the matrices J(k) commute, the eigenvectors of the matrix U₀ are the Kronecker products of the eigenvectors for the one-dimensional system:

$$\begin{equation*}{\mathbf{F}}_{{\alpha }_{1}{\alpha }_{2}\dots {\alpha }_{d}}={\mathbf{f}}_{{\alpha }_{1}}\otimes {\mathbf{f}}_{{\alpha }_{2}}\otimes \dots \otimes {\mathbf{f}}_{{\alpha }_{d}},\quad {\alpha }_{i}=1,2,\dots ,L,i=1,2,\dots ,d.\end{equation*}$$

Then again we reduce the eigenvalue problem ${\mathbf{U}}_{0}^{\left(d\right)}{\mathbf{F}}_{{\alpha }_{1}{\alpha }_{2}\dots {\alpha }_{d}}={\mu }_{{\alpha }_{1}{\alpha }_{2}\dots {\alpha }_{d}}{\mathbf{F}}_{{\alpha }_{1}{\alpha }_{2}\dots {\alpha }_{d}}$ to calculation of the values ${\mu }_{{\alpha }_{1}{\alpha }_{2}\cdot {\ldotp\ldotp}{\alpha }_{d}}={\mathbf{F}}_{{\alpha }_{1}{\alpha }_{2}\dots {\alpha }_{d}}^{+}{\mathbf{U}}_{0}^{\left(d\right)}{\mathbf{F}}_{{\alpha }_{1}{\alpha }_{2}\dots {\alpha }_{d}}$. With account for equation (3) we obtain

$$\begin{equation}{\mu }_{{\alpha }_{1}{\alpha }_{2}\cdot \cdot .{\alpha }_{d}}=\sum _{{k}_{1}=0}^{l}\sum _{{k}_{2}=0}^{l}\dots \sum _{{k}_{d}=0}^{l}w\left({k}_{1},{k}_{2},\dots ,{k}_{d}\right)\prod _{i=1}^{d}{\lambda }_{{\alpha }_{i}}\left({k}_{i}\right).\end{equation} \tag{ 14 }$$

Concluding this section, let us make some remarks. The set of the eigenvectors $\left\{{\mathbf{F}}_{m}\right\}$ and the eigenvalues $\left\{{\mu }_{m}\right\}$ we obtained above (for simplicity we use the indices m = 1, 2, ..., N, N = L^d in place of the more complex indices α₁, α₂, ..., α_L = 1, 2, ..., L) forms the eigenbasis of the Hamiltonian operator for the Ising model:

$$\begin{equation*}\hat{H}=\sum _{m=1}^{N}{\mu }_{m}\left\vert {\mathbf{F}}_{m}\right\rangle \left\langle {\mathbf{F}}_{m}\right\vert ,\hat{H}{\mathbf{F}}_{m}={\mu }_{m}{\mathbf{F}}_{m}.\end{equation*}$$

We can present each spin configuration S⁺ = (s₁, s₂, ..., s_N ), s_i = ±1 as an expansion in these basis eigenvectors: $\mathbf{S}={\sum }_{1}^{N}{a}_{m}{\mathbf{F}}_{m}$. Our analysis shows that account for long-range interactions does not change the basic eigenvectors $\left\{{\mathbf{F}}_{m}\right\}$ as well as the expansion coefficients a_m . The situation differs when we examine the energy of the state, which we can write as

$$\begin{equation*}E\left(\mathbf{S}\right)={\mathbf{S}}^{+}\hat{H}\mathbf{S}=\sum _{m=1}^{N}{a}_{m}^{2}{\mu }_{m}.\end{equation*}$$

A change of the long-range interaction does not influence the coefficients ${a}_{m}^{2}$, but the eigenvalues μ_m change. In other words, the long-range interaction does not change the expansion of the configuration in the basis eigenvectors; however, it leads to a shift of the energy of the state E(S).

The authors of paper [7] developed an algorithm for calculation of the free energy of a finite spin system using eigenvalues and eigenvectors of its connection matrix. In this paper, they analyzed the one- and two-dimensional Ising systems accounting for interactions with the nearest neighbors only. In particular, they showed that their approach allowed us to approximate accurately the known exact expressions for the free energy even when the size of the spin system is not very large. We find it interesting to estimate the influence of the long-range interaction on the behavior of the free energy using the obtained expressions (2), (8), and (12).

Let us examine the one-dimensional system. Let w(1) = 1 and w(2) ≡ b be the constants of interaction with the nearest and the next nearest spins, respectively. Then from equation (2) it follows that μ_α = 2 cos φ_α + 2b cos 2φ_α , where φ_α = 2π(α − 1)/L and α = 1, 2, ..., L. We obtain the density P(μ) by integrating the delta function $\delta \left(\mu -2\enspace \mathrm{cos}\enspace {\varphi }_{\alpha }-2b\enspace \mathrm{cos}\enspace 2{\varphi }_{\alpha }\right)$ with respect to the variable ${\varphi }_{\alpha }\in \left[0,2\pi \right]$. In what follows, we show that the spectrum changes drastically when b = 1/4. This is the reason why we examine separately the spectrum to the left and to the right from this point.

We begin with the simplest case b = 0. After rather simple calculations we obtain

$$\begin{equation*}P\left(\mu \right)=\frac{L}{\pi \sqrt{4-{\mu }^{2}}},\quad -2{< }\mu {< }2.\end{equation*}$$

We see that the spectrum is limited to the interval μ ∈ (−2, 2) and there is a divergence at the spectrum ends (figure 1).

• (a)  
  Let examine the values of b inside the interval 0 < b < 1/4. In this case the spectral density is nonzero only inside the interval μ ∈ (−2 + 2b, 2 + 2b), where

  $$\begin{equation}P\left(\mu \right)=\frac{L}{2\pi }\frac{{Q}_{0}}{\sqrt{\left(2+2b-\mu \right)\left(\mu +2-2b\right)}},\quad -2+2b{\leqslant}\mu {\leqslant}2+2b.\end{equation} \tag{ 15 }$$

  In equation (15) Q₀ = Q₀(μ) is a slow function without singularities on the interval in question:

  $$\begin{equation*}{Q}_{0}={\left\vert \frac{\left(1+4b+\sqrt{1+4b\left(\mu +2b\right)}\right)\left(1-4b+\sqrt{1+4b\left(\mu +2b\right)}\right)}{1+4b\left(\mu +2b\right)}\right\vert }^{1/2}.\end{equation*}$$

  We see that account for the interaction with the next nearest neighbors shifts the spectrum in the positive direction by a value equal to 2b. In the same time, divergences at the ends of the spectrum (μ → ±2 + 2b) remain. Note, when b → 1/4 the left end of the spectrum, where $P\left(\mu \right)\approx L/\pi \sqrt{\left(1-4b\right)\left(\mu +2-2b\right)}$, is much higher its right end where $P\left(\mu \right)\approx L/\pi \sqrt{\left(1+4b\right)\left(2+2b-\mu \right)}$.
• (b)  
  When b = 1/4 the spectral density is nonzero on the interval μ ∈ [−3/2, 5/2]. In this case equation (15) takes the form

  $$\begin{equation*}P\left(\mu \right)=\frac{L}{2\pi }\frac{{\left(2+\sqrt{3/2+\mu }\right)}^{1/2}}{{\left(3/2+\mu \right)}^{3/4}{\left(5/2-\mu \right)}^{1/2}},\quad -3/2{\leqslant}\mu {\leqslant}5/2.\end{equation*}$$

  As we see, in this case the divergence at the left end of the spectrum (μ → −3/2) is much stronger.
• (c)  
  Let us analyze the values of b > 1/4. We introduce the notations

  $$\begin{equation*}{\mu }_{\mathrm{min}}=-2b-\frac{1}{4b},\qquad {\mu }_{\mathrm{a}\mathrm{v}}=-2+2b,\qquad {\mu }_{\mathrm{max}}=2+2b.\end{equation*}$$

  Here μ_min and μ_max are the lower and the upper boundaries of the spectrum, respectively. Next, μ_av is a point between the spectrum boundaries: μ_min < μ_av < μ_max. The spectrum is a composite: in the different regions of the interval it is described by two different functions that does not match

  $$\begin{equation*}P\left(\mu \right)=\frac{L}{\pi \sqrt{\mu -{\mu }_{\mathrm{min}}}}\enspace \left(\frac{{Q}_{1}}{\sqrt{{\mu }_{\mathrm{a}\mathrm{v}}-\mu }}+{Q}_{2}\right)\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace {\mu }_{\mathrm{min}}{\leqslant}\mu {< }{\mu }_{\mathrm{a}\mathrm{v}},\end{equation*}$$

  $$\begin{equation*}P\left(\mu \right)=\frac{L}{\pi {Q}_{1}\sqrt{\left(\mu -{\mu }_{\mathrm{min}}\right)\left({\mu }_{\mathrm{max}}-\mu \right)}}\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace {\mu }_{\mathrm{a}\mathrm{v}}{\leqslant}\mu {\leqslant}{\mu }_{\mathrm{max}},\end{equation*}$$

  where Q₁ = Q₁(μ) and Q₂ = Q₂(μ) are slow functions without singularities at the abovementioned intervals

  $$\begin{equation*}{Q}_{1}={\left\vert \frac{4b-1+\sqrt{4b\left(\mu -{\mu }_{\mathrm{min}}\right)}}{4b+1+\sqrt{4b\left(\mu -{\mu }_{\mathrm{min}}\right)}}\right\vert }^{\enspace 1/2},\qquad {Q}_{2}=\frac{\sqrt{4b}}{{Q}_{1}{\left\vert {\left(4b+1\right)}^{2}-4b\left(\mu -{\mu }_{\mathrm{min}}\right)\right\vert }^{\enspace 1/2}}.\end{equation*}$$

  We see that at the point μ = μ_av the function P(μ) is discontinues: when μ → μ_av − 0, there is a singularity $P\left(\mu \right)\sim {\left(\mu -{\mu }_{\mathrm{min}}\right)}^{-1/2}$; when μ → μ_av + 0 the spectrum density is finite: $P\left(\mu \right)=2bL/\pi {\left(4b-1\right)}^{3/2}$. In figure 1, we show how the spectrum changes when b increases.

Zoom In Zoom Out Reset image size

The following picture arises from our analysis. If b ⩽ 1/4, the spectrum as a whole shifts by the value 2b; the lower and upper spectrum boundaries are μ_min = −2 + 2b and μ_max = 2 + 2b, respectively; the width of the spectrum remains constant. When b > 1/4, an increase of b leads to the spectrum broadening. Namely, the value of μ_min decreases: μ_min = −2b − 1/4b and the value of μ_max increases: μ_max = 2 + 2b. The spectrum width also increases as μ_max − μ_min = 2 + 4b + 1/4b. For an arbitrary b there are divergences both at the lower and upper spectrum boundaries. Moreover, if b > 1/4 an additional divergence at the point μ = μ_av ≡ −2 + 2b appears. In figure 2(a), we show this picture schematically.

Zoom In Zoom Out Reset image size

We analyze the two-dimensional Ising system in the simplest case supposing an isotropic interaction and accounting for interactions with the nearest and next nearest neighbors only. Let the constant of interaction with the nearest neighbors be equal to one and by b = w(1, 1) we denote the interaction with the next nearest neighbors. Then from equation (9), we obtain μ_αβ = 2 cos φ_α + 2 cos φ_β + 4b cos φ_α 2 cos φ_β , where φ_α = 2π(α − 1)/L, and φ_β = 2π(β − 1)/L. Consequently, to obtain the spectral density P(μ) we have to integrate the delta function $\delta \left(\mu -{\mu }_{\alpha \beta }\right)$ with respect to the variables ${\varphi }_{\alpha }\in \left[0,2\pi \right]$ and ${\varphi }_{\beta }\in \left[0,2\pi \right]$. We examine separately the behavior of the spectrum to the left and to the right from the point b = 1/2, since at this point it changes drastically.

• (a)  
  We begin with the case 0 ⩽ b < 1/2, where the spectral density is nonzero only when −4 + 4b ⩽ μ ⩽ 4 + 4b.When b = 0, we obtain

  $$\begin{equation}P\left(\mu \right)=\frac{{L}^{2}}{2{\pi }^{2}}K\left(m\right),\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\enspace m=1-\frac{{\mu }^{2}}{16},\quad \mu \in \left[-4,\enspace 4\right].\end{equation} \tag{ 16 }$$

  In equation (16), K(m) is a complete elliptic integral of the first kind. We see that the spectrum is symmetric with respect to the point μ = 0 and since $K\left(m\right)\approx -\mathrm{ln}\left\vert \mu \right\vert$ at m → 1, it has a logarithmic divergence when μ → 0 [see figure 3(a)].When b ≠ 0 (b < 1/2), we obtain a more general expression

  $$\begin{equation}P\left(\mu \right)=\frac{{L}^{2}}{2{\pi }^{2}\sqrt{1+b\mu }\enspace }K\left(m\right),\qquad m=\frac{1-{\left(b-\mu /4\right)}^{2}}{1+b\mu }\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace \mu \in \left[-4+4b,\enspace 4+4b\right].\end{equation} \tag{ 17 }$$

  As it follows from equation (17), when b increases, the spectrum as a whole shifts to the right by the value of 4b. In the same time, the maximum of the function P(μ) at μ = −4b, where there is the logarithmic singularity at m → 1, shifts to the left. In addition, at the right end of the spectrum (μ = 4 + 4b) the value of the spectral density decreases as P(μ) = L²/4π (1 + 2b) and at the left end (μ = −4 + 4b) it increases as P(μ) = = L²/4π(1 − 2b). We show the changes of P(μ) in figure 3(a).
• (b)  
  When b = 1/2, equation (16) takes the form

  $$\begin{equation*}P\left(\mu \right)=\frac{{L}^{2}}{2{\pi }^{2}\sqrt{1+\mu /2}\enspace }K\left(m\right),\quad m=\frac{6-\mu }{8}\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace \mu \in \left[-2,\enspace 6\right].\end{equation*}$$

  As we see in figure 3(b), the entire spectrum is to the right from the point μ = −2. Near this point when μ → −2 and, consequently, m → 1, we have $P\left(\mu \right)\sim -\mathrm{ln}\left\vert \mu +2\right\vert /\sqrt{\mu +2}$. Then the density has both the logarithmic and power divergences. When μ increases the value of P(μ) decreases gradually and reaches its minimal value P(μ) = L²/8π at the right end of the spectrum where μ = 6.
• (c)  
  The case b > 1/2. Now the spectral density is nonzero at the interval −4b ⩽ μ ⩽ 4 + 4b. There are three parts of this interval where the discontinuous function P(μ) is described by three different expressions:

  $$\begin{equation*}P\left(\mu \right)=\frac{{L}^{2}}{2{\pi }^{2}\enspace }R\enspace K\left(m\right),\end{equation*}$$

  where

  $$\begin{equation*}R=\frac{2}{\sqrt{{\left(b-\mu /4\right)}^{2}-1}\enspace },\qquad m=\frac{{\left(b+\mu /4\right)}^{2}}{{\left(b-\mu /4\right)}^{2}-1}\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace -4b{\leqslant}\mu {< }-1/b;\end{equation*}$$

  $$\begin{equation*}R=\frac{2}{b+\mu /4\enspace },\qquad m=\frac{{\left(b-\mu /4\right)}^{2}-1}{{\left(b+\mu /4\right)}^{2}}\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace -1/b{< }\mu {\leqslant}-4+4b;\end{equation*}$$

  $$\begin{equation*}R=\frac{1}{\sqrt{1+b\mu }\enspace },\qquad m=\frac{1-{\left(b-\mu /4\right)}^{2}}{1+b\mu }\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}\enspace -4+4b{\leqslant}\mu {\leqslant}4+4b.\end{equation*}$$

  At the left end of the first interval μ ∈ [−4b, −1/b], the density P(μ) has a finite value $P\left(\mu \right)={L}^{2}/2\pi \sqrt{4{b}^{2}-1}$; at the right end of this interval there is a logarithmic divergence since m → 1 and $P\left(\mu \right)\sim -\mathrm{ln}\left\vert \mu +1/b\right\vert /\left(4{b}^{2}-1\right)$ when μ → −1/b − 0. At the left end of the second interval there is also the same logarithmic discontinuity when μ → −1/b + 0 and at the right end of this interval the function P(μ) has a finite value P(μ) = L²/2π(2b − 1). The function P(μ) has no singularities on the third interval μ ∈ [−4 + 4b, 4 + 4b]. When μ increases, the density decreases smoothly from the value P(μ) = L²/4π(2b − 1) when μ = −4 + 4b up to the value P(μ) = L²/4π(2b + 1) when μ = 4 + 4b. In figure 3(b) we show all the transformations of the function P(μ) when b increases.

Let us summarize the analysis of this subsection. When b changes inside the interval 0 ⩽ b ⩽ 1/2, the entire spectrum shifts to the right by the value of 4b: the lower boundary of the spectrum is μ_min = −4 + 4b and its upper boundary is μ_max = 4 + 4b. Inside the interval there is a logarithmic discontinuity at the point μ = μ_av (μ_av = −4b when 0 ⩽ b ⩽ 1/2). When b > 1/2, the further increase in b leads to the broadening of the spectrum and change of its form. The lower spectrum boundary μ_min = −4b decreases and its upper boundary μ_max = 4 + 4b increases. The spectrum density has two singularities: the divergence at μ = μ_av and the finite jump at μ = μ₀, where μ_av = −1/b and μ₀ = −4 + 4b. Figure 2(b) presents this picture schematically.

In the case of Ising systems of higher dimensions, due to significant mathematical difficulties we were unable to obtain analytical expressions for the spectral density. However, it is possible to make some general conclusions. For example, when we account for the nearest neighbors only, we can present the spectral density P_d (μ) for the d-dimensional system as a convolution of the spectral densities of the systems of lower dimensions:

$$\begin{equation*}{P}_{d}\left(\mu \right)=\underset{-\infty }{\overset{\infty }{\int }}{P}_{{d}_{1}}\left(x\right){P}_{{d}_{2}}\left(\mu -x\right)\mathrm{d}x,\quad {d}_{1}+{d}_{2}=d.\end{equation*}$$

Due to the integration the divergence of P_d (μ) is weaker the divergences of ${P}_{{d}_{1}}\left(x\right)$ and ${P}_{{d}_{2}}\left(x\right)$. However, there are many cases when the density P_d (μ) has singularities. For example, suppose that ${P}_{{d}_{1}}\left(x\right)$ and ${P}_{{d}_{2}}\left(x\right)$ have rather strong divergences at the points x₁ and x₂, respectively: ${P}_{{d}_{1}}\left(x\to {x}_{1}\right)\sim {\left\vert x-{x}_{1}\right\vert }^{-{r}_{1}}$ and ${P}_{{d}_{2}}\left(x\to {x}_{2}\right)\sim {\left\vert x-{x}_{2}\right\vert }^{-{r}_{2}}$, where 0.5 ⩽ r₁ r₂ < 1. Then the density P_d (μ) (d = d₁ + d₂) also has a singularity at the point μ = x₁ + x₂. This may be a power singularity ${P}_{d}\left(\mu \to {x}_{1}+{x}_{2}\right)\sim {\left\vert \mu -{x}_{1}-{x}_{2}\right\vert }^{-r}$ when r > 0, or a logarithmic singularity ${P}_{d}\left(\mu \to {x}_{1}+{x}_{2}\right)\sim -\mathrm{ln}\left\vert \mu -{x}_{1}-{x}_{2}\right\vert$ when r = 0. In both formulas r = r₁ + r₂ − 1. Moreover, the singularity may be a combination of the logarithmic and power singularities. To avoid misunderstanding we note that the functions ${P}_{{d}_{1}}\left(x\right)$ and ${P}_{{d}_{2}}\left(x\right)$ are normalized functions and according definition may have only integrable singularities.

Zoom In Zoom Out Reset image size