Gibbs free energy and Helmholtz free energy for a three-dimensional Ising-like model

The critical behavior of a 3D Ising-like system is studied at the microscopic level of consideration. The free energy of ordering is calculated analytically as an explicit function of temperature, an external field and the initial parameters of the model. Within a unified approach, both Gibbs and Helmholtz free energies are obtained and the dependencies of them on the external field and the order parameter, respectively, are presented graphically. The regions of stability, metastability, and unstability are established on the order parameter-temperature plane. The way of implementation of the well-known Maxwell construction is proposed at microscopic level.


Introduction
Basic principles of phenomenological theory for second order phase transitions (PT) were formulated by Landau in late 30-40s [1,2], originally to describe superconductivity. The key assumption of the theory is that in the vicinity of the critical point, the free energy can be expanded in a power series in the order parameter, the equilibrium value of which is found to obey the minimum condition for free energy. As is now known [3], this expansion is well defined in the neighborhood of T c (critical value of temperature) except for a narrow interval determined by the Ginzburg criterion [4]. In this interval, the crucial role is played by the order parameter fluctuations while Landau's theory is essentially a mean-field theory. Nevertheless, it is of great use as a qualitative tool for understanding the nature of PT.
The above mentioned power series expansion is often called Landau's free energy. It turns out that similar quantity can also be derived (not constructed) in a theory considering PT at a microscopic level. Indeed, in the collective variables (CV) method proposed in [5,6] to describe the critical behavior of spin systems, such a quantity appears in a natural way. In particular, a microscopic analogue of Landau free energy was calculated in [7], where the Ising-like system was considered in zero external field.
The purpose of the present paper is to extend the results obtained earlier by the CV method for physical characteristics of 3D Ising-like model to the presence of an external field and look at what is going on when the field changes sign. As a result, we obtain the Gibbs free energy, Helmholtz free energy, and the order parameter and establish the regions of stability, metastability, and unstability for the system under investigation.

Method
In our research we consider the system of N Ising spins placed on sites of a simple cubic lattice of the spacing c. The Hamiltonian of such a system in an external field is well known Here, Φ(r i,j ) is a short-range interaction potential between spins located at the i-th and j-th sites, the spin variables σ i take on values ±1, H is the external field. We do not restrict the summation in (2.1) to the nearest neighbors. To this end, Φ(r) = const × exp (−r/b) can be chosen as the interaction potential with an effective range b.
It is also known that this problem has not yet been solved. Universal critical characteristics of three-dimensional (3D) systems are successfully described in a variety of approaches. Among them it is worth mentioning the field theoretical and renormalization group (RG) methods [8], which also enable one to calculate nonuniversal characteristics (critical amplitudes, critical temperatures), but without taking into account the dependency on initial parameters of the Hamiltonian. An alternative method for theoretical investigation of the model under consideration is the method of CV [5], which has an advantage of being a successive microscopic approach, which, in turn, makes it possible for thermodynamic functions and physical characteristics of the system to be explicitly obtained as functions of temperature, the field, and microscopic parameters of the model.
In the framework of the "ρ 4 -model" approximation the functional representation for the partition function in terms of CV ρ k is as follows: Here, the quantity d(k) contains the Fourier transform of the interaction potential The explicit expressions for a n can be found in [9] (equation (1.25)) , N 0 = N/s 3 0 where the quantity s 0 defines the region of validity for the parabolic approximation of Fourier transform of the interaction potential (for details, see [10]). The summation in (2.2) is performed over wave vectors of the first Brillouin zone corresponding to a reciprocal effective lattice with the lattice constant cs 0 .
In the work [11] the method for calculation of the partition function (2.2) of the model near the second order PT point was generalized to the case of the presence of the fixed external field of an arbitrary magnitude. Unlike the work [12], here no assumption concerning either the strength or the weakness of the applied field is made and, therefore,no perturbation series are used. The calculation procedure is based on Kadanoff's idea of constructing effective spin blocks [13]. We divide the phase space of the CV ρ k into layers, the RG parameter being s, and average the Fourier transform of the potential in each n-th layer. The first step in calculating equation (2.2) is to integrate over ρ k with k ∈ [k max , k max /s], k max = π/cs 0 . The result of integration proves to be represented in the form of the initial expression, equation (2.2), with renormalized coefficients a 4 . If we perform a step-by-step integration of the partition function over n p layers, we arrive at The partial partition functions Q n of the n-th level are characterized by a set of coefficients a (n) 1 , d n (0), a (n) 4 , for which the recurrence relations (RR) hold (the work [9] is devoted to this problem). The quantity n p is called the exit point from the critical regime of the order parameter fluctuations. It defines the number of iterations at which the system is still in the scaling region.
In [12], two regions in (h, τ )-plane were distinguished, which are of weak and of strong field. The calculations in the regions were performed using different forms of exit points. These quantities were subject to certain conditions which we do not discuss in the present work. If the field was considered to be strong, the formula n p = n h = − lnh/ ln E 1 − 1 was used, while in the case of a weak field another expression n p = m τ = − lnτ / ln E 2 − 1 was taken. Hence, the question immediately arises what should be done at intermediate values of a field. The problem is solved by constructing a new expression for the exit point, proposed in [11]: where some temperature fields h [14], the signs "+" and "-" are related to T > T c and T < T c respectively, and The quantities E 1 and E 2 are the eigenvalues of the matrix of the RG transformation linearized near the fixed point of RR.
In the limiting cases, as h → 0 or as τ → 0, equation (2.5) takes on the form of m τ or n h , respectively, and therefore can be applied to the case of an arbitrary field. Nevertheless, it is worth mentioning that the inequalityh h c defines a strong field, andh h c defines a weak field. Note also that there is an arbitrariness in choosing n p , which is discussed in [10,11] in more detail. What is important is that for n > n p one should be able to integrate the partition function with Gaussian measure density.
We will adhere to the approach adopted in [15], where the coefficients c k1 , f 0 , h 0 are also defined (see equations (4.4), (4.17), (4.23), and (4.49) in [15]), and take E 1 = 24.551, E 2 = 8.308. The information on the physical meaning of coefficient n 0 can be found in [16]. Numerical values of all coefficients used to obtain the final graphical results are presented in table 1. It is worth stressing that having fixed the value of RG parameter s = s * (s * = 3.5977 in the present calculation), there remains only one initial parameter, which is the ratio of the effective range of interaction b to the lattice constant c. All the other coefficients can be expressed via b/c [10,17].

Results and discussion
Based on the above mentioned method, in the work [11] as well as in [16] free energy of 3D Ising-like model was calculated as a logarithm of the partition function (2.4) multiplied by −kT and presented in the form of several contributions Here, the term F a is the analytical part of free energy and does not affect the critical behavior of the system. The last two terms in r.h.s. of (3.1) have non-analytical dependence on temperature and on the external field. The explicit expressions for F and F The quantity γ (±) s from (3.2) includes contributions from the critical regime of the order parameter fluctuations and from the limiting Gaussian regime (LGR)(for details, see [11,16]). The mentioned critical regime is characterized by the RG symmetry.
The quantity E 0 (σ ± ) is the contribution from the collective variable ρ 0 . As is known from the theory of CV [5], the mean value of ρ k=0 is connected with the order parameter and consequently the quantity F 0 is the free energy of ordering [7]. If there is no external field, F 0 is analogous to Landau's free energy in phenomenological theory of second order PT [7]. Since we work in the ensemble of N particles, with field h and temperature τ being independent variables, this analogy is now not direct. Following Stanley [14], and based on convincing arguments of the work [18], we associate the partition function of our model with Gibbs free energy. An appropriate thermodynamic potential is Gibbs free energy provided the magnetic field and temperature are independent ("natural") variables. If the role of independent variables is played by an order parameter (here magnetization per spin) and temperature, then the thermodynamic potential associated with the partition function is Helmholtz free energy. Hence, Landau's free energy is rather Helmholtz than Gibbs one. The expressions (3.1)-(3.3) in turn present Gibbs free energy for the considered model, although this was not specified in earlier works, particularly in [11,16]. Therefore, in order to obtain a microscopic analogue of Landau's energy one should perform the Legendre transformation. Let A 0 denote this analogue. Then, where M 0 is the magnetization of the system under consideration such that The expression for E 0 (σ ± ), which was found in [11] (see equation (4.4) there) for T > T c and in [16] The quantity σ ± was found from condition in the form This results in the cubic equation for σ (±) 0 where the denotations r n and u n are connected with d n (0) and a (n) 4 through d n (0) = s −2 r n , a (n) 4 = s −4 u n , and according to [11,16] are expressed as Note, that the coefficients p, q should be marked by the superscript ± (as should be done for n p from (2.5)), but where it does not cause confusion, we drop this superscript out. One should just remember about a difference in temperature scales between the cases of T > T c and of T < T c . Equation (3.9) can be solved by Cardano's method. The solutions are presented graphically in figure 1. From the plot we see that for a given value of field there exists some temperature τ 0 < 0 such that equation (3.9) has three real solutions for τ τ 0 and one real solution for τ > τ 0 . This quantity, τ 0 , corresponds to the condition Q = 0 where Q is the discriminant of the cubic equation (3.9): Q = (p/3) 3 + (q/2) 2 (3.12) with q and p from (3.10). When h = 0, then τ 0 = 0. If h = 0, then τ 0 is determined by setting r.h.s. of (3.12) equal to zero. In [16], τ 0 was found numerically and the plot of τ 0 versus h was drawn. As we will see in what follows, the solution τ 0 = τ 0 (h) to the equation Q = 0 defines a curve in "order parameter-temperature" plane which can be identified with the spinodal of a fluid. Taking into account (2.5), and (3.8), the quantity E 0 (σ) takes on where we have introduced the following notation (3.14) In the previous works [10,11,15,16] the system was considered in the external field h 0. Analytical results for (Gibbs) free energy [11,16], the order parameter [10,16], and the susceptibility [10] were obtained. The purpose of this paper is to extend the method to the region h < 0 and to take into account all solutions to the equation (3.9). We look at Gibbs free energy F 0 and at the contribution from it to the order parameter M 0 of equation (3.5). Thereupon, Helmholtz free energy is calculated.
The main results of this paper are presented in figures 2, 3, and 4. In  . Thick line is magnetization in zero external field and corresponds to the coexistence curve (binodal). In negative external field, the solution σ01 gives rise to Curves 1-5 on the magnetization-temperature plane. In positive external field, σ03 gives rise to Curves 6-10. Curves 11-20, which correspond to σ02, can be interpreted as non-physical. Thick dashed curve corresponds to the saturation curve (spinodal). In order to recover that, one should compute magnetization corresponding to σ01 at h < 0 and τ0 = τ0(|h|) (upper branch) and corresponding to σ03 at h > 0 and τ0 = τ0(h) (lower branch). The magnitudes of field, at which the results are drawn, are |h| = 0, 10 −7 , 10 −6 , 5 · 10 −6 , 10 −5 , 2 · 10 −5 .

43002-6
Free energy of an Ising-like model the system in different coordinates. Especially, points a, c, d, and f correspond to such a relation between the field h and temperature τ that the equality τ = τ 0 (h) holds. Thus, in each of these we have Q = 0. Points b and e are associated with first order phase transition and correspond to h = 0 and τ = −0.001. The isotherm in figure 2 (b) was drawn with the help of parametric Recall that Gibbs free energy is a concave function of field [14]. Therefore, from figure 2 (a) we can conclude that the solution σ (−) 02 gives rise to a non-physical result. It corresponds to the region of unstability (c − d part of the curve in figure 2 (a), (b)). This is not surprising if we note that σ In the case of liquid-gas PT it would be the region between the binodal and the spinodal curves. Hence, we regard these results important from the perspective of using this method for description of the critical behavior in simple fluids. Figure 2 (c) to some extent repeats figure 1 (b) and is placed here for convenience of representation. The quantity σ 0 is related to the scaling function of magnetization via formula (4.8) in work [16]. Finally, in figure 2 (d) the equation of state computed by means of (3.5) is presented at τ = −0.001.
A more detailed representation of the equation of state is shown in figure 3. The derivatives of Gibbs free energy with respect to field are presented as functions of temperature, different solutions to equation (3.9) being taken into account. As we can see, such an accounting allows us to establish spinodal and binodal curves of the model and, in this respect, to find the stable, metastable, and unstable regions on "order parameter-temperature" plane.
The order parameter dependence of free energy A 0 is presented in figure 4 for three different temperatures T > T c , T = T c and T < T c . Here again a, b, c, . . . denote the states as in figure 2. We can see that Helmholtz free energy has two minima below T c . The equilibrium values for the order parameter are defined by these minima. Knowing their positions at different temperatures, we can recover the coexistence curve. Points b and e coincide with the minima at τ = −0.001. Note that based on the principle of minimum Gibbs free energy and from figure 2 (a), we can conclude that with the field changing from h < 0 to h > 0 the system will tend to move along a − b − e − f part of the curve. Applied to Helmholtz free energy, it means that the system will tend to jump from state b to state e. Therefore, one may supplement figure 4 with a double-tangent construction. As 43002-7 a consequence, on M 0 − h plane we have a horizontal segment for M 0 at h = 0 (see figure 2 (d)), which corresponds to a Maxwell construction derived at microscopic level. On the other hand, based on the maxima of Gibbs free energy with respect to M 0 (see figure 2 (b)), we are able to construct a saturation curve. Let us just remember that dependence F 0 on the order parameter is formal due to M 0 being not "natural" variable for F 0 . However, the same results for the binodal and the spinodal can immediately be obtained if one appropriately chooses the solutions to (3.9).
The forms of isotherms in figure 4 are similar to those in Landau's theory [3]. However, the peculiar feature is that in the present work the free energy is an explicit function of temperature, external field and microscopic parameters (in the present case a ratio b/c) of the model, with a non-analytical dependence on its arguments. This gives rise to a non-classical critical behavior of the system, i.e., to the critical exponents taking on non-classical values. These values for the most important exponents are reported in table 2. There are also collected the values obtainable by CV method in higher, ρ 6 approximation as reported in Conclusions of the work [7]. Note that in our calculation we have neglected the critical exponent η responsible for the behavior of the pair correlation function at T = T c . Accounting for corrections to scaling is beyond the scope of the present paper as well. In table 2 the classic values for the critical exponents are presented as well as the results of other authors. We stick to the following notation. The temperature behavior of magnetization is governed by β, of the heat capacity by α, of susceptibility by γ, of the correlation length by ν. The field behavior of magnetization is governed by δ, M ∼ |h| 1/δ . We have calculated some of the exponents with the help of scaling laws. Based on earlier results [15], we are also able to estimate the critical exponents ϕ and ψ that describe the field behavior of the heat capacity, C ∼ |h| −ϕ , and of the entropy, S ∼ |h| ψ , respectively. Their numerical values are ϕ = 0.122 and ψ = 0.539, which differ from ϕ = 0 and ψ = 2/3 in mean-field theory [14]. It should also be mentioned that the scaling properties for physical characteristics of the considered model are discussed in earlier works. In section 3 of the work [10] the scaling functions for (Gibbs) free energy, for the order parameter, and for susceptibility were calculated.

Conclusions
In this work, an analytical expression for Gibbs free energy of Ising-like system is obtained and investigated near the second order phase transition. The emphasis is made on the presence of the external field and on the situation where the field changes its direction. Some well-defined concepts of phenomenological theory of phase transitions -e.g. Landau's energy, the Maxwell construction, the double-tangent construction -are derived and reexamined at microscopic level. We hope that the approach described and the results presented should be useful in further investigations of critical behavior in 3D Ising-like systems with analytical methods; in particular, in order to obtain the coexistence curve and the spinodal curve for a system at least in the vicinity of critical point. There exists a general belief [14,20,22] that simple fluids belong to the Ising universality class. Hence, we regard these results important from the perspective of using this method for description of the critical behavior in simple fluids.