Towards the description of the phase behavior of electrolyte solutions in slit-like pores. Density functional approach for the restricted primitive model

We develop a density functional approach for the phase behavior of the restricted primitive model for electrolyte solutions confined to slit-like pores. The theory permits to evaluate the effects of confinement on the ionic vapor – ionic liquid coexistence envelope. We have shown that due to confinement in pores with uncharged walls the critical temperature of the model decreases compared to the bulk. Also the coexistence envelope of the transition is narrower in comparison to the bulk model. The transition between dense and dilute phase represents capillary evaporation. We have analyzed changes of the density profiles of ions during transition. Possible extensions of this study are discussed.


Introduction
The so-called restricted primitive model of electrolytes (RPM) has been widely used to describe thermodynamic and structural properties of electrolyte solutions and molten salts.It consists of charged hard spheres with charges eZ i , i = 1, 2 (or +, −), e is the electron charge, all of the same diameter d, immersed into a structureless dielectric medium which is characterized only by the dielectric constant ε.Moreover, the neglect of the solvent structure explicitly results in the absence of ion-solvent effects.The model is a system satisfying the electroneutrality in general.In the case of symmetric RPM |Z 1 | = |Z 2 | = Z.In contrast to Debye-Hückel level of description of electrolyte solutions, the RPM is designed to account for the effects of the excluded volume and of electrostatic interactions between ions.The state of the RPM is fully described with just two parameters, i.e., reduced temperature -T * = kT εd/e 2 Z 2 and reduced density of ions -ρ * b = ρ * b1 + ρ * b2 , ρ * bi = ρ bi d 3 (the subscript b is used here to distinguish the bulk model).
The bulk, i.e. uniform, RPM has been successfully studied by means of Ornstein-Zernike equation supplemented by various closures.One of the most popular ones is the mean spherical approximation (MSA) [1,2].This approximation permits to obtain an analytical solution to the pair correlation functions by using Baxter factorization method.Correlation functions can thus be used to provide the equation of state.The MSA provides reasonable results for the pair distribution functions of ions, for activity coefficients and other thermodynamic properties in general sense.However, neither virial nor compressibility equations permit to obtain gas-liquid type transition in the RPM.This situation is in contradiction with computer simulation data.
It is well known that the problem lies in thermodynamic inconsistency of the mean spherical approximation.Only, the equation of state resulting from the energy route, predicts the first-order phase transition between vapor of ionic species and dense ionic liquid, with the critical temperature equal approximately to T * cb ≈ 0.0786 and the critical density equal to ρ * cb ≈ 0.0144.Unfortunately, there is a big discrepancy between the above values and the critical parameters from computer simulations, see e.g.table 1 of reference [3].A more sophisticated generalized mean-spherical approximation (GMSA) [4] has been designed to reconcile the thermodynamic inconsistency of the MSA.This approach is, however, more computationally demanding to be implemented in the study of structural and thermodynamic properties including possible phase equilibria.The critical parameters emerging from the GMSA treatment are similar to the ones obtained in the MSA.Better values for critical parameters with respect to simulation data have been obtained within the framework of recent theories which take into account the association of ionic species [3,5].
Our interest, however, is to proceed with the discussion of the phase behavior of inhomogeneous restricted primitive model.It is worth mentioning that the MSA solution to the bulk primitive model has provided essential progress in understanding the problem of electrolyte in contact with solid surfaces.Historically, the description of adsorption of ions on planar charged walls and in pores of different simple geometry has been attempted first, by using the so-called singlet integral equations for the density profiles of ions, see e.g.[6] for a comprehensive review.Later, the second-order integral equations, which allow for simultaneous evaluation of the local densities and of nonuniform pair correlation functions (but are much more demanding computationally compared with the singlet integral equations) have been solved [7].Most recently, similar to the developments for inhomogeneous simple fluids and their mixtures, density functional theories have been shown successful and powerful methods of describing inhomogeneous RPM [8][9][10][11][12][13][14][15][16].
Any density functional theory (DFT) determines thermodynamic properties of an inhomogeneous fluid from the Helmholtz free energy F and its functional dependence on the local densities of particles {ρ i (r)}.The free energy functional commonly is decomposed into the sum of three contributions, namely the ideal, F id , the hardsphere, F hs and the electrostatic, F el , terms.Various formulations of the DFT have been discussed in the literature.They differ in the formulation of expressions for the F hs and F el terms.
To our best knowledge, except for the approaches [17] and [18], all the DFTs of nonuniform RPM have followed the formulation of the electrostatic free energy functional introduced by Mier y Teran et al. [8], by Rosenfeld [9] and by Kierlik and Rosinberg [10].This formulation is based on the perturbation with respect to a bulk homogeneous reference fluid.However, not all the systems are amenable to such a treatment, since the effects invoked by nonuniformity are not necessarily small.Recently Gillespie et al. [19,20] have proposed a version of the electrostatic free energy functional that replaces a spatially uniform reference fluid with a suitably chosen position-dependent reference fluid.The inhomogeneous reference fluid densities are computed from the local densities by a self-consistent iteration procedure.
In spite of a great number of publications, in which DFT's theories have been applied in the studies of several aspects of adsorption of ions, the problem of phase transition in nonuniform RPM has been considered only by Groh et al. [17] and by Weiss and Schroer [18], who have investigated the gas-liquid interface for the RPM.This situation is in contrast with the case of nonuniform simple fluids, where the approaches of DFT type have been principally used in describing the surface-induced phase transitions.
Thus the aim of this work is thus to attempt to solve such a problem and to propose a density functional approach, which can be used in determining the effect of confinement on the gas-liquid phase diagram of the RPM.The proposed theory differs from the theories developed so far with respect to the choice of the F el term.
In the present approach we apply the electrostatic part of the free energy for the RPM resulting from the energy route within the MSA framework.The electrostatic functional is written down analogously to the hard-sphere contribution (or e.g. the associative free energy contribution in the case of uncharged but associating fluids, see e.g.[21]).In other words we use an expression for the bulk electrostatic free energy and replace the uniform densities by specifically designed averaged densities.

Theory
Just to recall, we restrict ourselves to binary mixture of ionic species.The symbols d i = d = 1, Z i and µ i denote, respectively, the hard-sphere diameter, valency of ions and the chemical potential of species (i = 1, 2).Moreover, in this work we consider 1: The fluid is confined in a slit-like pore of the width H.The walls of the pore are assumed to be infinitely thick.The distance normal to the walls is denoted by z.
The interaction between the ions is, where r is the distance between ions.We also assume that the dielectric constant ε is uniform throughout the entire system and is the same at any point.The interaction between the ions of species i and pore walls is given by where v i (z) and w i (z) are the nonelectrostatic and the electrostatic contributions to the external potential, respectively.The nonelectrostatic potential provides a confinement and is chosen in the simple form, Also we assume that the inner surface of each of the pore walls is charged with the same charge density, σ.Due to the ions of equal diameter and the equivalency of pore walls, the electrostatic interaction between each of the ions and a wall possesses symmetry with respect to the pore center at z = 0.The electrostatic interaction between an ion and the single surface is given by where σ is the charge density on a pore wall and z is the distance from the wall.The electrostatic potential Ψ(z) is determined by the Poisson-Boltzmann (PB) equation where ρ i (z) are the density profiles of ionic species.The integration of the PB equation is carried out taking into account the symmetry of density profiles with respect to the pore center such that dΨ(z) dz z=0 = 0.
Also, we assume the second boundary condition, the value of the electrostatic potential on each of the pore walls is fixed.After performing integration, one arrives at the equation A coupling between the density profiles and the charge density at each of the two pore walls, σ, is given by In DFT the grand potential of an inhomogeneous fluid is written in the form [22], where the free energy functional, F [{ρ i }], is decomposed into ideal (id), hard sphere (hs) and electrostatic (el) terms as follows: The ideal term is known exactly, For the hard-sphere term, however, we apply the expression resulting from a recent version [23] of the Fundamental Measure Theory [11], with the free energy density consisting of the terms dependent on scalar and vector weighted densities, n α (r) (α = 0, 1, 2, 3) and n α (r and where ξ(r) = |n V 2 (r)|/n 2 (r).In equations ( 10), ( 12) and ( 13) the spatial dependence of all variables has been suppressed, for the sake of simplicity.We have neglected the tensor contribution to the free energy [24].However, inaccuracy introduced by this approximation is negligible for the problem in question.The weighted densities, n α (r), α = 0, 1, 2, 3, V 1, V 2, are defined as sums of the spatial convolutions of the average densities and the corresponding weight functions where w αi (r) are given in [11].
Commonly, in the previously used DFT approaches, the electrostatic contribution to the free energy, F el [{ρ i }], has had the form of integrals over the short-ranged part of the direct correlation functions of ionic species.However, in this work we propose and test the following representation of this functional: where ρi (z) denote suitably defined inhomogeneous average densities corresponding to the reference fluid.One of the simplest possible choices of the form for f el [{ρ i (z)}] is to employ for this functional the expression resulting from the MSA equation of state evaluated via the energy route [1][2][3].Namely, where with κ denoting the inverse Debye screening length, The last three expressions as written above correspond to an electroneutral bulk fluid (the diameter of ions is our length unit, d = 1, as has been mentioned at the beginning of this section).A confined fluid is in equilibrium with the bulk fluid which is considered to be at density ρ b = ρ b1 + ρ b2 .The grand thermodynamic potential of the bulk RPM, Ω b , is then calculated according to equations ( 9)- (18), in which the density profiles are substituted everywhere by constant bulk densities, {ρ bi }.In calculating the properties of the bulk system, the electrostatic potential, Ψ, is set equal to zero, obviously.
In order to determine the electrostatic contribution to the free energy of inhomogeneous fluid, we need now to define a reference fluid and the corresponding averaged densities {ρ i (z)}.This reference fluid should be defined in such a way that the averaged densities {ρ i (z)} should not violate electroneutrality condition at any point r, because otherwise equation (16) could not be applied.This requirement does not mean that the nonuniform system in question, i.e. the system with the spatially varying densities {ρ i (z)} is locally electroneutral.
To proceed, we follow the development proposed by Gillespie et al. [19,20].Let us define weighted densities where W i (r) is a weight function.Gillespie et al. [19,20] provided the definition of the weight functions W i (r).The weight function is just a normalized step function The radius of the sphere over which averaging is performed is not a well-defined characteristic.In fact, it is an approximation for the electrostatic length scale or the "capacitance" radius (i.e. the ion radius plus the screening length) In addition, Gillespie et al. [19,20] have required that the reference fluid should have the same ionic strength as the system with weighted densities ρi (z).Finally, the averaged densities {ρ i (z)} resulting from the approach [19,20], applied to 1:1 electrolyte, are given as, Most general definition of the averaged densities {ρ i (z)} for multicomponent systems of ions with arbitrary valency can be found in [19,20].
Since equations ( 20), ( 21) and ( 22) are coupled, the evaluation of R f requires an iteration procedure.This iteration loop would be necessary, in addition to the principal iteration procedure, in order to obtain density profiles of species emerging from the minimization of the excess grand thermodynamic potential according to our discussion which is due slightly further.However, in this work we have used further simplification and have assumed that the second term of equation ( 21) is computed by using bulk value for Γ, i.e. by using {ρ bi } in spite of {ρ i (r)} in this equation.This approximation is along the lines of the approximation invoked previously by Groot [25].We are aware that a detailed structure of the confined fluid depends on the parameter R f , cf. the results of [19,20].However, the principal aim of this work is to study the phase behavior and to get a qualitative insight into the changes of the gas-liquid coexistence due to confinement.For this reason we have decided to apply the theory as simple as possible.Note also that the idea of using different weighted densities to calculate the hard-sphere and the electrostatic contributions to the free energy functional has been used by Patra et al. [13,14].The latter approach has followed Denton-Ashcroft formulation [26] of density functional theory for simple fluids.
The density profile equation is obtained by minimizing the excess grand potential functional It reads, where the one-particle direct correlation functions for the system in question, c (1)hs i (z; [{ρ i }]), and their bulk counterparts, c (1)hs i ([{ρ bi }]), are given explicitly in [11,19].In order to complete the theory, we write down the relation between the chemical potential, µ (µ = µ 1 = µ 2 ) and the bulk fluid density, where Γ = κ /(2 The knowledge of the density profiles of ions permits to calculate the excess adsorption isotherms of ionic species in the pore via common definition, and the total average density of ions in the pore, The normalization of the average density is performed over the region of the pore in which the density profiles are nonzero.
The principal objective of the present study lies in constructing and analysing the phase diagrams of the RPM under confinement.The phase diagrams have been obtained ordinarily.We analyze the dependence of the grand thermodynamic potential, ∆Ω, on the chemical potential of species, at several temperatures and look for crossing points between branches corresponding to different phases.For the sake of convenience the electric potential on the pore walls, V , has been chosen as a variable that specifies the system, as common in DFT calculations.Consequently, the charge density on the pore walls, σ, results from the calculations of the profiles.

Results and discussion
The formalism has been presented above in detail.It would be consistent to test the method by comparing theoretical results with computer simulation data.It is actually possible to evaluate the performance of the theory at temperatures higher than the critical temperature of the bulk RPM.Unfortunately, such a difficult problem as simulation of the phase equlibria for RPM in slit-like pores has not been solved so far.Therefore we are not able to perform a comprehensive test of the theory at present.Nevertheless, we are convinced that the results below are at least qualitatively correct.Moreover, our intention is to starightforwardly apply the method to the problem of phase equilibria in the RPM under confinement.It can be applied to the slit-like pores with either uncharged or charged walls.However, at this stage of our research, we restrict the presentation of the results only to the case of 1-1 electrolyte in the framework of restricted primitive model confined to slit-like pores with uncharged walls, i.e.V = 0.A more extended investigation of the effects of the value for electrostatic potential on the phase equilibria will be presented in future work which at present is in preparation.Here, we would like to determine how the phase behavior of the RPM is effected by purely geometric confinement.For illustrative purposes, in figure 1 we show an example of the dependence of Ω * ex ≡ Ω ex /kT on µ * i = µ i /kT (µ * = µ * 1 = µ * 2 ) at three temperatures, T * , in a rather narrow pore (H = 11).At each temperature there are two branches of the excess grand potential that cross at the equilibrium transition point, indicated by an arrow.The corresponding plots of the excess adsorption isotherms are given in figure 1 (b).Vertical jumps in this figure indicate transitions between adsorbed dilute and dense phases.Next figures, 1 (c) and 1 (d), show the "liquid-like" (figure 1 (c)) and the "gas-like" (figure 1 (d)) density profiles at coexistence.Under the assumption, V = 0, the profiles of two ionic species are identical.It is important to mention that the profiles, especially for dense phases (figure 1 (c)), exhibit strong depletion close to the pore walls.At lower temperature (T * = 0.065) the profiles attain almost constant value in the central part of the pore.However, at higher temperature (T * = 0.075), which is closer to the critical temperature of the electrolyte in this pore, depletion extends almost to the pore center.Similar calculations have been performed for wider pores (H = 21, H = 31 and H = 61).In figure 2, we demonstrate the density profiles of coexisting phases in a wide pore, H = 61, at low temperature (T * = 0.070) and at high temperature (T * = 0.078) which is close to the critical temperature of the fluid in this pore.Almost constant plateau values of density for both, dilute and dense phases can be observed at low temperature.Depletion of the profiles occurs only close to the pore walls.At high temperature, close to the critical one, strong depletion effects are observed in the dense phase at coexistence.Also, the profiles for dense and dilute phases at this temperature are almost equal close to the pore walls such that phase separation occurs in the central part of this wide pore.Summary of our calculations in the form of phase diagrams is displayed in figure 3. Namely, figure 3 (a) presents ρ * −T * projection, i.e. the relation between the average density of ions in the pore and reduced temperature.The effect of purely geometric confinement on the phase behavior of the RPM is in some sense similar to the effect of such confinement on Lennard-Jones fluid.We are convinced that the coexistence envelopes terminate at the critical point for all pores considered.Precise evaluation of the value for critical temperature and critical density would require very tedious calculations, however.
The effect of geometric confinement yields the lowering of the critical temperature.However, the effect is rather small.The wider is the pore the closer is the coexistence envelope to the bulk phase diagram.However, even in such a wide pore as H = 61, there is a well observed difference between the bulk and confined fluid characteristics with respect to the values of density of coexisting phases.In particular, the value for density of "liquid-like" confined phase is seen lower than the density of similar phase in the bulk.On the other hand, the values for density of dilute "vapor-like" phase are less affected by confinement.The observed similarity of the phase diagrams for a set of confined models and their bulk counterpart can be attributed to the dominant role of long-range electrostatic interaction between ions.
Less pronounced difference between the coexistence lines of the bulk and confined models is seen in figure 3 (b).In that figure we show the µ * − T * projections of the phase diagrams.However, it is clearly seen that the transition in the confined models occurs at higher values of the chemical potential compared to the bulk.In other words, all the discussed first -order transitions describe capillary evaporation of the RPM in slit-like pores with uncharged walls.
We would like to finish this study by some conclusive remarks.This work is only a first step in the investigation of the effects of confinement on the phase behavior of confined ionic fluids.From the very beginning we have focused on the principal issue unsolved previously.Namely, we have established that the presented theory describes first-order phase transition in the restricted primitive model confined by two hard walls that form slit-like pore.Also the theory captures the effects of confinement on the shape of the coexistence envelope and on critical parameters.
However, for the sake of consistency one should start probably from testing the theory with respect to available computer simulation data.A set of simulation data of which we are aware refers to the behavior of the RPM in contact with charged hard wall at temperatures higher than the critical temperature of the bulk model.These data are from the works of Torrie and Valleau [27,28] and from Boda et al. [29,30].We are not aware of simulations concerned with the RPM between hard uncharged plates.Therefore, it is urgent to establish the effects of surface charge density on the walls on the distribution of ionic species in pores as well as on the overall thermodynamic behavior and on the phase behavior of the confined model in particular.At present any kind of simulation data relevant to the coexistence envelope of the RPM under confinement are not available.
Our study permits several interesting and important extensions, however.It seems interesting to compare the behavior of the RPM in slit-like pores with hard walls and in pores with attractive uncharged walls.Solvation forces between plates separated by an electrolyte should be analyzed.Solvent effects at a primitive level can be also included.It is possible to consider various primitive models by taking into account the asymmetry of charges and the asymmetry of diameters of ions.However, an interesting line of research can be followed by including the association effects between ionic species as a supplement to Coulomb interactions between ions, see e.g.[3].These effects essentially improve the agreement of the theoretical bulk phase diagram with computer simulation results and consequently one can expect a more adequate thermodynamic description of adsorption of ionic fluids in pores and of phase behavior in particular.

Figure 1 .
Figure 1.(a) Examples of the dependence of excess grand thermodynamic potential of the system, Ω * ex on reduced chemical potential of one ionic species µ * (µ * = µ * 1 = µ * 2 ), at a given temperature, T * , for the RPM confined in narrow slit-like pores, H = 11.The lines from left to right correspond to T * = 0.65, T * = 0.70 and T * = 0.75.The arrow in each case points to the corresponding equilibrium phase transition.(b) Examples of the behavior of excess adsorption one ionic species on reduced chemical potential µ * .The nomenclature of lines is as in part a. (c) The density profiles of a dense "liquid-like" phase at coexistence at T * = 0.65 (solid line), T * = 0.70 (short-dashed line) and T * = 0.75 (long-dashed line).(d) The density profiles of dilute "gas-like" phase at coexistence with relevant dense phases shown in part c.The nomenclature of lines is as in part c.

Figure 3 .
Figure 3. Phase diagrams for the RPM confined in uncharged slit-like pores.(a): ρ * − T * projections.(b): µ * − T * projections.In both panels bulk phase coexistence is shown by solid line, H = 61 -filled circles, H = 31 -empty squares, H = 11 -empty triangles.The dashed and dotted lines connect symbols for the sake of better visualization.