Analytical Expressions for Ising Models on High Dimensional Lattices

Abstract

We use an m-vicinity method to examine Ising models on hypercube lattices of high dimensions $d\ge 3$. This method is applicable for both short-range and long-range interactions. We introduce a small parameter, which determines whether the method can be used when calculating the free energy. When we account for interaction with the nearest neighbors only, the value of this parameter depends on the dimension of the lattice $d$. We obtain an expression for the critical temperature in terms of the interaction constants that is in a good agreement with the results of computer simulations. For $d=5,6,7$, our theoretical estimates match the numerical results both qualitatively and quantitatively. For $d=3,4$, our method is sufficiently accurate for the calculation of the critical temperatures; however, it predicts a finite jump of the heat capacity at the critical point. In the case of the three-dimensional lattice ($d=3$), this contradicts the commonly accepted ideas of the type of the singularity at the critical point. For the four-dimensional lattice ($d=4$), the character of the singularity is under current discussion. For the dimensions $d=1, 2$ the m-vicinity method is not applicable.

1. Introduction

Statistical physics provides effective methods of analysis, allowing us to investigate large systems of elementary “agents” and to determine macroscopic characteristics—in particular, the free energy based on interactions between the “agents”. If we know the free energy—by means of computer simulations or theoretical calculations—we can calculate such properties of the system as the internal energy, magnetization, heat capacity, and susceptibility. The singularities of the temperature dependences of these characteristics define the critical temperatures at which the internal restructurings take place in the system and phase transitions occur.

Such results (even if they are not quite accurate) are important, and not only for physicists. In the second half of the 1980s, the statistical physics methods were applied to estimate the storage capacity of the Hopfield neural network. Later on, a lot of investigations in the field of neural science based on the statistical physics followed (see, for example, [1,2,3]). At the same time, the statistical physics methods became popular in the combinatorial optimization problems [4,5,6]. Starting from the mid-1990s, a new scientific branch named econophysics appeared. In econophysics, the statistical physics methods are the main instruments for analyzing economic models [7,8].

Physics provides a wide variety of methods for the calculation of the free energy, from computer simulations (where one uses the Metropolis or the Wang-Landau algorithms [9,10,11]) to cumbersome theoretical approaches of the type of the renormalization group or the transfer-matrix methods [12,13,14].

In the present paper, we sum up the results obtained when developing the m-vicinity method (at first, we called it “the n-vicinity method”, but then it became clear that the more appropriate term was “the m-vicinity method”) for analysis of the Ising systems. Our method allows us to calculate the free energy for an arbitrary connection matrix. We present a review of our results for the Ising models on hypercubic lattice of the dimensions $d=3,4,5,6$, and $7$. The dimensions $d=1$ and $d=2$ are absent in this list because this method is not applicable for such lattices.

There is an enormous number of papers of theoretical or numerical studies on the behavior of spin systems on specific types of lattices. Here, we cite only the papers that we used in the course of our work. Cubic lattices are studied, for example, in [9,15,16,17,18,19]. Four-dimensional lattices are discussed in [20,21,22]. Studies on higher-dimensional systems are very rare. Our single source on $d\ge 5$ was [23].

In the next section, we justify the main approximation of our method. It consists of the substitution of the Gaussian distribution in place of the unknown state distribution of the given Hamiltonian.

In Section 3, we show that the Gaussian approximation is the first-order term in the expansion of the density of states in a perturbation theory series in a small parameter ${\epsilon }_{\max }$. In the case of the planar lattice (d = 2), ${\epsilon }_{\max }\approx 0.7$, and this value is not sufficiently small. When the lattice dimension $d$ increases, the value of ${\epsilon }_{\max }$ decreases quickly, and the Gaussian approximation works very well.

In Section 4, we present in detail our results obtained with the aid of Gaussian approximation of the density of states. The mean and the variance of the Gaussian distribution that we use coincide with the first and the second moments of the density of states of the given system. We define the boundaries of the method applicability and obtain analytical expressions for the critical characteristics of the system. In particular, they are the critical value of the inverse temperature and the jump of the heat capacity. Our analytical results match quite accurately with the results of computer simulations. We find out that the higher the dimension of the lattice, the better our estimates; when $d\ge 5$, the relative error is of the order of the tenth or the hundredth percent.

In Section 5, we check whether the account of the second-order terms of the perturbation theory improves our results for the critical temperature. Here, we approximate the density of states with a distribution whose first three moments coincide with the moments of the state distribution of the system. (To do this a more advanced knowledge of the mathematical statistics is necessary [24,25].) We find that, while the role of second order parameters is negligible, we can achieve almost perfect agreement with computer simulations by introducing an adjustable parameter.

In Section 6, we present a detailed comparison of the theoretical results and computer simulations for the Ising model on the cubic lattice ($d=3$). For such a lattice, we examine the role of the long-range interaction; in particular, we discuss the interactions with the next-nearest neighbors and the next-next-nearest neighbors.

Finally, in Section 7, we sum up the strengths and weaknesses of the m-vicinity method. The details of the calculations are in the Appendix A.

2. Main Approximation of m-Vicinity Method

Let us examine the Ising model on a multidimensional cubic lattice, which is a system of $N$ spins $s_i=\left\{\pm 1\right\}$, $i=1,2,\dots N$ situated at the nods of a hypercubic lattice. In what follows, we assume the periodic boundary conditions.

The Hamiltonian of the system is

$${{\rm E}}_H=E-mH,\hspace{1em}E=-\frac{1}{2N}{\displaystyle \sum _{i,j=1}^N{J}_{ij}{s}_i{s}_j},\hspace{1em}m=\frac{1}{N}{\displaystyle \sum _{i=1}^N{s}_i},$$

where $J={\left(J_{ij}\right)}_1^N$ is a connection matrix, H is a magnetic field, and m is a magnetization of the state $s=(s_1,s_2,\dots ,s_N)$. The partition function of the system is

$$Z={\displaystyle {\sum }_{E}D(E)\exp \left(-N\beta E_H\right)},$$

where $\beta$ is the inverse temperature, the summation is carried out over all the values of the energy $E$, and $D\left(E\right)$ is the density of states (the degeneracy of the energy states).

In the general case, we do not know the energy distribution $D(E)$. It seems that we can define it with the aid of the central limit theorem. Indeed, the value of E is the sum of $N(N-1)/2$ weakly connected random variables, and the Gaussian distribution

$$2^N\sqrt{\frac{N}{2\pi {\sigma }_0}}\exp \left(-\frac{1}{2}N\frac{E^2}{{\sigma }_0^2}\right), \mathrm{where} {\sigma }_0^2=\frac{1}{2N}{\displaystyle \sum _{i,j=1}^N{J}_{ij}^2},$$

(1)

describes correctly the central part of the distribution $D(E)$. However, Equation (1) is not applicable at the tails of the true distribution, while the tails provide the main contribution to the formation of the phase transition. Many authors mentioned this fact (see [15]). The m-vicinity method allows us to overcome this difficulty. The essence of the method is as follows. We divide the set of $2^N$ configurations into $N+1$ subsets ${\mathsf{\Omega }}_m$, which will be called the m-vicinities. The m-vicinity ${\mathsf{\Omega }}_m$ contains all the configurations $s$ with the same magnetization $m$. These configurations $s$ differ from the configuration of the ground state $s_0=(1,1,\dots ,1)$ by opposite signs of $n$ spins, where $n=N(1-m)/2$ and the number of such configurations is equal to $\left(\begin{array}{c}N\\ n\end{array}\right)$.

The density $D(E,m)$ is the energy distribution for a given m-vicinity, and the partition function of the system is

$$Z=\sum _{n=0}^N\sum _{E}D(E,m)\exp [-N\beta (E-mH)].$$

(2)

In the general case, we do not know the true distribution $D(E,m)$. However, we do know the exact values of its mean energy $E_m$ and its variance, which we denote as $N^{-1}{\sigma }_m^2$ (see Appendix A or [26,27,28]). In the case of the Ising model on the hypercube, these expressions in the limit $N\to \infty$ are sufficiently simple,

$$E_m=E_0{m}^2, N^{-1}{\sigma }_m^2=N^{-1}{\sigma }_0^2\cdot {(1-m^2)}^2, E_0=-\frac{1}{2N}{\displaystyle \sum _{i,j=1}^N{J}_{ij}},$$

(3)

where $E_0=E(s_0)$ is the energy of the ground state $s_0$ of the Ising Hamiltonian.

Let us explain the purpose of dividing the whole set of configurations into m-vicinities. In the vicinity ${\mathsf{\Omega }}_m$, the energy $E$ behaves as a random value and due to the central limit theorem, we can approximate accurately the central part of $D(E,m)$ by Gaussian distribution with the mean $E_m$ and the variance $N^{-1}{\sigma }_m^2$:

$$D(E,m)=\frac{\sqrt{N}}{\sqrt{2\pi }{\sigma }_m} \left(\begin{array}{l}N\\ n\end{array}\right) \exp \left[-\frac{1}{2}N{\left(\frac{E-E_m}{{\sigma }_m}\right)}^2\right], n=N\left(1-m\right)/2.$$

(4)

It is evident that the sum ${\displaystyle {\sum }_{m}D(E,m)}$ differs from the Gaussian distribution (1) and it better describes the tails of the true distribution.

Let us note here that in the limit ${\sigma }_m\to 0$, we can replace the exponent in Equation (4) by a delta function $\delta (E-E_m)$. Then Equation (2) takes the classical form known from the mean field theory, which provides the Bragg–Williams results [29]. However, the value of ${\sigma }_m$ differs from zero for all the types of the connection matrices (an exclusion is the case of the complete graph when all the matrix elements are equal—see Equation (A4)) and the replacement of the Gaussian (4) by the delta function is not correct. It is the account for the value ${\sigma }_m\ne 0$ that leads to a much better agreement of theoretical estimates and results of simulations. Of course, the distribution $D(E,m)$ is not purely Gaussian since its higher odd moments are not equal to zero (see Appendix A); however, their contribution is sufficiently small. In what follows, we analyze when the Gaussian approximation (4) is applicable and how a deviation of the density of states $D(E,m)$ from Equation (4) influences the results. In the next section, we show that expression (4) is the first-order term of the perturbation theory in a small parameter ${(2q)}^{-1/2}$, where $q$ is an effective number of neighbors (see Equation (16)).

3. Small Parameter in m-Vicinity Method

The basis of the m-vicinity method is the abovementioned approximate description of the density of states. To analyze the approximation, we briefly repeat the calculations of the paper [30]. The starting point is as follows. We do not know a true energy distribution $D(E,m)$ but we do know the first moments of this distribution. In particular, we know the mean and the variance (3). Let us define the small parameter allowing us to expand the function $D(E,m)$ in a perturbation theory series.

We present $D(E,m)$ in the form

$$D(E,m)=\left(\begin{array}{l}N\\ n\end{array}\right)\exp \left[-N\phi (m,E)\right],$$

where $\phi =\phi (m,E)$ is an unknown function and use the Stirling formula to replace the summation in (2) by integration. Then, up to an insignificant constant, we obtain the partition function of the form

$$Z~{\displaystyle \underset{0}{\overset{1}{\int }}dm{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}dE} e^{-NF(m,E)}},$$

(5)

where

$$\begin{array}{l}F(m,E)=S(m)+\beta (E-mH)+\phi (m,E),\\ S(m)=-\ln 2+\frac{1}{2}\left[\left(1+m\right)\ln \left(1+m\right)+\left(1-m\right)\ln \left(1-m\right)\right].\end{array}$$

(6)

Let us estimate the integral (5) using the saddle point method. The equations for the saddle point are

$$\frac{\partial F}{\partial m}=\frac{1}{2}\ln \left(\frac{1+m}{1-m}\right)+\frac{\partial \phi }{\partial m}-\beta H=0, \frac{\partial F}{\partial E}=\frac{\partial \phi }{\partial E}+\beta =0.$$

(7)

The solutions of these equations $m=M$ and $E=U$ are the spontaneous magnetization and the internal energy, respectively. Substituting these values in Equation (6), we obtain the free energy $f(\beta )=F(M,U)$.

Now, we turn to defining the small parameter of the m-vicinity method. Since the magnetic field H does not influence the distribution $D(E,m)$, we set here $H=0$. We write the function $\phi (m,E)$ as a perturbation theory series in the vicinity of point $E=E_m$:

$$\phi =\frac{1}{2}{\epsilon }^2+\frac{1}{3!}{\kappa }_3{\epsilon }^3+\frac{1}{4!}{\kappa }_4{\epsilon }^4+\dots , \mathrm{where} \epsilon =\frac{E-E_m}{{\sigma }_m}.$$

(8)

The quantities ${\kappa }_k$ up to a sign coincide with the semi-invariants of the distribution $D(E,m)$ (see Appendix A).

The main idea of the m-vicinity method is the possibility to restrict ourselves by accounting for just a few first terms of the series (8). We can do this only if $\left|\epsilon \right|<<1$ and in this case the m-vicinity method is sufficiently accurate. Let us clarify that we are not interested in the values of $\epsilon =\epsilon (m,E)$ over the whole region of definition of the parameters $m$ and $E$. We have to know the expansion (8) only in a small vicinity of the saddle point, which is close to the values $m=M$ and $E=U$. Consequently, the small parameter we are looking for is

$${\epsilon }_0=\frac{U-E_0{M}^2}{{\sigma }_0(1-M^2)}.$$

The smallness of the parameter ${\epsilon }_0$ is the condition for applicability of the m-vicinity method.

In [30], there is a detailed analysis of the values of this parameter for different models. Here, we restrict ourselves with an estimate of ${\epsilon }_0$ by means of the reverse-reasoning method. Suppose the parameter ${\epsilon }_0$ is small and it is sufficient to use only the first term of the series (8) and set $\phi ={\epsilon }^2/2$. Then the second of the Equation (7) takes a simple form $\epsilon =-\beta {\sigma }_m$ and we can rewrite the first of the Equation (7) as follows:

$$\epsilon \left(\epsilon -\frac{E_0}{{\sigma }_0}\right)=-\frac{1-m^2}{4m}\ln \left(\frac{1+m}{1-m}\right).$$

(9)

We are interested in the behavior of the quantity $\epsilon ={\epsilon }_0$, where ${\epsilon }_0$ is the root of Equation (9). The value of $\left|{\epsilon }_0\right|$ reaches its maximum at the critical point $\beta ={\beta }_c$ (see [30]). In the limit $m\to 0$ ($\beta \to {\beta }_c$), Equation (9) takes the form

$${\epsilon }^2-\epsilon \frac{E_0}{{\sigma }_0}+\frac{1}{2}=0.$$

The solution of this equation is

$${\epsilon }_{\max }=-\frac{\left|E_0\right|}{2{\sigma }_0}+\frac{1}{2}\sqrt{\frac{E_0^2}{{\sigma }_0^2}-2},$$

(10)

and $\max \left|{\epsilon }_0\right|=\left|{\epsilon }_{\max }\right|$. Consequently, the m-vicinity method is applicable when the small parameter $\left|{\epsilon }_0\right|\le \left|{\epsilon }_{\max }\right|$. According to Equation (10), it is necessary that $E_0^2/{\sigma }_0^2\ge 2$ or $\left|{\epsilon }_{\max }\right|\le 2^{-1/2}$.

A detailed analysis shows (see [27]) that a stricter inequality $E_0^2/{\sigma }_0^2\ge 8/3$ or $\left|{\epsilon }_{\max }\right|\le 6^{-1/2}$ defines the framework of the m-vicinity method. In another case, we obtain a jump of the spontaneous magnetization at the critical point, and this contradicts the known results.

We would like to mention here that when analyzing the dependence of our results on the lattice dimension $d$, we see that the ratio $E_0^2/{\sigma }_0^2$ depends on the effective number of the neighbors $q$ introduced below in Equation (16). This parameter defines the number of interactions that we take into account. We can rewrite Equation (10) in terms of $q$; then ${\epsilon }_{\max }=\sqrt{q}\left(-1+\sqrt{1-4/q}\right)/2\sqrt{2}$. From this equation, it follows that when $q>>1$ the value of $\left|{\epsilon }_{\max }\right|\sim 1/\sqrt{2q}$. Consequently, when q increases, the value of the small parameter $\left|{\epsilon }_{\max }\right|$ decreases quite rapidly and the accuracy of the m-vicinity method increases.

In Figure 1, we show the dependence of $\left|{\epsilon }_{\max }\right|$ on $d$ when we account for the nearest neighbors only (in this case $E_0^2/{\sigma }_0^2=d$). For a planar lattice, $\left|{\epsilon }_{\max }\right|\approx 0.7$. Obviously, it is problematic to use it as a small parameter. Moreover, for a planar lattice, the inequality $E_0^2/{\sigma }_0^2\ge 8/3$ is not satisfied. Next, for a cubic lattice, the last inequality is satisfied and the value of $\left|{\epsilon }_{\max }\right|\approx 0.36$ is more appropriate as a small parameter. Again, when $d$ increases, the value of $\left|{\epsilon }_{\max }\right|$ decreases, and the accuracy of the method increases accordingly.

Summing up, we can say that the Gaussian approximation is the first order of the perturbation theory in the small parameter $\left|{\epsilon }_{\max }\right|$.

4. First Order of Perturbation Theory: Gaussian Approximation

In this section, we use the Gaussian approximation (4) for the distribution $D(E,m)$. If in Equations (6)–(8) we set $\phi ={\epsilon }^2/2$, Equation (6) becomes

$$F(m,E)=S(m)+\beta (E-mH)+\frac{1}{2}{\left(\frac{E-E_m}{{\sigma }_m}\right)}^2;$$

(11)

the equations for the saddle point (7) are

$$\frac{1}{2}\ln \left(\frac{1+m}{1-m}\right)=2m\beta \left(\left|E_0\right|+\epsilon {\sigma }_0\right)+\beta H, \epsilon =-\beta {\sigma }_0(1-m^2).$$

(12)

Eliminating the variable $\epsilon$ from these equations, we obtain an equation of state

$$\frac{1}{2}\ln \left(\frac{1+m}{1-m}\right)=2m\beta \left[\left|E_0\right|-\beta {\sigma }_0^2(1-m^2)\right]+\beta H.$$

(13)

This equation differs from the well-known equation of Bragg and Williams [29] by a term proportional to ${\beta }^2$. After a transformation of Equation (11), with the account for Equations (12) and (13), we obtain the expression for the free energy:

$$f(\beta )=\frac{1}{2}\ln \left(\frac{1-m^2}{4}\right)+\beta \left|E_0\right|m^2-\frac{1}{2}{\beta }^2{\sigma }_0^2(1-m^2)(1+3m^2)$$

(14)

that defines $f$ as a function of $\beta$ and $m$. In Equation (14), the spontaneous magnetization $m=m(\beta ;H)$ is the solution of Equation (13). Again, by setting ${\sigma }_0^2=0$, we recover the known result from the mean field theory.

We would like to point out that the equality $E_0^2/{\sigma }_0^2=d$ is true only when we account for the interaction with the nearest neighbors. Then the inequality $E_0^2/{\sigma }_0^2\ge 8/3$ excludes the one- and two-dimensional Ising models from the consideration. Consequently, in this case, we can use the $m$-vicinity method only for the spatial dimensions $d\ge 3$.

4.1. Critical Point

In this subsection, we set $H=0$ and define the critical temperature under this assumption. For this purpose, we rewrite Equation (13) as

$$\frac{1}{2m}\ln \left(\frac{1+m}{1-m}\right)=qb\left[1-b(1-m^2)\right].$$

(15)

where we introduce dimensionless characteristics

$$b=\beta \frac{{\sigma }_0^2}{\left|E_0\right|}, q=\frac{2E_0^2}{{\sigma }_0^2}=\frac{{\left({\displaystyle \sum _{i,j=1}^N{J}_{ij}}\right)}^2}{N{\displaystyle \sum _{i,j=1}^N{J}_{ij}^2}}.$$

(16)

The defined variables (16) are convenient since in Equation (15), the only parameter that depends on the type of the lattice is $q$.

When $m\to 0$ ($\beta \to {\beta }_c$), Equation (15) takes the form

$$\frac{1}{q}=b(1-b).$$

(17)

Now, with the account of Equation (16), we obtain the expression for the critical temperature

$${\beta }_c=b_c\frac{\left|E_0\right|}{{\sigma }_0^2}, \mathrm{where} b_c=\frac{1}{2}\left(1-\sqrt{1-\frac{4}{q}}\right).$$

(18)

The critical value of $b_c$ depends only on the number $q$, and this parameter depends on the mean and the variance of the elements of the connection matrix. Since $q$ is a characteristic of the interactions that we take into account, we regard it as an effective number of the neighbors. In particular, if we account for an isotropic interaction with the nearest neighbors only, ${\beta }_c=b_c$, $q=2d$, and $d$ is the dimension of the lattice. Note that in this case, $q$ is exactly equal to the number of spins with which the given spin interacts. Then Equation (18) describes pretty well the results of computer simulations for all the dimensions, which we examined (see Figure 2 and

Table 1. Critical temperatures and relative errors for the lattice dimensions $3\le d\le 7$. Results of mean field theory are also shown.

	$\mathbf{Critical} \mathbf{Temperature} {\mathit{\beta }}_{\mathit{c}}$
	d = 3	d = 4	d = 5	d = 6	d = 7
Best numerical results [9,21,23]	0.22166	0.14970	0.11392	0.09230	0.07771
Account for third moment, Equation (40), $k=-0.8$	0.22155	0.14894	0.11372	0.09227	0.07772
Gaussian approximation (18)	0.21132	0.14645	0.11270	0.09175	0.07742
Account for third moment, Equation (40), $k=1$	0.20196	0.14366	0.11151	0.09113	0.07706
Mean field theory	0.16667	0.12500	0.10000	0.08333	0.07143
	Relative error$\mathit{E}\mathit{r}\mathit{r}$
Account for third moment, Equation (40), $k=-0.8$	0.02%	0.25%	0.09%	0.02%	0.01%
Gaussian approximation (18)	2.39%	1.10%	0.54%	0.30%	0.18%
Account for third moment, Equation (40), $k=1$	4.65%	2.06%	1.07%	0.64%	0.42%
Mean field theory	14.16%	8.99%	6.50%	5.10%	4.21%

).

4.2. Analytical Expressions

Let us list the basic thermodynamic characteristics that we obtained from Equations (11)–(14).

(1) The interval $\beta <{\beta }_c$. When $\beta <{\beta }_c$ and $m=0$, from the above-mentioned equations, we obtain that the free energy, the internal energy, and the heat capacity are

$$f=-\ln 2-\frac{1}{2}{\beta }^2{\sigma }_0^2, U=-\beta {\sigma }_0^2, \mathrm{and} C={\beta }^2{\sigma }_0^2,$$

(19)

respectively.

(2) The interval $\beta \ge {\beta }_c$. To obtain the analytical expressions for the values $f$, $U$, and $C$, we solve Equation (15) for $b$:

$$b=\frac{1}{2(1-m^2)} \left[1-\sqrt{1-\frac{2(1-m^2)}{qm} \ln \left(\frac{1+m}{1-m}\right) } \right]$$

(20)

and transform Expressions (12) and (14) to the forms

$$f=S(m)-\frac{1}{2}qb{m}^2-\frac{1}{4}q{b}^2{\left(1-m^2\right)}^2,$$

(21)

$$U=E_0\left[m^2+b{\left(1-m^2\right)}^2\right],$$

(22)

$${\sigma }_E^2={\sigma }_0^2{\left(1-m^2\right)}^2+4E_0^2\frac{m^2\left(1-m^2\right){\left[1-2b\left(1-m^2\right)\right]}^2}{1-qb\left(1-m^2\right)\left[1-b\left(1-3m^2\right)\right]}.$$

(23)

Here, ${\sigma }_E^2=-d^{2}f/d{\beta }^2$ is the energy variance that is related to the heat capacity via $C={\beta }^2{\sigma }_E^2$. The expression (23) is a result of differentiating the second of Equation (12) with respect to variable $\beta$.

If we regard the magnetization as an independent variable $m\in [0,1]$ and consider the values of $f$, $U$, and ${\sigma }_E^2$ as functions of m, then Equations (20)–(23) define implicitly the dependence of $f$, $U$, and ${\sigma }_E^2$ on the inverse temperature $\beta =b\left|E_0\right|/{\sigma }_0^2$ at the interval $\beta \ge {\beta }_c$.

4.3. Critical Parameters

To examine the behavior of the thermodynamic characteristics near the critical point, we introduce a relative inverse temperature

$$t=\frac{\beta -{\beta }_c}{{\beta }_c}.$$

Omitting intermediate calculations, we only present the most important critical dependences.

(1) From Equations (19)–(23), it follows that at $\beta ={\beta }_c$ the free energy and the internal energy are continuous functions, and the heat capacity has a jump. Indeed, when $\beta >{\beta }_c$ Equation (23) holds, and in the limit $m\to 0$ ($\beta \to {\beta }_c$) we obtain

$${\sigma }_E^2={\sigma }_0^2+6E_0^2\frac{{\left(1-2b_c\right)}^2}{1-3q{b}_c^2},$$

(24)

where $b_c$ is defined by Equation (18). Comparing this expression with Equation (19), we see that at $\beta ={\beta }_c$ the energy variance has a jump and consequently the heat capacity also has a jump:

$$\Delta C=\frac{3}{2}{q}^2{b}_c^2\frac{{\left(1-2b_c\right)}^2}{1-3q{b}_c^2}.$$

(25)

For the lattice dimensions d = 5, 6, and 7, the authors of [23] used computer simulations to estimate the jumps of the heat capacity. The comparison of their results with the values following from Equation (25) shows:

$$\begin{array}{l}d=5,\hspace{1em}\hspace{1em}\Delta C^{(\exp )}=1.8703,\hspace{1em}\hspace{1em}\Delta C^{(theor)}=1.8469\\ d=6,\hspace{1em}\hspace{1em}\Delta C^{(\exp )}=1.7403,\hspace{1em}\hspace{1em}\Delta C^{(theor)}=1.7394\\ d=7,\hspace{1em}\hspace{1em}\Delta C^{(\exp )}=1.6860,\hspace{1em}\hspace{1em}\Delta C^{(theor)}=1.6824.\end{array}$$

.

We see that Equation (25) provides very good agreement with the computer simulations; the larger $d$ is, the better this agreement.

(2) When $\beta >{\beta }_c$, the value of the spontaneous magnetization near the critical point ($t\to 0$) obtained from Equation (15) is

$$M=A\sqrt{t}, A={\left(\frac{1-q{b}_c^2}{\frac{1}{3}-q{b}_c^2}\right)}^{\frac{1}{2}}.$$

(26)

This expression differs from the dependence $M\sim t^{\frac{1}{8}}$ that is valid for the two-dimensional Ising model (q = 4); however, it qualitatively coincides with the expressions $M=2\sqrt{t}$ and $M=\sqrt{3 t}$ obtained in the framework of the van der Waals theory [13] and the mean field theory, respectively. Note that our result tends to the result of the mean field theory in the limit $q>>1$. It is also worthwhile to note that Equation (26) predicts a larger magnetization than the mean field theory, $A>\sqrt{3}$ for any $q\ge 4$. If $q\le 11$ ($A>2$), the van der Waals magnetization is less than the value (26) and when $q>11$ ($A<2$), it is larger than the value (26).

(3) From Equation (13), it follows that at the critical point, the susceptibility has a jump

$${\chi }^{-1}=\{\begin{array}{c}\begin{array}{l}-q t \sqrt{1-4/q}\hspace{1em},\hspace{1em}t<0\\ \end{array}\\ \hspace{1em}2q t \sqrt{1-4/q}\hspace{1em},\hspace{1em}t>0.\end{array}$$

(27)

Comparing with the analogous expression of the mean field theory, we see that an extra factor $\sqrt{1-4/q}$ appears in Equation (27), and it tends to 1 when $q>>1$. As opposed to the mean field theory [13], in Equation (27), $q$ is the effective number of the neighbors (16).

(4) It is easy to see that our model satisfies the similarity hypothesis. Indeed, when we expand expression (13) in small parameters, $m$ and $t$, we obtain the dependence $H=H(m, t)$ that can be rewritten in the classical form ${\beta }_{c}H=m{\left|m\right|}^{\delta -1}{h}_s(t{m}^{-2})$ with the critical exponent $\delta =3$ and the scaling function

$$h_s(x)=\left(\frac{1}{3}-q{b}_c^2\right)-x\left(1-q{b}_c^2\right).$$

(28)

4.4. Magnetization Distribution

The integral

$$P(m)=Z^{-1}{\displaystyle \underset{E_0}{\overset{\left|E_0\right|}{\int }}D(E,m)dE}$$

(29)

defines the probability of finding the system in a state with the magnetization m. In the Gaussian approximation we use here, $D(E,m)$ is defined by Equation (4). We estimate this integral with the aid of the saddle point method. The value of $P(m)$ is accurate to a normalization constant $P_0$

$$P(m)=P_0 e^{-N\mathsf{\Phi }(m)}, \mathrm{where} \mathsf{\Phi }(m)=S(m)+\beta E_m-\frac{1}{2}{\beta }^2{\sigma }_m^2.$$

(30)

In Figure 3a, we show the typical behavior of the curves (30). As we might expect, after crossing the critical point, the bimodal distribution replaces the unimodal distribution. We can use Equation (30) when analyzing the Binder cumulant $Q=1-\langle m^4\rangle /3{\langle m^2\rangle }^2$ (see [31]). In Figure 3b, we show the curves $Q=Q(B)$ for cubic lattices whose linear sizes are $L=$ 8, 10, and 12. We are interested in the value of the cumulant $Q_c=Q({\beta }_c)$ at the critical point. To calculate $Q_c$ we use the Taylor expansion of the function $\mathsf{\Phi }(m)$ in Equation (30):

$$\mathsf{\Phi }(m)\approx -\frac{1}{2}{\beta }^2{\sigma }_0^2+\frac{1}{2!}{a}_2{m}^2+\frac{1}{2!}{a}_4{m}^4,$$

(31)

where

$$a_2={\frac{d^2\mathsf{\Phi }}{d{m}^2}|}_{m=0}=1+2\beta E_0+2{\beta }^2{\sigma }_0^2, a_4={\frac{d^4\mathsf{\Phi }}{d{m}^4}|}_{m=0}=1-12{\beta }^2{\sigma }_0^2.$$

This expansion is useful when $\beta \le {\beta }_c+O\left(N^{-1/2}\right)$ and the value of $P(m)$ is noticeably nonzero only when $m<<1$. In this case, the distribution (30) takes the form

$$P(m)=P_0\exp \left[-\frac{1}{2!}N{a}_2{m}^2-\frac{1}{4!}N{a}_4{m}^4\right]$$

which allows us to calculate the normalization constant $P_0$ easily as well as the mean values $\langle m^2\rangle$ and $\langle m^4\rangle$. It is not difficult to see that when $\left|\beta -{\beta }_c\right|>O\left(N^{-1/2}\right)$, we can consider the distribution $P(m)$ as purely Gaussian in full agreement with the results of [31]. However, when $\beta \to {\beta }_c$ we have $a_2\to 0$ and at the critical point, the distribution takes the form $P(m)=P_0\exp \left[-N{a}_4{m}^4/4!\right]$. Using this expression to calculate the critical value of the Binder cumulant [31] in the limit $N\to \infty$ we obtain

$$Q_c=1-\frac{\mathsf{\Gamma }(5/4)\cdot \mathsf{\Gamma }(1/4)}{3\mathsf{\Gamma }{(3/4)}^2}\approx 0.2705.$$

(32)

It is interesting that in the limit $N\to \infty$ the value $Q_c$ does not depend on any parameters of the model (the lattice type, the character of the long-range interaction, and so on). Although we obtained this result in the framework of the Gaussian approximation (4), it has a general character. Indeed, let $D(E,m)$ be an unknown function, and integration (29) results in expression (30) where $\mathsf{\Phi }(m)$ is also an unknown function. By a symmetry argument, it follows that only even powers of m are present in the Taylor expansion of this function. At the critical point, the second derivative of function $\mathsf{\Phi }(m)$ equals to zero (the bimodal distribution replaces the unimodal) and the Taylor series starts from the term $~m^4$. We assume that the function $\mathsf{\Phi }(m)$ has no singularities and its derivatives are finite. Then, we can leave only the term $~m^4$ in the exponent of the distribution $P(m)$. The reason is that when integrating over $m$ the account for the terms of the higher orders leads to corrections of the order $N^{-1/2}$. In other words, we can present the magnetization distribution as $P(m)=P_0\exp \left[-N{a}_4{m}^4/4!\right]$ where $a_4$ is an unknown, which is canceled out in the calculation of $Q_c$ and does not influence the final form of expression (32).

4.5. Density of States

In Section 4.2, we derived Equations (19)–(23), which allowed us to obtain implicitly the logarithmic density of states

$$\mathsf{\Psi }(E)=\beta E-f(\beta ); E=df/d\beta$$

using the Legendre relations. When $\beta$ changes from 0 to $\infty$, the value of $E$ changes from 0 to $E_0$ and for each $\beta$ we obtain a pair of values $E$ and $\mathsf{\Psi }(E)$. In such a way, we generate the function $\mathsf{\Psi }(E)$, which we suppose to be symmetric: $\mathsf{\Psi }(-E)=\mathsf{\Psi }(E)$. In Figure 4a, we present the comparison of our results with computer simulations.

Let us determine an explicit form of the dependence $\mathsf{\Psi }=\mathsf{\Psi }(E)$. The integral

$$D(E)={\displaystyle \underset{-1}{\overset{1}{\int }}D(E,m)dm}\approx c_1\cdot {\displaystyle \underset{-1}{\overset{1}{\int }}{e}^{-N\mathsf{\Lambda }(m,E)}dm}=c_2{e}^{N\mathsf{\Psi }\left(E\right)}$$

defines the density of states or, in other words, the number of states with the energy $E$. Here $\mathsf{\Lambda }(m,E)=S(m)+{(E-E_m)}^2/2{\sigma }_m^2$; $c_1$ and $c_2$ are non-essential constants. Therefore, the logarithmic density of states is $\mathsf{\Psi }(E)=\mathsf{\Lambda }(m_E,E)$, where $m_E=m_E(E)$ is the saddle point, which is a solution of the equation $\partial \mathsf{\Lambda }(m,E)/\partial m=0$. After some transformations, we can write this equation as

$$\frac{{\left(1-m_E^2\right)}^3}{2m_E}\ln \left(\frac{1+m_E}{1-m_E}\right)=q\left(1-\frac{E}{E_0}\right)\left(\frac{E}{E_0}-m_E^2\right).$$

(33)

From this equation, it follows that

$$\mathsf{\Psi }(E)=\{\begin{array}{ccc}\ln 2-\frac{E^2}{2{\sigma }_0^2},& \mathrm{when}& E\ge E_0{b}_c\\ -S(m_E)-\frac{1}{2{\sigma }_0^2}{\left(\frac{E-E_0{m}_E^2}{1-m_E^2}\right)}^2,& \mathrm{when}& E<E_0{b}_c,\end{array}.$$

(34)

where $b_c$ is defined by Equation (18) and

$$m_E^2\approx 1-2r\cdot \cos {\phi }_E, r=\sqrt{\frac{q}{3}\left(1-\frac{E}{E_0}\right)}, \mathrm{and} {\phi }_E=\frac{1}{3}\left[\pi +\arccos \left(\frac{9r}{2q}\right)\right].$$

(35)

Equation (35) is an approximate result obtained taking into account that the left-hand side of Equation (33) contributes significantly to the solution of this equation only when $m_E<<1$. Under this condition $\ln \left[(1+m_E)/(1-m_E)\right]/2m_E\approx 1$ Equation (33) reduces to a cubic equation for the quantity $(1-m_E^2)$. The solution of this equation has the form (35). The approximate solution (34) differs by fractions of a percent from the exact solution obtained by means of the Legendre relations (see Figure 4b).

5. Second Order of Perturbation Theory: Account for Third Moment

In this section, we analyze what happens if, when approximating the distribution $D(m,E)$, we account for the third moment. We examine only the simplest case, supposing that all the nonzero elements of the connection matrix are equal ($J_{ij}=J$). Otherwise, the obtained expressions are too cumbersome, and it is very difficult to analyze them. We restrict ourselves only by first two terms of the series expansion (8) and set

$$\phi =\frac{1}{2}{\epsilon }^2+\frac{1}{3!}{\kappa }_3{\epsilon }^3, \mathrm{where} \epsilon =\frac{E-E_m}{{\sigma }_m}.$$

(36)

We define the coefficient ${\kappa }_3$ from the following considerations. It is necessary that the third moment of our approximation $\exp \left[-N\cdot \phi (m,E)\right]$ coincides with the third moment of the true distribution $D(m,E)$: ${\displaystyle \int {(E-E_m)}^3\exp \left[-N\cdot \phi (m,E)\right]dE={\mu }_3(m)}$. In the Appendix A, we show that for this, it is necessary that the coefficient ${\kappa }_3$ is equal to the third semi-invariant of the distribution $D(m,E)$: ${\kappa }_3=-{\mu }_3(m)/{\sigma }_m^3$. In the Appendix A, we also obtain the expression for the third moment of the distribution $D(m,E)$, which is equal to ${\mu }_3(m)=-2q{m}^2{\left(1-m^2\right)}^2$. Here $q$ is the effective number of the neighbors (see Equation (16)).

Then

$${\kappa }_3=\frac{2q{m}^2}{{\sigma }_0^3(1-m^2)},$$

and expression (6) for the function $F(m,E)$ takes the form ($H=0$):

$$F(m,E)=S(m)+\beta E+\frac{1}{2}{\epsilon }^2+\frac{1}{3!}{\kappa }_3{\epsilon }^3.$$

We recall that $E_m=E_0{m}^2$, ${\sigma }_m={\sigma }_0(1-m^2)$, $E_0=-q/2$, and ${\sigma }_0=\sqrt{q/2}$.

The system of equations that define the saddle point has the form

$$\begin{array}{l}\frac{\mathrm{\partial }F}{\mathrm{\partial }E}=\beta +\frac{\epsilon +{\kappa }_3{\epsilon }^2/2}{{\sigma }_m}=0,\\ \frac{\mathrm{\partial }F}{\mathrm{\partial }m}=\frac{1}{2}\ln \frac{1+m}{1-m}+\left(\epsilon +\frac{1}{2}{\kappa }_3{\epsilon }^2\right)\dot{\epsilon }+\frac{{\dot{\kappa }}_3}{3!}{\epsilon }^3=0,\end{array}$$

(37)

where $\dot{\epsilon }=\partial \epsilon /\partial m$. The first of the equations of (37) provides the relation

$$\epsilon +{\kappa }_3\frac{{\epsilon }^2}{2}=-\beta {\sigma }_m.$$

(38)

Moreover, by direct calculations, we obtain

$$\dot{\epsilon }=\frac{2m}{1-m^2}\left(\epsilon +\frac{\left|E_0\right|}{{\sigma }_0}\right), {\dot{\kappa }}_3=\frac{4qm}{{\sigma }_0{\sigma }_m^2},$$

and substituting these equalities in the second Equation (37), we finally have

$$\frac{1}{2m}\ln \frac{1+m}{1-m}=2\beta \left(\epsilon {\sigma }_0+\left|E_0\right|\right)-\frac{2q}{3{\sigma }_0{\sigma }_m^2}{\epsilon }^3.$$

(39)

Generally speaking, for deriving an equation relating $m$ and $\beta$ it is necessary to solve Equation (38) and to determine $\epsilon =\epsilon (m)$ and substitute it into Equation (39). Analyzing this equation, it would be possible to define the region of applicability of the m-vicinity approximation with account for the third moment and to obtain an expression for the critical temperature. It is rather difficult to solve this problem analytically. This is the reason why we restrict ourselves by analyzing the influence of the third moment ${\mu }_3(m)$ on the value of the critical temperature defined by Equation (39) when $m=0$.

When $m\to 0$, from Equation (38) it follows that $\epsilon \to -\beta {\sigma }_0$. Then the second term in the right-hand side of Equation (39) tends to $-2q{\beta }^3/3$ and this equation itself takes the form $1=2\beta \left(-\beta {\sigma }_0^2+\left|E_0\right|\right)+2q{\beta }^3/3$. By virtue of the expressions for $E_0$ and ${\sigma }_0$ we obtain the equation for the critical value of the inverse temperature ${\beta }_c$:

$$\frac{1}{q}=\beta -{\beta }^2+\frac{2k}{3}{\beta }^3, \mathrm{where} k=1.$$

(40)

If in this equation we set $k=0$, we obtain Equation (17) that corresponds to the above-discussed case of the Gaussian approximation. We introduce the parameter $k$ since in what follows, we will use it as a fitting parameter.

Solving Equation (40) for $d\in [3,7]$ numerically, we see the worse agreement of the obtained results with the computer simulations (see Table 1). Previously (see [24], Section 17.6) it was mentioned that an account for higher moments does not always lead to an increase in the quality of an approximation. You can expect such a result only for a narrow class of distributions. This question was under analysis when studying the Gram–Charlier and the Edgeworth series expansions [24,25]. In particular, it turned out that when we omit the higher order terms of the infinite series (8), the error is of the order of the first omitted term. To compensate this error, we use the fitting parameter $k$. The author of [24] (s.17.6) suggests using this receipt if the optimal value of k gives a better agreement with all the numerical values. The questions of convergence of an infinite series are of secondary importance when solving a specific problem.

With those arguments in mind, we found that an agreement with the best numerical results is better when $k<0$. In this case, the solution of Equation (40) takes the form

$${\beta }_c=B\cos \varphi -\frac{1}{2\left|k\right|},$$

(41)

where

$$B=\frac{\sqrt{1+2|k|}}{|k|}, \varphi =\frac{1}{3}\left(2\pi -\arccos R\right),R=-\frac{1}{(1+2|k|{)}^{3/2}}\left(1+3|k|+\frac{6k^2}{q}\right).$$

For d = 3, the best agreement is reached when $k=-0.8$. It turned out that this value of k is optimal for all examined dimensions: the solutions of Equation (41) with the same fitting parameter provide a very good agreement with the computational data for the lattices of all the dimensions $3\le d\le 7$ ($q=2d$). As we see from Table 1, accounting for the third moment and introduction of the fitting parameter $k=-0.8$ reduce the relative error

$$Err=100\%\cdot \left(\frac{{\beta }_c^{(\exp )}-{\beta }_c^{(theor)}}{{\beta }_c^{(\exp )}+{\beta }_c^{(theor)}}\right)$$

by 1 to 2 orders of magnitude, and it becomes comparable with the computational error. When we account for the interactions with the second and third neighbors, the introduction of the third moment also improves the agreement with the computer simulations significantly (see the next section).

6. Comparison with Computer Simulations: Three-Dimensional Ising Model

To estimate the accuracy and the correctness of the obtained equations, we performed computer simulations using the Metropolis algorithm and the algorithm of Wang and Landau [11]. We restricted ourselves to the examination of the three-dimensional Ising model supposing that for lattices of higher dimensions ($d\ge 4$), the agreement with the numerical results would be only better. In the course of our simulations, we explored the functions $f=f(\beta )$, $U=U(\beta )$, and $C=C(\beta )$. This allowed us to define the dependence of the critical temperature ${\beta }_c={\beta }_c(q)$ on $q$, to calculate the logarithmic density of states $\mathsf{\Psi }(E)$, and to analyze how the magnetization distribution changed when the inverse temperature increased.

Recently, there is a strong interest in an account for interactions with the second and the third neighbors (see, for example, [32,33]). However, we do not know any analytical estimates for the critical temperatures. Our Equation (18) is quite accurate for such problems too.

We use the Metropolis algorithm to examine the three-dimensional lattices of the linear sizes $L$ = 20, 32, and 64 with periodic boundary conditions. We look for the dependences of the critical parameters on the effective number of neighbors $q$. We suppose that all the interaction constants with the 6 nearest neighbors are equal to one ($J_{NN}=1$). We start with varying the value of the constant of interaction with the 12 next-nearest neighbors $J_{NN}$ from 0 to 1, and this means that $q$ changes from 6 to 18. Then we fix the value $J_{NNN}=1$ and vary the interaction constant with the 8 next-next-nearest neighbors $J_{NNNN}$ from 0 to 1 so that $q$ changes from 18 to 26. The initial state of the system is random. To bring the system to a state close to equilibrium for a given temperature, at the first ${10}^4\cdot L^3$ Monte Carlo steps, we do not accumulate the statistical data. We flip spins according to the Metropolis algorithm in the order of their sequence in the lattice. After each flip, we measure the values of the energy and the magnetization. The total number of the Monte Carlo steps is $4\cdot {10}^5\cdot L^3$. We vary the inverse temperature with the step size $\Delta \beta \approx 1.8\cdot {10}^{-4}$ for $L=20$, with the step size $\Delta \beta \approx 7.5\cdot {10}^{-6}$ for $L=32$, and with the step size $\Delta \beta \approx {10}^{-4}$ for $L=64$. We define the critical point as the point of maximum of the energy variance ${\sigma }_E^2={\sigma }_E^2(\beta )$. At this point we fix the values of $U_{\max }=U({\beta }_{\max })$, ${\sigma }_{\max }^2={\sigma }_E^2({\beta }_{\max })$ as well as other values. We also use the Metropolis algorithm to calculate the energy and magnetization moments near the critical temperatures. In

Table 2. Values measured at $\beta ={\beta }_{\max }$.

${\mathit{J}}_{\mathit{N}\mathit{N}\mathit{N}}$	${\mathit{J}}_{\mathit{N}\mathit{N}\mathit{N}\mathit{N}}$	$\mathit{q}$	${\mathit{\beta }}_{\mathbf{max}}$	${\mathit{\sigma }}_{\mathbf{max}}^2$	$\left|\langle \mathit{E}\rangle \right|$	$\langle {\mathit{E}}^2\rangle$	$\langle {\mathit{E}}^4\rangle$	$\left|\langle \mathit{m}\rangle \right|$	$\langle {\mathit{m}}^2\rangle$	$\langle {\mathit{m}}^4\rangle$	${\mathit{Q}}_{\mathbf{max}}$
0.0000	0	6.00	0.22208	71.02	1.0247	1.0502	1.1041	0.06448	0.04742	0.00258	0.6173
0.0042	0	6.10	0.21978	74.45	1.0270	1.0550	1.1143	0.10751	0.04862	0.00271	0.6173
0.0410	0	7.00	0.20058	92.80	1.0423	1.0866	1.1823	0.11810	0.05241	0.00310	0.6243
0.0811	0	8.00	0.18328	107.06	1.0517	1.1065	1.2260	0.14503	0.04908	0.00273	0.6223
0.1213	0	9.00	0.16888	132.16	1.0623	1.1290	1.2770	0.08342	0.04356	0.00223	0.6084
0.1623	0	10.00	0.15648	129.91	1.0834	1.1743	1.3813	0.17143	0.04511	0.00229	0.6252
0.2048	0	11.00	0.14548	154.83	1.1020	1.2151	1.4792	0.15042	0.04237	0.00206	0.6184
0.2500	0	12.00	0.13548	189.87	1.1300	1.2776	1.6360	0.10335	0.04287	0.00212	0.6154
0.2991	0	13.00	0.12608	199.76	1.1546	1.3339	1.7833	0.01969	0.03976	0.00183	0.6145
0.3543	0	14.00	0.11709	243.49	1.1972	1.4342	2.0623	0.07870	0.04215	0.00203	0.6194
0.4189	0	15.00	0.10799	287.88	1.2179	1.4844	2.2101	0.06442	0.03106	0.00123	0.5765
0.5000	0	16.00	0.09869	325.40	1.3223	1.7497	3.0701	0.08284	0.04788	0.00249	0.6373
0.6169	0	17.00	0.08769	426.66	1.4011	1.9647	3.8730	0.03163	0.03256	0.00142	0.5536
0.7106	0	17.50	0.08059	482.31	1.4871	2.2134	4.9152	0.12453	0.04236	0.00200	0.6290
0.8566	0	17.90	0.07149	586.79	1.5548	2.4198	5.8771	0.02487	0.02444	0.00080	0.5519
0.8958	0	17.95	0.06949	600.78	1.6359	2.6783	7.1979	0.03046	0.03645	0.00151	0.6216
0.9324	0	17.98	0.06769	655.33	1.6793	2.8225	7.9946	0.03795	0.03277	0.00129	0.5992
1.0000	0.0000	18.00	0.06448	736.77	1.6900	2.8587	8.2047	0.01656	0.02573	0.00087	0.5643
1	0.0063	18.10	0.06429	768.80	1.6916	2.8644	8.2385	0.06596	0.03199	0.00123	0.5993
1	0.0637	19.00	0.06248	770.88	1.7294	2.9938	8.9982	0.07366	0.03418	0.00136	0.6133
1	0.1307	20.00	0.06049	848.33	1.7475	3.0569	9.3844	0.03053	0.03480	0.00140	0.6135
1	0.2023	21.00	0.05848	904.70	1.7442	3.0458	9.3187	0.08452	0.03413	0.00135	0.6132
1	0.2806	22.00	0.05649	926.04	1.7881	3.2008	10.2901	0.18149	0.03500	0.00150	0.5918
1	0.3693	23.00	0.05439	989.95	1.8190	3.3125	11.0228	0.08570	0.04307	0.00201	0.6390
1	0.4755	24.00	0.05199	1060.91	1.7446	3.0478	9.3385	0.00448	0.02501	0.00080	0.5762
1	0.6177	25.00	0.04919	1157.27	1.7578	3.0944	9.6306	0.01029	0.02170	0.00063	0.5537
1	0.7225	25.50	0.04739	1324.75	1.8598	3.4641	12.0698	0.03932	0.03389	0.00131	0.6203
1	0.8707	25.90	0.04499	1418.32	1.8758	3.5239	12.4942	0.09043	0.02928	0.00101	0.6080
1	0.9076	25.95	0.04449	1414.33	1.9694	3.8839	15.1680	0.20027	0.04111	0.00182	0.6402
1	0.9410	25.98	0.04399	1455.97	1.9562	3.8324	14.7722	0.15120	0.03737	0.00154	0.6334
1	1.0000	26.00	0.04309	1505.02	1.8808	3.5433	12.6366	0.01037	0.02359	0.00071	0.5765

, we present the numerical data for $L=64$.

This calculation had two goals. The first was to analyze the role of the next terms of the expansion series (8). The second goal was to examine the influence of the account for the next-nearest and the next-next-nearest neighbors on the character of the dependence $C=C(\beta )$ and find if the type of the singularity at the critical point changed to a finite jump.

6.1. Dependence ${\beta }_c={\beta }_c(q)$

In Figure 5a, we show the dependence of ${\beta }_c$ on $q$. We compare the numerical values (the upper solid line) with those that follow from Equation (18) for the Gaussian approximation (the lower solid line). As we see, the solid curves have similar shapes, and when we scale the theoretical curve by a factor of $~1.06$, they practically merge. Such a coincidence cannot be accidental. From Equation (18) it follows that in the Gaussian approximation, the curve ${\beta }_c(q)$ has singularities at $q=18$ and $q=26$:

$${\frac{d{\beta }_c}{dq}|}_{ q\to 18-0}=-\frac{{\beta }_c}{6\sqrt{ 18-q}} \mathrm{and} {\frac{d{\beta }_c}{dq}|}_{ q\to 26-0}=-\frac{{\beta }_c}{3\sqrt{ 26 (26-q)}}.$$

The curve obtained by means of computer simulations shows the same singularities. In our simulations, we check this statement very accurately, changing $q$ in the vicinities of these points with a very small step size ($\Delta q~{10}^{-3}$).

Agreement of the numerical results with our theory becomes much better when we account for the third moment in Equation (36) and solve numerically the equation for the saddle point (37) (accounting for Equations (A4) and (A11) of the Appendix A). In this case, the relative error does not exceed a fraction of a percent. When we account for interactions beyond the nearest neighbors, the expression for the third moment ${\mu }_3(m)$ becomes so cumbersome that it does not allow us to obtain analytical expressions for the critical temperature analogous to Equations (40) and (41). On the other hand, we can use a simple empirical formula

$${\beta }_c=\frac{{\beta }_c^{(0)}}{1+0.207\cdot (12J_{NNN}+8J_{NNNN})},$$

(42)

where ${\beta }_c^{(0)}\approx 0.22165$ is the critical value given by formula (41), accounting for the nearest neighbors only ($J_{NNN}=J_{NNNN}=0$). Comparing the solid and dashed lines in Figure 5b, we see that expression (42) describes the computer simulations much better than formula (18)—in comparison with the Gaussian approximation, the relative error decreases by an order of magnitude; it decreases nearly by two orders of magnitude comparing with the mean field theory.

6.2. Dependence $C=C(\beta ,q)$

Our analysis shows that in the case of a three-dimensional Ising model, our account for the next-nearest and the next-next-nearest neighbors does not qualitatively change the behavior of the heat capacity at the critical point. Next, when the number of spins N increases, the peak of the curve $C=C(\beta ,q)$ only increases and the character of the singularity remains the same for all $6\le q\le 26$.

The dependences of the heat capacity on the critical temperature calculated numerically differ significantly from the ones defined by Equations (19) and (23). This means that we cannot use Equations (21)–(23) to describe the dependences $C=C(\beta ,q)$. Equation (25) that defines the finite jump of the heat capacity and works well when $d>4$ is not applicable in the three-dimensional case.

In Figure 6, we show the dependences of the maximal energy variance ${\sigma }_{\max }^2={\sigma }_{\max }^2(q)$ on the effective number of neighbors. From the figure it is obvious that, first, the energy variance at the critical point increases when $q$ increases. Second, we see that with an increase in the lattice size $L$ the differences between the computer simulations and the theoretical curve grow quickly. Nevertheless, the dependences of ${\sigma }_{\max }^2={\sigma }_{\max }^2(q)$ on $q$ repeat qualitatively the theoretical curve (compare solid and dashed lines in Figure 6).

The situation with the lattices of larger dimensions $d$ is quite the opposite. Let us discuss it here shortly. The larger the dimension of the lattice, the better the theoretical expressions (21)–(23) describe the numerical results for the heat capacity in the vicinity of the critical point. When $d\ge 5$, the agreement is good, both qualitatively and quantitatively.

In Figure 7, the dashed lines are the theoretical curves $C=C(\beta )$ (Equations (21)–(23)); the solid lines are the results of computer simulations. For $d=4$, the value of $L=80$ (see [21]), and for the dimensions $d=5,6$ and $7$, we use the data from [23]. We see that our formulas describe the numerical results pretty well. When $d$ increases, the agreement of the theory and the simulations improves notably; when $d=7$ the theoretical and the numerical curves almost merge. In the same figure, we also present the graph for $d=4$; however, to date, there is no real understanding of what happens in the four-dimensional case. We do not know for sure if there is a jump of the heat capacity or an infinite singularity. In the same figures, we show the curves obtained in the framework of the mean field theory. We see that our approximation describes the results of the numerical experiments much better.

6.3. Magnetization Distribution $P(m)$ and Logarithmic Density of States

With the aid of the Metropolis algorithm, we calculated the magnetization distribution near the critical temperature for the three-dimensional Ising model, taking into account the nearest neighbors only ($L$ = 8, 10, and 12). To increase the accuracy of the Monte Carlo method, we performed $3\cdot {10}^7\cdot L^3$ steps. To fix the inverse temperature where the bimodal distribution replaces the unimodal, we used the step size $\Delta \beta ={10}^{-3}$. For $L=12$, we present the curves of the magnetization distribution in Figure 3a. The obtained distributions allow us to calculate the Binder cumulants $Q(\beta )$ (see [10,31]). In Figure 3b, the curves $Q=Q(\beta )$ intersect at a point $\beta =0.222, Q=0.488$. This value of $\beta$ defines the critical temperature for the infinite lattice. At the curves in Figure 3b, the square markers are the points where the bimodal distributions replace the unimodal. At these points, the cumulants are approximately equal to 0.29, and this value is sufficiently close to the value predicted by Equation (32). When $L$ increases, these points shift to the point of the curves intersection that corresponds to the critical value ${\beta }_c$ in the asymptotic limit $L\to \infty$.

To obtain the density of states, we used the algorithm of Wang and Landau [11]. We performed the calculations for the cubic lattice of the size $L=40$ without parallel computing and only accounting for the interaction with the nearest neighbors. As a criterion of the flatness of the histogram of visited states, we adopted the condition that all the values of the histogram had to be larger than 80% of its mean value. When this condition was satisfied, the algorithm reduced the modification factor according to the formula $f_{i+1}^{(\mathrm{mod})}=\sqrt{f_i^{(\mathrm{mod})}}$. The simulations stopped when the modification factor became less than $f_{final}^{(\mathrm{mod})}=\exp ({10}^{-10})$.

In Figure 4a, we present the graphs of the logarithmic density of states calculated numerically as well as the ones based on the equations of Section 2. The density of states calculated using Equations (19)–(23) and the Legendre relations was in a good agreement with the numerical data. We see that the maximal obtained error, which is about 0.7%, corresponds to the critical energy $E=U_c$. The results of the approximate formula (34) are somewhat less accurate; however, in these calculations, the deviation from the numerical results is also less than 0.8% (see Figure 4b).

7. Discussion

We used the m-vicinity method to investigate the Ising model on $d$-dimensional hypercube lattices for $3\le d\le 7$. (We recall that this method is not applicable for the lattices of lower dimensions.) The m-vicinity ${\mathsf{\Omega }}_m$ consists of all the configurations of the same magnetization m. We based our method on the series expansion of the logarithm density of states in the m-vicinities (Section 2). In the main part of the paper, we discuss the interaction with the nearest neighbors. Only in Section 6, for a cubic lattice, we analyzed both the short-range and long-range interactions.

When we account only for an isotropic interaction of the nearest neighbors, the small parameter is $\left|{\epsilon }_{\max }\right|=\sqrt{d}\left(1-\sqrt{1-2/d}\right)/2$ (see Equation (10)). The value of this parameter decreases from 0.366 when $d=3$ to 0.159 when $d=7$ (see Figure 1). Not surprisingly, the agreement between the theoretical results and computer simulations becomes better when the dimension of the lattice increases.

Summing up, we would like to list the main features of the m-vicinity method.

(1) The first order of the perturbation theory is equivalent to the Gaussian approximation of the true density of states in the m-vicinities supposing that the first two moments of the density of states and its Gaussian approximation coincide.

In this case, we obtained the simple analytical expression (18) for the critical value of the inverse temperature, which described quite accurately the results of computer simulations for different lattices (see Figure 2). When $d=3$, the relative error between the theoretical estimate ${\beta }_c$ and the computer simulations is 2.39%. When $d$ increases, the relative error decreases to 0.18% for $d=7$ (see Table 1).

In the framework of the same approximation, we obtained analytical expressions (19)–(23) that define implicitly the dependence of the free energy and its derivatives on the temperature. Based on these results, we calculated the dependence of the heat capacity $C(\beta )$ on the inverse temperature $\beta$. As it follows from Figure 7, for $d\ge 5$, the obtained graphs match the curves obtained numerically [23]. The values of the critical exponents that follow from Equations (25)–(28) coincide with the known mean field results. For $d\ge 5$, Equation (25) predicts a finite jump of the heat capacity whose value is only $0.6\%$ less than the value obtained by computer simulations [23].

(2) The second order of the perturbation theory requires going beyond the Gaussian approximation and accounting for the third moment in the expansion (36). We suppose that the first three moments of the approximate and the true distributions coincide. However, this does not automatically improve the agreement of the theoretical estimates of the critical temperature and its numerical values. The reason is an irregular convergence of the Gram–Charlier and the Edgeworth series expansions [24]. The agreement improves significantly when we introduce the fitting parameter that is the same for all the dimensions $d$ (see Section 5). After that, the second order perturbation Formula (41) describes the computer simulations within a fraction of the percent (see Table 1). From Figure 5, it follows that the same is true when we also account for the interactions with the next-nearest and the next-next-nearest neighbors.

(3) The results of our analysis show that when the dimension of the Ising model $d\ge 5$, the m-vicinity method describes the properties of the system pretty well, both qualitatively and quantitatively. When $d=4$, the type of the singularity of heat capacity still remains unclear [20,22]. Consequently, the question about the applicability of the m-vicinity method in this case remains open.

In the case $d=3$, the approach discussed here is incorrect because it predicts a finite jump of the heat capacity at the critical point and the classical values of critical exponents. This result contradicts the widely accepted concept of the singularity at the critical point in the modern phase-transitions theory. Nevertheless, this method allows one to calculate the critical temperature quite accurately (see Table 1) as well as to describe its dependence on the number of the neighbors (see Figure 5 and Equations (41) and (42)). We conclude that when $d=3$ our theory provides good results for the dependences of the free energy and the logarithmic density of states but not for their derivatives. Indeed, when for $d=3$ we use Equations (34) and (35) to calculate the logarithmic density of states, the result is in a good agreement with the computer simulations data. From Figure 4 we see a notable deviation from the numerical curve ($~0.7\%$) only in a narrow vicinity of the point $E=U_c$.

(4) Finally, we would like to note that a good agreement of the theoretical results for the density of states $\mathsf{\Psi }=\mathsf{\Psi }(E)$ with the data of computer simulations allows us to use the obtained expressions as an initial approximation for the Wang-Landau algorithm. We hope that our results will allow one to speed the algorithm up and to increase its accuracy.