HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 43(2) pp. 55-66 (2015) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2015-0010 SOME ANALYTIC EXPRESSIONS FOR THE CAPACITANCE AND PROFILES OF THE ELECTRIC DOUBLE LAYER FORMED BY IONS NEAR AN ELECTRODE DOUGLAS HENDERSON ∗ Department of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602, USA The electric double layer, which is of practical importance, is described. Two theories that yield analytic results, the venerable Poisson-Boltzmann or Gouy-Chapman-Stern theory and the more recent mean spherical approximation, are discussed. The Gouy-Chapman-Stern theory fails to account for the size of the ions nor for correlations amoung the ions. The mean spherical approximation overcomes, to some extent, these deficiencies but is applicable only for small electrode charge. A hybrid description that overcomes some of these problems is presented. While not perfect, it gives results for the differential capacitance that are typical of those of an ionic liquid. In particular, the differential capacitance can pass from having a double hump at low concentrations to a single hump at high concentrations. Keywords: Electric double layer, capacitance, Gouy-Chapman-Stern theory, mean spherical approx- imation, density functional theory, computer simulation 1. Introduction A double layer (DL) or an electric double layer (EDL) is formed when charged particles are attracted to a charged surface. The most obvious case is an electrolyte near a charged electrode (as in a battery). However, DNA can play a role that is analogous to the electrode. Ions can be attracted to membranes. A membrane can be thought of as a pseudo electrode. Ions are absorbed (often selec- tively) into physiological channels in membranes. Such channels permit the transport of nutrients into the cell and the removal of waste from the cell and are essential to the functioning of cells and life. The reader’s attention is drawn to some recent reviews of EDLs [1–3]. It is the case of an electrolyte near a charged flat sur- face that is considered here. This is the simplest case; it is an interesting and important application of statistical mechanical theory. The theory of the DL is important to our understanding of batteries. It can be used in the analy- sis of experimental electrochemical data and in analytical chemistry. In the model DL that is presented here, the electrode is approximated as a smooth flat charged surface located at x = 0. This surface is impenetrable and the ions are confined to the region x > 0. The charge of the elec- trode is located on the surface. There is no charge inside the electrode (x < 0). The electric field does not pene- trate the surface. The electrode is a classical metal. Ob- viously, this is an approximation but there has been very ∗Author for correspondence: doug@chem.byu.edu little work that takes into account the electronic structure of the electrode. The electrode charge is presumed to be uniform; the charge density of the electrode is σ and has the units of C/m2. Ions in the electrolyte near an elec- trode that have a charge opposite to that of the electrode are attracted to the electrode and form a layer whose net charge is equal in magnitude, but opposite in electric sign, to the charge of the electrode. The electrode and the at- tracted charge are together called an EDL. The charge in the EDL of the electrolyte can be spread over an extended region, usually called the diffuse layer, and need not con- sist solely of counterions whose charge is opposite to the electrode charge. The counterions can bring some coions with them. There may be regions of alternating charge where the coions predominate. However, the net charge of the attracted charged region in the electrolyte is equal in magnitude but opposite in sign to that of the electrode. Otherwise, the electric field would not vanish far from the electrode. For simplicity, the model electrolyte that is employed here is a fluid of charged hard spheres of diameter d. In this study, the electrolyte is assumed to be binary. For ad- ditional simplicity, the ions are assumed in this article to be symmetric both in the magnitude of their charge and diameter. The value of the charge of an ion of species i is zie, where zi is the ion valence and has the sign of the ion charge. The magnitude of the elementary charge is e. Because the ions are symmetric, |zi| = z. In the bulk, the density of the ions of species i is ρi = Ni/v, where Ni is the number of ions of species i in the bulk 56 HENDERSON and v is the volume of the system. Electrical neutrality requires that N1 = N2 or ρ1 = ρ2 or ∑ ziρi = 0. The solvent (usually water) of the electrolyte is characterized by a dielectric constant, �. Any change of the dielectric constant with a change of ion concentration is ignored. This model electrolyte is appropriately called the primi- tive model (PM). In the particular case considered here, where the ions all have the same diameter, this model is called the restricted primitive model (RPM). In this model, the interaction between a pair of ions, whose centers are separated by the distance r, is given by uij(r) = ∞ for r < d zizje 2 4π�0�r for r ≥ d , (1) where �0 is the permittivity of free space, and the interac- tion of an ion with the surface is given by uwi(x) = { ∞ for x < d/2 − σziex �0� for x ≥ d/2 , (2) where x is the distance between the center of the ion and the surface. Our task is to determine the density profile, ρi(x), of the ions, or equivalently, gi(x) = ρi(x)/ρi. Note that ρi(∞) = ρi, so that gi(∞) = 1. Once, the gi(x) are known, the charge profile (C/m2), for x > d/2, is given by q(x) = e ∑ i ziρihi(x), (3) where hi(x) = gi(x) − 1. In writing Eq. 3, the global charge neutrality condition ∑ ziρi = 0 has been invoked. The charge density on the electrode is given by σ = −e ∑ i ziρi ∫ ∞ d/2 hi(t)dt. (4) There is no point including the region 0 < t < d/2 in the integral since ∑ hi(t) = 0 in this region. The potential profile (in Volts) is given by φ(x) = − e ��0 ∑ i ziρi ∫ ∞ x (t−x)hi(t)dt. (5) In particular, the potential (Volts) of the electrode is given by V = φ(0) = − e ��0 ∑ i ziρi ∫ ∞ 0 thi(t)dt. (6) Note that these equations satisfy Poisson’s equation d2φ(x) dx2 = − q(x) ��0 . (7) Indeed, Eqs. 3 and 5 are obtained by integrating Poisson’s equation. An alternative procedure for computing the po- tential profile has been proposed by Boda and Gillespie [4] for simulation purposes. It is often convenient to use dimensionless, or re- duced, values that are denoted by an asterisk. For a sys- tem whose temperature (K) is T , the reduced temper- ature is T∗ = 4π��0dkT/z2e2. The reduced density is ρ∗i = ρid 3, the reduced electrode charge density is σ∗ = σd2/e, and the reduced potential is φ∗ = βeφ, where β = 1/kT , with k being the Boltzmann constant (the gas constant per particle). 2. Poisson–Boltzmann or Gouy–Chapman–Stern theory: comparison with simulations The classic theory of the EDL was developed by Gouy [5], Chapman [6], and Stern [7] (GCS) a century ago. The theory is based on Poisson’s equation together with the Boltzmann formula, gi(x) = { 0 x < d/2 exp[−βzieφ(x)] x ≥ d/2 . (8) In electrostatics, Poisson’s equation is exact and is equiv- alent to one of Maxwell’s equations. The Boltzmann for- mula is approximate and neglects ion size and correla- tions between the ions. Eq. 8 states that gi(x) for the coions is the reciprocal of gi(x) for the counterions. This is not true, in general [8]. Equation 8, when inserted into Poisson’s equation, yields what may be called the Poisson-Boltzmann (PB) or GCS approximation. This approximation is also em- ployed in the Debye-Hückel (DH) theory for bulk elec- trolytes that was developed some years later. However, because of the three dimensional geometry of the DH the- ory, the nonlinear PB equation cannot be solved analyti- cally and the PB equations in the DH theory are usually linearized. In the GCS theory, the resultant PB equation is a nonlinear second order differential equation. As has been pointed out, such equations generally do not yield analytic solutions. However, for the one dimensional ge- ometry of the planar DL that is considered here, an ana- lytic solution is possible in the case of the GCS theory. The resulting PB/GCS potential is βzeφ(x) 2 = ln { 1 + b/2 1 + √ 1 + b2/4 exp[−κy] } − ln { 1 − b/2 1 + √ 1 + b2/4 exp[−κy] } , (9) where y = x−d/2 > 0 and b = βzeσ ��0κ , (10) where κ is the Debye screening parameter that is given by κ = √ βz2e2ρ ��0 (11) Hungarian Journal of Industry and Chemistry DOUBLE LAYERS NEAR AN ELECTRODE 57 with ρ = ∑ ρi. The parameter b is another dimension- less measure of the electrode charge density. However, it is not as fundamental a quantity as σ∗ since it arises from a theory. The parameter κ is a screening parameter; it is an inverse measure of the distance over which the pro- files reach their asymptotic values within the GCS and DH theories. In the GCS theory the relationship between the poten- tial difference and electrode charge density is given by sinh [ βzeφd/2 2 ] = b 2 , (12) where φd/2 = φ(d/2) is often called the diffuse layer potential. Some relations that are equivalent to Eq. 12 are cosh [ βzeφd/2 2 ] = √ 1 + b2/4, (13) tanh [ βzeφd/2 2 ] = b/2√ 1 + b2/4 , (14) and tanh [ βzeφd/2 4 ] = b/2 1 + √ 1 + b2/4 . (15) The equivalence of Eqs. 12–15 is a result of identities among the hyperbolic functions. Thus, Eq. 9 can be written as βzeφ(x) 2 = ln { 1 + tanh [ βzeφd/2 4 ] exp(−κy) } − ln { 1 − tanh [ βzeφd/2 4 ] exp(−κy) } . (16) Alternative forms of Eqs. 9 and 16 are tanh [ βzeφ(x) 4 ] = tanh [ βzeφd/2 4 ] exp(−κy) (17) or tanh [ βzeφ(x) 4 ] = b/2 1 + √ 1 + b2/4 exp(−κy). (18) In the GCS theory, the potential difference across the EDL is V = − ze ��0 ∑ i ρi ∫ ∞ 0 thi(t)dt = σd 2��0 + φd/2, (19) where φd/2 is given by Eq. 12. Thus, the capacitance, C = σ/V , of the EDL is 1 C = d 2��0 + 2 sinh−1(b/2) ��0κb (20) and the differential capacitance, Cd = ∂σ/∂V , of the EDL is given by 1 Cd = d 2��0 + 1 ��0κ √ 1 + b2/4 . (21) Equations 20 and 21 are formally identical to a diffuse layer capacitor with capacitance, Cdl = ��0κ b/2 sinh−1(b/2) , (22) or differential capacitance Cdld = ��0κ √ 1 + b2/4 (23) in series with an inner–layer parallel plate capacitor with capacitance (or differential capacitance), Cil = Cild = 2��0 d . (24) At contact, gi(d/2) = exp[βzieφd/2] = 1 + b2 2 − zi z b √ 1 + b2 4 , (25) so that gsum(d/2) = 1 2 [g1(d/2) +g2(d/2)] is, in the GCS theory, given by gsum(d/2) = 1 + b2 2 . (26) This is to be compared with the exact result (for the re- stricted PM) due to Henderson and Blum [9] and Hender- son, Blum, and Lebowitz [10], gsum(d/2) = p ρkT + b2 2 , (27) where p is the osmotic pressure of the electrolyte. The second term in the above equation is just the Maxwell electrostatic stress. Thus, Eq. 27 is just a force balance condition where the momentum transfer to the electrode is equal to the sum of the osmotic term and the Maxwell stress. The GCS theory deals with the electrostatic term correctly but replaces the osmotic pressure with the ideal gas result p = ρkT because of the neglect of the ion diameters. For comparison with the mean spherical approxima- tion (MSA), which is a linear response theory that will be considered in the next section, it is worthwhile to give the linearized GCS theory results, obtained for the case of small electrode charge. In this case, gi(x) = { 0 for x < d/2 1 −βzieφ(x) for x ≥ d/2 , (28) or gi(x) = 0 for x < d/2 1 + βzieσ ��0κ exp (−κy) for x ≥ d/2 . (29) The potential profile in the diffuse layer is given by φ(x) = φd/2 exp(−κy) (30) with φd/2 given by φd/2 = σ ��0κ . (31) 43(2) pp. 55-66 (2015) DOI: 10.1515/hjic-2015-0010 58 HENDERSON 0.00 0.01 0.02 0.03 0.04 (C d dl ) -1 [cm 2 /µF] 0.02 0.03 0.04 0.05 0.06 0.07 C d -1 [ c m 2 /µ F ] -8 -4 -2 0 2 4 8 Figure 1: Experimental values of the inverse differential capacitance, Cd, of an aqueous solution of NaH2PO4 at 25 ◦C near a hanging drop mercury electrode as a function of the inverse diffuse layer capacitance, Cdld , obtained from Eq. 23. The points are the experimental results of Parsons and Zobel. The light straight lines of unit slope give the results of the GCS theory but with the experimental in- ner layer capacitance obtained empiricially. The numbers at the low concentration end of the lines give the elec- trode charges in units of µCcm−2. The heavy solid curve gives the results of the MSA using a dipolar hard sphere model for the solvent together with an estimate of the con- tribution of the electronic structure of the electrode and is intended only as an aid to the eye. This figure has been reproduced, with permission, from Ref. [1]. The potential difference across the EDL is V = σd 2��0 + σ ��0κ (32) and the capacitance (and differential capacitance) of the EDL is given by 1 C = d 2��0 + 1 ��0κ . (33) In the GCS theory, in the limit of large κ (high concentra- tion) or large b (high electrode charge), the diffuse layer capacitance is large and, as a result, the inner layer capac- itance dominates, due to the reciprocal or series additivity of Eqs. 20 and 21. Hence, in the GCS theory, C = 2��0/d is the limiting (maximum) value of C or Cd at large con- centrations or large electrode charges. Further, the differ- ential capacitance at low concentrations looks something like a parabola but flattens out at large σ. At high concen- trations, the differential capacitance is constant. Any ad- ditional shape in the experimental differential capacitance is added by an empirical fit of the inner layer capacitance to the experimental results. However, the diffuse layer ca- pacitance is presumed to be given adequately by the GCS theory. Parsons and Zobel (PZ) [11] have plotted their exper- imental results for the inverse of the differential capaci- tance as a function of the inverse of the diffuse layer dif- ferential capacitance, given by Eq. 23. Such a plot is of- ten called a Parsons-Zobel plot. If the GCS theory were correct, this should result in a straight line. The extrapola- tion of the straight line to 1/Cdld = 0 (high concentration and/or high electrode charge density) should, if the GCS theory were correct, yield the reciprocal of the inner layer capacitance. As is seen in Fig. 1, at first sight the exper- imental results of PZ (the points) do seem to follow a straight line and, conventionally, are presumed to provide an experimental verification of the GCS theory. In Fig. 1, the light straight lines are the GCS results. The solid curve is the result of the MSA that has not yet been dis- cussed. For the moment the solid line can be considered to be an aid to the eye in following the trend of the experi- mental results. In the conventional GCS picture, the inner layer capacitance might not be given by Eq. 24 but might differ because of the presumed effect of the presence and nature of the solvent molecules and the electronic struc- ture of the electrode that are beyond the GCS theory. Pos- sible solvent effects might be a lower dielectric constant due to the alignment of the solvent molecules because of the strength of the electrode charge. The important point is that, in the GCS theory, such solvent effects are pre- sumed to be confined only to the inner layer. The GCS theory is conventionally considered to provide an ade- quate description of the diffuse layer where the ions are present. Also, it is thought to provide a description of the EDL when combined with some treatment of the solvent molecules in the inner layer, or even an empirical fit. In- deed, the extrapolation of the straight lines to 1/Cdld = 0 is one method of obtaining presumed “experimental” val- ues of Cild . However, a careful examination of the PZ experimen- tal results in Fig. 1 indicates that the experimental dif- ferential capacitance does not follow Eq. 21 at high con- centrations (the left side of the figure) but rises above the extrapolated intercept, possibly without limit. Until re- cently, most experimentalists have ignored this point and did not concern themselves with this issue because ions are not soluble in water when their concentration is large and it is difficult to obtain results with other solvents. Additionally, experimental results are difficult to obtain at high electrode charges for conventional electrolytes. However, as we shall see, DLs in ionic liquids can be formed at high concentrations and the deficiencies of the GCS theory become quite apparent. Given that experiments on aqueous systems are dif- ficult in regimes where problems with the GCS theory become apparent, it is useful to consider computer sim- ulations. One simulation technique is the Monte Carlo (MC) method. Until recently, it has been the most com- mon simulation tool in DL studies. In MC simulations the ions undergo a random walk and the profiles and other properties of interest are obtained by averages over this random walk. A simple random walk would take forever before useful results could be obtained. However, mean- ingful results can be obtained by means of a biased ran- dom walk that confines the ions to regions in which they Hungarian Journal of Industry and Chemistry DOUBLE LAYERS NEAR AN ELECTRODE 59 0 5 10 15 b 0 2 4 6 βe φ( d/ 2) a 0.0 0.1 0.2 σd 2 /e 0 2 4 6 βe φ( d/ 2) b Figure 2: Diffuse layer potential of a 1:1 electrolyte (d = 3 Å) at room temperature as a function of b (part a) and σ (part b). The curves are, from top to bottom, for 0.01, 0.1, and 1 M solutions. The symbols give the simulation results. The solid curve in part a gives the GCS results. The dotted lines connect the MC results for easier visual- ization. The lines in part b give the GCS results. Part a is reproduced, with permission, from Ref. [1]. have a high probability of residing. The simulation cell consists of a parallelopiped with a charged wall (the elec- trode) at x = 0 and another wall (charged or uncharged) at x = L, where L is so large that the two walls do not interfere. Periodic boundary conditions are used in the other two directions. The size of the cell is chosen to be large enough that electrostatic screening eliminates the effects of the periodic image cells. The number of ions of each species is chosen so that the system is electroneu- tral. The first use of MC simulations for the study of the EDL was that of Torrie and Valleau [12, 13]. After their seminal studies, there was a hiatus in simulation studies of the EDL. However in recent years, there has been a renewed interest in simulations of the EDL that includes the work of Bhuiyan et al. [8], Boda et al. [14, 15], and Lamperski et al. [16, 17]. Another simulation technique is the molecular dy- namics (MD) method in which the equations of motion are solved and the properties of the system of interest are obtained by averaging over the positions and velocities of the ions. A simulation cell that is similar to that used 0 1 2 3 4 5 g i (x ) a 0 1 2 3 4 5 x/d 0 1 2 3 4 5 g i (x ) b Figure 3: Normalized density profiles, gi(x), of a 1:1 elec- trolyte (d = 3 Å) at 1 M and room temperature for the state for which the MC values are σd2/e = 0.1685 and βeφ(d/2) = 2.6. The points give the simulation results and the curves give the GCS results. The comparison is made at the same charge density (part a) and the same dif- fuse layer potential (part b). in MC simulations is employed. In recent years, there has been an interest in MD simulations of the EDL, especially for ionic liquids. Some representative studies are those of Vatamanu et al. [18, 19], Hu et al. [20], and Feng [21] et al.. A comparison with simulation gives an unambigu- ous test of the GCS theory since uncertainties resulting from empirical fits of the diffuse layer capacitance can- not arise. The GCS theory and simulations both use the same model and interaction parameters that are defined in Eqs. 1 and 2. Additionally, simulations and theory give results for the density profiles, gi(x), that cannot be obtained by present experimental methods. The simula- tions plotted in Figs. 2–5 are those of Boda et al. [22]. In Fig. 2a, the electrostatic potential βeφ(d/2) for a 1:1 electrolyte is plotted as a function of b for three con- centrations (0.01M, 0.1M, and 1M). If the GCS theory were correct, these curves would be identical and inde- pendent of concentration. Hence, there can be only one GCS curve in Fig. 2a. As is seen, φ(d/2) as a function of b actually decreases with increasing concentration. In Fig. 2b, φ(d/2) is plotted as a function of σd2/e. The CGS curves are greater than the simulation results, espe- 43(2) pp. 55-66 (2015) DOI: 10.1515/hjic-2015-0010 60 HENDERSON -0.2 -0.1 0.0 0.1 0.2 σd 2 /e -4 -2 0 2 4 6 βe φ( d/ 2) Figure 4: Diffuse layer potential of a 2:1 electrolyte (d = 3 Å) at room temperature as a function of σ. The curves are, from top to bottom, for 0.01, 0.1, and 1 M for positive σ and the reverse for negative σ. The symbols give the simulation results. The curves give the GCS results. cially as the concentration increases. The density profiles for a 1:1 electrolyte are plotted in Fig. 3. The compari- son is made at the same value of σ in part a and the same value of φ(d/2) in part b. In principle, there is no reason to choose whether the comparison should be made at the same σ, the same φ(d/2), or the same φ(0). It was natural for Torrie and Valleau to make their comparisons at the same σ because σ is the input variable in their method. However, φ is the natural variable in the GCS theory. In any case, the GCS theory looks best when σ is used as the input variable. Torrie and Valleau overstated things when they said that the GCS theory was reasonable for a 1:1 electrolyte. Their statement is most applicable if the comparison is made at the same value of σ. The value of the counterion profile would be in poor agreement at x = d/2 if φ(d/2) or φ(0) were used as the input vari- able. A similar comparison is made for a 2:1 electrolyte in Figs. 4 and 5. The agreement of the GCS theory with simulations is much poorer. The electrostatic interactions are stronger because of the presence of the divalent ions. When the divalent ions are the counterions, the potential, φ(d/2) has a maximum and then decreases with increas- ing electrode charge. This is not seen in the GCS results which are monotonic. Further, the simulation profiles are not monotonic whereas the GCS results are monotonic. The simulation profiles have oscillations. The EDL can consist of regions where counterions or coions predomi- nate. When the coions predominate, this phenomenon is known as charge inversion. Of course, the net charge in the diffuse layer is still equal in magnitude, but opposite in sign, to that of the electrode charge. This is required to screen the electrode charge and potential far from the electrode. Generally, experimentalists have been content to ig- nore the discrepancies in the results of the GCS theory and state that these differences are unimportant since they occur at high electrode charges or high concentrations 0 1 2 3 4 5 g i (x ) a 0 1 2 3 4 5 x/d 0 1 2 3 g i (x ) b Figure 5: Normalized density profiles, gi(x) of a 2:1 elec- trolyte (d = 3 Å) at 1 M and room temperature for the state for which the MC values are σd2/e = −0.1685 and βeφ(d/2) = −0.15. The points give the simulation re- sults and curves give the GCS results. The comparison is made at the same charge density (part a) and the same dif- fuse layer potential (part b). or for high valence electrolytes or nonaqueous systems, where experimental results are difficult to obtain. How- ever, this is short-sighted. As scientists, one of our goals is to understand what is happening. This cannot be done with an inaccurate theory even with curve fitting. In the remainder of this article, attention is directed to more ac- curate, but still analytic, theories. 3. Mean spherical approximation The mean spherical approximation (MSA) is a natural extension of the linearized GCS theory in which the size of the ions is taken into account. It was first ap- plied to the EDL by Blum [23]. The GCS is usually ob- tained by means of the solution of a differential equa- tion whereas the MSA is obtained from the solution of an integral equation. At first sight, the connection between the GCS and MSA theories is unclear. However, Hen- derson and Blum [24] demonstrated that the GCS the- ory could also be obtained from an integral equation. In fact, the linearized GCS integral equation is just the MSA integral equation with the effect of ion size ignored. Ac- tually, Henderson and Blum proved a more general re- sult. They showed that the GCS theory followed from the Hungarian Journal of Industry and Chemistry DOUBLE LAYERS NEAR AN ELECTRODE 61 hypernetted–chain approximation (HNCA) when ion size was neglected. The MSA can be regarded as a linearized version of the HNCA and, because of this, the stated rela- tion of the linearized GCS theory to the MSA follows. In this article, the HNCA is not considered because it does not yield analytic results and has severe problems when σ is large [25]. Also, the emphasis in this article is upon analytic, or at least explicit, results that can be valuable in practical calculations. The MSA result that is analogous to Eq. 29 is gi(x) = 0 for x < d/2 g0(x) − βzieσ ��0κ f(y) for x ≥ d/2 , (34) where g0(x) is the Percus-Yevick (PY) profile for hard spheres near a hard surface. Earlier, Henderson, Abra- ham, and Barker (HAB) [26] obtained an integral equa- tion for g0(x). The second term gives the electrostatic part of the profile for charged hard spheres near a charged hard surface. Blum [23] did not obtain a result for f(y) but he did obtain an analytic result for the Laplace trans- form of f(y),∫ ∞ 0 exp(−sy)f(y)dy = = s s2 + 2(Γσ)s + 2(Γσ)2(1 − exp[−s]) , (35) where 2Γ is a renormalized screening parameter that is related to κ by κ = 2Γ(1 + Γσ) or 2Γσ = √ 1 + 2κσ- 1. Note that for small κ (small concentrations), 2Γσ = κσ−(κσ)2/2+···. Thus, the MSA screening parameter is smaller than the GCS screening parameter. This suggests that the MSA EDL is wider than that of the GCS theory. This agrees with the simulation results. The notation of Blum has been followed. However, it might have been preferable if he had incorporated the factor of 2 into the definition of Γ so that Γ became κ at low concentrations. Note that at low concentrations, the Laplace transform of f(y) becomes 1/s(1 + κs). This means that f(y) = exp(−κy), (36) in the limit of low concentrations. Also, g0(x) = 1 in this limit. Thus, at low concentrations, the MSA becomes the GCS theory. Blum did not invert the Laplace transform of f(y). However, he did obtain the contact value of f(y) by examining the Laplace transform of f(y) at large s. He showed that f(0) = 1. Using the earlier result of HAB for g0(d/2), the contact value of gi(x) is gi(d/2) = 1 + 2η (1 −η)2 − βzieσ ��0 , (37) where η = πρd3/6. The MSA contact value is an im- provement over the GCS result that contains only the ideal gas term. However, the osmotic pressure should have both a hard sphere term and an electrostatic term. Additionally, the hard sphere term is accurate only for low values of ρ. Equation 37 does not contain the quadratic term b2 of Eq. 27. This is because the MSA is a linearized theory. A better expression for the osmotic term is p ρkT = 1 + η + η2 −η3 (1 −η)3 − Γ3 3πρ , (38) where Γ is the renormalized screening parameter that has been defined above. This result is obtained from the ap- plication of the MSA to bulk electrolytes. The MSA, as is the case for most theories, is not fully self-consistent. Henderson et al. [27] have compared this expression with their simulations (see their Fig. 1) and found it to be very accurate. Despite these problems, the MSA contact value given in Eq. 37 does represent an advance. By expansion of the Laplace transform of f(y), it is easy to show that the MSA EDL satisfies electroneutral- ity. That is, the charge in the EDL is equal in magnitude, but opposite in sign, to the electrode charge. Again, by expanding the Laplace transform, the MSA expressions for the total and diffuse layer potentials of the EDL are found to be V = σ ��0(2Γ) (39) and φd/2 = σ ��0κ [1 − (Γd)2]. (40) Thus, in the MSA, the capacitance (and differential ca- pacitance) of the EDL is 1 C = 1 ��0(2Γ) . (41) Expanding the expression that defines Γ, 2Γ = κ−κ2d/2 + κ3d2/2 + · · ·. (42) Therefore, 1 C = 1 ��0κ + d 2��0 − κd2 4��0 + · · · . (43) The MSA capacitance does not reach a maximum at 2��0/d but continues to increase, as is indicated in Fig. 1. The MSA (solid) curve in Fig. 1 was not calculated from Eq. 41 but from a more sophisticated version of the MSA, that is not discussed in detail here, which includes the contribution resulting from explicit solvent molecules and the electronic structure of the metal [28–30]. How- ever, the results of Eq. 41 are qualitatively similar to the more sophisticated results. In this paper, the solid curve serves to guide the eye. The inner layer capacitance con- tinues to be 2��0/d but it is simply the electrode charge divided by potential difference across the inner layer and not a ‘catch all’ for the deficiencies of the GCS theory. The diffuse capacitance is the electrode charge divided by the potential difference from the distance of closest approach to the bulk electrolyte and contains correction terms to the GCS theory. This is the reverse of the usual interpretation of the GCS theory where the GCS expres- sions are assumed to be accurate for the diffuse layer 43(2) pp. 55-66 (2015) DOI: 10.1515/hjic-2015-0010 62 HENDERSON 0.14 0.16 0.18 0.2 0.22 0.24 ρ* 1.4 1.6 1.8 2 2.2 2.4 C /F m -2 MC MSA GCS Figure 6: Double layer capacitance, C, as a function of the reduced density, ρ∗ = ρd3, at the reduced temperature T∗ = 0.08. The circles are the MC data of Henderson et al. [27] and the lines are the MSA and GCS results. The line through the circles is given as a guide to the eye. Because the MSA is a linearized theory, the MC, MSA, and GCS capacitances are for σ = 0. capacitance and everything else is ‘lumped’ into the in- ner layer capacitance. In the more sophisticated version MSA, the contributions due to the solvent molecules ap- pear in both the diffuse and inner layer potentials. The solvent molecule profile is as diffuse as that of the ions. The effect of the molecular nature of the solvent is not confined to the inner layer. The initial slope of the φ(d/2) vs. σ curves in Figs. 2 and 4 is just the inverse of the diffuse portion of the ca- pacitance. It is seen that the initial slope of the MC curves is well described by the GCS theory at low concentrations but increasingly falls below the GCS initial slope (the in- verse of the differential capacitance) with increasing con- centration as was seen in the experimental results in Fig. 1. Henderson et al. [27] compared the GCS and MSA differential capacitances with their simulation results for a broad range of densities and at a temperature that was meant to be qualitatively representative of an ionic liquid. As is seen in Figs. 6 and 7, the MSA results are consid- erably improved over the GCS results. The comparison with the MC results is made for a small value of σ be- cause the MSA is a linearized theory that is applicable only for small σ. Although an analytic expression for f(y) is not avail- able, Henderson and Smith [31] were able to obtain a zonal expansion for f(y). They showed that f(x) = ∞∑ n=1 fn(z)u(z), (44) where z = t−n + 1, t = x/d, u(z) is the Heaviside step function that is zero for z < 0 and one for z ≥ 0 and fn(z) = exp(−µ) µn (n− 1)! [jn−2(µ) − jn−1(µ)] (45) with µ = (Γd)z. The function jm(µ) is the spherical x / d 0.5 1.0 1.5 2.0 2.5 0 1 2 3 4 g (x /d ) 0 1 2 3 4 5 counterions co-ions counterions co-ions (a) (b) Figure 7: The electrode-ion normalized density profiles, gi(x/d), at the reduced density ρ∗ = ρd3 = 0.5 for the reduced temperature T∗ = 0.8 and surface charge den- sity, σ = 0.05 C/m2, is small enough that the MSA is applicable. The circles are the MC results of Henderson et al. [26] and the dashed line gives the MSA result. The line through the circles is given as a guide to the eye. This figure is reproduced, with permission, from Ref. [26]. Bessel function that is easily calculated using the recur- rence formula for this function. Hundreds of jm(µ) can be calculated without difficulty, even with a laptop com- puter. Henderson and Smith [31] also obtained a zonal ex- pansion for g0(x). Their result is g0(x) = ∞∑ n=1 gn0 (z)u(z), (46) where z is again given by z = t − n + 1 with t = x/d. The expressions for the gn0 (x) are rather complex. How- ever, Henderson and Smith gave results for n ≤ 5. The formulae for g0(x) and f(x) are not quite analytic since they involve infinite series. However, these results are ex- plicit and easily used. Fortran programs to obtain g0(x) and f(x) are given in Supplementary Material. Note that the program for g0(x) consists of two parts. One part cal- culates those parameters that depend only on the state of the electrolyte and the other subroutine in each code cal- culates profiles for a given x. The user should resist the temptation to combine the two parts into one. McQuarrie [32] did this in an appendix to his excellent book and pro- duced an inefficient, and probably incorrect, code that he referred to as ‘Henderson’s code’. Fortunately, his code Hungarian Journal of Industry and Chemistry DOUBLE LAYERS NEAR AN ELECTRODE 63 is illegible in the later printing of his book. If the reader does combine the codes, the reader is on his/her own and should not refer to the combined, or otherwise modified, code as ‘Henderson’s code’. The functions g0(x) and f(x) are oscillatory, in ac- cord with the simulations. They are improvements, qual- itative and quantitative, to the monotonic functions of the GCS theory. 4. A useful hybrid description As has been mentioned, the deficiencies of the GCS theory could, until recently, be dismissed as appearing mainly under conditions that are of limited experimental interest. However, there has been considerable recent in- terest in EDLs formed by ionic liquids. Ionic liquids can be thought of as room–temperature molten salts. Because there is no solvent, the ions do not become insoluble in some solvent and experimental results can be obtained at high concentrations. The fact that they exist at room tem- perature is a great experimental convenience. Kornyshev [33] has drawn attention to these electrochemical systems and aptly suggested that they provide a paradigm change in electrochemistry. He modestly ends the title of his im- portant paper with a question mark. An exclamation mark might have been more appropriate. As well as exposing the deficiencies of the GCS theory, EDLs in ionic liq- uids are important in green technologies, the design of novel energy storage devices, such as high-tech batter- ies and super-capacitors [34]. Ionic liquid DLs have at- tracted recent experimental [35, 36] and theoretical inter- est [18–21, 37–39]. The differential capacitance, as de- termined by MC simulations, of a simple model [38] of an ionic salt in which T∗ = 0.8 and d = 8 Å is given in Fig. 8 for ρ∗ = 0.04, 0.14 and 0.24. At low concen- trations, the differential capacitance is parabolic-like, as the GCS theory suggests. However, Cd does not become flat at large electrode charges. At higher concentrations, Cd at small electrode charges continues to increase with increasing concentration. This has been seen in Fig. 6. At high electrode charges, the capacitance decreases. The nature of this decrease seems to be independent of the concentration. This is similar to the GCS theory except that the capacitance is not flat at high electrode charges but decreases. The decrease is due to the fact that the ions are not point charges but occupy space and cannot sit on top of each other. The diffuse layer must become thicker and the capacitance decreases as the electrode charge in- creases. This is not because the distance of closest ap- proach of the ions increases. Strong secondary peaks in the counterion profile appear [37]. The beginnings of this trend were first observed by Torrie and Valleau [12] and seem to be quite universal. The GCS theory satisfies Eq. 27 at high electrode charges but fails at low electrode charges whereas the MSA gives reasonable results at small electrode charges. This implies that a repair of the GCS so that it gives the MSA results in the regime of the low electrode charges -2 -1 0 1 2 σ∗ 3 4 5 6 7 8 9 C df * 0.24 0.04 0.14 Figure 8: Differential capacitance, C∗df = Cdd/4π�0, ob- tained from MC simulation for the EDL of a model ionic liquid with d = 4 Å and T∗ = 0.8. The curves are, from bottom to top, for ρ∗ = 0.04, 0.14, and 0.24. This figure has been reproduced, with permission, from Ref. [40]. but leaves the high electrode charge part unchanged might be useful. Likely, there is no way to accomplish this in a fundamental way. Additionally, there are probably several semi-empirical ways in which this could be done. Henderson and Lamperski [40] have presented one pro- cedure. Because it is not based on any fundamental ideas, it is not a theory. It would be more appropriate to refer to their procedure as a description. They proposed that the differential capacitance for a symmetric salt could use- fully be written as 1 Cd = d′ 2��0 + d′ 2��0 √ 1 + b2/4 ( 1 Γd′ − 1 ) . (47) The parameter d′ is an adjustable parameter and repre- sents the effective thickness of the diffuse layer. At small b (small electrode charge), Eq. 47 yields 1 Cd = 1 2Γ��0 , (48) which is the MSA result. At large b (large electrode charge), Eq. 47 yields 1 Cd = d′ 2��0 . (49) The results of Eq. 47, using d′ = 2d for the system that Lamperski et al. simulated, were given by Henderson and Lamperski. Qualitatively, the results are very similar to the simulation results shown in Fig. 8. Better agreement could be obtained by making d′ increase with electrode charge. Figure 2 of Henderson and Lamperski suggests that Cd is proportional to 1/σ∗ at large σ∗ (electrode charge) with the proportionality constant being indepen- dent of concentration. This behavior was first predicted by Kornyshev [33] on the basis of a lattice theory and seems to be universal. As well as the simulations of Hen- derson and Lamperski, it has been seen experimentally 43(2) pp. 55-66 (2015) DOI: 10.1515/hjic-2015-0010 64 HENDERSON -1.5 -1 -0.5 0 0.5 1 1.5 σ* 4 6 8 10 C * d 0.24 0.14 0.04 Figure 9: Differential capacitance, C∗d = Cdd/4π�0, ob- tained from the hybrid description of the EDL of a model ionic liquid with d = 4 Å and T∗ = 0.8. The curves are, from bottom to top, for ρ∗ = 0.04, 0.14, and 0.24. The solid and broken curves give the results of the hybrid approach and GCS theory, respectively. by Islam et al. [35]. Something of the nature of d′ = d1 + d2|σ∗| (50) would give the desired decrease of Cd at large σ∗. Equa- tion 50 is sensible because it is consistent with the diffuse layer becoming thicker as the electrode charge increases. The results of this ansatz with d1 = 2d and d2 = d are given in Fig. 9. The results are similar to the simulation results in Fig. 8. This hybrid approach is capable of yield- ing a capacitance with a double hump at low concentra- tions and a single hump at high concentrations. This be- havior is predicted by simulations and all the good theo- ries of the DL of ionic liquids. A hybrid treatment of the profiles, gi(x), is possible. One could start with the MSA expression for gi(x) and add to this gGCSi (x; b) − g GCS i (x; b = 0). However, it must be realized that the MSA expressions for the pro- files are less accurate than the MSA expressions for the potential and capacitance. The potential and capacitance are integrals and tend to average out any inaccuracies in the profiles. 5. Conclusion The study of the electric DL is an important application of statistical mechanics that is of experimental and ap- plied interest. The GCS theory is popular with exper- imentalists because it is intuitively simple and easy to use in the routine analysis of experiments. However, the GCS theory has deficiencies. Its use leads to the idea that any problems with the GCS theory can be ‘swept under the carpet’ by placing all of these problems into an em- pirical treatment of the inner layer. In reality, the defi- ciencies of the GCS theory lie with the GCS treatment of the diffuse layer. Admittedly, it is difficult to observe this in aqueous systems. However, it is not impossible. The departure from linearity in the Parsons-Zobel plot (Fig. 1) is real and should not be ignored. The important field of ionic liquid electrochemistry requires something more adequate than the GCS theory. The best theories of the EDL are the modified Poisson-Boltzmann theory [41] and the density functional theory [42]. However, both the- ories are numerical and require an iterative numerical so- lution of a fairly large set of equations and may not be appealing in an experimental analysis. A hybrid descrip- tion, such as that explored here, preserves the advantage of an analytic treatment of the capacitance and is no more cumbersome than the GCS theory. Supplementary Information Fortran programs to calculate g0(x) and f(x−d/2) using MSA. The codes can be downloaded free of charge from http://tinyurl.com/hjic-2015-0010-suppl. Acknowledgement Professor Dezső Boda assisted with the preparation of this paper. Professor Lutful Bari Bhuiyan read the manuscript prior to submission and suggested several im- portant modifications. The author is grateful to both col- leagues for their continuing wise advice. REFERENCES [1] Henderson, D., Boda, D.: Insights from theory and simulation on the electrical double layer, Phys. Chem. Chem. Phys., 2009 11(20), 3822–3830 10.1039/b815946g [2] Cherstvy, A.G.: Electrostatic interactions in bio- logical DNA-related systems, Phys. Chem. Chem. Phys., 2011 13, 9942–9968 10.1039/C0CP02796K [3] Merlet, C., Rotenberg, B., Madden, P.A., Salanne, M.: Computer simulations of ionic liquids at elec- trochemical interfaces, Phys. Chem. Chem. Phys., 2013 15, 15781–15792 10.1039/c3cp52088a [4] Boda, D., Gillespie, D.: Calculating the electro- static potential profiles of double layers from sim- ulation ion density profiles, Hung. J. Ind. Chem., 2013 41(2), 125–132 ISSN: 0133-0276 [5] Gouy, G.: Sur la constitution de la charge electrique a la surface d’un electrolyte, J. de Phys., 1910 9(1), 457–468 10.1051/jphystap:019100090045700 [6] Chapman, D.L.: A contribution to the theory of electrocapillarity, Phil. Mag. Ser. 6, 1913 25(148), 475–481 10.1080/14786440408634187 [7] Stern, O.: Zur Theorie der elektrolytischen Dop- pelschicht, Zeit. Elektrochem., 1924 30(21–22), 508–516 [8] Bhuiyan, L.B., Outhwaite, C.W., Henderson, D.: Some simulation and modified Poisson-Boltzmann theory results for the contact values of an electrolyte near a charged electrode, J. Electroanal. Chem., 2007 607(1–2), 54–60 10.1016/j.jelechem.2006.10.010 Hungarian Journal of Industry and Chemistry DOUBLE LAYERS NEAR AN ELECTRODE 65 [9] Henderson, D., Blum, L.: Some exact results and the application of the mean spherical approxima- tion to charged hard spheres near a charged hard wall, J. Chem. Phys., 1978 69(12), 5441–5449 10.1063/1.436535 [10] Henderson, D., Blum, L., Lebowitz, J.L.: Exact for- mula for the contact value of the density profile of a system of charged hard-spheres near a charged wall, J. Electroanal. Chem., 1979 102(3), 315–319 10.1016/S0022-0728(79)80459-3 [11] Parsons, R., Zobel, F.: The interphase between mercury and aqueous sodium dihydrogen phos- phate, J. Electroanal. Chem., 1965 9(5–6), 333–348 10.1016/0022-0728(65)85029-X [12] Torrie, G.M., Valleau, J.P.: Electrical double-layers 1. Monte Carlo study of a uniformly charged sur- face, J. Chem. Phys., 1980 73(11), 5807–5816 10.1063/1.440065 [13] Torrie, G.M., Valleau, J.P.: Electrical double- layers 4. Limitations of the Gouy-Chapman the- ory, J. Phys. Chem., 1982 86(16), 3251–3257 10.1021/j100213a035 [14] Boda, D., Henderson, D., Chan, K.Y.: Monte Carlo study of the capacitance of the double layer in a model molten salt, J. Chem. Phys., 1999 110(11), 5346–5350 10.1063/1.478429 [15] Boda, D., Henderson, D., Chan, K.Y., Wasan, D.T.: Low temperature anomalies in the properties of the electrochemical interface, Chem. Phys. Lett., 1999 308(5-6), 473–478 10.1016/S0009-2614(99)00643-0 [16] Lamperski, S., Outhwaite, C.W.: Exclusion volume term in the inhomogeneous Poisson–Boltzmann theory for high surface charge, Langmuir, 2002 18(9), 3423–3424 10.1021/la011852v [17] Lamperski, S., Bhuiyan, L.B.: Counterion layer- ing at high surface charge in an electric dou- ble layer. Effect of local concentration approxi- mation, J. Electroanal. Chem., 2003 540, 79–87 10.1016/S0022-0728(02)01278-0 [18] Vatamanu, J., Borodin, O., Smith, G.D.: Molecular insights into the potential and temperature depen- dences of the differential capacitance of a room- temperature ionic liquid at graphite electrodes, J. Am. Chem. Soc., 2010 132(42), 14825–14833 10.1021/ja104273r [19] Vatamanu, J., Borodin, O., Bedrov, D., Smith, G.D.: Molecular dynamics simulation study of the interfacial structure and differen- tial capacitance of Alkylimidazolium Bis- (trifluoromethanesulfonyl)imide [Cnmim][TFSI] ionic liquids at graphite electrodes, J. Phys. Chem. C, 2012 116(14), 7940–7951 10.1021/jp301399b [20] Hu, Z., Vatamanu, J., Borodin, O., Bedrov, D.: A molecular dynamics simulation study of the electric double layer and capacitance of [BMIM][PF6] and [BMIM][BF4] room temperature ionic liquids near charged surfaces, Phys. Chem. Chem. Phys., 2013 15(34), 14234–14247 10.1039/c3cp51218e [21] Feng, G., Jiang, D., Cummings, P.T.: Curvature effect on the capacitance of electric double lay- ers at ionic liquid/onion-like carbon interfaces, J. Chem. Theor. Comp., 2012 8(3), 1058–1063 10.1021/ct200914j [22] Boda, D., Fawcett, W.R., Henderson, D., Sokołowski, S.: Monte Carlo, density func- tional theory, and Poisson-Boltzmann theory study of the structure of an electrolyte near an elec- trode, J. Chem. Phys., 2002 116(16), 7170–7176 10.1063/1.1464826 [23] Blum, L.: Theory of electrified interfaces, J. Phys. Chem., 1977 81(2), 136–147 10.1021/j100517a009 [24] Henderson, D., Blum, L.: The Gouy-Chapman the- ory as a special case of the hypernetted chain ap- proximation, J. Electroanal. Chem., 1978 93(2), 151–154 10.1016/S0022-0728(78)80228-9 [25] Woelki, S., Henderson, D.: Application of the SRISM approach to the study of the capacitance of the double layer of a high density primitive model electrolyte, Cond. Matt. Phys., 2011 14(4), 43801 10.5488/CMP.14.43801 [26] Henderson, D., Abraham, F.F., Barker, J.A.: The Ornstein-Zernike equation for a fluid in contact with a surface, Mol. Phys., 1976 31(4), 1291–1295 10.1080/00268977600101021 [27] Henderson, D., Lamperski, S., Outhwaite, C.W., Bhuiyan, L.B.: A mean spherical approximation study of the capacitance of an electric double layer formed by a high density electrolyte, Coll. Czechoslovak Chem. Comm., 2010 75(3), 303–312 10.1135/cccc2009094 [28] Carnie, S.L., Chan, D.Y.C.: The structure of electrolytes at charged surfaces: Ion–dipole mix- tures, J. Chem. Phys., 1980 73(6), 2949–2957 10.1063/1.440468 [29] Blum, L., Henderson, D.: Mixtures of hard ions and dipoles against a charged wall - the Ornstein- Zernike equation, some exact results, and the mean spherical approximation, J. Chem. Phys., 1981 74(3), 1902–1910 10.1063/1.441282 [30] Schmickler, W., Henderson, D.: The interphase be- tween jellium and a hard sphere electrolyte. A model for the electric double layer, J. Chem. Phys., 1984 80(7), 3381–3386 10.1063/1.447092 [31] Henderson, D., Smith, W.R.: Exact analytical for- mulas for the distribution functions of charged hard spheres in the mean spherical approximation, J. Stat. Phys., 1978 19(2), 191–200 10.1007/BF01012511 [32] McQuarrie, D.A.: Statistical mechanics (Univer- sity Science Books, Mill Valley), 2000 ISBN-13: 978- 1891389153 [33] Kornyshev, A.A.: Double-Layer in Ionic Liq- uids: Paradigm Change?, J. Phys. Chem. B, 2007 111(20), 5545–5557 10.1021/jp067857o [34] Winter, M., Brodd, R.J.: What Are Batteries, Fuel Cells, and Supercapacitors?, Chem. Rev., 2004 104(10), 4245–4270 10.1021/cr020730k 43(2) pp. 55-66 (2015) DOI: 10.1515/hjic-2015-0010 66 HENDERSON [35] Islam, M.M., Alam, M.T., Ohsaka, T.: Electrical double-layer structure in ionic liquids: A corrob- oration of the theoretical model by experimental results, J. Phys. Chem. C, 2008 112(42), 16568– 16574 10.1021/jp8058849 [36] Lockett, V., Horne, M., Sedev, R., Rodopoulos, T., Ralston, J.: Differential capacitance of the double layer at the electrode/ionic liquids interface, Phys. Chem. Chem. Phys., 2010 12(39), 12499–12512 10.1039/C0CP00170H [37] Wu, J., Jiang, T., Jiang, D., Jin, Z., Henderson, D.: A classical density functional theory for interfacial layering of ionic liquids, Soft Matter, 2011 7(23), 11222–11231 10.1039/c1sm06089a [38] Lamperski, S., Henderson, D.: Simulation study of capacitance of the electrical double layer of an elec- trolyte near a highly charged electrode, Mol. Sim., 2011 37(4), 264–268 10.1080/08927022.2010.501973 [39] Lamperski, S., Sosnowska, J., Bhuiyan, L.B., Hen- derson, D.: Size asymmetric hard spheres as a con- venient model for the capacitance of the electrical double layer of an ionic liquid, J. Chem. Phys., 2014 140(1), 014704 10.1063/1.4851456 [40] Henderson, D., Lamperski, S.: Simple Description of the Capacitance of the Double Layer of a High Concentration Electrolyte, J. Chem. Eng. Data, 2011 56(4), 1204–1208 10.1021/je101106z [41] Outhwaite, C.W., Bhuiyan, L.B.: An improved modified Poisson-Boltzmann equation in electric- double-layer theory, J. Chem. Soc. Faraday. Trans. II., 1983 79, 707–718 10.1039/F29837900707 [42] Jiang, J., Cao, D., Henderson, D., Wu, J.: A contact- corrected density functional theory for electrolytes at an interface, Phys. Chem. Chem. Phys., 2014 16(9), 3934–3938 10.1039/C3CP55130J Hungarian Journal of Industry and Chemistry