On the contact values of the density profiles in an electric double layer using density functional theory

A recently proposed local second contact value theorem [Henderson D., Boda D., J. Electroanal. Chem., 2005, 582, 16] for the charge profile of an electric double layer is used in conjunction with the existing Monte Carlo data from the literature to assess the contact behavior of the electrode-ion distributions predicted by the density functional theory. The results for the contact values of the co- and counterion distributions and their product are obtained for the symmetric valency, restricted primitive model planar double layer for a range of electrolyte concentrations and temperatures. Overall, the theoretical results satisfy the second contact value theorem reasonably well, the agreement with the simulations being semi-quantitative or better. The product of the co- and counterion contact values as a function of the electrode surface charge density is qualitative with the simulations with increasing deviations at higher concentrations.


Introduction
One of the more interesting recent developments in the electric double layer research has been the advancement of contact value theorems involving the charge profile in a primitive model (PM) planar double layer (charged hard spheres moving in a dielectric continuum next to a planar electrode) (see, for example, references [1][2][3][4]). Such exact conditions, or sum rules as they are often called, are important per se in statistical mechanics since they permit unambiguous assessment of various approximate theories and hence aid in theoretical development.
The most famous contact theorem in the double layer literature is the one formulated by Henderson and Blum [5], and Henderson, Blum, and Lebowitz (HBL) [6] over thirty years ago. It is a condition on the contact value of the total density profile in a planar double layer, and for a symmetric valency restricted primitive model (RPM) (equisized ions in the PM) planar double layer -the model system of interest in this paper, the HBL relation reads g sum (d/2) = [g co (d/2) + g ctr (d/2)]/2 = a + b 2 2 . (1) Here g co , g ctr are the co-and counterion singlet distribution functions, d is the common ionic diameter, and a = p/(ρk B T ) is the bulk osmotic coefficient with p being the bulk pressure, k B the Boltzmann constant, and T the absolute temperature. The quantity b = zeσ/(ǫ 0 ǫ r k B T κ) is a dimensionless parameter where z, e, and σ are, respectively, the absolute value of the ionic valency, the magnitude of the elementary charge, and the uniform surface charge density on the electrode, ǫ 0 is the vacuum permittivity, and ǫ r is the relative permittivity of the continuum solvent. Also, κ = z 2 e 2 ρ/(ǫ 0 ǫ r k B T ) is the Debye-Hückel parameter (inverse Debye screening length) with ρ = i ρ i where ρ i is the mean number density of the i th ionic species. In the RPM case the contact distance, that is, the distance of closest approach of an ion to the electrode, occurs at d/2, where d co = d ctr = d. Note that g sum (x) (x is the perpendicular distance from the electrode into the solution) is related to the total density profile where ρ i (x) is the singlet density profile of the i th species. It is of interest that the second term in equation (1) is just the Maxwell stress. Although equation (1) was obtained from statistical mechanics it is consistent with Maxwell's equations.
Equation (1) is a local expression and the consequent ease of its use has made the HBL contact condition very appealing in double layer research over the years. For example, the classical Gouy-Chapman-Stern (GCS) [7][8][9] theory of the double layer satisfies equation (1) but with a = 1, the ideal gas value. Thus, for an electrolyte with osmotic coefficient substantially different from unity, the GCS theory can lead to appreciable errors especially at low surface charges.
Sum rules such as equation (1) are useful not only for assessing theories but also for the insight they provide. For example, because the coion contact value becomes small at a large surface charge and the counterion contact value becomes large, according to this equation, the latter contact value increases as square of the surface charge density. On the other hand, the local electroneutrality condition (also a sumrule) requires that the area of the charge profile be equal but opposite in sign to the electrode charge density. As a result this area increases linearly with the electrode charge density. Consequently, at large electrode charge, the oscillations and charge inversions in the charge profile should diminish but cannot disappear altogether relative to the contact value.
Analogous relations for the contact value of the total charge profile in the double layer -the theme of the present paper, have been relatively recent. A formal, rigorous relation was derived by Holovko et al. [1,10] and Holovko and di Caprio [2] using the Bogoliuobov-Born-Green-Yvon hierarchy. For symmetric valency RPM planar double layer their expression is as follows: where β = 1/(k B T ), ψ(x) is the mean electrostatic potential, and g diff (x) = [g ctr (x) − g co (x)]/2. Again, g diff (x) is now related to the total charge profile, viz., with z i being the valency of the ionic species i and z = z co = −z ctr . The definition of g diff (x) is, for convenience only, designed to make this quantity positive in general. Since the use of this equation implies a knowledge of ψ(x), and g i (x) throughout the double layer, the expression is non-local.
Independently, Henderson and Boda (HB) [3] have proposed an approximate, local expression for g diff (d/2) at low electrode charges from empirical considerations, viz., To date a formal, analytic connection between equations (3) and (5) remains obscure, although in a later paper Holovko et al. [11] have outlined a very approximate connection. In a series of papers Henderson and Bhuiyan [4] and Bhuiyan and co-workers [12][13][14][15] have tested equation (3) against exact Monte Carlo (MC) simulation data for a spectrum of physical states including asymmetric electrolytes [13][14][15] and found the equation to be remarkably consistent with the simulations. Theoretical support came from an application of the modified Poisson-Boltzmann (MPB) equation, which was found to satisfy equation (3) to a very good degree [12,13]. Bhuiyan and Henderson [16] have also compared the two relations (equations (2) and (3)) numerically using simulations and the conclusion was that although exact, equation (2) is difficult to implement numerically because of its non-local nature. We note here that an approximate, non-local relation for g diff (d/2) has also been suggested by Henderson and Bhuiyan [17]. Following convention, we will call equation (1) the first contact value theorem, while equations (3) and (5) represent two versions of the second contact value theorem.
Another interesting recent result that also concerns us in this paper is the behavior of the product of the co-and counterion contact values f = g co (d/2)g ctr (d/2) in the RPM planar double layer. The classical GCS result for this quantity is strictly unity under all circumstances and thus constitutes a basis for the classical theory. For example, the value of the counterion contact g ctr (d/2) is as high as the reciprocal of the coion contact g co (d/2). However, the corresponding simulation results [18] dramatically alter the classical picture. The product f is seen to be not only different from unity, but also that its characteristics as a function of the electrode charge change with the salt concentration. At low concentrations there is a maximum before f becomes vanishingly small at high electrode charge. As the concentration increases, the height of the maximum decreases and at sufficiently high concentrations the maximum disappears completely with f decreasing monotonously. Again, theoretical support for such a behavior of f came from the MPB [18] and although the hypernetted chain/mean spherical approximation theory does not show a maximum, the product f does become very small when the electrode charge is large [19].
In this study we propose to utilize the HB second contact value theorem and the existing MC simulation results from the literature for f to assess the density functional theory (DFT) of the planar double layer. The DFT has been one of the more successful theories of the electric double layer phenomenon and compares favorably with the MPB across planar, cylindrical, and spherical geometries (see for example, references [20][21][22]). Early applications of the DFT to the planar double layer were made by Tang et al. [23] and Mier y Teran et al. [24]. Later Rosenfeld's [25] techniques were utilized by Mier y Teran et al. [26] and Boda et al. [27][28][29][30]. For even recent publications on application of the DFT to the planar double layer, we refer the interested reader to the works by Gillespie et al. [31,32], Valiskó et al. [33], Wang et al. [34], Yu et al. [35], and Pizio et al. [36]. Since there is more than one version of the DFT for the planar double layer, in the next section we will briefly outline the DFT method used in this paper. Results will be shown in section 3, and some conclusions drawn in section 4.

Molecular model
As indicated in the previous section, the model double layer system consists of a binary, symmetric valency RPM next to a non-penetrable, non-polarizable, uniformly charged planar electrode with a surface charge density of σ. Since for a given salt concentration, solvent dielectric constant, and temperature, b has a linear dependence on σ, it is often convenient to specify σ in terms of b.
The ion-ion interaction potential in the Hamiltonian is thus where r is the distance between a pair of ions. We also assume that the dielectric constant, ǫ r , is uniform throughout the entire system. The bare interaction between an ion of species i and the wall is given by where v i (x) and w i (x) are the non-electrostatic and electrostatic (Coulombic) parts of the ion-wall potential. The non-electrostatic contribution is a hard-wall potential The electrostatic part w i (x) is given by (9)

Density functional theory
The essence of the density functional theory (see for example, reference [37]) is that an expression for the grand potential, Ω, as a functional of the singlet density profiles, ρ i (x), of each of the species i , is initially constructed. At equilibrium the grand potential is minimal with respect to variations in the density profiles, viz., This condition is then used to calculate the density profiles and other relevant quantities like the free energy.
In the density functional theory the grand potential of an inhomogeneous fluid can be written in the form [30,36] where µ i denotes the chemical potential of species i . The free energy functions F ({ρ i }) is decomposed into ideal (id), hard-sphere (hs), and electrostatic (el) terms as follows F ( The ideal term is known exactly For the hard-sphere term, however, we apply the expression resulting from a recent version of the Fundamental Measure Theory [38], with the free energy consisting of the terms dependent on scalar and vector weighted densities, for details see reference [36]. Following Pizio et al. [36] electrostatic contribution to the free energy, F el ({ρ i }), is represented by where {ρ i (x)} denotes a set of suitably defined inhomogeneous average densities of a reference fluid. One of the simplest possible choices of f el ({ρ i (x)}) is to apply the expression resulting from the MSA equation of state evaluated via the energy route, namely [39] For a symmetric valency situation as in the present case the reduced temperature is T * = 4πk B T ǫ 0 ǫ r d e 2 z 2 Moreover, The inverse Debye screening length κ can be cast in terms of T * The last three expressions above correspond to an electroneutral fluid, so that the construction of the averaged densitiesρ i (x) at the electroneutrality condition is satisfied. In our approach we follow the development proposed by Gillespie et al. [31,32] described briefly below.
Let us define the weighted densitiesρ i (x) as where W (|r − r ′ |) is a weight function. Gillespie et al. made the assumption, viz., where θ(|r − r ′ |) is the step-function. The radius of the sphere over which averaging is performed, R f , is approximated by the "capacitance" radius, that is, by the ion radius plus the screening length In addition, Gillespie et al. [31,32] required that the fluid with the densities {ρ i (x)} have the same ionic strength as the system with weighted densities, {ρ i (x)}. Consequently, in the case of a symmetric 1:1 electrolyte the averaged densities {ρ i (x)} are given bȳ Because equations (15) and (17)- (20) are coupled, the evaluation of R f requires an iteration procedure.
This iteration loop has to be carried out in addition to the main iteration procedure for evaluating the density profiles.
The mean electrostatic potential ψ(x) is determined by the Poisson equation The integration of the Poisson equation is carried out subject to the boundary conditions lim z→∞ ψ(x) = 0 and lim x→∞ ψ ′ (x) = 0.
Having specified all the contributions to the free energy functional, the requisite density profiles can be obtained by minimizing the grand potential (cf. equation (11)).
All the details of our approach can be found in reference [36].

Results and discussion
The DFT equations have been solved numerically using the established methods (see for example, references [23,36,40]. We will also present the classical GCS results for comparison purposes, which for the RPM case can be obtained analytically. It is convenient to discuss the results in terms of universal reduced parameters such as the reduced density ρ * = i ρ i d 3 i and the reduced temperature T * defined earlier. Calculations were done at two different reduced temperatures, T * = 0.150 and 0.595, respectively, and at each reduced temperature a number of physical states were treated. The value of the ionic diameter was kept at d = 4.25 ×10 −10 m throughout. Although a 1:1 valency system was used in the actual calculations, in view of universality of T * this becomes a moot point since for a given T * a 1:1 system at T is equivalent to a 2:2 system at 4T . For example, in the specific case of T * = 0.15, a 1:1 valency case corresponds to ∼75 K, while a 2:2 valency case corresponds to ∼300 K. In implementing the HB contact condition, Henderson and Bhuiyan [4] found it convenient to recast equation (5) in the form for symmetrical valency electrolytes. We have followed the procedure here. We note though that a straightforward linear plot of equation (5) with a as the slope has also been done [16]. In    of symbols and notation as in figure 1. MC data from reference [18]. the trend that the DFT g diff (d/2, b)/b (upper panels of the figures) follows the corresponding simulation results very closely for not too high b for the range of concentration treated. Only a very slight discrepancy is seen at b = 0, which is a consequence of the fact that the DFT does not satisfy the HBL first contact theorem exactly. The classical GCS theory satisfies equation (5) but with a = 1, the ideal gas value. This is clear from the figures and for ρ * = 0.02, 0.03, 0.05, and 0.10 ( figures 1-4), where a is somewhat less than unity, the GCS theory leads to deviations from the MC data. The results at a different temperature T * = 0.595 and at ρ 3 = 0.00925 (c = 0.1 mol/dm 3 ) and ρ 3 = 0.0925 (c = 1.0 mol/dm 3 ) are shown in figures 7 and 8, respectively. Here too the trends shown by the DFT g diff (d/2, b)/b and their agreement with the corresponding simulations are similar to that seen in figures 1-6. We note that the HNC satisfies equation (1) with the first term being a function of the hard sphere compressibility, and equation (5) with a = 1 in the first term. For contact values, it is little better than the GCS theory.
(see for example, equation (18) of reference [18]). This equation is exact in the limit b →0. It is easy to see at b = 0 that the value of f depends on the value of a. Furthermore, the initial slope of f is negative for a < 1 (figures 1-4), the initial slope is approximately zero and the plots are initially flat when a ∼ 1 (figures 5 and 7), and the initial slope is negative when a > 1 (figures 6 and 8). Note again that since in the GCS theory a = 1, the right hand side of equation (19) is unity and the initial slope is zero being consistent with the observations.   Figure 1. MC data from reference [12].
An important property of the simulation data is that the contact product f tends to very small values at large values of b. As the surface charge increases the coion population near the electrode is depleted, while the counterion population increases. However, the latter also induces packing problems that inhibit distant counterions from migrating too close to the electrode surface. All these lead to the observed behavior of f . In the GCS theory however, the decrease in g co is always proportional to the increase in g ctr so that classically f = 1 consistently. The DFT f generally follows the MC trend in figures 2-8. Although in figure 1 the lack of DFT data beyond b = 8 implies that one cannot be definitive, in view of the results in the rest of the figures, it is a fair conjecture that here also the contact product will assume small values at still higher values of b.

Conclusions
In this paper we have examined the predictions of a density functional theory of the planar electric double layer with regards to (i) the HB second contact value theorem, and (ii) the behavior of the product of the DFT contact values of the co-and counterion distributions vis-a-vis exact MC simulation data from the literature. The principal finding regarding (i) is that generally the DFT follows the MC results very closely over the range of concentrations and temperatures studied. There is only just a hint of discrepancy at b = 0, which is probably tied to the approximation used for the hard-sphere term in the free energy functional used to construct the grand potential. By contrast, the GCS results show greater deviations from the simulations, especially at lower concentrations when the MC a is less than unity. It is of interest to note that the degree to which the DFT satisfies the HB contact condition is very similar to what some of us have observed with the MPB theory [12,13] with both of the approaches showing slight deviations at b = 0. This is not surprising since none of the theories satisfies the HBL first contact value theorem exactly.
With respect to (ii) above, on the other hand, our calculations indicate that the DFT is broadly in qualitative agreement with the characteristics of the simulations, with the product f of the contact values tending to small values with increasing surface charge on the electrode. A maximum in f as a function of b seen at lower concentrations although this occurs as, what might be termed, a delayed maximum . Importantly though, the behavior of the the initial slope of f as the electrolyte concentration increases follows the MC trend. Admittedly, however, there is a quantitative discrepancy between the DFT results and the MC data beyond c ∼ 1 mol/dm 3 . This is understandable in view of the fact that the product of the contact values of the distributions can be a rather more sensitive quantity than their difference so that a slight error in either of the contact values tends to become magnified in the contact product.
The version of density functional that is employed in this paper gives good results for the contact values but is less satisfactory in predicting oscillatory profiles. In contrast, other versions of density functional theory [35,40] are better at producing oscillatory profiles but are less successful for contact values. There is more to be done in the development of a fully satisfactory density functional theory.