Ion size effects in a primitive level model of the diffuse double layer

An analytical expression is developed for the potential drop across the diffuse layer φ in terms of a cubic polynomial in the corresponding estimate in the Gouy-Chapman approximation φ GC . The coefficients of this polynomial are defined in terms of the hard sphere volume fraction η and the MSA dimensionless reciprocal distance parameter Γ . The resulting expression is shown to describe the Monte-Carlo estimates of φ obtained in a primitive level simulation of diffuse layer properties.


Introduction
We have made considerable effort in recent years to develop an analytical model for the diffuse double layer which takes into consideration the effects of finite ion size [1][2][3].Our approach is based on the equations from the hypernetted chain approximation (HNCA) as a solution of the Ornstein-Zernike integral equation for a charged wall in an electrolyte solution [4].The HNCA keeps the non-linear character of the double layer problem but it does not have an analytical solution.Thus, the estimates of the wall-ion correlation functions must be obtained by iteration.Henderson and Blum [5] described an approximate way of solving the HNC equations in which the first estimates of the wall-ion correlation functions are assumed to be equal to those given by the Gouy-Chapman (GC) theory.As is well known, this theory ignores the effects of ion size.The integrals needed to evaluate the same correlation functions in the HNCA are then evaluated using the GC estimates.These results are then assumed to be sufficiently good estimates of the wall-ion correlation functions with consideration of ion size.
The above approach which is referred to as the Henderson-Blum (HB) approach [5] is clearly approximate but a good beginning in a more complete solution of the problem.In order to assess the results of the HB approach one makes use of Monte Carlo (MC) data also obtained for a primitive level system.In the present problem, "primitive" refers to the fact that the molecular nature of the solvent is ignored and it is treated as a dielectric continuum with the permittivity of the pure solvent.In addition the dielectric properties of the wall are assumed to be the same as those of the solvent.This means that the formation of images in the wall, which occurs when it is a medium with a different permittivity such as a metal, are ignored.MC data obtained at the primitive level have been published recently for 1-1 electrolytes assuming realistic sizes for the component ions [6,7].
In the present paper our previous work [1][2][3] is extended on the basis of the equations of the HNCA, and the appropriate MC data are used to obtain an analytical expression for the potential drop across the diffuse layer for a 1-1 electrolyte.

Theory
In the GC theory for the diffuse double layer, the Poisson-Boltzmann equation is solved assuming a dielectric continuum and ions which are point charges with no volume.The potential drop across the diffuse layer ϕ d is calculated from the following equation in the case of a 1-1 electrolyte: where ϕ d is the dimensionless value of the potential drop, σ m , the charge density on the wall, Θ GC , the GC constant, and c e , the electrolyte concentration.The GC constant is given by where ε s is the relative permittivity of the pure solvent, ε 0 , the permittivity of free space, and the electrolyte concentration is expressed in M. The dimensionless potential is related to the potential in volts by the equation The value of ϕ d is now used to calculate the correlation functions for the ions with the wall.At the outer Helmholz plane (oHp), the correlation function for the cation is given by and that for the anion, It follows that the wall-sum correlation function is Combining equations ( 1) and ( 6), it follows that Here, is the dimensionless electrical field at the oHp.Now, these equations are written in the HNCA.The wall-cation correlation function at the oHp is and that for the anion, Ξ 0 is a function which depends on the volume fraction η and is calculated as an integral of the wall-sum correlation function.The ionic volume fraction for a 1-1 electrolyte with a concentration c e in M is given by N L is the Avogadro constant, and σ, the ion diameter in meters.On the other hand, H 0 is a function which depends on the reciprocal distance parameter Γ defined in the mean spherical approximation (MSA).The dimensionless expression for Γ is where κ is the Debye-Hückel reciprocal distance given by Here ε s is the relative permittivity of the solvent at 25˚C (78.46 for water) and ε 0 , the permittivity of free space.The function H 0 is calculated as an integral of the wall-difference correlation function.Exact expressions for Ξ 0 and H 0 were given in earlier papers [1,2].The wall-sum correlation function in the HNCA is When the ionic diameter goes to zero, Ξ 0 goes to unity and H 0 , to zero.Then the HNC expression for g ws reduces to that in the GC theory (equation ( 7)).Explicit expressions for ϕ d are now derived.In the GC theory, and in the HNCA, where the value obtained in the GC theory is now designated as ϕ d GC .The strategy used here is to expand the hyperbolic functions as power series in E.
In order to proceed further, the dependence of the functions Ξ 0 and H 0 on E must be established.These functions were estimated following the method suggested by Henderson and Blum [5].Accordingly, the integrals defining Ξ 0 and H 0 were calculated using the wall-ion correlation functions given by the GC theory, but limiting the field E to values in the range −3 E 3. The results for Ξ 0 could be fit to the equation where and a 1 is also a function of η.In the case of H 0 , the dependence on E is given by where and b 3 is also a function of Γ.More details about these functions are given later in this paper.

Results and discussion
MC simulations of the diffuse double layer were carried out for the restricted primitive model for five electrolyte concentrations and three ion sizes.The ion diameters chosen (200, 300 and 400 pm) cover the range found for simple monoatomic ions.The electrolyte concentrations, namely, 0.1, 0.2, 0.5, 1.0, and 2.0 M cover the range where significant departure from the GC theory is expected.Finally, data were obtained for charge densities between 5 and 40 µC cm −2 at increments of 5 µC cm −2 .Details about the procedure followed for the MC calculations are given in earlier papers [6,7].
MC data obtained for an electrolyte concentration of 0.1 M are shown in figure 1.For smaller ion diameters (200 and 300 pm) the value of ϕ d MC is always less than the GC estimate, ϕ d GC .However, for a diameter of 400 pm, the MC data fall below the GC estimates for smaller electrode charge densities but then begin to increase with respect to ϕ d GC , eventually becoming greater than ϕ d GC .In all cases, the MC data could be fit to the equation Results at a concentration of 0.5 M are shown in figure 2. The MC data are qualitatively similar to those at 0.1 M except that the values of ϕ d MC and ϕ d GC are quantitatively smaller for a given value of electrode charge density σ m .Finally, results at 2.0 M are shown in figure 3.In this case the value of ϕ d MC is always less than that  of ϕ d GC .However, the curvature of the plot of ϕ d MC against ϕ d GC remains positive for an ionic diameter of 400 pm.
In order to relate the MC results to the calculations based on the integral equation approach, equations ( 15) and ( 16) must be expanded as an infinite series in the dimensionless field E. This leads to the result It is also possible to express ϕ d as an infinite series in ϕ d GC which gives the equation The first two coefficients in this series are and If the present approach is to agree with the MC data then only the first two coefficients in equation ( 23) are significant, and equation (23) becomes the same as equation ( 21).Furthermore, we assume that the results of the HB approach give sufficiently good estimates of a 0 and b 1 and accept the values given by equations ( 18) and (20).This assumption is easily checked by comparing the values of d 1 determined in a least squares fit of equation ( 21) to the MC data with values estimated on the basis of equation ( 24).This comparison is shown in figure 4 where the value of d 1 estimated from the MC data is plotted against a −1/2 0 − Γ 2 for five different concentrations, and for three values of the ionic diameter.There is a tendency for the estimates of d 1 obtained from MC results to be too high but the maximum deviation between the value obtained from the MC data and equation ( 24) is never greater than 10 percent.Thus, it is reasonable to assume that d 1 is given by equation (24) to a very good approximation.
The MC data were then refitted to equation (21) in a one parameter least squares fit in which the values of d 1 were forced to be those given by equation (24).The resulting values of d 3 vary with both c e and σ but no clear trend is apparent.This is due to the fact that d 3 is a complex function of the four parameters, a 0 , a 1 , c 1 , and c 3 .Two of these parameters, a 0 and c 1 , determine the linear coefficient d 1 and can be accepted as correct.However, the contributions of a 1 and c 3 must be reevaluated if the curvature coefficient d 3 is to be correctly estimated.Therefore, it is reasonable to evaluate a new curvature coefficient e 3 defined as Values of e 3 are plotted against electrolyte concentration c e in figure 5. Clearly, e 3 increases with increase in electrolyte concentration but the change for σ = 200 pm is very small.In addition e 3 becomes positive for higher values of c e and σ.This means that the contribution from the term in a 1 must dominate under these conditions.In order to optimize the values of the curvature parameter the MC data for e 3 were fit to the equation where f 1 (η) is a function of the volume fraction η, f 2 (Γ), a function of the reciprocal distance parameter Γ, and α 1 and α 2 are numerical constants.Guided by the HB results various simple functions involving both integer and fractional powers of η and Γ were tested.The best fit to the MC results was obtained for the relationship Thus, we find that and These results are very different from those obtained from the analysis of the functions Ξ 0 and H 0 estimated using the HB approach.This means that the HB approach can only be used in the limit of extremely small fields.Finally, the equation used to estimate ϕ d is where the parameters d 1 and d 3 are given by and The results obtained using equations (31)-(33) are also shown in figures 1 to 3. It is clear that the revised results fit the MC data very well.Moreover, the revised estimates of ϕ d are easily estimated from the GC value ϕ d GC once the volume fraction η and the reciprocal distance parameter Γ have been calculated.Thus, the present approach gives simple analytical equations which can be used to estimate the diffuse layer potential with consideration of ion size effects.In the case of 1-1 electrolytes, the correction to the GC estimates of ϕ d are important at high concentrations and electrode charge densities.The failure of the GC theory is much more important for 2-1, 1-2, and 2-2 electrolytes.Application of the present approach to these systems will be considered in a future paper.

Figure 1 .
Figure 1.Plots of the potential drop across the diffuse layer ϕ d according to the MC simulation (filled symbols) and the present model (open symbols) for an electrolyte concentration of 0.1 M and for three ion sizes: 200 pm ( , ), 300 pm (•, •) and 400 pm ( , ).The results at 300 and 400 pm are shifted vertically by 2 and 4 units, respectively, for the sake of clarity.

Figure 2 .
Figure 2. As in figure 1 but for an electrolyte concentration of 0.5 M.

Figure 3 .
Figure 3.As in figure 1 but for an electrolyte concentration of 2.0 M.