Analytical Expressions for Ising Models on High Dimensional Lattices

We use an m-vicinity method to examine Ising models on hypercube lattices of high dimensions d≥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.


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 ≥ 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 ε max . In the case of the planar lattice (d = 2), ε max ≈ 0.7, and this value is not sufficiently small. When the lattice dimension d increases, the value of ε 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 ≥ 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.

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 = {±1}, i = 1, 2, . . . N situated at the nods of a hypercubic lattice. In what follows, we assume the periodic boundary conditions. The Hamiltonian of the system is where β is the inverse temperature, the summation is carried out over all the values of the energy E, and D(E) 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 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 Ω m , which will be called the m-vicinities. The m-vicinity Ω 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, . . . , 1) by opposite signs of n spins, where n = N(1 − m)/2 and the number of such configurations is equal to N n .
The density D(E, m) is the energy distribution for a given m-vicinity, and the partition function of the system is 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 σ 2 m (see Appendix A or [26][27][28]). In the case of the Ising model on the hypercube, these expressions in the limit N → ∞ are sufficiently simple, 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 Ω 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 σ 2 m : It is evident that the 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 σ m → 0 , we can replace the exponent in Equation (4) by a delta function δ(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 σ 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 σ m = 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 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)).

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 where ϕ = ϕ(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 where Let us estimate the integral (5) using the saddle point method. The equations for the saddle point are ∂F ∂m 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 (β) = 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 ϕ(m, E) as a perturbation theory series in the vicinity of point E = E m : The quantities κ 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 |ε| << 1 and in this case the m-vicinity method is sufficiently accurate. Let us clarify that we are not interested in the values of ε = ε(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 The smallness of the parameter ε 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 ε 0 by means of the reverse-reasoning method. Suppose the parameter ε 0 is small and it is sufficient to use only the first term of the series (8) and set ϕ = ε 2 /2. Then the second of the Equation (7) takes a simple form ε = −βσ m and we can rewrite the first of the Equation (7) as follows: We are interested in the behavior of the quantity ε = ε 0 , where ε 0 is the root of Equation (9). The value of |ε 0 | reaches its maximum at the critical point β = β c (see [30]). In the limit m → 0 ( β → β c ), Equation (9) takes the form The solution of this equation is and max|ε 0 | = |ε max |. Consequently, the m-vicinity method is applicable when the small parameter |ε 0 | ≤ |ε max |. According to Equation (10), it is necessary that A detailed analysis shows (see [27]) that a stricter inequality E 2 0 /σ 2 0 ≥ 8/3 or |ε max | ≤ 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 2 0 /σ 2 0 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 ε max = √ q −1 + 1 − 4/q /2 √ 2. From this equation, it follows that when q >> 1 the value of ε max ∼ 1/ 2q . Consequently, when q increases, the value of the small parameter |ε max | decreases quite rapidly and the accuracy of the m-vicinity method increases.
In Figure 1, we show the dependence of |ε max | on d when we account for the nearest neighbors only (in this case E 2 0 /σ 2 0 = d). For a planar lattice, |ε max | ≈ 0.7. Obviously, it is problematic to use it as a small parameter. Moreover, for a planar lattice, the inequality E 2 0 /σ 2 0 ≥ 8/3 is not satisfied. Next, for a cubic lattice, the last inequality is satisfied and the value of |ε max | ≈ 0.36 is more appropriate as a small parameter. Again, when d increases, the value of |ε max | 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 |ε max |.  Summing up, we can say that the Gaussian approximation is the first order perturbation theory in the small parameter max ε .

First Order of Perturbation Theory: Gaussian Approximation
In this section, we use the Gaussian approximation (4) for the distribution D (6)

If in Equations
Eliminating the variable ε from these equations, we obtain an equation of st This equation differs from the well-known equation of Bragg and Williams [2 term proportional to 2 β . After a transformation of Equation (11), with the acco Equations (12) and (13), we obtain the expression for the free energy: quently, in this case, we can use the m -vicinity method only for the spatial dime

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 ϕ = ε 2 /2, Equation (6) becomes the equations for the saddle point (7) are Eliminating the variable ε from these equations, we obtain an equation of state This equation differs from the well-known equation of Bragg and Williams [29] by a term proportional to β 2 . After a transformation of Equation (11), with the account for Equations (12) and (13), we obtain the expression for the free energy: that defines f as a function of β and m. In Equation (14), the spontaneous magnetization m = m(β; H) is the solution of Equation (13). Again, by setting σ 2 0 = 0, we recover the known result from the mean field theory.
We would like to point out that the equality E 2 0 /σ 2 0 = d is true only when we account for the interaction with the nearest neighbors. Then the inequality E 2 0 /σ 2 0 ≥ 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 ≥ 3.

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 where we introduce dimensionless characteristics 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 → 0 ( β → β c ), Equation (15) takes the form Now, with the account of Equation (16), we obtain the expression for the critical temperature 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, β 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).
where we introduce dimensionless characteristics The defined variables (16) are convenient since in Equation (15), the only para that depends on the type of the lattice is q .
Now, with the account of Equation (16), we obtain the expression for the critica perature The critical value of c b depends only on the number q , and this paramet pends on the mean and the variance of the elements of the connection matrix. Since a characteristic of the interactions that we take into account, we regard it as an ef number of the neighbors. In particular, if we account for an isotropic interaction w nearest neighbors only, с , and d is the dimension of the lattice. No in this case, q is exactly equal to the number of spins with which the given spin int Then Equation (18) describes pretty well the results of computer simulations for dimensions, which we examined (see Figure 2 and Table 1.)

Analytical Expressions
Let us list the basic thermodynamic characteristics that we obtained from Equations (11)-(14).
(1) The interval β < β c . When β < β c and m = 0, from the above-mentioned equations, we obtain that the free energy, the internal energy, and the heat capacity are respectively.
(2) The interval β ≥ β c . To obtain the analytical expressions for the values f , U, and C, we solve Equation (15) and transform Expressions (12) and (14) to the forms Here, σ 2 E = −d 2 f /dβ 2 is the energy variance that is related to the heat capacity via C = β 2 σ 2 E . The expression (23) is a result of differentiating the second of Equation (12) with respect to variable β.
If we regard the magnetization as an independent variable m ∈ [0, 1] and consider the values of f , U, and σ 2 E as functions of m, then Equations (20)-(23) define implicitly the dependence of f , U, and σ 2 E on the inverse temperature β = b|E 0 |/σ 2 0 at the interval β ≥ β c .

Critical Parameters
To examine the behavior of the thermodynamic characteristics near the critical point, we introduce a relative inverse temperature Omitting intermediate calculations, we only present the most important critical dependences.
(1) From Equations (19)- (23), it follows that at β = β c the free energy and the internal energy are continuous functions, and the heat capacity has a jump. Indeed, when β > β c Equation (23) holds, and in the limit m → 0 ( β → β c ) we obtain where b c is defined by Equation (18). Comparing this expression with Equation (19), we see that at β = β c the energy variance has a jump and consequently the heat capacity also has a jump: 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: We see that Equation (25) provides very good agreement with the computer simulations; the larger d is, the better this agreement.
(2) When β > β c , the value of the spontaneous magnetization near the critical point ( t → 0 ) obtained from Equation (15) is This expression differs from the dependence M ∼ t 1 8 that is valid for the twodimensional Ising model (q = 4); however, it qualitatively coincides with the expressions M = 2 √ t and M = √ 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 > √ 3 for any q ≥ 4. If q ≤ 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 Comparing with the analogous expression of the mean field theory, we see that an extra factor 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 β c H = m|m| δ−1 h s (tm −2 ) with the critical exponent δ = 3 and the scaling function

Magnetization Distribution
The integral 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 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 − m 4 /3 m 2 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(β c ) at the critical point. To calculate Q c we use the Taylor expansion of the function Φ(m) in Equation (30): where This expansion is useful when β ≤ β c + O N −1/2 and the value of P(m) is noticeably nonzero only when m << 1. In this case, the distribution (30) takes the form which allows us to calculate the normalization constant P 0 easily as well as the mean values m 2 and m 4 . It is not difficult to see that when |β − β c | > O N −1/2 , we can consider the distribution P(m) as purely Gaussian in full agreement with the results of [31]. However, when β → β c we have a 2 → 0 and at the critical point, the distribution takes the form P(m) = P 0 exp −Na 4 m 4 /4! . Using this expression to calculate the critical value of the Binder cumulant [31] in the limit N → ∞ we obtain takes the form Using this expression to calculate the critical value of the Binder cumulant [31] in the limit N →∞ we obtain

Density of States
In Section 4.2, we derived Equations (19)- (23), which allowed us to obtain implicitly the logarithmic density of states It is interesting that in the limit N → ∞ 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 Φ(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 Φ(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 Φ(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 −Na 4 m 4 /4! 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).

Density of States
In Section 4.2, we derived Equations (19)- (23), which allowed us to obtain implicitly the logarithmic density of states using the Legendre relations. When β changes from 0 to ∞, the value of E changes from 0 to E 0 and for each β we obtain a pair of values E and Ψ(E). In such a way, we generate the function Ψ(E), which we suppose to be symmetric: Ψ(−E) = Ψ(E). In Figure 4a, we present the comparison of our results with computer simulations. , which we suppose to be symmetric: ( ) ( ) E E Ψ − = Ψ . In Figure 4a, we present the comparison of our results with computer simulations. Let us determine an explicit form of the dependence defines the density of states or, in other words, the number of states with the energy E .
From this equation, it follows that 2 0 2 0 where c b is defined by Equation (18) and From this equation, it follows that where b c is defined by Equation (18) and

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 We define the coefficient κ 3 from the following considerations. It is necessary that the third moment of our approximation exp[−N · ϕ(m, E)] coincides with the third moment of the true distribution D(m, E): In the Appendix A, we show that for this, it is necessary that the coefficient κ 3 is equal to the third semi-invariant of the distribution D(m, E): κ 3 = −µ 3 (m)/σ 3 m . In the Appendix A, we also obtain the expression for the third moment of the distribution D(m, E), which is equal to µ 3 (m) = −2qm 2 1 − m 2 2 . Here q is the effective number of the neighbors (see Equation (16)). Then , and expression (6) for the function F(m, E) takes the form (H = 0): We recall that E m = E 0 m 2 , σ m = σ 0 (1 − m 2 ), E 0 = −q/2, and σ 0 = q/2. The system of equations that define the saddle point has the form where . ε = ∂ε/∂m. The first of the equations of (37) provides the relation Moreover, by direct calculations, we obtain and substituting these equalities in the second Equation (37), we finally have Generally speaking, for deriving an equation relating m and β it is necessary to solve Equation (38) and to determine ε = ε(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 µ 3 (m) on the value of the critical temperature defined by Equation (39) when m = 0.
When m → 0 , from Equation (38) it follows that ε → −βσ 0 . Then the second term in the right-hand side of Equation (39) tends to −2qβ 3 /3 and this equation itself takes the form 1 = 2β −βσ 2 0 + |E 0 | + 2qβ 3 /3. By virtue of the expressions for E 0 and σ 0 we obtain the equation for the critical value of the inverse temperature β c : 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 ∈ [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 where 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 ≤ d ≤ 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 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).

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 ≥ 4), the agreement with the numerical results would be only better. In the course of our simulations, we explored the functions f = f (β), U = U(β), and C = C(β). This allowed us to define the dependence of the critical temperature β c = β c (q) on q, to calculate the logarithmic density of states Ψ(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 · 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 · 10 5 · L 3 . We vary the inverse temperature with the step size ∆β ≈ 1.8 · 10 −4 for L = 20, with the step size ∆β ≈ 7.5 · 10 −6 for L = 32, and with the step size ∆β ≈ 10 −4 for L = 64. We define the critical point as the point of maximum of the energy variance σ 2 E = σ 2 E (β). At this point we fix the values of U max = U(β max ), σ 2 max = σ 2 E (β 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, 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(β) and find if the type of the singularity at the critical point changed to a finite jump.

Dependence β c = β c (q)
In Figure 5a, we show the dependence of β 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 β c (q) has singularities at q = 18 and q = 26: 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 ( ∆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 µ 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 where β (0) c ≈ 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.
(0) 1 0.207 (12 8 ) where ( 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.

Dependence
( , ) C C 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  (42)); lower solid curve corresponds to Gaussian approximation (see Equation (18)); dotted line is β c that follows from mean field theory. (b) Relative error: lower dashed line corresponds to account for third moment (42); solid line corresponds to Gaussian approximation (18); upper dotted line corresponds to mean field theory.

Dependence C = C(β, 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(β, q) only increases and the character of the singularity remains the same for all 6 ≤ q ≤ 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(β, 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 σ 2 max = σ 2 max (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 σ 2 max = σ 2 max (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 ≥ 5, the agreement is good, both qualitatively and quantitatively.
In Figure 6, we show the dependences of the maximal energy variance 2 2 max max ( ) 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 2 2 max max ( ) q σ σ = on q repeat qualitatively the theoretical curve (compare solid and dashed lines in Figure 6). (see [21]), and for the dimensions 5, 6 d = 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 7 d = the theoretical and the numerical curves almost merge. In the same figure, we also present the graph for 4 d = ; however, to date, there is no real understanding of what happens in the fourdimensional 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. In Figure 7, the dashed lines are the theoretical curves C = C(β) (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. are from [21]; curves for d = 5, 6, and 7 are numerical results [23]. Dotted line in figures for d = 5, 6, 7 represents mean field theory.

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 Figure 7. Heat capacity C vs. β for d = 4, 5, 6, 7. Computer simulations-solid lines; theory (Equations (19) and (23))-dashed lines. Data for d = 4 and L = 80 are from [21]; curves for d = 5, 6, and 7 are numerical results [23]. Dotted line in figures for d = 5, 6, 7 represents mean field theory.

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 · 10 7 · L 3 steps. To fix the inverse temperature where the bimodal distribution replaces the unimodal, we used the step size ∆β = 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(β) (see [10,31]). In Figure 3b, the curves Q = Q(β) intersect at a point β = 0.222, Q = 0.488. This value of β 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 β c in the asymptotic limit L → ∞ .
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 . The simulations stopped when the modification factor became less than f (mod) f inal = 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).

Discussion
We used the m-vicinity method to investigate the Ising model on d-dimensional hypercube lattices for 3 ≤ d ≤ 7. (We recall that this method is not applicable for the lattices of lower dimensions.) The m-vicinity Ω 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 |ε max | = (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 β 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(β) on the inverse temperature β. As it follows from Figure 7, for d ≥ 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 ≥ 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 nextnearest and the next-next-nearest neighbors.
(3) The results of our analysis show that when the dimension of the Ising model d ≥ 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 Ψ = Ψ(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.

Appendix A. Calculation of the First Moments of Energy Distribution
Previously, we described in detail how to calculate the moments of the distribution of the energy states belonging to the vicinity Ω m [26][27][28]. In these papers, we obtained the exact expressions for the first two moments of the energy distribution that are the mean E m and the variance σ 2 m . For a finite system and a general connection matrix J = J ij N i,j=1 , we presented these expressions in finitely combinatorial forms. In the case of the Ising connection matrix for a d-dimensional hypercube, we, in addition, found E m and σ 2 m in an asymptotic limit when the number of spins N tends to infinity. Next, we show how to obtain the combinatorial and the asymptotic formulas for the third moment µ 3 .

Appendix A.1. Notations
In what follows, for the sake of simplicity when averaging over the m-vicinities Ω m , we sometimes will omit the subscript m. For example, if a state s = (s 1 , s 2 , . . . , s N ) belongs to an m-vicinity Ω m and E(s) is the energy of the state s, we will write the first three moments as We are interested in calculating the second and third central moments-the variance σ 2 and the semi-invariant µ 3 : Usually, the energy E(s) = a · ∑ N i,j=1 J ij s i s j of the state s includes a normalization coefficient a. In what follows, we for simplicity omit this coefficient and write the energy as E(s) = sJs + ≡ N ∑ i,j=1 J ij s i s j . (A1) After obtaining the final expressions for the moments, we have to multiply E , σ 2 , and µ 3 by a, a 2 , and a 3 , respectively.

Appendix A.2. Connection Matrix of General Form
Suppose that the connection matrix has a general form and the configuration s 0 = (s 01 , s 02 , . . . , s 0N ) is the ground state whose energy is E 0 = s 0 Js + 0 . To calculate the energy moments E r , r = 1, 2, . . ., it is necessary to present E r as a sum of terms which contain only the products of spin variables with non-repeating indices. For example, we can write the expressions for E and E 2 as where χ ij and χ ijkr are tensors, whose components are nonzero and equal to one only if they do not contain the repeating indices.
The same as when calculating < E > in Equation (A2), we can obtain the second and the third energy moments and show that the averaging over Ω m leads to the following equalities: where h 0i = ∑ N j=1 J ij s 0j is a local field acting on the i-th spin belonging to state s 0 . The matrix J (3) in the expression for the third moment is J (3) = J 3 ij . The expressions (A3) and (A4) provide us with the most general formulas for the first, second, and third moments of the distribution of energies of the states from the mvicinity Ω m in terms of the elements of the connection matrix J. (If it is necessary to use the normalized expression for the energy E(s), see the note after Equation (A1).)

Appendix A.3. Ising Model on Hypercube
In this subsection, we analyze the Ising model with the short-range interaction, when only the nearest spins interact, and other elements of the connection matrix are equal to zero. Let the nonzero matrix elements be equal to a constant J. We factor out this constant in all our calculations. Then, there are q = 2d ones in each row of the matrix J, all the other matrix elements are equal to zero, and the ground state is a configuration s 0 = (1, 1, 1, . . . , 1) ∈ R N .
It is easy to see that the following equations E 0 = qN, Tr J 2 = qN, s 0 J 2 s + 0 = q 2 N, Tr J 3 = 0, hold. Then the expressions (A4) take the form To obtain asymptotic expressions for the moments ( N → ∞ ), at first it is necessary to obtain the asymptotic expressions for the coefficients K 2 , K 4 , and K 6 , expanding them with respect to the small parameter δ = 1/N: here, m = (N − 2n)/N ∈ [−1, +1] is the magnetization of the state s ∈ Ω m .