Calculation Method for the Three-Dimensional Ising Ferromagnet Thermodynamics within the Frames of $\rho^6$ Model

Calculation of thermodynamic functions of the three-dimensional Ising ferromagnet above and below critical temperature is performed in the approximation of sixfold basis distribution ($\rho^6$ model). Comparison with the results for the $\rho^4$ model indicates that dependence of the thermodynamic functions on the renormalization group parameter $s$ becomes weaker. The optimal interval of the renormalization group parameter values is determined.


Introduction
Significant results in the description of the system thermodynamic properties in the vicinity of the transition point have been obtained by means of the collective variables (CV) approach. The method to deriving explicit expressions for the thermodynamic and correlation functions of the three-dimensional Ising model at temperatures both above and below critical temperature T c has been suggested within this approach. The calculations are performed with a non-Gaussian measure density. The measure is represented as an exponential function of the CV, the argument of which contains, along with the quadratic term, higher powers of the variable with the corresponding interaction constants. The simplest non-Gaussian measure density is the quartic one (ρ 4 model) with the second and the fourth powers of the variable in the exponent. Then the sixfold measure goes containing the sixth power of the variable (ρ 6 model), and so forth.
The results of the theory depend on the renormalization group (RG) parameter s due to an approximation of the Ising model partition function calculation using the non-Gaussian measure densities. This dependence decreases essentially if the non-Gaussian measure density becomes more complicated. Calculations of the correlation length critical exponent ν within the ρ 2m models with m = 2, 3, 4, 5 confirm this statement [1][2][3]. It has been established that the ρ 6 model provides an adequate description of the Ising model critical behaviour, in particular, the critical exponents, at the RG parameter values in the interval 2 ≤ s ≤ 4.
Investigation of the ρ 6 model within the numerical realization of the CV method has been performed in [4]. Analytical derivation of the explicit expressions for the ρ 6 model thermodynamic functions is the subject of the present paper. The foundations for such kind of investigations have been developed in [5][6][7][8][9], where the quartic distribution was used as a basis measure.
6 are related to the coefficients of the n + 1-th layer by the recurrent relations (RR) [11][12][13]. The solutions of these relations [13] are used in the calculation of the system thermodynamic characteristics.
2 Thermodynamic functions of the ρ 6 model in the regions of critical and limit Gaussian regimes (CR and LGR) above T c It is convenient to rewrite the model partition function as [14] Z = 2 N Z CR Z LGR . (2.1) Let us consider Z CR given by It should be mentioned that in (2.2) η −1 ≡ η ′ , ξ −1 ≡ ξ ′ at n = 0. We represent the right-hand side (RHS) of (2.2) in the form of an explicit dependence on the phase layer number n in order to calculate Z CR . In the CR region, the basic h n , α n and intermediate η n , ξ n arguments are close to their values at the fixed point. Therefore, functions of these arguments can be written as power series of deviations of basic arguments from their values at the fixed point (see [15,16]). Using the obtained representations for I 0 (h n , α n ), I 0 (η n−1 , ξ n−1 ), C(η n−1 , ξ n−1 ), we determine from (2.2) the partial free energy corresponding to the n-th phase layer: Considering (2.4), we rewrite the partial energy of the n-th phase layer as The quantities u (0) , w (0) il were determined in [13]. Let us note that in the expressions for h n and α n we can neglect a qualitatively new term proportional to E n 3 which arises within the ρ 6 model considered (E 3 is not essential as compared to E 1 or E 2 ).
Summing up expressions for F n (2.6) over layers of the phase space from n = 0 to n = m τ , we obtain To calculate these expressions, we have used the formulas: where m 0 ,c 1 , c 20 , c 30 and f 0 , ϕ 0 were defined in [14]. The dependence of m 0 on s for the ρ 6 model under consideration is plotted in figure 1 (solid curve).
Here the dashed curve corresponds to the ρ 4 model. It is easy to notice that, when s = s * (s * is the value at which h n turns to zero at the fixed point; s * =3.5862 for the ρ 4 model and s * =2.7349 for the ρ 6 model), the values of m 0 for ρ 4 and ρ 6 models coincide. Further, we put F (2) CR = 0 at τ ≪ 1, since E 2 < 1, E 3 < 1, and m τ is large. Taking E 2 , E 3 into account gives rise to terms characterising corrections to scaling. As a result, the free energy of the CR region takes the form: Note that γ 0 , δ 0 are the functions of temperature, since they are expressed in terms ofc 1 , c 20 , c 30 and Q(d). Let us extract the temperature dependence in these quantities.
Near T c we have forc 1 for c 20 and c 30 we get, respectively, The values of β cΦ (0), the correlation length critical exponent ν = ln s/ ln E 1 , exponents of the corrections to scaling ∆ 1 = − ln E 2 / ln E 1 , ∆ 2 = − ln E 3 / ln E 1 , and coefficients of the expressions forc 1 , c 20 , c 30 are given in tables 1,2. In the present paper, the numerical calculations are performed at b/c = 1 and arithmetically averaged Fourier transform of the potential. Table 1: 1 for different s.
lj occurring in (2.12) -(2.14) are presented in [13,14]. The coefficient γ 0 can be written as For δ 0 we obtain Coefficients of the expressions for γ 0 (2.15), δ 0 (2.16) at different values of the RG parameter s are given in table 3. Table 3: Values of coefficients in equations for γ 0 (2.15) and δ 0 (2.16). Hence, the free energy of the CR region reads The numerical values of coefficients γ are given in table 4. Knowledge of F CR allows one to calculate other thermodynamic functions of the system in the CR region at T > T c . For the entropy S CR , internal energy U CR and specific heat C CR we get , γ 1 , γ 2 and quantities c ν , γ + , γ − contained in γ (CR) ± 3 for different values of the parameter s.
In the region of LGR, the expression for the part Z LGR of the partition function (2.1) reads To calculate Z LGR it is convenient to select two regions of the wave vector values [14]. The first, transition region, corresponds to the values of k close to B mτ ; the second, Gaussian region, corresponds to small values of the wave vector (k → 0). Hence, we have LGR Z LGR . (2.20) The contribution to free energy from the phase space layers following the point of exit from the CR region is [14]. The plots ofm  Introducing an infinitely small external magnetic field h = µ B H (µ B is the Bohr magneton), we can write the part of the free energy corresponding to Z (2) LGR as Here, A general expression describing the contribution of long-wave fluctuations to the free energy (LGR region) reads The values off T R ,f ′ ,γ + 4 are given in table 5. The entropy, internal energy, specific heat corresponding to the LGR region are defined by the relations 3 Contributions to the thermodynamic functions of the model from the critical and inverse Gaussian regime (CR and IGR) regions below T c The free energy at T < T c can be written as [6,10] where F 0 = −kT N ln 2 is the free energy of the system of N non-interacting spins, F CR is the contribution to the free energy from the short-wave fluctuation phases of the spin moment density (CR region), and F IGR is the contribution from the long-wave phases of the fluctuations (IGR region). The number µ τ of the CV phase space layer, separating the short-wave and long-wave phases of fluctuations, is an important characteristic of the system. It is determined from the equation Here r (0) corresponds to the fixed point of the RR, r µτ +1 is determined from the solutions of the RR equations (see, for example, [13,14]), δ is a constant (δ ≤ 1). In the present paper we put δ = 1 (see [8]). Let us write an the solution of which is We need to sum the partial free energies over the layers of the CV phase space to calculate F CR . Extracting an explicit dependence on the layer number, using relations (3.3) and we obtain The value of γ − is given in table 4. The entropy, internal energy and specific heat of the system corresponding to the CR region read Let us calculate now the contribution to the free energy from the IGR region Consider the first term on the RHS of (3.9). Making use of the relations The first term on the RHS of (3.9) is equal to To find the second term on the RHS of (3.9), we need to calculate Z µτ +1 . The coefficients occurring in it equal d µτ +1 (k) = r µτ +1 s −2(µτ +1) +qk 2 ,q = 2βΦ(0)b 2 , a (µτ +1) 4 = u µτ +1 s −4(µτ +1) , a (µτ +1) 6 = w µτ +1 s −6(µτ +1) , r µτ +1 = −r µτ +1 βΦ(0),r µτ +1 = 2f 0 , (3.13)

Let us perform the change of variables
in the expression for Z µτ +1 in order to extract the free energy related to the ordering that has appeared in the system. Here <σ > is being determined from the extremum condition [17,10] <σ > 2 = 10 a Simultaneously, we include in the treatment a constant external field h = µ B H, which sustains the separated average moment, and separate from the sums over k the terms with k = 0. We obtain The prime on the sum over k means that k = 0, We have for the integrand coefficients of (3.16) The quantities p l are being determined by the equations The dependence of the quantity B 6 = 4r µτ +1 − 10 w µτ +1 (−1 + b 2 ) on the parameter s is plotted in figure 3. The dashed line shows the dependence of the analogous quantity 4f 0 on s for the ρ 4 model. One can see from figure Expanding exp{p 0 + p 1 ρ 0 + p 2 ρ 2 0 + p 3 ρ 3 0 + p 4 ρ 4 0 } in series and restricting ourselves to the terms of second order, we integrate (3.16) over ρ k with k = 0 using the Gaussian basis distribution. Gathering up the series over the averages with respect to the Gaussian distribution in the exponential, we obtain [15] ))],

21)
Expressions for I m occurring in (3.21) can be found via the transition to the spherical Brillouin zone and integration over k ∈ (0, B µτ +1 ] [17,18].
The quantity I 2 can be represented in the form: where The other I l (l = 3, 4, 5, 6) are determined by analogous relations The values of α m (m = 1, 2, 3, 4, 5, 6) are given in table 6.
which cancels terms proportional to odd powers of ρ 0 in the exponent of the integrand. We obtain The quantitiesB,Ḡ,D are given in table 7. Using the method of the steepest descent for Z µτ +1 (3.29), we find Hereρ is the extremum point of the expression E 0 (ρ) = Dρ 6 + Gρ 4 −Bρ 2 − βhρ (3.32) arising in the exponential of the integrand of (3.29) at the substitution Having (3.12) and (3.31), we can calculate the contribution to the system free energy at T < T c from the long-wave phases of the spin moment density fluctuations (3.9): The quantity γ µτ 3 determines the free energy after the exit from the CR region, and γ <σ> integrand of (3.29) has an extremum. Performing the substitution (3.33) in that integrand, we obtain where E 0 (ρ) is given in (3.32). Owing to the factor N in the exponent in (4.1), the integrand has a sharp maximum at the pointρ corresponding to the equilibrium value of the order parameter. The value ofρ can be found from the extremum condition ∂E 0 (ρ) ∂ρ = 0 or In the case h = 0 we obtain a biquadratic equation, which is reduced by means of substitutionρ 2 = y Extracting the temperature dependence, we obtain the equation for the mean spin moment < σ >=ρ = √ y: The value of <σ > (0) is given in table 7.
The susceptibility per particle χ can be found from equation (4.2) by differentiating it with respect to H and using the relation χ = µ B ∂<σ> ∂H : Separating the temperature dependence in the coefficients D, G,B (see (3.30)), we obtain a final expression for the susceptibility.

Thermodynamics of the system in the vicinity of the phase transition point
Having calculated the contributions to the system free energy from the shortwave and long-wave modes of the spin moment density oscillations both above and below T c , we can compute the total free energy the entropy, internal energy and specific heat. We obtain the total free energy of the system at h = 0 taking (2.17), (2.23) and (3.6), (3.35) into account: where (see table 9) The coefficients γ ± 3 can be written as a product of the universal partγ ± 3 (ta- ) and non-universal factor c 3 ν depending on the microscopic parameters of the Hamiltonian (the lattice constant c, the effective interaction radius b and the valueΦ(0) of the interaction potential Fourier transform at k = 0): We have for the entropy, internal energy and specific heat of the system with the coefficients given by the relations The formula for the specific heat can also be rewritten as The plus and minus signs correspond to T > T c and T < T c , respectively. The system susceptibility per particle at infinitely small values of the external field H at T > T c (χ = − 1 N ∂ 2 F LGR ∂H 2 , see (2.23)) and T < T c (see (4.6)) is given by Here (see table 7), Plots of the temperature dependence of the system free energy F/N , entropy S/kN , specific heat C/kN , mean spin moment < σ > (4.5), susceptibility χ (5.7) at s = 2, 2.5, 3 are shown in figures 4-8.
The corresponding plots of the thermodynamic characteristics calculated within the ρ 4 model, with the confluent corrections [9] being taken into account, are also given here (dashed curves). Comparison of these plots shows that the dependence of the thermodynamic functions on the parameter s for the ρ 6 model is weaker than for the ρ 4 model. The dependence of F/N on s for these two models is represented in figure 9.  Let us note that the calculations performed are best suited for the intermediate values of s, close to the quantity s * , at which h n turns to zero at the fixed point (s * =2.7349 for the ρ 6 model and s * =3.5862 for the ρ 4 model). The use of the difference form of the RR based on a non-Gaussian measure density works especially well for this region of s. At small values of s (s → 1), some complications arise when the unit element is extracted. In this limit, the RR should be represented as the perturbation series with respect to the Gaussian distribution (h n is large at s → 1, and expansions in h −2 n , α n can be used [19,20]). There also exists an upper limit for s. At large s, one must take into account the correction due to the potential averaging [10,21], which increases with s.

Conclusions
The method for the calculation of the three-dimensional Ising model thermodynamics is developed in the sixfold distribution approximation. Both temperature regions above and below the critical value of T c are considered. The main distinguishing feature of the approach is the separate allowance for the contributions from the short-and long-wave fluctuation phases of the spin moment density to the free energy of the system near T c . Within the framework of the ρ 6 model, we obtained the explicit expressions for the critical amplitudes of the thermodynamic functions of the three-dimensional Ising ferromagnet and calculated the coefficients of the free energy, the universal characteristics (the critical exponents, the ratios of the critical amplitudes). Calculation of the free energy, entropy, specific heat, mean spin moment, susceptibility was performed for different values of s. Comparison of the results obtained within ρ 4 and ρ 6 models indicates that the dependence of thermodynamic functions on the RG parameter s is weaker for the ρ 6 model.
The values of s close to s * are optimal for the method presented. Obtaining analytical expressions for critical amplitudes and system thermodynamic Table 10: Values of the critical exponents, ratios of the critical amplitudes and their combinations for s = s * (s * =2.7349 for the ρ 6 model and s * =3.5862 for the ρ 4 model) obtained by means of the CV method. Data calculated within the field theory approach (FTA) [22][23][24][25][26][27] and high-temperature series (HTS) [28][29][30][31][32].

Model
Ref characteristics as functions of the Hamiltonian microscopic parameters is the advantage of the method developed. The leading critical amplitudes for the specific heat and other thermodynamic characteristics are represented as a product of the universal part, independent of microscopic parameters, and the non-universal factor, which depends on these parameters.