Free Energy of Three-Dimensional Uniaxial Magnet in the Higher Non-Gaussian Approximation and in the Presence of an External Field

A three-dimensional Ising-like system in a homogeneous external field is studied on the basis of the higher non-Gaussian measure density (the ρ6 model). The presented solutions of recurrence relations for the coefficients of the effective measure densities and the generalized point of exit of the system from the critical regime are used for calculating the free energy of the system at temperatures T > Tc (Tc is the phase transition temperature in the absence of an external field). A calculation technique is based on the first principles of statistical physics and is naturally realized without any general assumptions and without any adjustable parameters. The obtained expression for the free energy does not involve series expansions in the scaling variable and is valid near the critical point not only in the regions of the so-called weak and strong external fields, but also in the crossover region between these fields, where power series in the scaling variable are not effective.


Introduction
The Ising model is one of the most studied models in the theory of the phase transitions, not only because it is considered as the prototype of statistical systems showing a nontrivial power-law critical behaviour, but also because it describes several physical systems [1].Many systems characterized by short-range interactions and a scalar order parameter undergo a transition belonging to the Ising universality class.Despite the simplicity of the Ising model and great success in the research by means of various methods (see, e.g., [1]), the problem of analytic description of the Ising-like magnets in the three-dimensional (3D) space near the critical point is still unsolved exactly even in the case of the absence of an external magnetic field.The presence of a magnetic field complicates this problem.
The description of the phase transitions in the 3D magnets, usually, is associated with the absence of exact solutions and with many approximate approaches for obtaining different system characteristics.In this paper, the behaviour of a 3D Ising-like system near the critical point in a homogeneous external field is studied using the collective variables (CVs) method [2][3][4].The main peculiarity of this method is the integration of short-wave spin-density oscillation modes, which is generally done without using perturbation theory.The CV method is similar to the Wilson nonperturbative renormalization group (RG) approach (integration on fast modes and construction of an effective theory for slow modes) [5][6][7].The term collective variables is a common name for a special class of variables that are specific for each individual physical system [2,3].The CV set contains variables associated with order parameters.Because of this, the phase space of CV is most natural for describing a phase transition.For magnetic systems, the CV ρ k are the variables associated with modes of spin-moment density oscillations, while the order parameter is related to the variable ρ 0 , in which the subscript "0" corresponds to the peak of the Fourier transform of the interaction potential.
The free energy of a 3D Ising-like system in an external field at temperatures above T c is calculated using the non-Gaussian spin-density fluctuations, namely, the sextic measure density.The latter is represented as an exponential function of the CV whose argument includes the powers with the corresponding coupling constants up to the sixth power of the variable (the ρ 6 model).
The present paper supplements the earlier works [8][9][10][11], in which the ρ 6 model was used for calculating the free energy and other thermodynamic functions of the system in the absence of an external field.The ρ 6 model provides a better quantitative description of the critical behaviour of a 3D Ising-like magnet than the ρ 4 model [10].For each of the ρ 2m models, there exists a preferred value of the RG parameter s = s * (s * = 3.5862 for the ρ 4 model, s * = 2.7349 for the ρ 6 model, s * = 2.6511 for the ρ 8 model, and s * = 2.6108 for the ρ 10 model) nullifying the average value of the coefficient in the term with the second power in the effective density of measure at the fixed point.The values of s close to s * are optimal for the given method of calculations.The difference form of the recurrence relations (RRs) between the coefficients of effective non-Gaussian densities of measures operates successfully just in this region of s.It was established (see, e.g., [10,11]) that as the form of the density of measure becomes more complicated, the dependence of the critical exponent of the correlation length ν on the RG parameter s becomes weaker gradually, and, starting from the sextic density of measure, the value of the exponent ν, having a tendency to saturate with increasing m (which characterizes the order of the ρ 2m model, m = 2, 3, 4, 5), changes insignificantly (see Figures 1 and 2).The point s ≈ s * in Figure 1 corresponds to the beginning of the ν(s) curve stabilization for each of the ρ 2m models.The value of the exponent ν in Figure 2 is calculated for s = s * .The ρ 2 model (Gaussian approximation) leads to the classical value ν = 0.500.In the case when s = s * , we have ν = 0.605 for the ρ 4 model and ν = 0.637 for the ρ 6 model.The value of the critical exponent ν for the ρ 6 model agrees more closely with the other authors' data for the 3D Ising model than the estimate in the ρ 4 model approximation, for example, with the values determined using the fixed-dimension perturbative RG (ν = 0.6304(13) [12]), high-temperature series (ν = 0.63002 (23) [13]), and Monte Carlo simulations (ν = 0.6296 (7) [14]).The Ising model corresponds to the ρ 2m model approximation, where the order of the model 2m ≥ 4. The ρ 4 model allows us to go beyond the classical analysis and to describe all qualitative aspects of the second-order phase transition.As is seen from Figures 1 and 2, the critical behaviour of a 3D Ising-like system within the CV method can be described quantitatively at 2m ≥ 6, and, in particular, at 2m = 6.It was shown in [10,15] that the graphs of the temperature dependences of the order parameter (the spontaneous magnetization) and specific heat for the ρ 6 model agree more closely with the Liu and Fisher's results [16] than the corresponding plots for the ρ 4 model.The correctness of the choice of the ρ 6 model for investigations is also confirmed in [17,18], where Tsypin proved that the term with the sixth power of the variable in the effective potential plays an important role.
The methods existing at present make it possible to calculate universal quantities to a quite high degree of accuracy (see, e.g., [1]).The advantage of the CV method lies in the possibility of obtaining and analysing thermodynamic characteristics as functions of the microscopic parameters of the initial system [10,11,15,19,20].The results of calculations for a 3D Ising system on the basis of the ρ 4 and ρ 6 models are in accord with the results obtained by other authors.Comparison of the critical exponents and universal ratios of critical amplitudes with the data calculated within the field-theory approach [21][22][23] and high-temperature expansions [24][25][26][27][28] can be found in our articles [10,15].In [29], the scaling functions of the order parameter and susceptibility, calculated on the basis of the free energy for the ρ 4 model, were graphically compared with other authors' data.Our results accord with the results obtained within the framework of the parametric representation of the equation of state [30] and Monte Carlo simulations [31].
The expressions for the thermodynamic characteristics of the system in the presence of an external field have already been obtained on the basis of the simplest non-Gaussian measure density (the ρ 4 model) in [32][33][34][35] using the point of exit of the system from the critical regime as a function of the temperature (the weak-field region) or of the field (the strong-field region).In [32,33], the thermodynamic characteristics are presented in the form of series expansions in the variables, which are combinations of the temperature and field.Our calculations in the ρ 4 model approximation were also performed for temperatures T > T c [34] and T < T c [35] without using similar expansions for the roots of cubic equations appearing in the theoretical analysis.In this paper, the free energy of a 3D uniaxial magnet within the framework of the more complicated ρ 6 model is found introducing the generalized point of exit of the system from the critical regime.The expression for this point takes into account the temperature and field variables simultaneously.In our earlier article [29], the point of exit of the system from the critical regime was found in the simpler non-Gaussian approximation (the ρ 4 model) using the numerical calculations.In contrast to [29], the point of exit of the system in the present paper is explicitly defined as a function of the temperature and field.This allows one to solve our problem, which consists in obtaining the free energy of a 3D Ising-like system in the higher non-Gaussian approximation and in the presence of an external field without involving numerical calculations and without using series expansions in the scaling variable.

General Relations
We consider a 3D Ising-like system on a simple cubic lattice with N sites and period c in a homogeneous external field h.Such a system is described by the Hamiltonian where r jl is the distance between particles at sites j and l, and σ j is the operator of the z component of spin at the jth site, having two eigenvalues +1 and −1.The interaction potential has the form of an exponentially decreasing function Here A is a constant and b is the radius of effective interaction.For the Fourier transform of the interaction potential, we use the following approximation [2,10,11]: where B is the boundary of the Brillouin half-zone 3 .In the CV representation for the partition function of the system, we have [2,36] Here, the summation over the wave vectors k is carried out within the first Brillouin zone, β = 1/(kT) is the inverse temperature, the CV ρ k are introduced by means of the functional representation for operators of spin-density oscillation modes is the Jacobian of transition from the set of N spin variables σ l to the set of CV ρ k , and δ k1+•••+k2n is the Kronecker symbol.
Proceeding from ( 4) and ( 5), we obtain the following initial expression for the partition function of the system in the ρ 6 model approximation: Here The expressions for the remaining coefficients are given in [8][9][10][11].These coefficients are functions of s 0 , that is, of the ratio of microscopic parameters b and c.We shall use the method of "layerby-layer" integration of ( 6) with respect to variables ρ k described in [2][3][4].The integration begins from the variables ρ k with a large value of the wave vector k (of the order of the Brillouin half-zone boundary) and terminates at ρ k with k → 0. For this purpose, we divide the phase space of the CV ρ k into layers with the division parameter s.In each nth layer, the Fourier transform of the interaction potential is replaced by its average value (the arithmetic mean in the given case).
The integration over the zeroth, first, second,. .., nth layers of the CV phase space leads to the representation of the partition function in the form of a product of the partial partition functions Q n of individual layers and the integral of the "smoothed" effective measure density The expressions for Q n , Q(P n ) are presented in [8][9][10][11], and n+1) .The sextic measure density of the (n+1)th block structure W (n+1)

6
(ρ) has the form where 1 and a (n+1) 2l are the renormalized values of the coefficients a 1 and a 2l after integration over n + 1 layers of the phase space of CV.The coefficients a 4 = s −4n u n , and a (n) 6 = s −6n w n are connected with the coefficients of the (n + 1)th layer through the RR whose solutions in the region of the critical regime are used for calculating the free energy of the system.Here The quantity q = qβ Φ(0) determines the average value of the Fourier transform of the potential β Φ(B n+1 , B n ) = β Φ(0) − q/s 2n in the nth layer (in this paper, q = (1 + s −2 )/2 corresponds to the arithmetic mean value of k 2 on the interval (1/s, 1]).The basic arguments h n and α n are determined by the coefficients of the sextic measure density of the nth block structure.The intermediate variables η n and ξ n are functions of h n and α n .The expressions for both basic and intermediate arguments as well as the special functions appearing in (11) are the same as in the absence of an external field (see [8][9][10][11]).The quantities E l in (10) are the eigenvalues of the matrix of the RG linear transformation We have E 1 = R 11 = s (d+2)/2 .Other nonzero matrix elements R i j (i = 2, 3, 4; j = 2, 3, 4) and the eigenvalues E 2 , E 3 , E 4 coincide, respectively, with the quantities R i1 j1 (i 1 = i − 1; The quantities f 0 , ϕ 0 , and ψ 0 characterizing the fixed-point coordinates as well as the remaining coefficients in (10), are also defined on the basis of expressions corresponding to a zero external field.The temperature-independent quantities f 0 , ϕ 0 , and ψ 0 as well as the renormalized quantities f 0 = f 0 /q, ϕ 0 = ϕ 0 /q 2 , and ψ 0 = ψ 0 /q 3 independent of the potential averaging are presented in Table 1 for the optimal RG parameter s = s * = 2.7349.

Contributions to the Free Energy of the System in the Presence of an External Field
Let us calculate the free energy F = −kT ln Z of a 3D Isinglike system above the critical temperature T c .The basic idea of such a calculation on the microscopic level consists in the separate inclusion of the contributions from short-wave (F CR , the region of the critical regime) and long-wave (F LGR , the region of the limiting Gaussian regime) modes of spinmoment density oscillations [2][3][4]: Here F 0 = −kTN ln 2 is the free energy of N noninteracting spins.Each of three components in ( 14) corresponds to individual factors in the convenient representation for the partition function given by (7).The contributions from short-and long-wave modes to the free energy of the system in the presence of an external field are calculated in the ρ 6 model approximation according to the scheme proposed in [8][9][10][11].Short-wave modes are characterized by an RG symmetry and are described by the non-Gaussian measure density.The calculation of the contribution from long-wave modes is based on using the Gaussian measure density as the basis one.Here, we have developed a direct method of calculations with the results obtained by taking into account the short-wave modes as initial parameters.The main results obtained in the course of deriving the complete expression for the free energy of the system are presented below.

Region of the Critical Regime. A calculation technique
based on the ρ 6 model for the contribution F CR is similar to that elaborated in the absence of an external field (see, e.g., [4,9,10]).Carrying out the summation of partial free energies F n over the layers of the phase space of CV, we can calculate F CR : An explicit dependence of F n on the layer number n is obtained using solutions (10) of RR and series expansions of special functions in small deviations of the basic arguments from their values at the fixed point.The main peculiar feature of the present calculations lies in using the generalized point of exit of the system from the critical regime of order-parameter fluctuations.The inclusion of the more complicated expression for the exit point (as a function of both the temperature and field variables) [37] leads to the distinction between formula ( 16) for F CR and the analogous relation at h = 0 [9,10].The quantity h = h / f 0 is determined by the dimensionless field h , while the quantity h c = τ p0 is a function of the reduced temperature τ = (T − T c )/T c .Here 1 characterizes the coefficient c 1 in solutions (10) of RR, ν = ln s/ ln E 2 is the critical exponent of the correlation length.At h = 0, n p becomes m τ = − ln τ/ ln E 2 − 1 (see [4,9,10]).At T = T c (τ = 0), the quantity n p coincides with the exit point n h = − ln h/ ln E 1 − 1 [38].The limiting value of the field h c is obtained by the equality of the exit points defined by the temperature and by the field (m τ = n h ).
Having expression (17) for n p , we arrive at the relations [39] where Δ 1 = − ln E 3 / ln E 2 and Δ 2 = − ln E 4 / ln E 2 are the exponents, which determine the first and second confluent corrections, respectively.Numerical values of the quantities E l (l = 1, 2, 3, 4), ν, Δ 1 , and Δ 2 for the RG parameter s = s * = 2.7349 are given in Table 2.In the weak-field region ( h h c ), quantities (18) can be calculated with the help of the following expansions: In the strong-field region ( h h c ), these quantities satisfy the expressions It should be noted that the variables h/ h c (the weak fields) and ( h c / h) 1/p0 (the strong fields) coincide with the accepted choice of the arguments for scaling functions in accordance with the scaling theory.In the particular case of h = 0 and τ / = 0, (19) are defined as We shall perform the further calculations on the basis of (18), which are valid in the general case for the regions of small, intermediate (the crossover region), and large field values.The inclusion of E np+1 3 (or H 3 ) leads to the formation of the first confluent corrections in the expressions for thermodynamic characteristics of the system.The quantity E np+1 4 (or H 4 ) is responsible for the emergence of the second confluent corrections.The cases of the weak or strong fields can be obtained from general expressions by using (19) or (20).We disregard the second confluent correction in our calculations.This is due to the fact that the contribution from the first confluent correction to thermodynamic functions near the critical point (τ = 0, h = 0) is more significant than the small contribution from the second correction ( h2 + h 2 c 1, Δ 1 ≈ 0.5, and Δ 2 is of the order of 3, see Table 2).
Proceeding from an explicit dependence of F n on the layer number n [4,8,9] and taking into account (18), we can now write the final expression for F CR (16): Here c (0) 20 characterizes c 2 in solutions (10) of RR, and the coefficients are determined by the components of the quantities The components δ (i) 0 (i = 0, 1, 2) satisfy the earlier relations [4,8,9] obtained in the case of a zero external field.The components γ (i)  0 are given by the corresponding expressions at h = 0 under condition that the eigenvalues E 1 , E 2 , and E 3 should be replaced by E 2 , E 3 , and E 4 , respectively.
Let us now calculate the contribution to the free energy of the system from the layers of the CV phase space beyond the point of exit from the critical regime region.The calculations are performed according to the scheme proposed in [2,4,10,11].As in the previous study, while calculating the partition function component Z LGR from (15), it is convenient to single out two regions of values of wave vectors.The first is the transition region (Z (1)  LGR ) corresponding to values of k close to B np , while the second is the Gaussian region (Z (2)  LGR ) corresponding to small values of wave vector (k → 0).Thus, we have (1)  LGR Z (2)  LGR . (25)

Transition Region.
This region corresponds to m 0 layers of the phase space of CV.The lower boundary of the transition region is determined by the point of exit of the system from the critical regime region (n = n p + 1).The upper boundary corresponds to the layer n p + m 0 + 1.We use for m 0 the integer closest to m 0 .The condition for obtaining m 0 is the equality [9,10] where A 0 is a large number (A 0 ≥ 10).
The free energy contribution + ln I 0 η np+m , ξ np+m (27) corresponding to Z (1)  LGR from ( 25) is calculated by using the solutions of RR.
The basic arguments in the (n p + m)th layer w np+m u 3/2 np+m (28) can be presented using the relations obtained on the basis of ( 10) and (18).We arrive at the following expressions: In contrast to H c , the quantity H 3 in expressions (30) for h np+m and α np+m as well as in expression (21) for F s takes on small values with the variation of the field h (see Figure 3).The quantity H c at h → 0 and near h c is close to unity, and series expansions in H c are not effective here.Power series in small deviations (h np+m − h (0) np+m ) and (α np+m − α (0) np+m ) for the special functions appearing in the expressions for the intermediate arguments allow us to find the relations The quantities η (0) np+m , η , and ξ (0) np+m , ξ , where Proceeding from expression (27) for f LGR1 (m), we can now write the following relation accurate to within H 3 : LGR1 (m)c (0) 20 H 3 , LGR1 (m) = ϕ The quantities b The final result for F (1)  LGR (see (27) and ( 34)) assumes the form F (1)  LGR = −kTN s −3(np+1) LGR1 (m), LGR1 (m). ( On the basis of ( 26) and (30), it is possible to obtain the quantity m 0 determining the summation limit m 0 in formulas (36): Let us now calculate the contribution to the free energy of the system from long-wave modes in the range of wave vectors using the Gaussian measure density.

Region of Small
Values of Wave Vector (k → 0).The free energy component LGR = corresponding to Z (2)  LGR from ( 25) is similar to that presented in [4,9,10].The calculations of the first and second terms in (39) are associated with the calculations of the quantities where and (1)  n p −1 c (0) 20 H 3 ) satisfy the corresponding expressions from ( 29) and ( 30) at m = m 0 + 1.
Introducing the designation and presenting it in the form we obtain the following relations for the coefficients: The quantities (45) determine the function (1) Here F * (n p −1)

0
, where Taking into account (41) and (43), we rewrite formulas (40) as The second term in (39) is defined by the expression 1 2 Relations ( 48) and (49) make it possible to find the component F (2)  LGR in the form On the basis of ( 36) and (50), we can write the following expression for the general contribution F LGR = F (1)  LGR + F (2)   LGR to the free energy of the system from long-wave modes of spin-moment density oscillations: LGR + f LGR c (0) 20

Total Free Energy of the System at T > T c
The total free energy of the system is calculated taking into account ( 14), (21), and (51).Collecting the contributions to the free energy from all regimes of fluctuations at T > T c in the presence of an external field and using the relation for s −(np+1) from ( 18), we obtain LGR , l = 0, 1. (52) The coefficients γ (CR) 0 , γ 1 , and γ 2 are defined by ( 23), g 1 is presented in (48), and γ + 4 is given in (50).The coefficients of the nonanalytic component of the free energy F [see (52)] depend on H c .The terms proportional to H 3 determine the confluent corrections by the temperature and field.As is seen from the expression for F, the free energy of the system at h = 0 and τ = 0, in addition to terms proportional to τ 3ν (or h 6/5 c ) and h 6/5 , contains the terms proportional to τ 3ν+Δ1 and h 6/5+Δ1/p0 , respectively.At h / = 0 and τ / = 0, the terms of both types are present.It should be noted that Δ 1 > Δ 1 / p 0 .At h = h c , we have τ 3ν+Δ1 = h 6/5+Δ1/p0 and the contributions to the thermodynamic characteristics of the system from both types of corrections become of the same order.
The advantage of the method presented in this paper is the possibility of deriving analytic expressions for the freeenergy coefficients as functions of the microscopic parameters of the system (the lattice constant c and parameters of the interaction potential, that is, the effective radius b of the potential, the Fourier transform Φ(0) of the potential for k = 0).

Conclusions
A 3D Ising-like system (a 3D uniaxial magnet) in a nonzero external field is investigated in the higher non-Gaussian approximation based on the sextic distribution for modes of spin-moment density oscillations (the ρ 6 model).The simultaneous effect of the temperature and field on the critical behaviour of the system is taken into account.An external field is introduced in the Hamiltonian of the system from the outset.In contrast to previous studies on the basis of the asymmetric ρ 4 model [32,33,40], the field in the initial process of calculating the partition function of the system is not included in the Jacobian of transition from the set of spin variables to the set of CV.Such an approach leads to the appearance of the first, second, fourth, and sixth powers of CV in the expression for the partition function and allows us to simplify the mathematical description because the odd part is represented only by the linear term.When the field is included in the transition Jacobian, the measure density involves the odd powers of CV in addition to the even powers.
An initial expression for the partition function of the system is constructed in the form of a functional with explicitly known coefficient functions [see (6)].The partition function is integrated over the layers of the CV phase space.The RR and their solutions near the critical point are written for the ρ 6 model containing the field term.The presence of the field caused appearance of the additional relation for t n in (10).The fixed-point coordinates and the eigenvalues of the RG linear transformation matrix as well as the exponents of the confluent corrections obtained in the higher non-Gaussian approximation are given.
The main distinctive feature of the presented method for calculating the total free energy of the system is the separate inclusion of the contributions to the free energy from the short-and long-wave spin-density oscillation modes.The expression for the generalized point of exit of the system from the critical regime contains both the temperature and field variables.The form of the temperature and field dependences for the free energy of the system is determined by solutions of RR near the fixed point.
The expression for the free energy F (52) obtained at temperatures T > T c without using power series in the scaling variable and without any adjustable parameters can be employed in the field region near h c (the crossover region).The limiting field h c satisfies the condition of the equality of sizes of the critical regime region by the temperature and field (the effect of the temperature and field on the system in the vicinity of the critical point is equivalent) [32,33,38,40].In the vicinity of h c , the scaling variable is of the order of unity, and power series in this variable are not effective.Proceeding from the expression for the free energy, which involves the leading terms and terms determining the temperature and field confluent corrections, we can find other thermodynamic characteristics (the average spin moment, susceptibility, entropy, and specific heat) by direct differentiation of F with respect to field or temperature.

Figure 1 :Figure 2 :
Figure1: Evolution of the critical exponent of the correlation length ν with increasing parameter of division of the CV phase space into layers s.Curves 1, 2, 3, and 4 correspond to the ρ 4 , ρ 6 , ρ 8 , and ρ 10 models, respectively.

3 Figure 3 :
Figure 3: Dependence of quantities H c and H 3 on the ratio h/ h c for the RG parameter s = s * = 2.7349 and the reduced temperature τ = 10 −4 .

Table 1 :
Quantities, which characterize the coordinates of the fixed point.

Table 2 :
The eigenvalues E l and the exponents ν, Δ 1 , and Δ 2 for the ρ 6 model.