Phase equilibria of size- and charge-asymmetric primitive models of ionic fluids: the method of collective variables

We study the phase diagrams of z:1 size-asymmetric primitive models (PMs) using a theory that exploits the method of collective variables. Based on the explicit expression for the relevant chemical potential (conjugate to the order parameter) which includes the correlation effects up to the third order we consider some particular cases of asymmetric PMs (1:1, 2:1 and 3:1 systems). Coexistence curves and critical parameters are calculated as functions of size ratio λ = σ + /σ −. Similar to simulations, our theory predicts the reduction of coexistence regions as well as a decrease of critical parameters T * c and ρ * c with an increase of size asymmetry.


Introduction
The study of phase diagrams of ionic systems in which the phase separation is mainly driven by electrostatic forces is of great fundamental interest.Electrolyte solutions, molten salts and ionic liquids are examples of systems with strong Coulomb interactions.The simplest representation of the systems dominated by Coulomb interactions is provided by the so-called primitive models (PMs).In these models, the ionic fluid is described as an electroneutral mixture of charged hard spheres immersed in a structureless dielectric continuum.A special case of the two-component PM consisting of monovalent equisize hard spheres is called a restricted primitive model (RPM).Because of its simplicity, the RPM has been extensively studied both theoretically and by simulations.Early studies established [1] that the system has a low-temperature gas-liquid like phase transition.A reasonable theoretical description of gas-liquid critical point of RPM was accomplished at a mean-field (MF) level using integral equation methods [2] and Debye-Hückel (DH) theory [3].Over the years numerous simulation studies have been devoted to the location of critical point but reliable estimations have been obtained only recently.Critical temperature is now known to an excellent accuracy while critical density only to a moderate accuracy [4,5].
From practical and fundamental points of view, asymmetric PMs are much more important than perfectly symmetric systems.The investigation of phase behavior of PMs with size and charge asymmetry was started more recently.The key findings from simulation studies of asymmetric models are as follows: the suitably normalized critical temperatures decrease with size and charge asymmetry while the critical densities increase with charge asymmetry but decrease with size asymmetry [6][7][8][9][10][11][12][13][14].Comparison of simulated critical parameters and theoretical predictions for asymmetric models has revealed that several established theories, such as the mean spherical approximation (MSA) and the original DH theory are not capable of predicting the trends observed in simulations [15,19].Moreover, both the original DH theory and the MSA predict no dependence on charge asymmetry in the equisize case.Only a direct inclusion of additional effects such as association or geometrical constraints in the case of size asymmetry allowed one to qualitatively reproduce the trends found in simulations [18,19].It should be noted that most simulations have concentrated on critical parameters and only a few studies [7,8,11] have addressed the phase coexistence in asymmetric PMs.
Recently [20], we have proposed a method for the study of gas-liquid phase diagram of PMs using the statistical field theory that exploits the method of collective variables (CVs) (see [21][22][23][24]).The approach is based on determining the chemical potential conjugate to the order parameter and enables one to take into account the effects of higher-order correlations.In [25], for the general case of an asymmetric PM, we derived an analytical expression for the relevant chemical potential which included the effects of correlations up to the third order.The results obtained for 1:1 and 2:1 systems qualitatively reproduced the dependence of critical temperature and critical density on the size asymmetry.In this paper which is an extension of [25], we mainly focus on the effects of size asymmetry on gas-liquid coexistence in asymmetric PMs.
The layout of the paper is as follows.In section 2 we outline the theory.We present the functional of grand partition function of a size and charge-asymmetric PM in terms of CVs and determine the order parameter.The method of calculating the relevant chemical potential is also briefly described.In section 3 we calculate the gas-liquid coexistence data and the critical parameters for some particular cases (1:1, 2:1 and 3:1 PMs).We conclude in section 4.

Model
We consider a classical two-component system consisting of N + cations carrying a charge q + = zq with diameter σ + and N − anions carrying a charge q − = −q with diameter σ − .The ions are immersed in a structureless dielectric medium.The system is electrically neutral: α=+,− q α ρ α = 0 and ρ α = N α /V is the number density of the αth species.
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 σ β .We call the two-component hard sphere system a reference system (RS).Thermodynamic and structural properties of RS are assumed to be known.φ C αβ (r) is the Coulomb potential: φ C αβ (r) = q α q β /(Dr), where D is dielectric constant, and hereafter we put D = 1.The model is characterized by the size and charge asymmetry parameters: The fluid is at equilibrium in the grand canonical ensemble.The grand partition function of the model (1) can be written as follows: where the following notations are used: is the inverse de Broglie thermal wavelength; (dΓ) is the element of configurational space of the particles.
It is worth noting that regularization of the potential φ C αβ (r) inside the hard core is arbitrary to some extent.For example, different regularizations for the Coulomb potential were considered in [26,27].Within the framework of Gaussian approximation of the grand partition function, the best estimation for critical temperature is achieved for optimized regularization [28] that leads to an optimized random phase approximation.However, this approximation does not work properly in higher orders of perturbation theory [27].Here we use the Weeks-Chandler-Andersen (WCA) regularization scheme for φ C αβ (r) [29].

Functional integral. The Gaussian approximation
Using the CV based theory, developed in [24] for a multicomponent continuous system with short-and long-range interactions in the grand canonical ensemble, we can write the functional of the grand partition function of the PM (3) with the interaction potential (1) in the form of a functional integral [25]: where Ξ MF [ν α ] is the grand partition function of the model in the MF approximation and να is the renormalized chemical potential Here is the CV which describes the value of the k-th fluctuation mode of the number density of the αth species, the indices c and s denote real and imaginary parts of ρ k,α and each of ρ c k,α (ρ s k,α ) takes all the real values from −∞ to +∞. ω k,α is conjugate to the CV ρ k,α .(dρ) and (dω) are volume elements of the CV phase space where the product over k is performed in the upper semi-space In the case of the WCA regularization we obtain for β φC αβ (k) where the following notations are introduced: The nth cumulant M α1...αn coincides with the Fourier transform of the nth connected partial correlation function of the RS [24].M α1...αn depends on both the wave vectors k i and the partial chemical potentials να .
Let us consider the Gaussian approximation of Ξ[ν α ] setting M α1...αn ≡ 0 for n 3.Then, after integration in (4) over ω k,α we obtain where Cαβ (k) is the Fourier transform of the two-particle direct correlation function in the Gaussian approximation The diagonalization of the square form in (7) leads us to the equation Eigenvalues ε1 (k) and ε2 (k) in the long-wavelength limit are found to be [25] ε1 where T * and δ are given by ( 5)-( 6) and is a reduced total number density.Equation ( 9) leads to GQQ (k = 0) = 0, where GQQ (k = 0) is the Fourier transformation of the charge-charge connected correlation function.It reflects the fact that the first moment Stillinger-Lovett rule is satisfied.
It is worth noting that equation ε2 (k = 0) = 0 (see (10)) leads us to the same expression for gas-liquid spinodal as that obtained in [30] but for another type of regularization of the Coulomb potential inside the hard core.
Eigenvectors ξ k,1 and ξ k,2 in the long-wavelength limit have the form [25]: where CVs ρ 0,N = δρ 0,+ +δρ 0,− and ρ 0,Q = zδρ 0,+ −δρ 0,− describe long-wavelength fluctuations of the total number density and charge density, respectively.As is seen, CV ξ 0,1 describes fluctuations of the charge density.In the general case z = 1, ξ 0,2 is a linear combination of CVs ρ 0,N and ρ 0,Q with the z-dependent coefficients.At z = 1, CV ξ 0,2 solely describes fluctuations of the total number density.Thus, we suggest that CV ξ 0,2 is connected with the order parameter of the gas-liquid critical point.Finally, after integration in (8) we arrive at the logarithm of the grand partition function in the Gaussian approximation where Φ C and M 2 are matrices of elements β φC αβ (k) and M αβ (ν α ; k), respectively.

Beyond the Gaussian approximation
We study the gas-liquid phase diagrams of asymmetric PMs using the method proposed in [20].First, we pass from the initial chemical potentials ν + and ν − to their linear combinations Chemical potentials ν 1 and ν 2 are conjugate to CVs ξ 0,1 and ξ 0,2 , respectively.Since we suggest that CV ξ 0,2 (see ( 12)) is connected with the order parameter, ν 2 appears to be of special interest in our study.It can easily be shown that ν 2 is an electrochemical potential.
We start with the grand partition function in the Gaussian approximation (13) replacing the cumulants M αβ (ν α ; k) by their values in the long-wavelength limit Taking into account that ln Ξ HS and M α1α2...αn are functions of full chemical potentials ν 1 and ν 2 , we present ν 1 and ν 2 as 1 and ν 0 2 being the MF values of ν 1 and ν 2 , respectively.Then, we self-consistently solve equations for the relevant chemical potential ∆ν 2 by means of successive approximations.The procedure of searching for a solution of equations ( 14)-( 15) is described in [20].As a result, in the first nontrivial approximation corresponding to ∆ν 1 = 0 we find the following expression for ∆ν 2 [25] where Apart from M α1α2 formulas ( 16)-( 17) include the third order cumulants M α1α2α3 or equivalently the third order connected correlation functions of the RS.In the case of a two-component hardsphere system, the analytical expressions for a second order cumulant can be obtained in the Percus-Yevick approximation using the Lebowitz' solution [31,32].The corresponding formulas for M α1α2 (k = 0) can be found in [25].In order to derive the expressions for the third order cumulants one can use the recurrent relation where ν 0 αi is the MF value of chemical potential ν αi that due to the electroneutrality condition coincides with the hard-sphere chemical potential of the α i th species.
Finally, the expression for the full chemical potential ν 2 can be written as follows where with ν HS + (ν HS − ) being the hard-sphere chemical potential of the αth species and ν S 2 being the combination of the self-energy parts of chemical potentials ν + and ν − ν The explicit expressions for ν HS 2 and ν S 2 can be obtained using the results of [32] supplemented by the electroneutrality condition and are presented in [25].
Below we use formulas ( 16)-( 21) for the study of the gas-liquid phase equilibria in asymmetric PMs.

Monovalent PMs with size asymmetry
First we consider a monovalent PM with size asymmetry corresponding to z = 1 and λ = 1.Because of symmetry with respect to the exchange of + and − ions, only λ < 1 (or λ > 1) need be considered in this case.Based on the expressions ( 16)-( 21) (at z = 1) supplemented by the Maxwell construction we calculate the coexistence data.Figure 1 demonstrates the calculated coexistence curves of the monovalent size-asymmetric PM at λ = 1.0, 0.75, 0.5, 0.25.As is seen, the region of coexistence narrows when the size asymmetry increases.This effect becomes increasingly pronounced as λ becomes smaller that is in qualitative agreement with the simulation findings [7,8].
Estimates of the critical temperature and the critical density are given by their values for which the maxima and minima of the van der Waals loops coalesce.The calculated critical parameters are reported in table 1.Similar to the simulation results, both the critical temperature T * c and the critical total number density ρ * c decrease when the size asymmetry increases.As was found from simulations [13], the critical volume fraction η c of a monovalent PM demonstrates an interesting nonmonotonic behaviour: when the size asymmetry increases, the volume fraction is seen to increase before it begins to rapidly decrease.For an arbitrary z, the volume fraction defined as  (11)) by the relation Our results for the dependence of η c on size asymmetry at z = 1 are shown in figure 2. The appearance of the first maximum (at λ = 1) is due to the fact that the expression for ∆ν 2 given by ( 16) tends to the random phase approximation when λ → 1.As was shown in [25], equation ( 16) reduces to the RPA in the limiting case of the RPM (z = 1 and λ = 1).For λ 2, the calculated trend of η c is similar to that found in simulations.Moreover, our results indicate that η c changes its trend at λ * ≈ 4 while the simulations yield λ * ≈ 4.26 (|δ * | ≈ 0.62).

Size-and charge-asymmetric PMs
Now we use formulas ( 16)- (21) for the study of the size-and charge-asymmetric PMs with z = 1 and λ = 1.In particular, we consider 2:1 and 3:1 models.As before, in order to calculate the coexistence curves and the corresponding critical parameters we apply the Maxwell construction.The phase coexistence envelopes for 2:1 systems with size asymmetries are presented in figure 3 for λ < 1 and in figure 4 for λ > 1.Similar to the monovalent case, the coexistence regions become     narrower with the increase of size asymmetry.The corresponding critical parameters decrease when size asymmetry increases that qualitatively agrees with simulation results [10].The calculated data for critical temperatures and critical densities are given in table 2. As far as we know, no simulation results are available regarding the effect of size asymmetry on the coexistence region for z:1 PMs.
Figures 5 and 6 demonstrate the effects of size asymmetry on the phase diagrams of 3:1 systems.The estimated values of the critical parameters are presented in table 3. Again, we can observe both the decrease of the critical parameters and the reduction of the coexistence regions with the increase of size asymmetry.As is seen from figure 5, our results demonstrate the phase transition at λ = 0.2.This contradicts the prediction made in [11] on the basis of an extrapolation of the simulation data.As was suggested in [11], the phase transition vanishes at λ 0.2 (or equivalently at δ −0.67) in the case of 3:1 PM.Thus, this issue deserves further investigation.

Summary
We have studied the effects of size and charge asymmetry on the gas-liquid phase diagrams of two-component PMs.For this purpose we have used the CVs based theory that enables us to derive an exact expression for the functional of grand partition function of the model and on this basis to develop the perturbation theory.As was shown in [20], the well-known approximations for the free energy, in particular Debye-Hückel limiting law and MSA, can be reproduced within the framework of this theory.
In this paper, extending our previous studies [25], we have calculated the coexistence curves for the 1:1, 2:1 and 3:1 PMs at different values of the size asymmetry parameter.In all cases considered our results have demonstrated the coexistence region narrowing with the increase of size asymmetry.Since such a behaviour qualitatively agrees with the simulation results for the monovalent PMs we expect that our predictions for the 2:1 and 3:1 PMs are also qualitatively correct.Moreover, we have calculated the critical temperatures and critical densities for the 3:1 PM at λ = 0.2, 1, 1.5, 2, 3.As before, our theory in this case predicts a decrease of the both critical parameters with the increase of size asymmetry.
Summarizing, we can state that the theory developed for the determination of the relevant chemical potential has for the first time enabled us to get a qualitatively correct description of the effect of size asymmetry on the phase diagrams of asymmetric PMs without directly including the ionic association.

Figure 1 .
Figure 1.Coexistence curves of the monovalent PM at different values of λ.

Figure 2 .
Figure 2. The critical volume fraction of the monovalent PM as a function of size asymmetry.

Table 2 .
Critical parameters T * c and ρ * c of the 2 : 1 PM for different values of λ.

Table 3 .
Critical parameters T * c and ρ * c of the 3:1 PM for different values of λ.