J. Phys. Chem. B 2007, 111, 1443-1448
1443
Series Approach to Modeling Ion Size Effects for Symmetric Electrolytes in the Diffuse Double Layer Thomas G. Smagala and W. Ronald Fawcett* Department of Chemistry, UniVersity of California, DaVis, California 95616 ReceiVed: October 26, 2006; In Final Form: December 14, 2006
An analytical model is developed for the potential drop and differential capacity across the diffuse layer which considers the effects of ion size on these properties. For symmetric electrolytes, this potential drop is expressed in terms of a cubic polynomial in the corresponding estimate in the Gouy-Chapman theory. Optimal polynomial coefficients and model validation for 1:1 and 2:2 electrolytes are provided by fits of Monte Carlo data obtained for a restricted electrolyte in a primitive solvent. Simple relationships between these coefficients and parameters commonly associated with the mean spherical approximation are obtained. It is shown that the series approach accurately describes potential drops and differential capacities of the diffuse layer for 1:1 and 2:2 electrolytes for the chosen assumptions.
Introduction In our recent work,1-4 we have emphasized improvements to the classical Gouy-Chapman (GC) theory5 that are easily understood and applied by experimentalists. In one approach, an expression for the potential drop across the diffuse layer was developed as a series in the corresponding GC estimate.1,2 The coefficients of the series can be related to parameters associated with ion size effects in the mean spherical approximation (MSA). This approach can also be used to estimate the differential capacity of the diffuse layer, but it does not give expressions for the ionic surface excesses or the potential distribution in the diffuse layer. In an alternative approach,3,4 an empirical model for all important diffuse layer properties was developed on the basis of empirical expressions for the ion profiles in the diffuse layer. However, our previous work was limited to 1:1 electrolytes and cannot deal with the very interesting effects seen in MC data when ions of higher charge are involved.6,7 In the present paper, the series approach is improved for 1:1 electrolytes and extended to 2:2 electrolytes. In earlier work1,2,8 the functions that account for the effects of ion size were estimated on the basis of the Henderson-Blum (HB) approximation.9,10 More specifically, they were estimated assuming that the ion profiles in the diffuse layer are those given by GC theory. In the present paper, the HB approximation is avoided, and the ion profiles are taken directly from the MC simulations. In this way, we have been able to extend the series approach to electrolytes with multivalent ions. The details of our recent work are given below. A variety of different and more accurate double layer theories exist.11-18 The series approach used here is derived from the hypernetted chain approximation integral equation,14,15,19 as applied to symmetric electrolytes. The term symmetric refers to the equal and opposite charges on the cations and anions. A method to estimate the diffuse layer potential drop and differential capacity for restricted electrolytes in a primitive solvent * Corresponding author. Phone: (530) 752-1105. Fax: (530) 752-8995. E-mail:
[email protected].
is presented. Validation of the new model for 1:1 and 2:2 electrolytes is provided by MC simulations. Series Expressions Obtained from the Integral Equation Approach. Details of the integral equation approach (HNC/ MSA theory) are found in the works of Carnie et al.14 and Lozada-Cassou et al.15 Numerical solution of these equations coupled with the Poisson equation requires iteration to solve for the ion-wall correlation functions, gi(x), and the dimensionless potential profile at position x in the diffuse layer, φ(x). The dimensionless potential is related to the actual potential, φ(x), by the equation
φ(x) ) Fφ(x)/RT
(1)
The origin of the system, x ) 0, is defined at the edge of the diffuse layer, also known as the outer Helmholtz plane. Thus, the potential at x ) 0 is the potential drop across the diffuse layer, φd, and is our main quantity of interest. In this work, the focus is on symmetric electrolytes. Species 1 is defined to be the cation with a charge of +z, and species 2 is the anion with charge -z. For convenience, the wall-sum and wall-difference correlation functions are introduced:
gws(x) ) [g1(x) + g2(x)]/2
(2)
gwd(x) ) [g2(x) - g1(x)]/2
(3)
and
Solution of the HNC/MSA equations is numerically tedious, and thus, simpler but more approximate methods have been devised. The series approach takes advantage of the condition for the contact value of the density profile for charged hard spheres near a charged wall, which for a symmetric electrolyte has the form
gws(0) ) a0 + E2/2
(4)
where a0 is the compressibility factor and E is the dimensionless field given by
10.1021/jp067039w CCC: $37.00 © 2007 American Chemical Society Published on Web 01/25/2007
1444 J. Phys. Chem. B, Vol. 111, No. 6, 2007
E ) σm/AGC
Smagala and Fawcett
(5)
The quantity σm is the charge density of the wall, not to be confused with σ, the ionic diameter. The Gouy-Chapman constant is defined as
AGC ) (2RT0sce)1/2
(6)
where 0 is the vacuum permittivity, s is the relative permittivity of the pure solvent, and ce is the bulk concentration in SI units. The currently accepted expression for a0 is
1 + η + η 2 - η3 Γ3 18η (1 - η)3
a0 )
(7)
which contains the MSA hard-sphere result of Carnahan and Starling20 and an electrostatic term.18 The MSA volume fraction, η, and reciprocal distance parameter, Γ, are defined as
η ) πFσ3/6
(8)
Γ ) (1 + 2κσ)1/2/2 - 1/2
(9)
and
for φd demonstrate that higher-order polynomials in E are required to match the curvature found from the MC work. This is the main reason why we1,2 introduced the substitution
E ) 2 sinh(zφdGC/2)
(17)
where φdGC is the GC estimate of the potential drop across the diffuse layer. The resulting series expression has the form
φd(φdGC) ) d1φdGC + d3(φdGC)3 + ...
(18)
and was found to be in excellent agreement with values of φd obtained from MC simulations using only the first two terms of the series. This occurs because the substitution removes most of the strongly nonlinear curvature and allows MC data for φd to be fitted by a simple polynomial equation. The success of expressions of the form of eq 18 clearly indicates that φd for 1:1 electrolytes is accurately modeled by an odd series, usually truncated to an odd cubic polynomial. The procedure is readily extended to arbitrary symmetric electrolytes. Rather than directly using eq 11, we instead introduced the approximate expression commonly used for gwd(x) to obtain
where κ is the Debye-Hu¨ckel reciprocal length. Fσ3 is the dimensionless ion density for a symmetric electrolyte, in which
Ξ0 cosh(zφd + Η0) ) Ξ0 + E2/2
F ) NL(c1 + c2)
Both eqs 11 and 19 yield odd series when expanded as a Taylor series of the form φd(φdGC), but the series generated from eq 19 is far simpler with nearly equal accuracy. The series coefficients are related to the parameters for Ξ0 and Η0, defined in eqs 14 and 15 by the expressions
(10)
NL being the Avogadro constant. When eq 4 is combined with HNC/MSA theory, the result is
Ξ0 cosh(zφd + Η0) ) a0 + E2/2
(11) d1 )
The quantities Ξ0 and Η0 are defined as
Ξ0 ) exp[2πFσ3
∫-11 (gws(t) - 1)Chs(t) dt]
(12)
and
Η0 ) 2z πFσ 2
3
∫0
1
and
d3 ) gwd(t)Ces(t) dt
(13)
The short-range direct correlation function, which is decomposed into a hard-sphere contribution, Chs(x), and an electrostatic contribution, Ces(x), is assumed to be given by the bulk MSA. The full mathematical definitions of these quantities are shown in previous work.10,14,15,19 We found earlier1,2 that the quantities Ξ0 and Η0 are well-represented by simple analytical functions,
Ξ0 ) a0 sech(a1E)
(14)
Η0 ) b1E + b3E3
(15)
and
It will be shown later that these same expressions work equally well in the present series approach. Using eqs 11, 14, and 15, it is possible to expand φd as a Taylor series to yield a rather simple odd series in E.
φd(E) ) c1E + c3E3 + ...
1 a01/2
- b1
(
6a12 1 1 1 + - b1 - 24b3 24 a 1/2 a 1/2 a 3/2 0 0 0
(20)
)
(21)
The remaining challenge is to find expressions for the unknown parameters, a1, b1, and b3. In our original work, a0 was also treated as an adjustable parameter, and the optimal value was determined to equal the hard-sphere compressibility factor of Carnahan and Starling.20 As pointed out above, Ξ0 and Η0 were previously estimated using ion density profiles from the GC model. However, these profiles are available from MC simulations. Thus, Ξ0 and Η0 can be more accurately calculated by numerical integration using the MC results. This is the method used throughout this work. The revised fits of Ξ0 and Η0 can then be used to relate the generalized polynomial coefficients, di, to the MSA parameters, η and Γ. Diffuse Layer Differential Capacity. The analytical expressions derived in the previous sections for the potential drop across the diffuse layer provide a rapid means of calculating the differential capacity of the diffuse layer, Cd. Once again, dimensionless variables are introduced for simplification. The dimensionless diffuse layer capacity is defined as
(16)
Truncation of this series after the cubic term yields a very simple result that can mimic the change in slope exhibited by φd as a function of E at high E. However, comparisons with MC data
(19)
C ˜d )
RT Cd F AGC
(22)
For theoretical calculations, C ˜ d is computed by differentiation.
Series Approach to Modeling Ion Size Effects
C ˜d )
dE d
J. Phys. Chem. B, Vol. 111, No. 6, 2007 1445
(23)
dφ
Using this definition, the GC differential capacity for a general symmetric electrolyte is
C ˜ d,GC ) z cosh(zφdGC/2)
(24)
At the point of zero charge (pzc), C ˜ d,GC (1:1) and C ˜ d,GC (2:2) are 1 and 2, respectively. The full forms of the Gouy-Chapman diffuse layer differential capacities are necessary for the calculation of C ˜ d in the series approach. For the case of symmetric electrolytes, the estimate of C ˜ d by the series approach follows directly from eq 18.
dφd 1 1 ) (d + 3d3(φdGC)2) ) C ˜ d dE C ˜ d,GC 1
(25)
For any electrolyte, the series approach provides a very simple result for calculating Cd at the pzc.
C ˜ d (pzc) ) C ˜ d,GC (pzc)/d1
(26)
Because of the scatter and small sampling of our MC data, accurate numerical differentiation is impossible. Thus, a comparison with the differential capacity of the diffuse layer estimated from the MC data is not shown in the subsequent section. Results and Comparisons with Monte Carlo Simulations. The details of the Monte Carlo simulations are discussed in our previous work.2 For each electrolyte, the MC calculations were performed for five different concentrations, eight different positive charge densities, and two or three different ion sizes. By symmetry, the potential profiles for the negative side are equal in magnitude, but opposite in sign that of the positive side potentials. This is why an antisymmetric odd series works well for describing φd. All calculations are carried out at room temperature, T ) 298 K, and assume a primitive solvent that has the same relative permittivity as water, s ) 78.46. Previous comparisons of diffuse layer models15,16 with MC data6,7 used the rather large ion size of 425 pm. This is greater than the diameter of most unsolvated alkali metal and halide ions.5 Instead, we choose a diameter of 300 pm for most calculations. This seems to be a more physically realistic choice to represent the diameter of an average unsolvated monatomic ion; however, calculations at 200 and 400 pm are also shown for comparison. The potential drop across the diffuse layer as a function of charge density was tabulated from MC simulations. The MC values of φd were plotted as a function of φdGC and fitted to polynomials of the form of eq 18. The resulting best fit curves for 1:1 and 2:2 electrolytes are shown in Figure 1. The best fit coefficients are given in Table 1. The broken line with unit slope shows the GC results for reference. The MC data and the series approach deviate from this line by an amount that depends upon concentration and ion size. In all cases, GC theory becomes progressively less valid at higher charge densities. Figure 1 clearly shows how GC theory fails to predict the change in slope and the extremum that occurs in the case of 2:2 electrolytes. Fluctuations in the MC data occurred as a result of the finite number of particles used in the simulations. Aside from these fluctuations, the series approach matches the behavior of φdMC quite successfully. The differential capacity of the diffuse layer is easily computed from the polynomial fits of Figure 1. Note that C ˜d
Figure 1. Dimensionless diffuse layer potential drops versus their Gouy-Chapman estimates for 1:1 and 2:2 electrolytes. Solid curves show odd cubic fits for concentrations of 0.1, 0.2, 0.5, 1, and 2 M from top to bottom. Dashed line shows GC theory. Points indicate MC data for concentrations of 0.1 (O), 0.2 (3), 0.5(0), 1(]) and 2 M (∆).
approaches infinity whenever φd exhibits an extremum. Thus, it is easier to plot 1/C ˜ d as a function of φdGC, the results for 1:1 and 2:2 electrolytes being shown in Figure 2. GC values are computed from eq 24. As stated previously, comparisons with MC data are not feasible. Once again, GC theory predicts a single curve, whereas the series approach shows branching dependent upon concentration and ion size. More sophisticated theories of the diffuse layer14-19 and experimental data5 suggest that C ˜ d at the pzc varies with both concentration and ion size. Table 2 shows how C ˜ d (pzc) varies over five different concentrations in the series approach, calculated by means of eq 26. Monte Carlo simulations are tedious, and thus, it is desirable to derive explicit analytical expressions for the coefficients di that do not require the use of MC data. Previously, we1,2,4 obtained these expressions from fits of the HB estimates of Ξ0 and Η0. In this work, however, the more accurate MC estimates of these quantities were used. The integrands of eqs 12 and 13 are gi(x) profiles from scattered MC data, which include fluctuations due to a finite number of particles used in the simulations. Consequently, there will be more scatter in the results for Ξ0 and Η0. Fortunately, unlike numerical differentiation, numerical integration leads to a reduction in random error associated with the integrand. The same simple analytical
1446 J. Phys. Chem. B, Vol. 111, No. 6, 2007
Smagala and Fawcett
TABLE 1: Best-Fit Coefficients from Least-Squares Regression of MC Diffuse Layer Potential Drop to the Corresponding GC Value (eq 18)a 1:1 electrolytes 200 pm 0.1 M 0.2 M 0.5 M 1M 2M a
2:2 electrolytes
300 pm
400 pm
300 pm
400 pm
d1
d3
d1
d3
d1
d3
d1
d3
d1
d3
0.978 0.971 0.978 0.993 0.961
-0.005 -0.006 -0.009 -0.014 -0.02
0.984 0.959 0.929 0.868 0.790
-0.0048 -0.0053 -0.0063 -0.0065 -0.0070
0.918 0.919 0.839 0.758 0.566
-0.0010 -0.0015 0.0011 0.0052 0.0180
1.135 1.127 1.111 0.926 0.859
-0.095 -0.125 -0.183 -0.226 -0.360
0.878 0.868 0.829 0.545 0.174
-0.072 -0.098 -0.157 -0.170 -0.190
Each set of data was fitted to an odd cubic polynomial.
Figure 2. Reciprocal dimensionless differential capacity of the diffuse layer versus φdGC for 1:1 and 2:2 electrolytes. Dashed curve shows GC theory. Solid curves are derived from the cubic fits of Figure 1 and are for concentrations of 0.1, 0.2, 0.5, 1, and 2 M from top to bottom.
TABLE 2: Dimensionless Differential Capacity of the Diffuse Layer at the Point of Zero Charge for 1:1 and 2:2 Electrolytes as Predicted by the Series Approach for Five Different Concentrations and an Ion Size of 300 pm 1:1 2:2
0.1 M
0.2 M
0.5 M
1M
2M
1.02 1.76
1.04 1.77
1.08 1.80
1.15 2.16
1.27 2.33
expressions given by eqs 14 and 15 work equally well for fitting the MC estimates of Ξ0 and Η0 in both 1:1 and 2:2 electrolytes. The resulting plots of the MC estimates of Ξ0 and Η0 along with their fits to eqs 14 and 15 are shown in Figure 3 for the case of a 1:1 electrolyte with an ion size of 300 pm. All ion
Figure 3. Plots of Ξ0 and Η0 as a function of E for 1:1 electrolytes with an ion size of 300 pm. MC estimates are computed using eqs 12 and 13 and are shown as points with the same meaning as in Figure 1. Solid curves show best fits of these MC estimates to eqs 14 and 15. Plots for other ion sizes and for 2:2 electrolytes have very similar shapes.
sizes considered for both 1:1 and 2:2 electrolytes exhibited very similar results. In our previous work,1,2,4 the coefficients a0, a1, b1, and b3 obtained from fits of Ξ0 and Η0 to eqs 14and 15 were further refined. Simple representations of each coefficient were obtained in terms of the MSA parameters, η and Γ. The same approach was attempted in this work. It was found that expressions obtained in this manner for a0, a1, b1, and b3 from very accurate fits of both Ξ0 and Η0 resulted in unsatisfactory representations of d1 and d3. However, those expressions served as a useful guideline for the optimal functional form of the di coefficients
Series Approach to Modeling Ion Size Effects
J. Phys. Chem. B, Vol. 111, No. 6, 2007 1447
Figure 4. Series coefficients obtained from regression of MC data as a function of their theoretical estimates from MSA parameters using eqs 20 and 21 and 27-32. Data are shown for both the linear (d1) and cubic (d3) coefficients of the series for both 1:1 and 2:2 electrolytes. The solid lines have a slope of unity. Points indicate coefficient data for ion sizes of 200 (2), 300 (b), and 400 pm (1).
in terms of the MSA parameters. It must be emphasized that accurate analytical representations of Ξ0 and Η0 are not the primary objective. From a practical sense, accurate expressions for d1 and d3 are far more important for the analytical computation of φdGC. With a0 given by eq 7 and using the definition of d1 given by eq 20, the following expressions for b1 were obtained:
b1(1:1) ) 0.5Γ
(27)
and
b1(2:2) ) Γ0.75
(28)
These expressions are valid for all the different concentrations and ion sizes used. Different expressions for the compressibility were considered, but the optimal a0 is given by eq 7. Similarly, expressions for a1 and b3 were obtained by forcing them to accurately fit the coefficient data for d3 via eq 21:
a1(1:1) ) 1.25 η0.5
(29)
a1(2:2) ) 1.1 η0.5
(30)
b3(1:1) ) (1/7)Γ1.75 - (4/25)η0.5
(31)
b3(2:2) ) 0.8Γ4 - 75η2
(32)
and
Combining eq 18 with eqs 20, 21, 25, and 27-32 yields a
complete analytical model for the computation of potential drops and differential capacities of the diffuse layer. Interested users can now apply the analytical model directly without performing MC simulations. Comparisons of each coefficient obtained from regression of the MC data with its theoretical estimate from the MSA parameters are shown in Figure 4. Discussion and Conclusion There have been many advances in recent years in modeling the diffuse double layer. The increased mathematical and computational complexity of these theories hinders their application to experimental electrochemical systems. However, the simple analytical models available are usually not sophisticated enough to accurately describe reality. Gouy-Chapman theory is completely inappropriate for describing diffuse layer potential profiles and differential capacities in higher valence electrolytes. The series approach described here provides a simple and effective way of obtaining improved estimates of φd and Cd for a restricted electrolyte. The capacity of the diffuse layer has often been used to estimate the capacity of the inner layer.5 Another application first described by Parsons and Zobel21 is to test for the presence of ionic specific adsorption. In such a test, the reciprocal of the experimental capacity Cex measured for various electrolyte concentrations at constant charge density is plotted against the corresponding reciprocal of the diffuse layer capacity. Using the relationship for capacitors in series, if ionic specific adsorption is absent,
1448 J. Phys. Chem. B, Vol. 111, No. 6, 2007
Smagala and Fawcett require long calculation times.6,7 The series approach presented in this work is the first completely analytical model to accurately describe nonmonotonic potential drops in the diffuse layer. Phenomena such as ion pairing that are certainly present in 2:2 electrolytes are not taken into account in the present model. To deal with ion pairing, one additional parameter is required: namely, the ion pairing equilibrium constant. In addition, the effects of unequal ion size, change in solvent permittivity with electrolyte concentration, and the ionic images formed in the conducting electrode should also be considered in developing a more realistic picture of the diffuse double layer. We hope the series approach presented here can provide a useful compromise between numerically tedious theories and the overly simple GC theory. Application to asymmetric electrolytes will be considered in a future paper.
Figure 5. Plot of the reciprocal differential capacity of the diffuse layer obtained from the series approach versus its corresponding estimate in GC theory for 1:1 electrolytes. Points were calculated using eqs 20, 21, 25, and 27-32 for charge densities of 0 (2), 0.04 (b), and 0.08 C m-2 (1) at concentrations from 0.01 M at the upper right to 1 M at the lower left. The latter two sets of data are shifted to the right by 2 and 4 m2 F-1, respectively, for clarity.
1 1 1 ) + Cex Cin Cd
(33)
It follows that such a plot should have unit slope and that an estimate of Cin can be obtained from the intercept on the y-axis. Thus, an important question to ask is whether Parsons-Zobel plots are valid when the value used for Cd is that from GC theory. To answer this question, plots of 1/Cd estimated using the present model were plotted against 1/Cd,GC for charge densities in the neighborhood of the pzc where the ParsonsZobel analysis is applied (see Figure 5). These plots are excellent straight lines with unit slope; however, they do not go through the origin. Thus, the consequence of using the GC estimate of Cd would be that an error is made in the estimate of Cin. The conclusion that ionic specific adsorption is absent would not be affected. Although the series approach does not provide expressions for gi(x) and φ(x), it is superior to the empirical model in matching the nonmonotonic behavior of φd and Cd at high charge densities. Previously, this complicated behavior could be modeled only by more sophisticated theories that require numerical solution11-16,19 or by Monte Carlo methods that
Acknowledgment. The authors are grateful to Dr. Dezso¨ Boda for making available his Monte Carlo program. Financial support was received from the National Science Foundation, Washington through a research grant (CHE 0451103). References and Notes (1) Fawcett, W. R.; Smagala, T. G. Condens. Matt. Phys. 2004, 7, 709. (2) Fawcett, W. R.; Smagala, T. G. J. Phys. Chem. B 2005, 109, 1930. (3) Smagala, T. G.; Fawcett, W. R. Z. Phys. Chem. 2006, 220, 427. (4) Fawcett, W. R.; Smagala, T. G. Langmuir 2006, 22, 10635. (5) Fawcett, W. R. Liquids, Solutions, and Interfaces; Oxford University Press: New York, 2004; Chapters 3 and 10. (6) Torrie, G. M.; Valleau, J. P. J. Chem. Phys. 1980, 73, 5807. (7) Torrie, G. M.; Valleau, J. P. J. Phys. Chem. 1982, 86, 3251. (8) Fawcett, W. R.; Henderson, D. J. Phys. Chem. B 2005, 109, 22608. (9) Henderson, D.; Blum, L. J. Electroanal. Chem. 1980, 111, 217. (10) Henderson, D.; Blum, L. Can. J. Chem. 1981, 59, 1906. (11) Levine, S.; Outhwaite, C. W. J. Chem. Soc., Faraday Trans. 2 1978, 74, 1670. (12) Outhwaite, C. W.; Bhuiyan, L. B.; Levine, S. J. Chem. Soc., Faraday Trans. 2 1980, 76, 1388. (13) Levine, S.; Outhwaite, C. W.; Bhuiyan, L. B. J. Electroanal. Chem. 1981, 123, 105. (14) Carnie, S. L.; Chan, D. Y. C.; Mitchell, D. J.; Ninham, B. W. J. Chem. Phys. 1981, 74, 1472. (15) Lozada-Cassou, M.; Saavedra-Barrera, R.; Henderson, D. J. Chem. Phys. 1982, 77, 5150. (16) Lozada-Cassou, M.; Henderson, D. J. Phys. Chem. 1983, 87, 2821. (17) Boda, D.; Fawcett, W. R.; Henderson, D.; Sokołowski, S. J. Chem. Phys. 2002, 116, 7170. (18) Boda, D.; Henderson, D.; Miery Teran, L.; Sokołowski, S. J. Phys. Condens. Matter 2002, 14, 11945. (19) Smagala, T. G.; Fawcett, W. R. J. Electroanal. Chem. 2006, 597, 43. (20) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969, 51, 635. (21) Parsons, R.; Zobel, F. G. R. J. Electroanal. Chem. 1963, 9, 333.
