Phase transition in a cell fluid model

We propose a method of describing a phase transition in a cell fluid model with pair interaction potential that includes repulsive and attractive parts. An exact representation of the grand partition function of this model is obtained in the collective variables set. The behavior of the system at temperatures below and above the critical one is explored in the approximation of a mean-field type. An explicit analytic form of the equation of state which is applicable in a wide range of temperatures is derived, taking into account an equation between chemical potential and density. The coexistence curve, the surface of the equation of state and the phase diagram of the cell Morse fluid are plotted.


Introduction
Description of phase transitions in fluids remains a ponderable area of investigations in statistical physics both in macroscopic and microscopic scales. During many years of research in the field due to hard work of scientists it became known that the phase transition is possible either at the infinite volume or at the thermodynamic limit and is characterized by a density jump on the isotherm below the critical temperature [1]. Plenty of successful studies where dedicated to explain the nature of a first-order phase transitions on the macroscopic level, but resemblant theoretical description on the microscopic stage is still a question. Most approaches to description of phase transitions and critical phenomena in fluids are based on complete scaling ideas [2,3]; methods of integral equations, and in particular selfconsistent Ornstein-Zernike approximation (SCOZA) [4,5]; perturbation series expansion, for example hierarchical reference theory [6]; non-perturbative renormalization group approach [7]; method based on the study of behavior of the virial equation of state with extrapolated coefficients [8] numerical methods and computer simulations. There are not so many recent works devoted to looking for the possible ways to solve this problem in the framework of grand canonical ensemble. Among them, for example, the collective variables method [9] used by Yukhnovskii to completely integrate the grand partition function of a system of interacting particles in the phase-space of collective variables and therefore investigate its behavior in a vicinity of the critical point, or the method developed by Tang [10] were he combined grand canonical ensemble with density functional to explain the occurrence of a first-order phase transition in homogeneous fluids.
The objective of our investigation is microscopic description of a first-order phase transition in a cell fluid model. We do this only on the base of grand canonical ensemble. The model we used as an approximation of a continuous system is quite similar to the cell gas model [11] Thermodynamic functions, correlation functions and free energy of the latter appeared to be close to the corresponding values of continuous systems if the cells size tends to zero [12] Somewhat similar model was effective for estimating the cluster size with hydrogen bounds [13].
To obtain numerical results the values of parameters of the interaction potential are required. For this purpose we chose the Morse potential as a potential of interaction. The Morse fluid has already been well studied, e.g., within the integral equation approach [14], by molecular dynamics simulations of simple fluids [15], by Monte Carlo simulations using both N pT plus test particle method [16] and the grand-canonical transition matrix method [17]. Appliance of such a potential appeared to be sufficient in description of a liquid-vapor coexistence in liquid metals [17]. This paper is laid as follows: in Section 2 the Jacobian of transition from individual to collective variables is calculated. A rigorous representation of the GPF of the cell fluid in the form of multiple integral over collective variables is obtained in Section 3. In Section 4 this expression is restricted to ρ 4 -model and calculated in the mean-field approximation. Section 4 and 5 are dedicated to exploration of behavior of the system in wide temperature region except a vicinity of the critical point, an equation of state of the cell fluid applicable in this region is derived. Conclusions are presented in Section 6.

A cell model for description of simple fluid behavior
To follow our goal we calculate the grand partition function (GPF) of the cell fluid model.
As well as in case of a cell gas model, the idea of the cell fluid consists in fixed partition of the volume V of a system with N particles on N v congruent cubic cells, each of volume v = V /N v [18]. The GPF of this system can be written in the form Here z = e βµ is the chemical activity, β is the inverse temperature, µ is the chemical potential. In this expression integration is made over coordinates of all the particles in the system is the potential energy of the model, difference between two cell vectors is noted as l 12 = | l 1 − l 2 |. Each l i takes values from a set Λ, defined as Here c is the side of a cell, N 1 is the number of cells along each axis, ρ l (η) is the occupation number of a cell The characteristic functions (indicators) I ∆ l (x) identify particles in cells and their contribution to interaction of the model. This is the way to plunk all the particles and sort them to different cells. The following expression is obvious Particles with different coordinates can get into the same cell. Interaction between constituents of one cell and particles hosted in another cell is expressed as functionŨ l 12 dependent on the distance between cells. This function consists of repulsive and attractive parts respectivelyŨ l 12 = Ψ l 12 − U l 12 . Both Ψ l 12 and U l 12 are positive. As an example we chooseŨ l 12 in the form of the Morse potential where α is the effective interaction radius. The parameter R 0 corresponds to the minimum of the functionŨ l 12 (Ũ (R 0 ) = −D determines the depth of potential well).
On this stage the GPF (1) of N interacting particles will be written in the form of N integrals over their coordinates and N v integrals over the collective variables (CV) ρ l [19] Herewith Let us transform the initial potential of interaction as follows Here the effective interaction potential W l 12 is expressed as Hence χ is the real positive parameter (χ > 0), β c = 1/k B T c , k B is the Boltzman constant, T c is some fixed temperature which will be defined later. As a result a part of the repulsive interaction β c χΨ l 12 δ l 1 ,l 2 > 0 is transferred to the Jacobian of transition from individual coordinates to collective variables. As we will see subsequently it will form a new measure of distribution of density fluctuations. As a result of the above-mentioned transformations the GPF takes the form here e βµN = exp βµ l∈Λ ρ l (η) and Ψ 0 = De 2R 0 /α To integrate over the coordinates x 1 , . . . , x N it's enough to make the following transformation of a part of (7) The parameter p is real and positive. It is expressed by the following formula After integrating over the coordinates x 1 , . . . , x n (see Appendix A) we obtain a functional representation of the GPF of the cell fluid model in similar form to the one we had before [20] Now the difference between the cell gas model [11] and the cell fluid model becomes clear. Molecules of a cell gas move in space in such a way, that each cell can house no more than one particle. This is somehow likely to the lattice-gas model, which is frequently used for description of simple fluids. As it is known in the lattice-gas only one particle can reside on each site, otherwise the site is vacant. Regarding to the cell fluid model we introduce the occupation numbers of cells ρ l (η) which can take values N = 0, 1, 2, . . .. This means that each cell of a cell fluid can host arbitrary number of particles. Due to the term e −pN 2 the more N increases the less N -th term contributes to the sum in (9). Therefore the probability to find many particles in one elementary cube is very small.
3 Functional representation of the grand partition function of the cell fluid model The GPF in the CV set ρ l is quiet compact but not perspective due to the sum of non-diagonal terms in the exponent with effective interaction. Let us use the Fourier representation in (9) Therefore to make any estimations it is convenient to choose the potential of interaction possessing the Fourier transform. The variable ρ k is the representation of the collective variable ρ l in reciprocal space Vector k takes values from the set B c corresponding to one cell The Fourier transform of the effective interaction potential W (k) has the form In case of the Morse potential (4) one has the following On this stage we propose two ways to deal with calculation of the GPF (9). The former is represented in our recent manuscript [20]. The latter, which on our opinion gives more transparent results, will be considered here. Firstly, let us make the Stratonovich-Hubbard transformation are real and imaginary parts respectively. For an element of a phase volume we have The stress above the product means that k > 0.
Note that the mandatory condition for transformation (13) is W (k) > 0. Taking into account the form of interaction potential (4), it is easy to see, that the effective potential W (k) (10) is positive when and χ > 0. The left hand side of inequality (14) is peculiar for the Morse potential, which was used [21,22,23] for description of metals (particularly, Cu, F e, Al, Au, N a and so on).
Taking into account transformation (13) the GPF (9) can be written in the form The Jacobian of transition J(t l ) factors in l Obviously, integrals over ρ l in (16) are delta-functions, which so as to say remove integration over ν l .
In consequence of this we come to the following expression of the GPF The series under the product of l in (18) can be expressed as a cumulant series Coefficients a n are explicit analytical functions of the parameters p and v (see Appendix B) in contradiction to our previous work [20], where such coefficients were expressed in a complicated way via the integrals containing restricted cumulant series. Therefore they could be only numerically calculated. Consequently, we obtained an approximated expression of the GPF [20]. Evidently p and v appeared to be instrumental parameters in the present approach. Moreover there is a monosemantic dependence between them (for more detail see Appendix C).
It is worthily to use a space oft k to make the following calculation of the GPF. In this case the Jacobian of transition fails to be a function of chemical potential, which actually makes calculation much easier. So let us transit in (18) to the variablest k using (17), and rewrite it in the form W (0) is the effective potential of interaction at | k| = 0 and Note that the functional representation (20) is exact for the cell fluid model, in comparison with visually similar expression obtained in our previous manuscript [20], which is approximate. To derive an equation of state and physical characteristics of the system one should calculate the GPF. Computation of (20) would be an exact solution of the problem, but this task is rather complicated, whereas expression (20) contains an infinite number of terms in the exponent. Thereby we use the ρ 4 -model approximation, which means that the cumulant series in (20) are restricted to n ≤ 4. We made sure [27,24,25,26], that this approximation gives qualitatively satisfactory description of a second-order phase transition with non-classical values of critical exponents for the Ising model. So we hope that the same approximation would work for the cell fluid model as well. 4 Calculating an explicit form of the GPF in the approximation of a mean-field type The expression, which was obtained in the previous section is likely to the partition function of a three-dimensional Ising model in an external field [28]. In our case µ plays the role of the field. To calculate the GPF (20) make a change of variables which is aimed to destroy terms proportional to the third power of the variablet k . As a result we obtain the following expression where The coefficient d(k) has the form Currently we are interested in the behavior of the cell fluid in a wide range of temperature. The results, which were obtained for a magnet in an external field [29,30] and also for fluids [31,9] point out, that the variable ρ 0 with k = 0 ( k is the wave vector of a fluctuation mode) gives a main contribution to a free energy in the region of temperatures far from the vicinity of the critical point.
Terms, which include ρ k with k = 0 fail to play a part to the average density of the system since integration over this variables in (21) is not related to the chemical potential, namely a derivative of these terms with respect to the chemical potential is equal to zero. Obviously such terms do not make the density dependent contribution to the pressure So let us make such a type of mean-field approximation and consider the collective variables ρ k with k = 0 both with neglecting all the ρ k with k = 0 in expression (21). In this approximation the GPF has the form The temperature of transition in the mean-field approximation [32,33] can be determined from the following condition Thereof we obtain the expression for the critical temperature where W (0, T c ) = W (0)| T =Tc is the value of the effective interaction potential W (0) at the critical temperature.
Taking into account the form of the effective potential (6) it is easy to make sure, that temperature dependence of d(0) can be expressed as follows is the Fourier transform of the interaction potential between the constituents of the system. Since the value N v in the exponent of (26) tends to infinity in the thermodynamic limit, we can use the Laplace method to make integration in this expression (see [20]) and obtain the GPF as follows where the value ρ 0 =ρ 0 can be found from the equation ρ 0 should give a maximum of E(ρ 0 ). Having an explicit expression of the GPF (31) one can write the average density of the system using (24) where It is important, that we can find the chemical potential as a function of average density, but it appeared to be not an easy thing to express M = f (n, τ ) using (33). However we can act in a different way to linkn with M by writing an equation that meets the condition of maximum of E(ρ 0 ) (27) where according to (33) Using the denotation m = M +n − n g , we can reduce equation (35) to the following form where Since p b < 0 the value of D b can be either positive or negative (see Figure. 1). This fact is an influence on the number of real solutions of equation (38). If D b < 0 equation (38) has three real solutions. If D b = 0 there are two values of densityn. In case of D b > 0 the equation has a single solution describing a decrease ofn with an increase of the chemical potential and that is not a physical behavior.
In the ranges of density 0 ≤n ≤n max D b is negative, as it is plotted on Figure. 1. Going back to calculation of the GPF (31), the case T > T c means a single real solution of (32) Then the equation of state or in other words the pressure as a function of temperature and density has the form Now considering T = T c in (43), it is easy to write an explicit dependence of pressure on density at the critical temperature The following formula expresses the chemical potential M 0 as a function of density Index 0 denotes, that M 0 corresponds to T = T c . The pressure P = P (n) expressed by (43) and P | T =Tc = P | T =Tc (n) expressed by (44) are plotted on Figure 4.

Behavior of the system at temperatures T < T c
The equation of state in case of T < T c can be written in the form where E µ is defined in (22), and valuesρ 0i i = 1, 2, 3 are the three solutions of equation (32) [20] ρ 01 = 2ρ 0r cos α m 3 ,ρ 02 = −2ρ 0r cos α m + π 3 ,ρ 03 = −2ρ 0r cos M q is the value of the chemical potential M when the discriminant Q of equation (32) is equal to zero Note that the region of values of M is restricted to the interval |M | ≤ M q . For all |M | ≤ M q we have Q < 0, so there exist three real solutions of equation (32) in this interval of M . In case of |M | > M q we have a single root of equation (32), as well as in case of T > T c (42). At M = −M q we have Q = 0, so using (42) we obtain and at M = M q using (42) we can find So the equation of state at T < T c can be written in the form of several terms using the Heaviside functions θ in accordance with values of the chemical potential M This information is enough to plot the function P = P (τ, M ). Corresponding curves are shown on Figure 5. There we can see two surfaces touching to each other on the line M = 0. However, plot of the pressure as a function of density and temperature is more informative. To get this dependence we have to express the chemical potential as a function of density according to the interrelation (33) and use the obtained expression in (55). For this purpose let us use (41) and compare the behavior ofρ 0 as a function of chemical potential on the one hand and on the other hand as a function of density. Consider each of the regions depicted on Figure 6.
Region I, M ≤ −M q . At M = −M q we have expression (49) forρ 0 , where ρ 0r is a function of temperature (47). On the other hand, there is equality (33) which links the average densityn with the chemical potential (in the present case M = −M q ) and the valueρ 0 = −2ρ 0r . Thereby according to (33) we can find the densityn 12 (the point of transition between Region I and Region II), which corresponds to M = −M q n 12 = n g − 2κã 2 ρ 0r + a 4 3 ρ 3 0r .
Region III, 0 ≤ M ≤ M q . This region starts from M = 0, whereρ 01 = √ 3ρ 0r and respectivelȳ and ends up at M = M q , which takes the consequences withρ 01 = 2ρ 0r , thereforē Region IV, M > M q . This region starts fromn 03 (59), which decreases until some fixed valuē n lim . The valuen lim is bounded. Plot ofn 20 andn 03 as functions of temperature is shown on Figure 7 (the binodal curve). Maximal value of M corresponds ton lim = 10.01. In case ofn > n lim the solution m 1 in (41) has to be replaced by m 3 . However this solution is non-physical, because it gives a decrease of the chemical potential with an increase of the density.
Having a temperature dependence of the boundary densitiesn 12 ,n 20 ,n 03 andn 34 one can rewrite the equation of state (55) in terms of (τ,n). For this purpose it is enough to apply the formula (41) in equation (55). This formula gives dependence of the chemical potential on temperature and average density. The boundary densities described above have to be substituted for corresponding values of the chemical potential in Heaviside functions Θ(M i − M q ) in (55). As a result the equation of state is expressed as follows The pressure P as a function of density and temperature is plotted on Figure 8. Easy to see, that in temperature region τ < 0 there is a jump of density depicting a coexistence curve. Numerical results show that one can observe a liquid phase at T < T c in the density regionn ≤ n lim .
On the diagram of state ( Figure 9) one can see the line of some limited pressure P lim = lim n→n lim P (n) connected to the maximal value of the fluid density n lim . This boundary value depends on the microscopic parameters of the system in other words it takes specific value for each system. The line of maximal density may appear as a consequence of using the ρ 4 -model (see (21) which is the approximation of (20)). This line sets the range of densities 0 ≤ n ≤ n lim applicable for description of a fluid behavior in this manuscript.  Black line between the gray and the white zones of the surface of pressure as a function of temperature and average density ( Figure 10) corresponds to the region of temperatures close to the critical T c , where the mean-field approximation is not applicable.

Conclusions
To describe a first-order phase transition in the cell fluid model we deal with calculation of the grand partition function. A rigorous functional representation of the GPF (20) is obtained, which is reduced to the ρ 4 -model (21). This expression has to be calculated in two stages. The former is described in this paper. It covers a wide range of temperatures except the critical point. The latter stage concerns a narrow vicinity of the critical point. It requires the consideration of contribution from all the variables ρ k with k = 0 within the field-theoretical approach.
In the present work we were restricted to use a simplest approximation, which is called the approximation of a mean-field type. Thereby determination of the GPF is led down to computation of a single integral using the Laplace method. As a result we obtained an explicit form of the GPF as a function of temperature and chemical potential in this approximation. The pressure as a function of temperature and average density (equations (24) and (25)) is found. Figure 6 is important in respect to this. It shows correspondence between the regions of chemical potential and density at temperatures below the critical one T c . The critical temperature is defined by (29). Referring to Figure 5 easy to see that transition of the effective chemical potential M across the value M = 0 changes the behavior of the pressure P at T < T c . Here we have lim P , however, derivatives of pressure on the left hand side and on the right hand side from M = 0 are different. The latter fact causes a jump of density starting from n 20 to n 03 (see Figure 6). Each of these values are functions of temperature. In case of T T c we do not observe such a jump (formally n 20 = n 03 ).
A special attention is paid to the equation of state. Both the coexistence curve and the spinodal are plotted in the case of particular parameter p (see Appendix C), which take specific value for each system. Analyzing the surface of equation of state (see Figure 10) we noticed that curves of pressure increase sharply in the region of high densities (which exceed the liquid phase density). Apparently, the ρ 4 -model has some boundary density n lim (see Figure 2) such that there is a range of physically existent densities 0 n n lim . To examine behavior of the system in the high density regionn n lim (see Figure 9) it might be essential to consider (20) in the approximation of a ρ 2m -model where m ≥ 3.
The forthcoming research is calculation of the GPF (21) using the renormalization group method [27] with correlation effects taken into account. Thereby it will be possible to describe the behavior of a simple fluid in wide temperature region including the critical point within a sole approach.
As it can be seen from the latter expression one can integrate over the coordinates of particles in the system. For this purpose let us consider an integral over dx 1 in detail A particle with the coordinate x 1 can get into one cell, for example into ∆ l , and cannot get into any other cell, so The rest of integrals over each of x i give the same contribution. Now the GPF can be written in the form Here √ 4πp e −pN 2 is the result of integrating over ϕ l : .

B
Let us find an "entropy" part of (18) (the Jacobian of transition). To ease the process apply where the functions T n (v, p) are rapidly convergent series for all p > 0. As it can be seen from (8)  .
Therefore we obtain both the minimal and the maximal density = 0, allowing us to calculate respective value of the parameter v (volume of elementary cell) for each given parameter p (see Figure 11). Using Figure 12, it is easy to make sure, that not all the values of v in combination with the parameter p would give the correct picture.
p v Figure 11: Plot of v as a function of p. Curves 1, 2, 3 plotted on Figure 12 match p = 0.4 and three different values of v. Curve 2, which is pertinent to p = 0.4 and v = 20.8689, reaches the least value atn = 0. Increasing the value of v at fixed p the least m 1 is reached forn < 0 (curve 3 for v = 31). Decreasing the value of v the least m 1 is reached forn > 0 (curve 1 for v = 11) Thus there is a clear criterion for selection of values of the parameters p and v, which are the arguments of the special function T n (p, v) (see Appendix B). The former p is defined by the formula (8). This parameter varies upon a type of investigated matter and depends on the relation R 0 /α, where R 0 and α characterize the potential of interaction (4). The latter parameter v provides a description of particular density region of the systemn (0 ≤n ≤n lim , wheren lim is the limited density of a liquid phase) and is unambiguously defined by values of the parameter p.
Moreover p and v have to meet the condition a 4 > 0. This condition provides convergency of such integrals as (26).
This result is very important, since all parameters of the model are uniquely linked with each other. So making estimations for numeric values of the coefficients a n expressed by (B) we set the value of parameter p, then find the corresponding value of v using the equation n min = 0. Taking into account (29) expression (8) obtains the following form p = χΨ 0 /(2ã 2 W (0)(T c )).
Using the later one can find the value of χ, which meets the condition W (k) > 0. Being dependent on p the parameter χ becomes a function of microscopic parameters of an interaction potential, particularly of the relation R 0 /α. In case of the interaction potentialŪ l 12 = Ψ l 12 − U l 12 with the Fourier transform (11) at R 0 /α = 3 ln 2 we have the following set of parameters p = 0.4, v = 20.8689, χ = 3.5947. (C)