Application of a density functional approach to nonuniform ionic fluids: the effect of association

In the present paper we discuss a density functional approach for nonuni-form ionic fluids, which takes into account the existence of ion pairs. The theory is based on a fundamental measure theory of hard-spheres, the theory of Gillespie et al., which leads to a more accurate description of the electrostatic part of the grand potential as well as on Wertheim's association theory. The results of model calculations indicate that the inclusion of the associative term in the grand potential leads to the structure of the double layer, which differs from the structure evaluated by neglecting the association. These differences are important at low temperatures only.


Introduction
The theoretical study of electrolytes and in particular association phenomena have been of great interest over the recent years.Ion association has a significant effect on physical and chemical characteristic of electrolytes.In recent years several theoretical methods have been proposed and discussed to take the phenomena of ion association into account [1].
The properties of an electrochemical interface within the so-called restricted primitive model (RPM) and several theoretical approaches have included the RPM model into theories aimed at describing the uniform, as well as the nonuniform electrolytes and molten salts.In particular, these systems have been intensively investigated using computer simulations [2][3][4][5] and integral equations theories [6][7][8][9].
However, the most successful approaches for determining the thermodynamic and structural properties of nonuniform RPM fluids are those based on recent density functional formulations.Density functional theory (DF) determines the structure and thermodynamic properties of a system minimizing the grand potential functional.Usually, the construction of the grand potential is based on a perturbation expansion around a hard-sphere reference system.One of the most accurate DF theories for simple hard-sphere fluid, which is known in the literature as the fundamental measure theory, was developed by Rosenfeld [10].In subsequent years there were several papers published devoted to the application of DF approaches to the description of nonuniform RPM fluids at adsorbing surfaces [11][12][13][14][15][16].
In many cases the DF approach is quite accurate in predicting the structure of nonuniform electrolytes.According to common DF treatment the Helmholtz free energy, F , is divided into the ideal term, F id , the contributions resulting from hardsphere interactions, F HS , and the electrostatic contribution, F el .Of course, the ideal term is known exactly.However, in order to determine all remaining contributions we need to introduce additional approximations.Various formulations of the functional F HS exist [17], including Rosenfeld's [10] and the improved formulations based on it [18][19][20].However, the contributions to the free energy resulting from electrostatic interactions have been evaluated using a perturbational approach with respect to a bulk uniform fluid.It is obvious that not all the systems are capable of perturbing around a bulk fluid.If local densities vary in large amounts within the system, then the bulk-fluid perturbation ansatz is probably not the best choice.
Recently Gillespie et al. [21,22] have drawn attention to this problem.They presented an approximate electrostatic excess free energy functional for charged hard sphere fluids.Their method is based on a perturbation expansion concerning a suitably chosen position-dependent reference fluid.This method accurately reproduces the results of Monte Carlo simulations [21,22].
Real bulk electrolyte solutions consist of free ions and ionic pairs and even more complex ionic aggregates.Therefore the description of thermodynamic properties of electrolyte solutions calls for applying the associative concept of Bjerrum [23] and Ebeling and Grigo [24] to the RPM.In fact, the ideas of Ebeling and Grigo have already been used to improve the equation of state calculations of bulk fluids [1,[25][26][27][28][29].In accordance with this approach, the density of ions is divided into the density of free and bound parts, which are connected by the mass action law.
Associative concept of nonuniform ionic fluids can be realized by two different methods.According to the first approach only unbounded ions are taken into account, i.e. the total fluid density is modified by subtracting the concentration of associated ions, which is obtained from the mass action law [30].The second method starts from the associative mean spherical approximation [27] that is based on Wertheim's association theory [31] of inhomogeneous chemically associating fluids [32].The incorporation of association into integral equation theory of uniform, as well as nonuniform electrolytes has been reported in [33][34][35][36].
Recently, Yu and Wu have focused attention on the case of fluids interacting via short-range associative forces [32].They have proposed an extension of the fun-damental measure density functional theory [10] to the case of fluids interacting via short-range associative forces [32].Their theory has been successfully applied to the description of nonuniform non-ionic fluids, exhibiting close-loop demixing phase diagrams in bulk systems [37,38].
In this paper we report the study on the behavior of nonuniform associating ionic fluids in contact with a charged surface.We propose a density functional approach, which combines the ideas of Gillespie et al [21,22] and Yu and Wu [32].To our best knowledge, the association effects in nonuniform electrolytes have been considered so far within the framework of the associative mean spherical approximation only, cf.[34,35] and the present work is the first attempt to incorporate the association of ions into density-functional approach.This work is also the first step towards the development of a theory which would predict phase transitions in confined electrolytes.

Theory
The model system in this study consists of a two-component mixture of ions (a generalization to n-component mixture is straightforward).The electrolyte is assumed to be a fluid of charged hard spheres, all of the same diameter d and of valences z i , i = 1, 2. The fluid is in contact with a charged hard wall.The interaction between the ions of species i and the wall is given by where w i (z) and v i (z) are the electrostatic and the non-electrostatic (van der Waals) parts of the external potential, respectively.The van der Waals interaction potential is where z is the distance from the surface.The electrostatic interaction between an ion and the surface is given by In the above e is the electron charge, z i is the valency of ions of i-th type, q is the surface charge density and is the dielectric constant.
In the present paper we apply the theory, which adapts the approach of Gillespie et al. [21,22].The first step toward the development of the grand potential of the system as a functional of the local density, ρ i (r), is its decomposition into: ideal (id), excess hard-sphere (hs), electrostatic (el), and associative (as) terms: The term Ω id is where e is the electron charge and Ψ(r) is the electrostatic potential.Similarly to [21,22] the last equation also involves non-ideal contributions to the configurational chemical potential µ i .The electrostatic potential Ψ(z) is obtained by integrating [39] Poisson's equation: The hard-sphere component of the grand potential is calculated from the version of the fundamental measure theory [10], proposed by Roth et al. [40] and by Yu and Wu [41] Ω hs with In the above ξ(r) = |n 2 (r)|/n 2 (r) and n α , α = 0, 1, 2, 3, denote the weighted densities defined as sums of the spatial convolutions of the average densities and the corresponding weight functions with the weight functions w α (r) given in [10,[18][19][20][21][22].
The electrostatic part of the grand potential is approximated by [21,22]: z i e drΨ(r)ρ i (r), (10) where ∆ρ i (r) = ρ i (r) − ρ ref i (r) and ρ ref i (r) is a suitably chosen position-dependent reference system density.c ref ij (r, r ) is a part of the reference system two-particle direct correlation function, which results from electrostatic interactions.
According to Gillespie et al. [21,22] the reference system density is determined by defining an auxiliary system.This system locally has the same ionic strength as the system under study and locally satisfies the condition of electroneutrality.For an 1:1 electrolyte the density of this system equals to ρi (r) = [ρ 1 (r) + ρ 2 (r)]/2.The reference system density is evaluated by averaging the densities ρi (r) over a sphere of the radius R f with the weight function The radius R f is position-dependent and for the 1: and T * is the usual reduced temperature, T * = kT d/e 2 .The definitions ρi and R f for an arbitrary (multicomponent) system are given in [21,22].
The evaluation of the radius R f requires the knowledge of the reference system density and vice versa.Therefore, their determination must be carried out using a self-consistent procedure.Gillespie et al. have proposed an application of an additional iterational method, which is described in detail in [21,22].
The short-range electrostatic part of the two-particle direct correlation function, resulting from the MSA approximation kT c is the complement of the electrostatic part of grand potential term.In the above equation λ i = d/2 + 1/2Γ is the capacitance length of ion species (obviously, in the system in question λ i = λ j ).
The last term of the grand potential arises from association.The associative contribution is approximated by using Wertheim's treatment [31].This approach, which has been widely used [37,38,42] to describe association in simple nonuniform fluids leads to the expression [41] Ω as [{ρ(r In the above the dissociation constant is The density of ions, ρ b , is divided into the density of free ions, ρ b0 and bounded pairs, ρ bb , which are connected by the mass action law [30] 1 where α = ρ b0 /ρ b , ρ b = ρ b0 +ρ bb is the degree of dissociation and K is the dissociation constant.The dissociation constant depends on density and temperature and at a given temperature it can be separated into two terms: into the density-independent and the density dependent terms.From the Ebeling's theory [3,24,29] we obtain the density-independent term The density-dependent term, however, is given by with g hs (d) being the contact value of the hard-sphere radial distribution function between unlike ions.This part has been evaluated from the Carnahan-Starling equation of state (cf.[1]).The electrostatic part of the radial distribution function between unlike ions, g el , arises from the MSA approach In the theory presented above two sets of averaged densities are applied: the set {n α } -to calculate the hard-sphere contribution to the grand potential and the densities {ρ ref i } -to calculate the electrostatic contribution.The association is the result of electrostatic interactions.Thus, to be consistent, the contact values of the radial distribution function have been evaluated for the set of densities {ρ ref i }.The derivative of Ω as is The expressions for the functional derivatives of Ω hs and Ω el are given in [10,18,21].Requiring the grand potential to be at a minimum, δΩ/δρ i (r) = 0, we obtain equilibrium density profiles.2) density profiles at a charged wall resulting from the present theory (dashed and solid lines, respectively), from the theory of Boda et al. [44] (dotted and dash -dotted lines) and from Monte Carlo simulations (filled and empty circles) [44].(b) A comparison of theoretical predictions with simulation data of Torrie and Valleau [4,5].All the parameters are given in the figure .The knowledge of one-dimensional profiles ρ i (z) allow us to calculate the adsorption isotherms and the surface charge density z i e ρ i (z)dz.(22) Our calculations were proceeded on a grid of the size of 0.02d.To evaluate the density profiles we have employed Piccard iterational procedure [43], as described in [21,22].The input quantity to our calculations was the value of the electrostatic potential at the wall Ψ * 0 = Ψ(z = 0)(4πed/ ) and the value of the charge at the wall was calculated from equation (22).All the calculations were continued until the maximum difference between two consecutive profiles was smaller than 10 −7 percent.Our calculations are carried out for symmetric 1:1 nonuniform mixtures.

Results and Discussion
In order to test the method we first perform some comparisons of the structural properties of the system with computer simulation data.Figure 1a presents a com-parison of theoretical density profiles of ions at ρ b = 0.04 and T * = 0.15, in contact with a charged wall q * = qd 2 /e = 0.00765 with Monte Carlo simulation data [44] and with previous calculations based on an earlier version of the density functional theory, described in detail in [44].We observe that the present theory provides a better agreement with simulations for this thermodynamic state than the previous version of the DF theory [44].Next comparisons are given in figure 1b.Here we show the results for the density profiles of ions in contact with a single charged wall at T * = 0.59492.The Monte Carlo data were taken from the works by Torrie and Valleau [4,5].In general, the present theory describes the distributions of ions reasonably well.However, the temperature is quite high in this case and under such conditions the present approach reduces to that of Gillespie et al. [21,22].
Since the association of ions plays an important role at low temperature only we concentrate our discussion on the results obtained at the reduced temperatures lower than 0.5.Here we have plotted the density profiles of co-ions and counterions.The effect of the association on the density profiles of counterions is smaller than on the profiles of co-ions.Figures 2a and 2c present the profiles obtained at a relatively low reduced temperature T * = 0.09, whereas figures 2b and 2d -at higher temperature, T * = 0.3.An increase of the temperature changes the structure of the solid-fluid interface.At a higher temperature we can observe "the smoothing" of the shape of the obtained profiles.As we can see in figures 2a (and 2c, respectively) low temperature causes a high concentration of the ions of a given type and causes "a sucking" of the ions of an opposite sign.A competition between repulsive ion-wall force and a tendency to form ionic pairs results in the appearance of the local density peak of co-ions at a distance of z * = z/d ≈ 0.8 from the wall.Since the first maximum of the counterions concentration is at z/d = 0.5, the shift of the location of the first local density peaks for counter-and co-ions suggests that the ionic pairs do not assume orientation with their axes parallel to the surface, but rather the associates are tilted with respect to the surface.An increase of the value of the electrostatic potential leads to the growth of concentrations of counterions at the wall and causes the second distinct local density peak of counterions.This is observed in figures 2a and 2c.
Figures 3a and 3b present the fraction of unbounded ions, α, for the systems in figure 2a and figure 2b respectively.We can observe that a decrease of temperature causes an increase of the association.In both cases, directly at the wall the fraction of unbounded co-ions is high, whereas the fraction of unbounded counterions is rather small.This is the result of a high wall-ion repulsion, which causes "the absence" of co-ions in the nearest vicinity of the surface.To complete the description of the studied system we present examples of the adsorption isotherm.Figure 4a shows the isotherms for the same value of electro-static potential at the wall Ψ * 0 but obtained at different reduced temperatures.We can see only small differences for the isotherm obtained at temperature higher than T * = 0.15.The illustration of the effect of the value of electrostatic potential on the adsorption is plotted in figure 4b.As we expected, an increase of Ψ * 0 causes the increase of adsorption of the system investigated.In figure 4c we compare isotherms obtained for Ψ * 0 = 0.1 at low temperatures with those obtained for the system neglecting association.At lower temperature, T * = 0.09 the isotherm obtained with the association "switched on" and "switched off" differ significantly.The differences decrease with an increase of the temperature.
In this paper we have applied the density functional theory to the study of nonuniform associating ionic fluids.The obtained results indicate that at low temperatures the inclusion of the associative term in the grand potential leads to the structure of solid-fluid interface.This structure differs from the structure evaluated for the system neglecting association, but the differences become negligibly small at high temperatures.

Figure 1 .
Figure 1.(a) A comparison of the counterion (1) and co-ion (2) density profiles at a charged wall resulting from the present theory (dashed and solid lines, respectively), from the theory of Boda et al.[44] (dotted and dash -dotted lines) and from Monte Carlo simulations (filled and empty circles)[44].(b) A comparison of theoretical predictions with simulation data of Torrie and Valleau[4,5].All the parameters are given in the figure.

Figure 3 .
Figure 3. Dependence of α(z) for counterions (solid) and co-ions (dashed) on the distance from wall.The curves in part (b) are for the same bulk densities as the relevant curves in part (a).All the parameters are given in the figure.

Figure 2
Figure2illustrates how the temperature and the value of electrostatic potential at the wall Ψ * 0 effect the structure of the solid-fluid interface.Here we have plotted the density profiles of co-ions and counterions.The effect of the association on the density profiles of counterions is smaller than on the profiles of co-ions.Figures2a and 2cpresent the profiles obtained at a relatively low reduced temperature T * = 0.09, whereas figures 2b and 2d -at higher temperature, T * = 0.3.An increase of the temperature changes the structure of the solid-fluid interface.At a higher temperature we can observe "the smoothing" of the shape of the obtained profiles.As we can see in figures 2a (and 2c, respectively) low temperature causes a high concentration of the ions of a given type and causes "a sucking" of the ions of an opposite sign.A competition between repulsive ion-wall force and a tendency to form ionic pairs results in the appearance of the local density peak of co-ions at a distance of z * = z/d ≈ 0.8 from the wall.Since the first maximum of the counterions

Figure 4 .
Figure 4. Adsorption isotherms (a) obtained for Ψ * 0 = 0.1.Part (b) presents isotherms obtained at low temperature for different values of Ψ * 0 .Part (c) shows comparison results evaluated by "switching on" (solid) and "switching off" (dashed) association, respectively.All the parameters are given in the figure.