Phase behaviour of the restricted primitive model

The phase behaviour of the restricted primitive model (RPM) is studied using a microscopic approach recently developed for the description of phase transitions in binary fluid mixtures. For the model we obtain the explicit expression for the functional of the grand partition function. Based on the functional we calculate the phase diagram of the RPM in the high temperature approximation (HTA) and then we do this calculation taking into account the terms of the higher orders in the effective Hamiltonian. In both cases the phase diagrams demonstrate the gas-liquid (GL) and charge ordering phase instabilities. In the latter case, the obtained value for the GL critical temperature is in good agreement with the MC simulation data whereas the critical density is underestimated.


Introduction
In recent years, much attention has been focused on an issue of the critical and phase behaviour of ionic fluids.For reviews of the experimental and theoretical situation see [1][2][3][4][5][6][7].
The simplest model capable of capturing the main features of ionic systems is the restricted primitive model (RPM) [3,5].Early studies [8] established that the model has a gas-liquid (GL) phase transition.A reasonable theoretical description of the GL critical point in the RPM was accomplished at a mean-field (MF) level using integral equation methods [5,9] and Debye-Hückel theory [10].
The universality class of criticality in the RPM has also been a subject of active research [20][21][22][23][24]. Very recent simulation studies have found strong evidence for Ising universal class [19,25].
The one more relevant issue concerns the full phase diagram of the RPM.The stability analysis shows that another phase transition occurs in the continuum RPM between the disordered and charge-ordered phases along a λ-line [7,26,27].But no λ-line was found so far in computer simulations for this model.
In this paper we address the issue of the critical and phase behaviour of the RPM using the theoretical approach proposed in [29,30] for the binary mixture.The theory has its origin in the approach based on a functional representation of a partition function by means of the collective variables (CV) method [31,32].Its particular feature is a choice of the phase space in which the system is considered.This phase space is formed by a set of CV and contains a variable connected with the order parameter.The approach allows one to determine, on microscopic grounds, the explicit form of the effective Ginzburg-Landau-Wilson (GLW) Hamiltonian and then to integrate the partition function in the neighborhood of the phase transition point using the non-perturbative renormalization group method [33].As a result, nonclassical critical exponents and analytical expressions for thermodynamic functions were obtained [33,34].More recently this theory was developed for a binary fluid mixture [29,30,[35][36][37][38].
The paper consists of two parts.In the first part we obtain the functional of the grand partition function of the RPM given in terms of the CV c k (connected with charge density fluctuation modes) and in terms of the variables γ k and h k=0 conjugate to the CV c k and ρ k=0 , respectively (ρ k=0 is connected with the k = 0 mode of total number density fluctuations).Restricting our consideration to the second powers of γ k (taking into account the higher powers of h k=0 ) we derive the equation for the chemical potential.Based on the chemical potential obtained from the linearized equation we calculate the spinodal curve.Its run suggests that two types of phase instabilities can occur in the RPM.One part of the spinodal is of the gas-liquid (GL) type while another one looks like a λ-line.We obtain the following values of the GL critical point: T * c = 0.084 and η c = 0.005 which agrees with the other MF theories [9] and we discuss the approximation used at this stage of our study.
In order to study the nature of the criticality of the RPM as well as to get the best estimates for its GL critical point we go beyond the above mentioned approximation taking into account the terms of higher orders in the effective Hamiltonian.The second part of the paper is devoted to this end.First, we consider the Gaussian approximation of the functional of the grand partition function.It yields the equation for the boundary of stability with respect to the charge density fluctuations.Then, applying the procedure proposed in [29,30] we obtain the expression for the grand thermodynamic potential in the vicinity of the GL critical point as a power series in the field hk=0 (up to h4 ) conjugate to the order parameter.The expression obtained has the form of the Landau free energy.We also calculate the phase diagram demonstrating both GL and charge ordering phase instabilities.The data for the GL critical point are found to be T * c = 0.0502 and η c = 0.022.

Functional representation of the grand partition function of the RPM
The RPM consists of N = N + +N − hard spheres of diameter σ with N + carrying charges +q and N − (= N + ) charges −q, in a medium of dielectric constant D. The interaction potential of the RPM has the form We split the potential U γδ (r) into short-and long-range parts using the Weeks-Chandler-Andersen partition [39].As a result, we have Φ γδ (r) = (q γ q δ )/Dσ, r σ (q γ q δ )/Dr, r > σ .
This simple form for Φ γδ (r) inside the hard core changes the behaviour of the Fourier transform for large k from usual Coulombic k −2 to k −3 decay.As was shown [40], this choice of Φ γδ (r) for r < σ produces rapid convergence of the series of the perturbation theory for the free energy.The Fourier transform of Φ γδ (r) = q 2 /Dr = Φ C (r) has the form where We start with the grand partition function for a two-component system (γ, δ = +, −): where (dΓ) = γ dΓ Nγ , dΓ Nγ = dr γ 1 dr γ 2 . . .dr γ Nγ (γ = +, −) is an element of the configurational space of the γth species; z γ is the fugacity of the γth species: µ γ is the chemical potential of the γth species.The second term in this expression is obtained as a result of integration over the momenta.m γ is the mass of the γth species, h is the Planck constant.µ γ is determined from the equation where N γ is the average number of the γth species in the grand canonical ensemble.The correlation effects of different scales are connected with the potentials ψ γδ (r) and Φ γδ (r): the potential ψ γδ (r) describes the behaviour of the particles at very short distances and provides their mutual impenetrability.The potential Φ γδ (r), on the contrary, describes an attraction between the particles which takes place at long distances.These effects are proportional to different parameters.In order to study them simultaneously it is necessary to consider the system successively in two phase spaces.First we write the grand partition function and the correlation functions of the RS in the phase space of the Cartesian coordinates.Here the pair interaction is described by the potential ψ γδ (r).The thermodynamic and structural properties of the RS are assumed to be known.Then the grand partition function of the full system is constructed in the phase space of the collective variables (CV) by means of the functions of the RS.
We introduce the grand partition function of the RS where µ 0,γ is the chemical potential of the γth species in the RS.Now we rewrite the attractive potential in the Fourier space where Φγδ (k) is the Fourier transform of Φ γδ (r).ρNγ (k) is the Fourier transform of the operator of the particle number density for the γth species by means of the relations [32,36] ρc sin(kr γ j ).
Here the indices c and s denote the real part and the coefficient of the imaginary part of both ρNγ (k) and ρ k,γ ; δ(• • •) is the Dirac delta function.Now we present Ξ in the following form: where Ξ 0 is the grand partition function of the RS.Then Ξ 1 can be written in the form of the functional integral [36] Ξ Here ) describes the value of the k-th fluctuation mode of the number of γ-th species particles and (dρ) is a volume element of the CV phase space: The prime means that the product over k is performed in the upper semi-space.µ 1,γ is a part of the chemical potential of the γth species which is determined from the equation J(ρ) = J(ρ + , ρ − ) is the Jacobian of the transition to CV averaged on the RS: Substituting the explicit forms for delta functions into the expression for J(ρ), we obtain where the variables ν k,γ are conjugate to the CV ρ k,γ and have the form and (dν) is a volume element of the phase space of the variables ν k,γ J(ν) has the form: We present J(ν) in the form of the cumulant expansion Here M γ 1 ...γn is the nth cumulant which is determined from the formula and is connected with S γ 1 ...γn (k 1 , . . ., k n ), the n-particle partial structure factor of the RS, by means of the relation Now we pass in these formulas to new variables ρ k and c k (according to ω k and γ k ) by means of the orthogonal linear transformations The variables ρ k and c k are CV connected with total density fluctuation modes and charge density fluctuation modes, respectively.
As a result, we can present the functional of the grand partition function of the RPM in the form: In the case of the RPM, the RS is a one-component hard-sphere system with the diameter σ (potential ψ γδ (r) = ψ(r)).For µ 1 we have (µ 1,+ = µ 1,− ): and equation ( 5) is rewritten as For the RPM J(ρ, c) has the same form as that for the symmetrical binary fluid with µ a = µ b [30]: In ( 9) the cumulants M (in) n with i n = 0 are connected with the nth structure factors of the RS [36]: Structure factors S n (0) with n > 2 can be obtained from S 2 (0) by means of a chain of equations for correlation functions [41].Cumulants with i n = 0 can be expressed in terms of M (0) n (see also formulae (4.8) in [36]): First we integrate in (6) over CV ρ k .This integration leads to the delta-functions As the result of integration over ω k we obtain where and condition (8) has the form: Expressions ( 11)-( 12) do not include the "field" variable h k with k = 0.As one can see below, this fact will give rise to the Landau type free energy of the RPM in the vicinity of the GL critical point.Formulas ( 11)-( 14) are the initial formulas in our study of the phase behaviour of the RPM.

Phase diagram of the RPM in the high temperature approximation (HTA)
Restricting our consideration in (12) to the second power of γ k and assuming M (2) n+2 h n 0 n! > 0, we first integrate in ( 11)-( 12) over γ k and then over c k .As a result, we obtain for Ξ 1 where the following notations are introduced: and formulas (10) are used for M (in) n .Using ( 14) and (15) we can obtain the following equation for the chemical potential where Phase diagram in the HTA.We solve equation ( 17) in the simplest approximation.Neglecting terms h 2 0 , h 3 0 etc. in the right hand side of ( 17) and setting in the left hand side of ( 17) we obtain for µ 1,+ (= µ 1,− ) On the other hand, we can get the same result for µ 1,+ taking into account the terms proportional to γ 2 and γ 2 h 0 in the exponent of ( 12) and setting M (2) 3 is the coefficient of γ 2 h 0 ).To this end we neglect in ( 11)- (12) the terms proportional to . In this case Ξ 1 has the form: Integrating in (22) over γ k and then over c k we obtain Figure 1.The phase diagram of the RPM calculated from (25) (see the text for explanation).
The full chemical potential µ + is equal to (see ( 4)): where µ 0,+ (= µ 0,− ) is the chemical potential of a one-component hard sphere system.The equation βρ ∂µ + ∂ρ = 0, where ρ is the total number density, gives the spinodal of the RPM in the approximation considered.Using (2) for ΦC (k) and the Percus-Yevick approximation for the RS, the equation for the spinodal curve can be written as The phase diagram calculated from ( 25) is shown in figure 1.The obtained data for the GL critical point are T * c 0.084 (T * c = 1/β * c ) and η c 0.005.It is evident that in the considered approximation we get the overestimated value for the GL critical temperature and the underestimated value for the critical density.But, in contrast to the previous results, the spinodal curve changes its run (at η 0.047) and then it directs to the higher temperature.The second positive slope of the spinodal indicates another type of the phase instability appearing in the RPM (similarly to the λ-line in a symmetrical binary fluid).Below, we shall calculate the phase diagram of the RPM taking into account in (12) the terms of the higher orders.

Phase diagram of the RPM: beyond the HTA
As was shown above, the spinodal (25) can be obtained as the result of the trick, namely: the integration in (12) was accomplished taking into account the term proportional to γ 2 h 0 and in the final result the coefficient at γ 2 h 0 was set equal to zero.Below we will sequentially take into consideration in (12) the terms of the higher order.Now let us rewrite ( 11)-( 12) as and the notations are the same as those in previous section.
Charge ordering phase instability.First, we restrict ourselves to the Gaussian approximation which corresponds to neglecting the terms proportional to 0 , etc. in the exponent of (26).After integration in ( 26) over γ k we obtain As is seen from ( 27), the equality holds at some values of the wave-vector k, temperature and density.Equation ( 28) determines the boundary of stability connected with the charge fluctuations (the field variable γ k is conjugate to the CV c k ): or where x * is determined from the condition which yields x * 4.0783.Substituting x * in ( 27) we obtain the boundary of stability with respect to fluctuations of the local charge density A similar result (for another choice of interaction inside the hard core) was obtained in [7,26] within the framework of the field-theoretical approach.The possibility of the charge-ordering transition in the continuous-space RPM model was also discussed in [27,28].
It is worth noting that the RPM does not demonstrate the GL phase transition in the approximation given by (27).In order to obtain the GL spinodal curve we should take into consideration the terms of the order higher than the second one (γ 2 h, γ 2 h 2 , etc.).

GL phase instability.
We restrict our consideration in (26) to the terms of the fourth order.In this case Ξ has the form: where the summation in (16) over n is restricted to 4. Now we follow the programme proposed in [29,36] for a two-component fluid system.First, we separate the two types of variables: the essential variables (which include the variable connected with the order parameter) and the non-essential variables.Then, integrating over the non-essential variables with the Gaussian density measure, we construct the basic density measure (the GLW Hamiltonian) with respect to the essential variables.
For the RPM in the vicinity of the GL critical point the variable h 0 (conjugate to the CV ρ 0 ) turns out to be the essential variable [30].Thus, we can present (33) as where The right hand side in (38) has the form of the Landau free energy expressed in terms of the field h0 conjugate to the order parameter.From the equation M2 = 0 we obtain the equation for the GL spinodal curve where S n (0) is the nth structure factor of the one-component hard-sphere system at k = 0.The phase diagram of the RPM is shown in figure 2. The curve with the maximum is the GL spinodal calculated using (39).The Percus-Yevick approximation is used for S 2 (0) (the expressions for S 2 (0), S 3 (0) and S 4 (0) are given in Appendix).The straight line calculated by (32) (the Gaussian approximation) corresponds to the charge ordering phase transition.The GL critical point (the maximum of the GL spinodal) is located at T * c = 0.0502 and η c = 0.022.While the value for T * c is in good agreement with the recent data of computer simulations [15,18] (T * c 0.05), the critical density is underestimated (ρ * c 0.04).The same phase diagram for the RPM was obtained in [42] by means of the Hubburd-Schofield method [43].
We can also obtain from (38) (using ( 14)) the expression for µ +,1 If we set M (0) 3 = 0 the expression (40) reduces to (20).It should be pointed out that the above described scheme of integration in the vicinity of the GL critical point cannot be used in the region close to the line defined by (32) (dotted line in figure 2): the variables γ k appear to be the essential variables in this region.

Conclusions
We use the recently developed approach in the study of the phase behaviour of ionic fluids.For the RPM we obtain the functional of the grand partition function given in terms of the CV c k , which are connected with charge density fluctuation modes, and in terms of the variables γ k and h k=0 conjugate to the CV c k and ρ k=0 , respectively (CV ρ k=0 is connected with the k = 0 mode of total number density fluctuations).
We conclude that within the framework of the same approach the following results are obtained for the RPM: • the phase diagram in the HTA which consists of two parts: the part corresponding to the GL phase transition and another one looking like the λ-line, the parameters of the GL critical point correlate with those obtained by other MF theories; • the explicit expression for the grand thermodynamic potential in the vicinity of the GL critical point as an expansion in terms of the field h0 conjugate to the order parameter (up to h4 0 ).The expression suggests the classical critical behaviour which disagrees with the findings of the recent simulation studies [19,25].Thus within the framework of this approach the criticality in the RPM needs further investigations.
• the GL spinodal in the approximation when the terms of the higher order (h 0 , γ k γ −k h 0 , γ k γ −k h 2 0 and h 4 0 ) are taken into account in the Hamiltonian, the value obtained for the GL critical temperature agrees well with MC simulation data; • the boundary of stability with respect to the charge density fluctuations (in the Gaussian approximation of the functional integral).This result confirms that obtained in [26].As was discussed in [27,28], in continuous systems a fluctuation-induced first-order order-disorder transition can occur when the instability is associated with fluctuations characterized by k = 0.
It is worth noting that this approach can be applied to the case when both longrange and short-range interactions are involved into the model simultaneously (the RPM+SR model).This task will be considered elsewhere.

Figure 2 .
Figure 2. The phase diagram of the RPM calculated from (37) (see the text for explanation).