Phase diagram of the restricted primitive model: charge-ordering instability

We study the phase behaviour of the restricted primitive model (RPM) using a microscopic approach based on the method of collective variables with a reference system. Starting from the Hamiltonian of the RPM we derive the functional of the grand partition function given in terms of the two collective variables: the collective variables ρk and ck describing fluctuations of the total number density and charge density, respectively. Within the framework of the Gaussian approximation we found the boundary of stability with respect to fluctuations of the charge density. It is shown that due to the approximated character of the theory the boundary of stability is very sensitive to the particular choice of the long-range part of potential inside the hard core. This point is discussed in more detail.


Introduction
For the last decade, both the phase diagrams and the critical behaviour of ionic solutions have been intensively studied using experimental and theoretical methods.These studies were stimulated by controversial experimental results, demonstrating three types of critical behaviour in electrolyte solutions: classical (or mean-field) and Ising-like behaviour as well as crossover between the two [1][2][3].In accordance with these peculiarities, ionic solutions were conventionally divided into two classes, namely: "solvophobic" systems in which Coulomb forces are not supposed to play a major role (the solvent is generally characterized by high dielectric constant), their critical behaviour is Ising-like; "Coulombic" systems in which the phase separation is primarily driven by Coulomb interactions (the solvent is characterized by low dielectric constant).
The basic model for theoretical studies of ionic systems, including electrolyte solutions, molten salt and colloids, is the so-called "restricted primitive model" (RPM), in which the ionic fluid is modelled as an electroneutral binary mixture of charged hard spheres of equal diameter immersed in a structureless dielectric continuum.Early studies [4][5][6] established that the model has a gas-liquid (GL) phase transition.The following main problems are of great interest for the RPM: 1) the location of the gas-liquid (GL) critical point; 2) the study of thermodynamic behaviour in the vicinity of the GL critical point; 3) the investigation of the full phase diagram; 4) the study of transition to the charge-ordered phase.Concerning the first issue, qualitative estimations of the critical temperature and the critical density have been found in the framework of the Debye-Hückel theory [7] as well as by the integral equation methods [8,9].In recent years the GL critical point has been much studied using Monte Carlo (MC) simulations in the grand canonical ensemble [10,11].Concerning the second issue, strong evidence for an Ising universal class is found by computer simulations [12] as well as by a recent theoretical study [13].However, a description of a crossover region when the critical point is approached still remains an open question.
Recent theoretical [14][15][16][17] and simulation studies [18][19][20][21][22] of the RPM revealed its rich phase behaviour.Moreover, the lattice version of the RPM (LRPM) shows a phase diagram that is rather different from the continuous RPM and the phase diagram namely, its high-density branch appears to be highly sensitive to changes of lattice structure [23].When constructing a full diagram of the RPM, the stability analysis showed [14,15] that another phase transition can occur in the model between the charge-disordered and charge-ordered phases along the so-called λ-line.The λline has been also found in earlier theories for the continuous RPM ( see [15] and references herein), but no λ-line has been found so far for this model in computer simulations.Thus, the problem of charge-density instability is still open for the discussions and needs further investigations.
In this paper we address an issue of a phase instability connected with chargeordering fluctuations, using the functional representation of the grand partition function in terms of collective variables (CVs).Note that the method of CVs, proposed initially in the 1950s [25,26] for the description of the classical charged many particle systems and developed later for the needs of the phase transition theory [27][28][29][30] was in fact one of the first successful attempts to attack the problems of statistical physics using the functional integral representation.At the same time, other functional approaches being based on the Stratonovich-Hubbard transformation [31,32] were originated.Both groups of theories are in fact in close relation and we are going to discuss this problem more in detail elsewhere [33].

Functional representation: 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 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 determined from the equation For the RPM we have µ + = µ − .Now we separate the interaction potential U γδ (r) into short-and long-range parts: where ψ γδ (r) is a potential of a short-range repulsion and Φ γδ (r) is a long-range part of the potential.At this stage of consideration we hit the problem of a proper treatment of the potential Φ γδ (r) for physically inaccessible regions.Evidently, the properties of the system are independent of the value of the potential Φ γδ (r) at the intermolecular distances smaller than σ.However, since we cannot calculate the partition function exactly, the approximate results will be sensitive to the concrete choice of the long-range part of the potential inside the hard core.Several choices of a particular form for Φ γδ (r) in the region r σ will be analyzed in the later part of this paper.Within the framework of the approach considered, the short-range repulsion part of the potential ψ γδ (r) is described in the space of the Cartesian coordinates of the particles.We call the hard sphere system with the diameter σ a reference system (RS).The interaction connected with an attraction (potential Φ γδ (r) ) is considered in the CV space.The transformation from the Cartesian coordinates to the CVs is performed by means of the transition Jacobian.
Using the method of CVs, developed for a two-component continuous system [34,35] we can rewrite the grand partition function of the RPM in the following form: Here Ξ 0 is the grand partition function of the RS.Φ C (k) is the Fourier transform of Φ ++ (r) (= Φ −− (r)), and ρ k and c k are the CVs describing total density and charge density fluctuations, respectively: The indices c and s denote the real and imaginary parts of CVs ρ k and c k .Each of ρ c k (c c k ) and ρ s k (c s k ) takes all the real values from −∞ to +∞, and (dρ) and (dc) are volume elements of the CV phase space: J(ρ, c) is the Jacobian of the transition to the CVs averaged over the RS.For the RPM, J(ρ, c) is of the same form as that for the symmetrical binary fluid [35]: M(in and variable ω k (γ k ) is conjugate to CV ρ k (c k ).Note that all odd terms with respect to γ k are equal to zero in (5)  are expressed as linear combinations of the nparticle partial structure factors of the RS S γ 1 ...γn (k 1 , . . ., k n ) and are calculated for γ 1 , . . ., γ n = +, − and n 4 in [30] (see Appendix B in [30]).
In general, the dependence of M (in) n (k 1 , . . ., k n ) on wave vectors k 1 , . . ., k n is rather complicated.However, since we are interested in the critical behavior, the small-k expansion of the cumulant can be considered.Hereafter we shall replace M (in) n (k 1 , . . ., k n ) by their values in long wavelength limit M (in) n (0, . . ., 0), where they can be expressed via well defined thermodynamic quantities and their derivatives.
Let us present J(ρ, c) as In (7) the linear term proportional to M (0) 1 is eliminated by the shift ρ k = ρ k + M(0) 1 δ k (the prime on ρ k is omitted for clarity).Cumulants of the second order have the forms: where S 2 (k) is a two-particle structure factor of a one-component hard sphere system.In the Carnahan-Starling approximation one has for S 2 (k = 0): Since M(0) 2 (0) is the positive and smooth function in the region under consideration and M(2) 2 is equal to constant, we can integrate in (7) over ω k and γ k using the Gaussian density measures as basic ones.This integration can be performed by several ways.For example, using the Euler equations we can determine ω * k (and γ * k ) which provide the maximum for the functional in the exponent in (7): As a result, we can present Ξ in the form [13]: where and the Hamiltonian H is given exclusively in terms of the total density fluctuation modes and charge density fluctuation modes: Here superscript i n indicates the number of variables c k at a (in) n .Taking into account the first terms in (8), coefficients a (in) n can be written as follows: M( 4) Another way of integration in ( 7) is outlined in Appendix.

Charge-ordering phase instability
Expression (11) defines indeed the Landau-Ginzburg-Wilson Hamiltonian H for the RPM with two fluctuating collective variables {ρ k } and {c k } describing the number and charge densities.Based on this one can analyze both the liquid-gas and charge-ordering phase transitions using standard methods developed in the phase transition theory [36].
Let us consider the coefficients a (0) 2 (at the second power of CV ρ k ) and a (2) 2 (at the second power of CV c k ) in (11).As is seen from ( 12), the coefficient a (0) 2 never reduces to zero for physical values of the density.The fact that the RPM does not demonstrate the GL phase instability in this approximation, is attributed to the absence of direct pair interactions of density fluctuations in the model as well as to the neglect of the effects of non-direct correlations via a charge subsystem at this level.In order to obtain the GL spinodal curve for this system we should take into consideration the terms of the higher order [14,37].This task was accomplished in [37] using the approach described and the calculated GL spinodal curve is shown in figure 1 (the curve with the maximum).[37] The curve with the maximum is the GL spinodal, the maximum is located at T * c = 0.0502 and η c = 0.022.The straight line calculated by (22) corresponds to the boundary of stability of charge-disordered phase.

By contrast, the coefficient a
(2) 2 can be zero for the certain choices of the value of potential Φ γδ (r) inside the hard core.Namely, the equality can hold at some values of the wave-vector k, temperature and density.To demonstrate this fact we consider several different choices of Φ γδ (r) at r σ.
First we separate the initial potential U γδ (r) into short-and long-range parts according to the Weeks-Chandler-Andersen partition [38]: This simple form for Φ γδ (r) inside the hard core changes the behaviour of the transform for large k from usual Coulomb k −2 to k −3 decay.As was shown [39], this choice of Φ γδ (r) for r σ produces rapid convergence of the series of the perturbation theory for the free energy.The Fourier transform ΦC (k) has the form where Substituting ( 18) into (15) we obtain the equation for the boundary of stability (a spinodal) connected with the charge fluctuations: or where y * is determined from the condition which yields y * 4.0783.Substituting y * in (19) we obtain the boundary of stability with respect to fluctuations of the local charge density The stright line in figure 1 corresponds to equation (22).A similar result but for another choice of interaction inside the hard core was obtained in [14,15] within the framework of the field-theoretical approach.In that case the potential ψ γδ (r) has the form ( 16) and the potential Φ γδ (r) is equal to zero in the region r σ.
This choice gives βρ ΦC (y) = 24β * η(cos y/y 2 ) and the equation for the boundary of stability with respect to charge density fluctuations has the form: It is worth noting that the boundary of stability of the charge-disordered phase determined by ( 22) and ( 23) is associated with finite-wavelength fluctuations with k * = y * /σ.As was discussed in [15], the terms of the higher order in (11) can lead to a significant shift and even disappearance of a λ-line.For the LRPM, following the outline of Brazovskii theory [41] Ciach studied the effect of fluctuations on the charge-ordered-charge-disordered instability and showed that a fluctuation-induced first-order phase transition should be expected instead of continuous transition found in the mean-field approximation [42] .
Consider now another choice of Φ γδ (r) for r σ determined from the so-called optimized random phase approximation (ORPA).It is worth noting that the ORPA is equivalent to the mean spherical approximation (MSA) provided the reference system is approximated by the Percus-Yevick theory.Based on the exact solution of the MSA for the RPM [43] Anderson and Chandler [44] found the following form of Φ γδ (r) inside the hard core: where the following notations are introduced: and κ is inverse Debye length.In this case one can find for the Fourier transform of the full potential Φ γδ (r) where and y = kσ.Although for (25) the l.h.s. of equation (15) approaches zero very closely at some values of κ and y, it does not become equal to zero in this case.Thus, for the Anderson and Chandler partition (24) the charge-ordering instability does not appear in the Gaussian approximation.As a particular case of the Chandler-Anderson regularization (24) one may consider the form for Φ γδ (r) inside the hard core derived by Caillol [45] in the way of minimization of the mean field functional for free energy with respect to variations of the smearing functions.In that case the Fourier transform of the potential Φ γδ (r) is always positive, so that no real solution of equation ( 15), connected with the existence of charge-ordering instability can be found.
Based on the results obtained we can suggest that the charge-ordered-chargedisordered phase transition does not take place in the continuous RPM in fluid state.In order to confirm this suggestion for some regularized potential discussed above the fluctuation terms of the order higher than the second one should be taken into account.

Discussion
In this paper we have studied the phase instability of the RPM connected with charge density fluctuations.For this purpose the method of CVs with a reference system was used.First we present the functional of the grand partition function given in terms of the two sets of CVs: {ρ k } and {c k } describing fluctuations of the total number density and the charge density respectively.The resulting Landau-Ginzburg-Wilson Hamiltonian has the form of an infinite series in powers of ρ k and c k .Based on the study of the Gaussian approximation we analyse the phase instabilities of the model.For the RPM this approximation does not produce the GL phase instability.In order to describe the GL phase transition, the terms of the higher order should be taken into account in the effective Hamiltonian [14,37].Actually, the charge-charge correlations cause the effective attraction between the ions which, in turn, leads to the GL phase transition in the RPM.However, within the framework of this approximation we could estimate the boundary of stability with respect to fluctuations of the charge density.As was shown, the boundary of stability is very sensitive to the choice of the particular form of the long-part of potential inside the hard core.
We consider several different forms of the potential Φ γδ (r) in the physically inaccessible region.On the one hand, within the Gaussian approximation we confirm the results obtained within the framework of the field-theoretical approach [14,15].Although equations (22) and (23) give the λ-lines with very different slopes, these two simple forms for the long-range part of the potential inside the hard core produce spinodal lines connected with the charge-ordering instability.On the other hand, the result obtained with (25) confirms the well-known fact that the MSA does not yield such type of instability.Similar conclusion is made for the regularized potential, proposed in [45].However, it is worth noting that simulation results demonstrate the λ-line for the LRPM [18] in contrast to the MSA.
The possibility of the charge-ordering transition in the RPM was also discussed in [15,23].The authors, comparing the form of the schematic full phase diagram of the LRPM (with the certain degree of space discretization) obtained in their approach with that obtained by computer simulations [20], suggested that a fluctuationinduced first-order phase transition of order-disorder type can be realized in the system.In that case the λ-line describes the boundary of stability with respect to a formation of ionic crystal [23].

Figure 1 .
Figure1.The phase diagram of the RPM obtained in[37] The curve with the maximum is the GL spinodal, the maximum is located at T * c = 0.0502 and η c = 0.022.The straight line calculated by(22) corresponds to the boundary of stability of charge-disordered phase.