The liquid-vapor interface of the restricted primitive model of ionic fluids from a density functional approach

We investigate the liquid-vapor interface of the restricted primitive model for an ionic fluid using a density functional approach. The applied theory includes the electrostatic contribution to the free energy functional arising from the bulk energy equation of state and the mean spherical approximation for a restricted primitive model, as well as the associative contribution, due to the formation of pairs of ions. We compare the density profiles and the values of the surface tension with previous theoretical approaches.


Introduction
According to the simplest model of ionic fluids, called the restricted primitive model (RPM), the fluid consists of charged hard spheres (at total density ρ) of equal diameter, d, half of which carry a charge of +q, while the other half carry a charge of −q. These spheres are immersed in a dielectric medium of dielectric constant ǫ. Reduced thermodynamic quantities, temperature and density, appropriate to this model are defined as follows T * = kT ǫd/q 2 and ρ * = ρd 3 .
Perhaps, the first successful statistical-mechanical description of this model was proposed by Debye and Hückel [1] about ninety years ago. The Debye-Hückel approach predicted the existence of the first-order transition for RPM and was subsequently applied in several works, e.g., to describe phase behavior of three-component ionic fluids [2]. More recently, the bulk RPM was studied by means of the Ornstein-Zernike (OZ) equation supplemented by the mean spherical approximation (MSA) [3][4][5], which permits to obtain an analytical solution for the pair correlation functions. The MSA correlation functions can be used to evaluate an equation of state [6,7]. In spite of quite reasonable results obtained for the structure and thermodynamics in a general sense, neither virial nor compressibility equations of state permit to obtain gas-liquid transition for the RPM, which contradicts the computer simulations [8][9][10][11][12][13][14][15]. By contrast, the equation of state resulting from the energy route predicts the first-order phase transition between an ionic vapor and a dense ionic liquid [6,7]. However, there is a big discrepancy between the values of critical temperature and critical density of MSA theory and the critical parameters resulting from computer simulations. Moreover, a more sophisticated generalized mean-spherical approximation [16], which was designed to reconcile thermodynamic inconsistency of the MSA, does not lead to a better agreement of the critical parameters with simulations as well.
Since both MSA and generalized MSA theories underestimate the critical density and overestimate the critical temperature, attempts were made to include ad hoc the effect of ion pairing (association) into the RPM, which is most pronounced along the vapor branch of the coexistence envelope and near the critical point [13,17]. The inclusion of the association into the theory improves the agreement of theoretical predictions for the critical constants with simulations (cf.  [13]). However, the overall shape of the coexistence envelope coming from several theoretical approaches is not well described compared to simulation data. In particular, the liquid densities along coexistence are poorly described for different temperatures. We should bear in mind that the Hamiltonian of the model with association used in the theoretical developments differs from the genuine RPM used in simulations. Recent developments in the modelling of different types of associations and clustering in liquids and solutions have been comprehensively and critically reviewed by Holovko in [18].
The liquid-vapor interface is one of the simplest examples of nonuniform systems. However, compared to the case of Lennard-Jones fluids [19], much less is known about the liquid-vapor interface in the RPM. Theoretical work on the gas-liquid interface for RPM was pioneered by Telo da Gama et al. [20], that used the gradient expansion method [21]. The problem of describing an ionic liquid-ionic vapor interface has been also studied by Groh et al. [22], who proposed a density functional approach to evaluate the surface tension and the density profiles through the interface. Their approach involved a local density approximation for the hard-sphere part of the free-energy functional and a nonlocal treatment of Coulombic contributions. The latter term was evaluated approximating the inhomogeneous pair correlation functions by their bulk counterparts. There are significant differences between the theory used by Groh et al. [22] and the density functional approaches that were proved successful in the description of electrical double layers at charged walls [23][24][25][26][27][28][29][30][31]. Namely, the theories of [23][24][25][26][27] treat the electrostatic part of the free energy by means of a second-order density expansion about the density of a reference fluid, which was taken to be the homogeneous bulk fluid far from the surface. Although such an approach can be appropriate for many purposes, but it is problematical when it comes to a liquid-vapor interface, where two bulk fluids are involved. Moreover, it is known that in the case of Lennard-Jones fluids, the corresponding second-order expansion with respect to a homogeneous fluid fails to account for the liquid-vapor coexistence. More recently, the description of a liquid-vapor interface was the subject of investigations of Weiss and Schröer [32]. Using a square-gradient type theory they computed the density profiles and the interfacial tension at different temperatures using Debye-Hückel theory and its extension to ion-pair formation, as well as adding the free energy term describing correlations between ion pairs as entities and free ions, as developed by Fisher and Levin [17].
Any density functional theory determines thermodynamic properties of an inhomogeneous fluid from the Helmholtz free energy F and its functional dependence on the local densities {ρ i (r)}. The free energy functional is commonly decomposed into the sum of three contributions, namely into the ideal, the hard-sphere, and the electrostatic terms. Of course, one can also include here an additional term, due to possible association of ions. Various formulations of the DFT have been discussed in the literature. As we have already stressed, a majority of the DFTs have followed perturbational second-order expansion of the electrostatic free energy functional with respect to a bulk homogeneous fluid [23][24][25][26][27].
Recently, Gillespie et al. [33,34] proposed a version of the electrostatic free energy functional that replaces a uniform reference system with a suitably chosen position-dependent reference fluid. The inhomogeneous reference fluid densities are then computed from the local densities by a selfconsistent iteration procedure. Actually, Gillespie et al. [33,34] proposed a reference fluid density functional which permits to construct a reference model that locally satisfies electroneutrality condition and has the same ionic strength at every point as the inhomogeneous fluid in question. Such a construction of the reference fluid permits to apply the expression for the electrostatic contribution to the free energy, which results from the equation of state for a bulk ionic system and makes it unnecessary to employ the direct correlation functions as input quantities. This kind of approach was proposed by us to study liquid-vapor transitions in RPMs confined to slit-like pores [31,35]. The results of our approach reasonably well reproduce the structure of ionic fluids at a charged wall at sufficiently high temperatures and correctly predict [36] the temperature dependence of the double layer capacitance, in agreement with experiments and simulation data [37]. Quite recently, the theory was also extended to include the effects of association between unlike ions [38]. However, the applications of the approach described above would not be complete without explor-

13603-2
ing liquid-vapor interface of ionic fluid. Therefore, in this communication we intend to apply the theory of [35,36,38] to the study of liquid-vapor interface of the RPM. We calculate the density profiles and the interfacial tension. Moreover, we would like to investigate the effect of association leading to the formation of ion pairs and whether this type of effects improves the description of the interfacial properties similarly to the bulk structure and thermodynamics.

Theory
We consider a binary mixture of ionic species. The symbols d, Z i = q i /e and µ i denote, respectively, the hard-sphere diameter, valence of ions and the chemical potential of species i = 1, 2, and e is the electron charge. The interaction between the ions is where r is the distance between a pair of ions. We also assume that Z 1 = −Z 2 = 1 and that the dielectric constant ε is uniform throughout the entire system. The DFT we use in this work is identical to that described in [35,38] and therefore we present only its brief description. If there is no external potential field the grand potential of the system is written in the form, According to the usual density functional treatment, the free energy density functional, , whereas for the hard-sphere contribution we apply an expression resulting from a recent version [39] of the Fundamental Measure Theory [40,41], with the free energy density consisting of terms dependent on scalar and on vector weighted densities, n α (r) (α = 0, 1, 2, 3) and n α (r) where ξ(r) = |n V 2 (r)|/n 2 (r). The definitions of weighted densities, n α (r), α = 0, 1, 2, 3, V 1, V 2, are given in [40,41]. The electrostatic contribution is evaluated from the approach described in [35,38], according to which . The reference fluid densities {ρ i } are evaluated by employing the approach of Gillespie et al. [33,34]. In the case of gas-liquid interface, the spatial symmetry of the RPM implies that the density profiles of the two ionic species should be the same. Therefore, in the presently considered situation, the approach of Gillespie et al. [33,34] is simplified and the reference fluid densities are given byρ where the weight function W i (r) is just a normalized step function

13603-3
The radius of the sphere over which the averaging is performed, R f , is approximated by the "capacitance" radius, i.e., by the ion radius plus the screening length The evaluation of R f requires an iteration procedure. This iteration loop should be carried out in addition to the main iteration procedure used to evaluate the density profiles, see for details [33,34]. Finally, the associative contribution, f as , is formulated at the level of the first-order thermodynamic perturbation theory [38,42] where α is the degree of dissociation according to the mass action law, 2α whereas K γ follows from the so-called simple interpolation scheme [13,38,44] The theory presented above uses different weighted densities to evaluate the hard sphere and electrostatic free energy contributions. Note that the idea of using different weighted densities to evaluate these contributions has also been employed by Patra et al. [25,26].
The density profile is obtained by minimizing the excess grand potential functional ∆Ω = Ω−Ω b , where Ω b is the grand potential of a bulk uniform system of density {ρ b,i }.
The evaluation of the density profiles requires the knowledge of the densities of both coexisting phases. Therefore, we have precisely evaluated the bulk phase diagram prior to the density profile calculations. Next, the local density calculations were carried out assuming that for z −L z and for L z the density of the fluid is equal to the bulk densities of gaseous and liquid phases, respectively. The value of L z was evaluated at each temperature by testing the convergence of the density profiles towards the final solutions. We have started with L z = 30d and the consecutive runs were carried out doubling the value of L z for each new run. All the integrations were performed using Simpson method with the grid size of 0.01d. The convergence criterion definitly states that the maximum difference between two consecutive iterations must be smaller than 1E-8 percent. The knowledge of the density profiles of ions permits to calculate the excess grand thermodynamics potential per unit surface area, A, i.e., the surface tension γ = ∆Ω/A. Figure 1 shows the phase diagram in the density-temperature plane. The phase diagram resulting from the MSA energy equation of state has been presented in several previous works, cf. [13,22]. Nevertheless, we have decided to include it here for completeness of the study. The solid line corresponds to the system without association (i.e. the associative free energy term is neglected in the grand thermodynamics potential functional), whereas the dashed line describes phase behavior of the model with the association effects included. Inclusion of chemical association into the theory 13603-4 provides a mechanism for the formation of ionic pairs. The fraction of pairs depends on density, on temperature and on the association constant. From physical point of view, the formation of pairs of oppositely charged ions, alters the effective interactions between all the particles. In effect, the critical temperature of liquid-vapor transition, T * c ≈ 0.0745, becomes lower compared to the critical temperature of the model without association, T * c ≈ 0.0786. The liquid-vapor phase diagrams are strongly asymmetric. At low temperatures the vapor phase is very dilute. For example, for the system without association at T * = 0.05 the vapor density is ρ * b ≈ 1.094E − 5, whereas the density of the coexisting liquid phase is ρ * b ≈ 0.2135.   First, we consider the system without association. The evaluation of the liquid-vapor density profiles was carried out assuming that the limiting values ρ i (−L z ) and ρ i (L z ) are equal to the MSA densities of coexisting phases. Of course, the agreement of our results with the results published  previously [22,32] essentially depends on the agreement between the bulk phase diagrams. One should keep in mind the above remark while analyzing figure 2, where we show a comparison of our results with those of Groh et al. [22]. The latter profiles were obtained from a different bulk theory (called by them "the MSA1 approach" [22]), which yields the value of the critical temperature, T * c ≈ 0.0846, different from that obtained by MSA theory. Note that a comparison of the entire phase diagrams resulting from the MSA and MSA1 approaches in given in figure 1 of [22].   The way in which the density profile approaches the asymptotic bulk values has been discussed by Groh et al. [22]. We have checked that the decay of the function ρ l b,i − ρ i (z) (where ρ l b,i is the liquid-phase bulk density) with z → L z is exponential, i.e. ln[ρ i (z) − ρ l b,i ] ∝ z/ξ. This point is illustrated in figure 3, where we have plotted ln[ρ l b,i − ρ i (z)] versus z/d. At larger distances this dependence is fairly well approximated by a straight line with the correlation coefficient, R 2 , very close to unity.

Results and discussion
We now consider the effect of the association, cf. figure 4. Part a shows a comparison of the total density profiles for the systems with (symbols) and without association (lines) at different temperatures. At low temperatures the liquid density for the system with the association included is higher than for the system with the association effects switched off, cf. figure 1. The profiles for the model with association are more "smeared out", but this is the effect of lowering the critical temperature due to association effects. Similarly to the case of a system without association we have found exponential decay of the function ρ l b,i − ρ i (z) with z → L z (the corresponding plot was omitted for the sake of brevity). Figure 4b shows the profiles of α(z) (the phase diagram in the α-T plane has been already presented in [38]) We stress that the interfacial region of α(z) is shifted towards rarefied phase compared to the density profiles. It means that the position of the "pseudo-Gibbs" dividing surface, which can be introduced for the dissociation degree in analogy to the usual liquid-vapor dividing surface is also shifted compared to the position of the usual dividing surface (which is located at z = 0). Note that similar shift of the dissociation degree was also observed in the case of liquid-vapor interfaces of non-ionic fluids [45].
The surface tension obtained from the density functional theory and from the MSA1 and MSA approaches, developed by Groh et al. [22], as well as from the square gradient theory [20] are shown in figure 5. We have plotted here the ratio γd 2 /kT c versus the temperature reduced by the critical temperature, T /T c . Except for the temperatures very close to the critical temperature, the present

13603-6
The liquid-vapor interface of the restricted primitive model  [22], the curve labelled "square gradient theory" was evaluated by Telo da Gama et al. [20], the curve abbreviated as DFT was obtained for the model without association, whereas black circles -for the model with the association effects included. theory predicts lower values of the surface tension than all the rest approaches. The associative free energy terms still lower the value of γd 2 /kT c . Note that the MSA theory grossly underestimates the densities of the coexisting liquid phase. This can suggest that the values of the surface tension are also underestimated. On the other hand, the value of the critical temperature is overestimated by the MSA approach, and hence can lead to too high values of the interfacial tension.
Near critical temperature, both models (i.e., the model with and without association) yield the standard mean-field critical behavior, i.e., γd 2 /kT c ∝ (1 − T /T c ) (3/2) (see figure 6). As the temperature approaches the critical point, the theory yields a nearly linear variation of ln[γd 2 /kT c ] with ln[(1 − T /T c )] and the slope of the straight line approximating the obtained data is almost perfectly equal to (3/2).  The width of the interfacial zone can be characterized by the parameter W , defined as [46]

13603-7
where z 0 is given by ρ(z 0 ) = (1/2)[ρ s (z = L z )+ρ(z = −L z )] and ρ(z = L z ) and ρ(z = −L z ) are the total densities of the coexisting liquid and vapor phases. Obviously, the interface becomes wider as the temperature increases. The effects of association slightly increases the interface width. As the temperature approaches the critical temperature, the interfacial width diverges. the exponent (0.5) is characteristic of mean-field type theories [22]. The scaling of W is consistent with the scaling of the surface tension. Our discussion involved the results of different, but solely theoretical approaches. However, it must be clarified how accurate our approaches are (with and without association) compared to the available computer simulation data. In particular, Alejandre and co-workers investigated the liquidvapor interface for the RPM using molecular dynamics simulation method supplemented by certain specific technical peculiarities dealing specifically with the RPM [47,48]. These authors presented the density profiles, ρ * (z) at four reduced temperatures, T * =0.035, 0.038, 0.040 and 0.043 [47]. On the other hand, a single density profile at T * = 0.391 relevant to our study was given in [48] (this profile is for the soft primitive model but with parameters permitting comparison with the RPM). The authors have not performed estimates of the critical parameters coming from their method of study of the RPM. Thus, they somewhat arbitrarily used the critical temperature obtained in [11], T * c = 0.0492, for their purposes. We have rescaled the temperature and density with respect to the critical parameters of both theoretical approaches (with and without association), respectively, and present here a comparison of the density profiles with computer simulation data in figure 8. Note that the critical density given in [11] was used to rescale the simulation data. We observe that the rescaled liquid phase density at coexistence coming from the DFT for the model with ion pairing is in reasonable agreement with the simulation data. The DFT for the model without association greatly overestimates these liquid densities. However, the width of the interfacial region is not described well. Theoretical approaches yield narrower interface width for two reduced temperatures in question. It is difficult to attribute this specific inaccuracy to the particular term of the free energy functional. However, it seems that the mean field consideration requires improvement in order to better describe interparticle correlations and hopefully to obtain a more adequate interface width.
As concerns the relationship between the simulation data for the surface tension and the results of theoretical approaches we would like to comment on the following points. A set of computer simulation results for the surface tension were put together with the predictions of some theoretical approaches [20, 22,32] (that were discussed in the introductory part of our communication) on arbitrary temperature scale, see figure 2 of [47]. Apparently, the theories yield the surface tension of the same magnitude as computer simulation data. This permitted to conclude that the surface tension from this set of theories [20,22] is overestimated because the effect of ion association is not  taken into account and solely the theory that accounts for association [32] is in good agreement with simulations. Then, it would seem that our approach taking ion pairing into account is the best, see figure 5. Moreover, a comparison of the theories, simulations and experimental data was performed [48] using an arbitrary scale. Such an idyllic picture is, however, entirely destroyed, if one rescales the reduced temperature scale by the critical temperature of each theory and of simulations. This is actually an adequate comparison. According to [47] the simulations yield, e.g. γd 2 /kT c =0.3 at T /T c =0.6. Such a value is very much higher than the results of any theory, c.f. figure 5, at this particular temperature. Moreover, this means that the inclusion of association effects makes things worse than from any other theory. This is difficult to accept, however. According to the more recent simulation work of the same authors [48] it is difficult to draw a definite conclusion about the dependence of the surface tension on temperature for the model in question due to the statistical error associated with the value for surface tension. Moreover, it has been discussed [48] that the inflection point at T * = 0.04 (again it is a reduced temperature of simulation) can be present on the dependence of surface tension on temperature as a result of chemical association of ions in the RPM in close similarity to hydrogen bonded fluids. In our opinion, additional work, both theoretical and simulational, is necessary to obtain definite answers about the behavior of surface tension on temperature for the model in question and the related models. From theoretical viewpoint one needs a theory that yields a more adequate shape of the bulk coexistence envelope. Then, the free energy expression coming from this approach, if available, would contribute to the development of better density functional approaches for inhomogeneous ionic fluids.
To summarize briefly, in this work we have applied density functional theory to the study of the liquid-vapor interface of a RPM fluid. Of course, the theory we use here is ad hoc. However, our previous calculations [35,36] have shown that it reproduces reasonably well the structure of a fluid at a charged and uncharged wall, predicts the dependence of the capacitance of the double layer on the temperature and yields "capillary evaporation" phase diagrams for confined ionic systems [35,36,38]. We have shown that the inclusion of ion pairing effect leads to a better agreement of the density profiles of ions across the liquid-vapor interface. A more sophisticated approaches [13], compared to the present one, including the effects of ionic association can lead to even better agreement of the coexistence envelope with simulations and of the density profiles across the liquid-vapor interface. We plan to extend our study in this respect in the nearest future. However, the inaccuracies of the values for the liquid phase density at coexistence at different temeperatures prevent us from obtaining a reasonable agreement with the available simulation data for the surface tension. It seems, however, that a more ample set of simulation data using different techniques is necessary. Nevertheless, the theories of this study yield a decreasing surface tension with increasing temperature as one intuitively would expect. Still we are not able to conclude that the theoretical approach is entirely successful both for the fluid-solid and fluid-fluid interfaces and look forward to improving it in future work.