Gas-liquid phase equilibrium in ionic fluids: Coulomb versus non-Coulomb interactions

Using the collective variables theory, we study the effect of competition between Coulomb and dispersion forces on the gas-liquid phase behaviour of a model ionic fluid, i.e. a charge-asymmetric primitive model with additional short-range attractive interactions. Both the critical parameters and the coexistence envelope are calculated in a one-loop approximation as a function of the parameter $\alpha$ measuring the relative strength of the Coulomb to short-range interactions. We found the very narrow region of $\alpha$ bounded from the both sides by tricritical points which separates the models with"nonionic"and"Coulombic"phase behaviour. This is at variance with the result of available computer simulations where no tricritical point is found for the finely-discretized lattice version of the model.


I. INTRODUCTION
Critical and phase behaviour of ionic fluids has been intensively studied for the last decades. These studies were stimulated by earlier experiments on ionic solutions that yielded three types of critical behavior: mean-field and Ising-like behavior as well as crossover between the two [1,2]. In accordance with these peculiarities, ionic solutions were conventionally divided into two classes, namely: "solvophobic" systems with Ising-like critical behavior in which Coulomb forces are not supposed to play a major role (the solvent is generally characterized by high dielectric constant) and "Coulombic" systems in which the phase separation is primarily driven by Coulomb interactions (the solvent is characterized by low dielectric constant) [2,3]. The critical behaviour of the Coulombic systems has been a subject of great debates. At present, the analysis of the existing experimental data for Coulombic and solvophobic systems has shown that the asymptotic critical behaviour is Ising-like for both classes [4,5]. Experiments that supported the expectation of mean-field critical behaviour could not be reproduced in later works [6]. Strong evidence for the Ising universal class for the Coulombic fluids has been also found by recent simulations [7][8][9][10] and theoretical studies [11][12][13]. More recent accurate experiments indicate a crossover to the mean-field behaviour with crossover temperatures much closer to the critical temperature than it is observed in nonionic systems [5,14,15]. A microscopic explanation of such behaviour has been recently given in [16,17].
In this paper we continue our systematic study of the phase behaviour of ionic fluids.
Our aim is to study the effect of the strength of electrostatic interactions (the dielectric constant of the solvent) on the phase separation in ionic solutions. It should be noted that ionic and nonionic ("classical") fluids differ greatly in the strength of the interactions. In simple neutral fluids, the interaction energy is of the order of the thermal energy, i.e., the reduced temperature T * = k B T /φ min ≃ 1 (φ min is the depth of the interaction potential at its minimum). By contrast, in typical molten salts, the interionic interactions are more than one order of magnitude larger than k B T , i.e., T * ≪ 1. In electrolyte solutions, the strength of the Coulomb interactions depends on the dielectric constant of the solvent. In addition, the phase diagrams of Coulomb dominated fluids are quite asymmetric when compared to nonionic fluids [18][19][20].
The evolution of the gas-liquid phase diagram of the charged Yukawa fluid with an increase of ionic interactions was studied in [21] using both the mean-spherical approximation and the Gibbse ensemble Monte Carlo simulations. It has been found that the gas-liquid coexistence envelope changes from a "nonionic type" (such as that of the one-component Yukawa fluid) to a "Coulombic type "(such as that of the primitive model) when the strength of the Coulomb interactions is increased. However, the issue of the watershed between the nonionic and ionic phase behaviour (the "solvophobic" and "Coulombic" systems) has not been addressed so far.
Here, we address this issue using a simple model of ionic fluids, i.e., the charge-asymmetric primitive model (PM) supplemented by short-range attractive interactions where the shortrange attraction is considered as an approximation to van der Waals interactions [16,[21][22][23].
The model without the Coulomb interactions (a hard-sphere square-well model) exhibits a gas-liquid coexistence typical of nonionic ("solvophobic") fluids. Another limiting model, i.e., the PM, demonstrates the "Coulombic" phase diagram. A theoretical background for this study is the statistical field theory that exploits the method of collective variables (CVs) [24][25][26][27]. Using this theory we have recently derived a microscopic-based effective Hamiltonian for the model with short-and long-range interactions [16]. Here, we use this Hamiltonian to study the gas-liquid phase behaviour of the model with parameters ranging from the purely non-Coulombic regime to the purely Coulombic regime.
The paper is arranged as follows. A theoretical background is given in section 1. The results for the gas-liquid critical parameters and phase diagrams are presented in section 3.
Concluding remarks are made in section 4.

A. Model
We start with a two-component model of ionic fluids. The pair interaction potential is assumed to be of the following form: where φ HS αβ (r) is the interaction potential between the two additive hard spheres of diameters σ α and σ β . The potential φ SR αβ (r) describes the short-range (van-der-Waals-like) attraction. φ C αβ (r) is the Coulomb potential: φ C αβ (r) = q α q β φ C (r), where φ C (r) = 1/(ǫr) and ǫ is the dielectric constant of the solvent. The ions of the species α = 1 carry an electrostatic charge q + = q 0 and those of species α = 2 carry an opposite charge q − = −zq 0 (q 0 is an elementary charge and z is the parameter of charge asymmetry). Overall charge neutrality requires that q + N + + q − N − = 0, where N + and N − are, respectively, the number of positive and negative ions. In general, the binary system of hard spheres interacting via the potential φ SR αβ (r) can exhibit both the gas-liquid and demixion phase transitions. Thus, we simplify the model assuming that (i) the hard spheres are of the same diameter σ α = σ β = σ and . With these restrictions, the uncharged system can only exhibit a gas-liquid phase separation and a possible demixion is ruled out. Then, we specify φ SR (r) in the form of the square-well (SW) potential of depth ε and range λ. The system of hard spheres interacting through the SW potential with λ = 1.5σ can serve as a reasonable model for simple (nonionic) fluids.
It is worth noting that in the treatments of models with hard cores, the perturbation potential is not defined uniquely inside the hard core. Here, we use the Weeks-Chandler-Andersen regularization scheme for the both potentials φ C (r) and φ SR (r) [28]. In this case, the Fourier transforms of these potentials have the form: and where x = kσ .

B. Effective Hamiltonian near the gas-liquid critical point
We consider the model (1) : where ρ k,N is the collective variable (CV) which describes the value of the k-th fluctuation mode of the total number density. In a one-loop approximation, the coefficients of Hamiltonian (4) have the form : a n,0 = −ρ n−1 C n,HS − ρ n−1 C n,C , n ≥ 3.
Here, we introduce the following notations. The subscript HS refers to the hard-sphere system and the superscript (2) in equation (7) denotes the second-order derivative with the respect of the wave vector k: The addend ∆ν N in Eq. (5) is related to the chemical potentials whereν α is determined byν ν α is the dimensionless chemical potential, ν α = βµ α −3 ln Λ α and µ α is the chemical potential of the αth species; β = 1/k B T is the reciprocal temperature; Λ −1 α = (2πm α β −1 /h 2 ) 1/2 is the inverse de Broglie thermal wavelength.
C n,HS is the Fourier transform of the n-particle direct correlation function of a onecomponent hard-sphere system at k = 0, ρ = N /V = ρ 1 + ρ 2 is the total number density.
Explicit expressions for C n,HS and C 2,HS for n ≤ 4 in the Percus-Yevick (PY) approximation are given in reference [16] (see Appendix in [16]).
The last term on the right-hand side of equations (5)-(8) results from the charge-charge correlations being taken into account. It takes the form where a 1,0 is the excess part of the chemical potential ν N connected with the short-range attractive and long-range Coulomb interactions. Equation a 1,0 = 0 yields the expression for the chemical potential in a one-loop approximation (the random phase approximation). All the coefficients in equation (4) except a 2,2 describing the square-gradient term can be derived from the one-loop free energy as the derivatives with respect to the density [16].
It is worth noting that there is no difference between the charge-asymmetric PM and the restricted primitive model (RPM) having z = 1 at this level of approximation [29][30][31][32].
Hereafter, we briefly refer to the model (1)-(3) as a RPM-SW model.

III. GAS-LIQUID PHASE DIAGRAM
Here, we study the gas-liquid phase diagram of the the RPM-SW model for the model parameters ranging from the RPM limit to the SW limit. Using equations (2)-(3) and (9)- (12) we rewrite the expressions for the coefficients of the effective Hamiltonian (see (5)- (8)) as follows: 2,HS + 4 5 where [16] i κ = κ D σ with κ 2 D = 24η/(T C σ 2 ) being the Debye number. In equations (13)- (18), T C is the reduced temperature defined as the ratio between the thermal energy k B T and the Coulomb energy of the opposite charged hard spheres at contact, E C = zq 2 0 /ǫσ, α is the ratio of the Coulomb and square-well energies at contact and η = πρσ 3 /6 is the packing fraction of the ions. For ε being fixed, α measures the strength of the Coulomb interactions. In this case, either the increase of ion charges or the decrease of dielectric constant for fixed charges results in the increase of the parameter α. It is worth noting that a n,0 with n ≥ 3 is independent of α when the reduced temperature is given by equation (19). The approximation (13) figure 2, the dependence of the reduced critical temperature given by (19) on the parameter α −1 ∼ ǫ is shown. As is seen, T C c increases almost linearly with an increase of the dielectric constant. The trend of T C c agrees with the experimental observations when one expresses the critical temperature of real ionic solutions in the RPM values [14]. When the strength of the  (20)). The line is a guide to the eye.
In order to compare our results with the available results of computer simulations [23] we rewrite the reduced temperature in equation (19) as k B T /ε, the usual definition for the SW potential. Figure 3 shows that the critical temperature scaled by the SW depth ε rapidly decreases with a decrease of the Coulomb interactions and then slowly approaches the SW limit k B T c /ε = 1.2667 [16]. Such behaviour is in agreement with the results of simulations [23].
The dependence of the reduced critical density ρ * c = ρ c σ 3 on the parameter α −1 is shown in figure 4. For very small α −1 , ρ * c increases slowly from the RPM value ρ * c ≃ 0.009 and then at α −1 ≃ 0.048 changes sharply to the value about 20 times larger than the RPM limit. For α −1 > 0.048, the critical density approaches the SW limit ρ * c = 0.2457 passing through a shallow minimum located at α −1 ≃ 0.1. Just as standard mean-field theories, the one-loop approximation considered in this paper underestimates the critical number density, the deviation from the simulation data increases when the strength of the Coulomb interactions becomes sufficiently large [21,31,33]. At the same time, our results demonstrate the behaviour of ρ * c for small α −1 which qualitatively differs from the available results of computer simulations. In simulations, the critical density rather rapidly increases with α −1 in the region α −1 ≤ 0.1 but it does not exhibit a jump. It is worth noting that the behaviour of ρ * c for α −1 > 0.1 qualitatively agrees with the results of simulations. It is instructive to view the range of the effective density-density interactions R N = a 2,2 /σ (see equation (15)) as a function of α −1 . It should be emphasized that the expression for a 2,2 given by (15) produces the correct result for the density correlation length ξ in the limit of charged point particles (see [16] and references herein). Figure 5 shows the behaviour of R N,c = R N (T c , ρ c ) with the variation of α −1 . The dimensionless amplitude of the density correlation length as a function of α −1 is shown in the inset. Remarkably, the trends of R N,c and ρ * c are very similar, especially for α −1 > 0.048. Furthermore, the behaviour of R N,c clearly shows the two distinctive regions separated by α −1 ≃ 0.048: the region where R N,c 0.5 and the region where R N,c > 0.6. Similarly, the two distinctive branches are clearly seen in the behaviour of ξ * 0 . This implies that α −1 ≃ 0.048 separates the models which demonstrate non-Coulombic phase behaviour ("solvophobic systems") from the models demonstrating Coulombic phase behaviour ("Coulombic systems"). to α = 20 intersects the moderate-density dashed curve. For α = 21.13, the spinodal intersects the low-density curve and it is tangent to the curve placed at higher densities.
The opposite situation appears for α = 20.77: the spinodal is tangent to the low-density dashed curve and intersects another dashed curve located at higher densities. The dashdotted curve which presents the loci of equation (22) crosses the dashed curves at their extremum points (exactly at the points of tangency of the solid and dashed curves). These special points, at which equations (21) and equation (22) hold are tricritical points [34].
The coordinates of the tricritical points are presented in table 1. The two tricritical points can be considered as limiting points which separate the low-density family of spinodals from the spinodals located at moderate densities. It should be noted that both a critical point and a tricritical point were predicted in [22] for a model system with additional short-range interactions added to the RPM (see Fig. 6B in [22]). On the other hand, no tricritical point was found in computer simulations [23] for the lattice version of the RPM-SW model with the lattice discretization parameter ζ = 10.   Coulomb interaction with respect to the short-range SW interaction, we calculate the gasliquid phase diagrams for α varying from 0 (the purely nonionic system) to ∞ (the purely Culombic system). It is worth noting that the parameter α is proportional to the inverse dielectric constant of the solvent ǫ when the ion charges are fixed.
Both the coexistence envelopes and the critical parameters have been calculated in a oneloop approximation which is equivalent to the random phase approximation. We have found that the very narrow region of α (20.77 ≤ α ≤ 21.13) separates the models demonstrating a nonionic type of phase behaviour from the models which demonstrate a "Coulombic" type of phase diagrams. This region is bounded from the both sides by the tricritical points: one point located at T SR trc,1 = 1.9233 and ρ * trc,1 = 0.0269 and another point with the coordinates T SR trc,2 = 2.1701 and ρ * trc,2 = 0.2905. The dependence of the reduced critical temperature on α (α −1 ) indicates a continuous variation from a phase transition driven by Coulomb interactions to a transition determined by "solvophobic" interactions. Such behaviour agrees with the available data of computer simulations performed for the lattice version of the RPM-SW model with high values of the lattice discretization parameter ζ (ζ = σ/l = 10) [23]. On the other hand, the reduced critical density ρ * c sharply changes with α in the region between the two tricritical points. This is at variance with the results of computer simulations mentioned above. Furthermore, no tricritical points were found for the lattice version of the model in the high ζ case in the computer simulations. It should be noted that the existence of the tricritical point for a similar model was predicted within a mean-field approximation in [22].
We have studied the range of the effective density-density interaction R N,c with the variation of α (α −1 ). The results clearly indicate the two distinctive regions in the phase diagram (R N,c 0.5 and R N,c 0.6) which are separated by α ≃ 21. The amplitude of the density correlation length ξ * 0 also shows different trends for α < 21 and α > 21. It should be noted that the values of ξ * 0 are larger for the Coulomb-dominated systems than for the non-Coulombic systems. The opposite situation takes place for R N,c .
Thus, we have obtained the results for a phase behaviour of the RPM-SW model that qualitatively differ from the findings of the available computer simulations. Further investigations are needed in order to confirm (or rule out) the presence of the tricritical points in the model with α ≃ 21. To this end, a two-loop approximation for the coefficients of the effective Hamiltonian (4) should be considered. On the other hand, it would be useful to perform computer simulations for the model parameters for which the theory predicts tricritical points.