Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Collect Czechoslov Chem Commun. Author manuscript; available in PMC 2010 July 20.
Published in final edited form as:
Collect Czechoslov Chem Commun. 2010 January 1; 75(4): 425–446.
doi:  10.1135/cccc2009098
PMCID: PMC2906829



The properties of the singlet ion distributions at and around contact in a restricted primitive model double layer are characterized in the modified Poisson–Boltzmann theory. Comparisons are made with the corresponding exact Monte Carlo simulation data, the results from the Gouy–Chapman–Stern theory coupled to an exclusion volume term, and the mean spherical approximation. Particular emphasis is given to the behaviour of the theoretical predictions in relation to the contact value theorem involving the charge profile. The simultaneous behaviour of the coion and counterion contact values is also examined. The performance of the modified Poisson–Boltzmann theory in regard to the contact value theorems is very reasonable with the contact characteristics showing semi-quantitative or better agreement overall with the simulation results. The exclusion-volume-treated Gouy–Chapman–Stern theory reveals a fortuitous cancellation of errors, while the mean spherical approximation is poor.

Keywords: Electric double layer, Contact value theorems, Density profile, Charge profile, Monte Carlo

Alongside exact numerical simulations, exact analytical conditions are invaluable in statistical mechanics of Coulomb fluids in comparative assessment of different theories and/or further theoretical development. These exact conditions can also help provide additional, useful checks on the consistency of simulations. One such condition in the electric double layer phenomenon is our principal concern in this paper. An electric double layer is formed when a charged electrode is brought in contact with a charged fluid and an ionic atmosphere develops in the vicinity of the electrode. The phenomenon has practical significance for a spectrum of systems in biology and industrial chemical processes. For a very recent review of aspects of the electric double layer theory, see ref.1

The analytical conditions we are interested in involve the contact values of the electrode-ion distributions in a planar geometry, that is, a planar double layer. One condition, formulated by Henderson and Blum2, and Henderson, Blum and Lebowitz3 (HBL) over thirty years ago, gives an expression for the contact value of the density profile in such a situation. For a primitive model (PM) planar double layer (spherical rigid ions of arbitrary diameters and charge moving in a dielectric continuum (solvent) next to a uniformly charged flat electrode), the condition reads


Here ρs, ds and gs are, respectively, the average number density, the diameter and the electrode-ion distribution function of ionic species s, p is the bulk osmotic pressure and σ is the uniform surface charge density on the electrode with ε0 and εr, respectively, being the vacuum permittivity and relative permittivity of the solvent. Also, kB is the Boltzmann constant and T the absolute temperature. The relation was derived from force balance considerations at the electrode–electrolyte interface. It is local, easy to implement and, as a consequence, has been useful in checking theoretical descriptions of the electric double layer over the past three decades1.

The other, relatively recent condition concerns the contact value of the charge profile in the planar double layer. This was derived by Holovko, Badiali and di Caprio4,5 (HBC), and Holovko and di Caprio6 starting from the Bogoliubov–Born–Green–Yvon (BBGY) hierarchy of equations7. For a double layer containing a restricted primitive model (RPM) electrolyte (the ions of the PM are now restricted to having a common diameter), and for symmetric valency salts – our interest in this paper, their relation is particularly simple, viz.,


where gsum(x) = (1/2)(gctr(x) + gco(x)), gdiff(x) = (1/2)(gctr(x) − gco(x)), the subscripts ‘ctr’ and ‘co’ denoting counter- and coions, respectively. The quantity ψ(x) is the mean electrostatic potential at a perpendicular distance × from the planar electrode, e is the elementary charge, β = 1/(kBT), d is the common ion diameter, and z the absolute value of the ion valency. Note that gsum(x) and gdiff(x) are related to the total density and charge profiles through ρ(x) = ρgsum(x) and q(x) = −zgdiff(x) with ρ = Σsρs. In contrast to the HBL contact condition the above relation is non-local and requires, for its implementation, a full knowledge of the potential and density profiles across the double layer. It is worth mentioning here that a non-rigorous, non-local relation for gdiff(d/2) has also been proposed by Henderson and Bhuiyan8.

Simultaneously, a further empirical local contact condition for gdiff(d/2) has been advanced by Henderson and Boda9 (HB), viz.,


for small b. Here a = p/(ρkBT) is the bulk osmotic coefficient, b = zeσ/(ε0εrkBTκ) is a dimensionless parameter, and κ is the Debye–Hückel constant, κ=z2e2ρ/ε0εrkBT. This was obtained primarily through intuition while examining Monte Carlo (MC) simulation data10. Details of the steps leading to Eq. (3) have been given in refs8,11,12. Some of us have been involved in assessing both these local and non-local relations using MC simulation data8,1113. What emerged was a wealth of evidence encompassing a host of physical states, which suggest the Eq. (3) to be remarkably accurate and perhaps exact as b → 0. Theoretical support came from the modified Poisson–Boltzmann (MPB) theory, which predicted results for some 1:1 and 2:2 valency13 electrolytes. In a later analysis involving 2:1/1:2 valency electrolytes, both MC simulations and the MPB theory results further testified to the validity of the HB contact condition14.

The situation was somewhat different with the HBC contact formula (Eq. (2)). Because of the non-local nature of this relation and inherent scatter in the MC simulation data, it was difficult to obtain accurate results through numerical integrations8. Nonetheless, within the numerical limitations, the broad conclusion was that the simulations were consistent with the HBC expression. However, a direct relationship between Eqs (2) and (3) remains elusive.

The seminal classical Gouy–Chapman–Stern1517 (GCS) theory of the double layer satisfies all of the above contact value theorems but with the caveat a = 1. This simply states the fact that the point ions in this theory behave like an ideal gas with no correlations among them. As indicated previously, of the formal theories the MPB was seen to satisfy the HBL contact condition to a very good degree although a limited number of states was probed13. There is evidence from earlier MPB studies18 that the MPB a calculated through Eq. (1) at b = 0 shows a good consistency with the corresponding bulk osmotic coefficients.

In the light of these encouraging results, it is tempting to try and investigate further how the MPB theory fares with regard to the contact value theorems over a broader range of physical states. In particular, it would be useful and of interest to see how well the theory satisfies the HBC contact value theorem, which has not been attempted to date. This is the focus of this paper. Besides the contact theorems, Eqs (2) and (3), we will also look at other contact quantities, viz., gco(d/2) and gctr(d/2) to obtain a comprehensive picture. Previous works8,19 suggested a global view is necessary rather than focus on individual quantities. Two other theories, an exclusion-volume-corrected GCS theory and the mean spherical approximation will also be considered. We will also supplement the theoretical calculations with fresh MC simulations of more states.


We will be employing the restricted primitive model of planar double layer as the model system throughout this work. We will also confine our attention, in general, to binary, symmetric z:z valency electrolytes in water-like solvent. The relative permittivity of the medium of the electrode is taken to be the same as that of the solvent so that no polarization charges occur at the interface. In this model, the HBL contact formula assumes the simple form


Thus, gsum(d/2,b = 0) = a gives the osmotic coefficient of the bulk fluid, and hence is a measure of the internal consistency of a theory. In the earlier works Henderson and Bhuiyan12 and Bhuiyan et al.13 found it convenient to use the HB condition (Eq. (3)) by writing it in the form


We will follow the same convention here.

To treat the HBC contact condition, we consider the identity




Here zs = z+ or z is the valency of a cation or anion, with g+ and g the corresponding cation and anion singlet distribution functions. Equation (6) is equivalent to Eq. (9) of ref.4, while for the GCS closure we have that G = 0. Integrating Eq. (6) for a z:z salt, we have


so that the integral of G gives the deviation or correction to the Holovko et al. result for an approximate theory. From Eq. (6) we can deduce that the HBC contact condition always holds for the GCS theory.

The Modified Poisson–Boltzmann and the Exclusion-Volume-Corrected Gouy–Chapman–Stern Theories

The MPB formalism in the electric double layer theory is an electrostatic potential energy approach that endeavours to incorporate the principal features neglected in the traditional GCS theory, namely, the fluctuation potential and an ionic exclusion volume term (EVT). To date, the MPB approach remains one of the more successful statistical mechanical theories of the double layer in planar1,20, spherical21, and cylindrical22 geometries. The MPB theory utilized in this work is at the MPB5 level, which is presently the most accurate formulation of the theory. Its development has been detailed elsewhere in the literature (see, for example, ref.23) and will not be repeated here.

Another theory used here is the exclusion-volume-corrected GCS theory. A class of such exclusion-volume-treated mean field theories forms the basis of the lattice-gas Poisson–Boltzmann theories and can be useful in physical situations where steric effects are dominant (see, for example, ref.24). In the present context, if we neglect the fluctuation potential term in the MPB theory, what remains is the GCS theory but with an added exclusion volume term. For example, we have the Poisson equation for ψ(x),


where the electrode-ion distribution function gs(x) is now given by


The exclusion volume term ξs(x) used here was developed from the BBGY hierarchy of equations25, viz.,


where H(x) is the Heaviside function or unit step function, and ξ = ξs for a z:z RPM electrolyte. In the rest of the paper the theory specified by Eqs (9)(11) will be denoted by GCS+EVT. An alternative form of the exclusion volume term was developed by Outhwaite and Lamperski26, which in many instances gave similar results as that with Eq. (11)27,28. From Eq. (8), we have found that the GCS+EVT theory deviation from the HBC contact condition is


When b = 0, ξ(d/2) is the BBGY uncharged hard sphere–hard wall singlet distribution function gBBGY(d/2), so


The Mean Spherical Approximation

The mean spherical approximation (MSA) for the RPM planar double layer was developed by Blum29. Its chief advantage is that it affords analytical expressions for the contact quantities, the potential, and the density profiles, which makes the theory very easy and convenient to use. For example, for a binary charge and size symmetric electrolyte, the contact values of the electrode-ion distributions are given by



Here gPY(d/2) is the contact value of the uncharged hard wall–uncharged hard sphere distribution in the Percus–Yevick (PY) theory


with η = (π/6)ρd3 being the packing fraction of all ions. From Eqs (14) and (15) we have immediately



These equations imply that neither of gsum(d/2) and gdiff(d/2)/b, has any b dependence to linear order. It is further noted that Eq. (18) suggests that the MSA satisfies the HB condition (Eq. (5)) with a = 1, although from Eq. (17), gsum(d/2,b = 0) ≠1.

The disadvantage of the MSA is that it is a linear theory valid only in the neighbourhood of zero surface charge.

Monte Carlo Simulations

The MC simulations were done in the canonical ensemble following similar procedure as in earlier cases1114. The main features were (i) the use of parallel charged sheets method, pioneered by Torrie and Valleau30, to incorporate the long-range nature of the Coulomb interactions, and (ii) the adjustments in the length of the central rectangular parallelopiped MC cell perpendicular to the electrode until the target bulk concentration was reached. The number of configurations sampled were typically of the order 108 out of which the first 107 were used to equilibrate the system. The statistical uncertainty in reproducing the bulk target concentration was about ±5%.


All simulations and the theoretical calculations were done at the fixed ionic diameter d = 4.25 × 10−10 m and fixed relative permittivity εr = 78.5, which corresponds to a water-like solvent. The MPB and the GCS+EVT equations were solved numerically using a quasi-linearization iterative technique31. It is convenient to discuss results in terms of universal reduced parameters, viz., the reduced density ρ* = ρd3 and the reduced temperature T* = 4πε0εrdkBT/(z2e2). Note that the reduced temperature is the reciprocal of the commonly used plasma coupling parameter. While the value of z was 1 or 2 corresponding to 1:1 or 2:2 valency electrolytes, the ranges of ρ* and T* encompassed 0.00463 < ρ* < 0.1 (corresponding to 0.05 < c < 1.08, c being the electrolyte concentration in mol/dm3), and 0.15 < T* < 0.60. We note that although at the given d and εr, T* = 0.15 corresponds to T ~ 75 K, a low temperature for a 1:1 electrolyte, T* = 0.15 corresponds to T ~ 300 K, room temperature for a 2:2 electrolyte.

The results of this work are shown in Figs 110 for different physical states of our model system. In each figure, the quantities gdiff(d/2,b)/b, Gint/b, gco(d/2,b) and gctr(d/2,b) are plotted in the panels a, b, c and d, respectively. Gint is defined from Eq. (8) as Gint=(1/2)d/2Gdx. We have not shown the gsum(d/2,b) because (i) at large b it has the well-known parabolic dependence on b and any decent theory, including the GCS theory, shows this; (ii) at small b, gco(d/2,b) and gctr(d/2,b) are of comparable magnitudes so that a good idea of gsum(d/2,b) can be had simply by examining gco and gctr; (iii) some results presented elsewhere11,32. The MSA results have not been plotted. Our calculations have shown that for the range of physical parameters treated the MSA gco(d/2,b) < 0 beyond b ~ 1, leading to non-physical values of all the contact quantities. Within b < ~1, the MSA results are either linear (gco(d/2,b) and gctr(d/2,b)), or constants (gdiff(d/2,b)/b and Gint), and are therefore featureless and of limited usefulness. For the purposes of this discussion we will denote the left-hand side of Eq. (2) by HBC+MPB or HBC+(GCS+EVT) gdiff(d/2,b) when these theories are used to evaluate the right-hand side integral, and simply MPB or GCS+EVT gdiff(d/2,b), respectively, when computed directly from the solution of these equations.

FIG. 1
The gdiff(d/2,b)/b (a), Gint/b (b), gco(d/2) (c), and gctr(d/2) (d) as functions of b at reduced density ρ* = 0.00925 and reduced temperature T* = 0.60 (1:1 valency, c = 0.1 mol/dm3, T = 298 K). MC data, symbols; MPB results, solid line; GCS+EVT ...
FIG. 10
As in Fig. 1, but at reduced density ρ* = 0.0925 and reduced temperature T* = 0.40 (1:1 valency, c = 1 mol/dm3, T = 200.36 K)

The inexactness of MPB a leads to very small deviations seen in MPB gdiff(d/2,b) at and around b = 0 relative to the simulations in panels a of the figures. In most cases the deviations are imperceptible increasing slightly with increasing electrolyte density and/or decreasing temperature. The approximate value of MPB a can be traced to an approximate treatment of the fluctuation potential in the MPB formulation23. This in turn implies that although the HBC contact formula, Eq. (2), is exact, the same equation with the MPB gsum(x) and dψ (x)/dx in the integrand (on the right-hand side) is not. Hence, the differences between the MPB gdiff(d/2,b) and HBC+MPB gdiff(d/2,b) seen in panels a, with the latter being generally smaller than the former. Again these differences are small at low electrolyte densities around room temperature and increasing at high densities and/or low temperatures (Figs 24, ,6,6, ,8810). An alternative way of analyzing these differences is through the quantity Gint, which is shown in panels b of the figures. The MC values shown here were evaluated using the earlier simulation data8,11,12 and fresh simulation data in the integral of Eq. (7). Not unexpectedly, the MC data tend to cluster around the baseline Gint = 0, with the MPB |Gint| remaining small overall. The simulation results reinforce what Henderson and Bhuiyan8 observed earlier, that within the statistical and numerical constraints, these results are consistent with the HBC contact condition.

FIG. 2
As in Fig. 1, but at reduced density ρ* = 0.0925 and reduced temperature T* = 0.60 (1:1 valency, c = 1 mol/dm3, T = 298 K)
FIG. 4
As in Fig. 1, but at reduced density ρ* = 0.0462 and reduced temperature T* = 0.15 (2:2 valency, c = 0.5 mol/dm3, T = 298 K)
FIG. 6
As in Fig. 1, but at reduced density ρ* = 0.10 and reduced temperature T* = 0.15 (1:1 valency, c = 1.08 mol/dm3, T = 75.14 K)
FIG. 8
As in Fig. 1, but at reduced density ρ* = 0.00925 and reduced temperature T* = 0.20 (1:1 valency, c = 0.1 mol/dm3, T = 100.18 K)

A feature of the present results is the consistency between the GCS+EVT gdiff(d/2,b) and the HBC+(GCS+EVT) gdiff(d/2,b) and the good agreement with the simulations. This arises because of the closeness of the GCS+EVT theory, to that of GCS, as GCS satisfies the HBC condition exactly. For the densities and b considered in the figures, the largest value of ξ(d/2) is approximately 1.2 (Fig. 6, ρ* = 0.10), while ξ(x) rapidly tends to unity away from the electrode. Taken at face value this agreement with HBC and simulation may be misleading. The contact values gco(d/2,b) and gctr(d/2,b) at small b are poor, this being clearly visible in panels c for gco(d/2,b). Because of the scale used in panels d such differences get masked. Neglecting the fluctuation potential means that just having the exclusion volume term leads to an overestimation of the required correction to the GCS theory. As b increases, the volume term becomes increasingly important, and at large b the GCS+EVT gctr(d/2,b) can be in better agreement with simulation than the MPB gctr(d/2,b). The features shown in panels c and d are in contrast to the results of panels a and b, and clearly suggest some cancellation of errors. Indeed, similar error cancellations in applications of the mean field Poisson–Boltzmann theory in polyelectrolyte literature have long been known. We refer the reader to ref.33 for a recent account.

The MPB gco(d/2,b) and gctr(d/2,b) are seen to display reasonable consistency with the simulations in panels c and d of the figures. In the earlier reported analysis of the HBC condition using MC data, one conclusion that Henderson and Bhuiyan8 reached is the fact that it is important that both gdiff(d/2,b) and gco(d/2,b) be examined simultaneously. This is because at large b the gdiff(d/2,b) and gctr(d/2,b) are similar, both being parabolic in b, and can be seen in the figures. Thus, gco(d/2,b) = gsum(d/2,b) − gdiff(d/2,b), being the difference of two large numbers, provides a check of the MC results as it must tend to zero at large b. Both the GCS+EVT and MPB theories satisfy this property.


The emphasis in this paper has been on a systematic scrutiny of the contact characteristics of the MPB singlet distributions in a RPM planar double layer with special reference to the HBC contact value theorem. Earlier work has shown the MPB theory does not satisfy the HBL contact formula exactly11,18,32, in that the MPB a is not being fully consistent with the bulk osmotic coefficient, which is due to an approximate treatment of the fluctuation potential in the theory. This may well explain why the theory does not satisfy the empirical HB condition. These features notwithstanding, past13 and present calculations clearly indicate that the deviations in the MPB predictions from the HBL and HB relations are minimal for symmetric valency electrolytes in the electrolyte solution regime at room temperature even at a respectable 1 mol/dm3 concentration (cf. Fig. 1 of ref.13). The deviations are still small for asymmteric 2:1/1:2 salts at room temperature up to about ~1 mol/dm3 (cf. Figs 13 of ref.14).

FIG. 3
As in Fig. 1, but at reduced density ρ* = 0.00462 and reduced temperature T* = 0.15 (2:2 valency, c = 0.05 mol/dm3, T = 298 K)

With regard to the HBC contact theorem, the non-local nature of the relation makes its use seemingly less straightforward. Being an approximate theory the MPB is not expected to satisfy the relation exactly. However, the trends seen here are consistent with that seen earlier with the deviations increasing with concentration and/or lowering of the temperature. The general picture that emerges when one takes into account the gco(d/2) and gctr(d/2) shows the MPB contact behaviour to be very reasonable relative to the simulation data.

The good consistency shown by the GCS+EVT theory with the HBC and HB expressions is due to its close affinity to the GCS theory, which satisfies these conditions. At low b the GCS+EVT contact values of the singlet distribution functions are in relatively poor agreement with simulation, but improve at larger b. Also the GCS+EVT value of a = gBBGY(d/2) differs from the GCS value of unity. Thus the GCS+EVT agreement with these contact conditions indicates some cancellation of errors. The MSA, on the other hand, is limited to small values of b beyond which non-physical results occur for the contact quantities.

Detailed contact analysis of other statistical mechanical theories of the double layer is lacking. The present work suggests it would be of interest to see, for example, how well the hypernetted chain and density functional theories satisfy the HBC and HB contact conditions.

FIG. 5
As in Fig. 1, but at reduced density ρ* = 0.02 and reduced temperature T* = 0.15 (1:1 valency, c = 0.216 mol/dm3, T = 75.14 K)
FIG. 7
As in Fig. 1, but at reduced density ρ* = 0.00925 and reduced temperature T* = 0.20 (1:1 valency, c = 0.1 mol/dm3, T = 100.18 K)
FIG. 9
As in Fig. 1, but at reduced density ρ* = 0.0925 and reduced temperature T* = 0.30 (1:1 valency, c = 1 mol/dm3, T = 150.27 K)


L. B. B. wishes to acknowledge an institutional grant through Fondos Institucionales Para la Investigación (FIPI), University of Puerto Rico, while W. S. A. thanks a studentship from FIPI. D. H. acknowledges support from the National Institutes of Health through Grant No. 3R01GM076013-04S1.


Dedicated to Professor Ivo Nezbeda on the occasion of his 65th birthday.


1. Henderson D, Boda D. Phys. Chem. Chem. Phys. 2009;11:3822. [PubMed]
2. Henderson D, Blum L. J. Chem. Phys. 1978;69:5441.
3. Henderson D, Blum L, Lebowitz J. J. Electroanal. Chem. 1979;102:315.
4. Holovko M, Badiali J, di Caprio D. J. Chem. Phys. 2005;123:234705. [PubMed]
5. Holovko M, Badiali J, di Caprio D. J. Chem. Phys. 2007;127:014106. [PubMed]
6. Holovko M, di Caprio D. J. Chem. Phys. 2008;128:174702. [PubMed]
7. McQuarrie DA. Statistical Mechanics. Sausalito: University Science Books; 2000.
8. Henderson D, Bhuiyan LB. Collect. Czech. Chem. Commun. 2008;73:558.
9. Henderson D, Boda D. J. Electroanal. Chem. 2005;582:16.
10. Boda D, Henderson D, Chan KY, Wasan DT. Chem. Phys. Lett. 1999;308:473.
11. Bhuiyan LB, Outhwaite CW, Henderson D. J. Electroanal. Chem. 2007;607:54.
12. Henderson D, Bhuiyan LB. Mol. Simulation. 2007;33:953.
13. Bhuiyan LB, Outhwaite CW, Henderson D, Alawneh M. Bangladesh J. Phys. 2007;4:93.
14. Bhuiyan LB, Outhwaite CW, Henderson D. Mol. Phys. 2009;107:343.
15. Gouy M. J. Phys. 1910;9:457.
16. Chapman DL. Philos. Mag. 1913;25:475.
17. Stern O. Z. Elektrochem. 1924;30:508.
18. Bratko D, Bhuiyan LB, Outhwaite CW. J. Phys. Chem. 1986;90:6248.
19. Henderson D, Bhuiyan LB. J. Chem. Theory Comput. 2009;5:1985. [PubMed]
20. Bhuiyan LB, Outhwaite CW. Phys. Chem. Chem. Phys. 2004;6:3467.
21. Bhuiyan LB, Outhwaite CW. Condens. Matter Phys. 2005;8:287.
22. Patra CN, Bhuiyan LB. Condens. Matter Phys. 2005;8:425.
23. Outhwaite CW, Bhuiyan LB. J. Chem. Soc., Faraday Trans. 2. 1983;79:707.
24. Koryshnev AA. J. Phys. Chem. B. 2007;111:5545. [PubMed]
25. Outhwaite CW, Bhuiyan LB. J. Chem. Soc., Faraday Trans. 2. 1982;78:775.
26. Outhwaite CW, Lamperski S. Condens. Matter Phys. 2001;4:739.
27. Bhuiyan LB, Outhwaite CW. J. Colloid Interface Sci. 2009;331:543. [PubMed]
28. Lamperski S, Outhwaite CW, Bhuiyan LB. J. Phys. Chem. B. 2009;113:8925. [PubMed]
29. Blum L. J. Phys. Chem. 1977;81:136.
30. Torrie GM, Valleau JP. J. Chem. Phys. 1980;73:5807.
31. Bellman R, Kalaba R. Quasilinearization and Nonlinear Boundary Value Problems. New York: Elsevier; 1965.
32. Outhwaite CW, Bhuiyan LB. J. Chem. Phys. 1986;85:4206.
33. Piñero J, Bhuiyan LB, Reščič J, Vlachy V. Acta Chim. Slov. 2006;53:316.