A modified Poisson-Boltzmann approach to homogeneous ionic solutions

The mean electrostatic potential approach to ionic solutions was initiated by the mean field work of Gouy and Chapman for inhomogeneous systems, and Debye and Hückel for bulk solutions. Any successful extension of the mean field theories requires an adequate treatment of both the exclusion volume and fluctuation terms. One such development has been the modified Poisson-Boltzmann approach. Although the bulk modified PoissonBoltzmann theory was introduced 35 years ago, only a limited amount of work has been put into its development due to the successful application of liquid state theories to ionic systems. We outline here the bulk modified Poisson-Boltzmann equation, comment on some of its successes, and mention some outstanding problems.


Introduction
Electrolyte solutions arise in a wide variety of chemical, physical and biological systems.Any true understanding of the structure and thermodynamics of these systems requires an appropriate molecular description based on statistical mechanics.The first successful approach was the classical work of Debye and Hückel (DH) [1] using the primitive model electrolyte.Subsequent progress was limited until the application of liquid state theories [2] and computer development.Increased computer power and speed now enables ever more models to be accurately simulated and more complex theoretical analysis to be evaluated numerically.Many excellent reviews have appeared over the years and we list a few of the most recent [3][4][5][6][7][8][9].
In this article we outline and comment upon the modified Poisson-Boltzmann (MPB) theory.An early review of this theory is given in reference [10].The basic idea of the MPB approach is to extend the classical mean electrostatic approach c C.W.Outhwaite of DH by expressing the distribution functions in the Kirkwood, Bogolubov, Born, Green, Yvon (KBBGY) hierarchies [11] in terms of mean electrostatic potentials.The resulting reformulation of the hierarchies allows both a clear identification of the mean field approximations in the DH theory and a natural framework for higher order closures.A further advantage is that any numerical work involving mean electrostatic potentials deals with essentially damped exponential functions.

Theory
We assume that the potential of mean force U N for a neutral system of N ions can be expressed as where e i e j 4πε 0 ε r r = e i e j ν ij (2) and u s ij is any short ranged spherically symmetric interaction.Here e i , e j are the charges on i and j respectively, r = r ij their separation, ε 0 the vacuum permittivity and ε r the solvent relative permittivity.When dealing with the primitive model (PM) electrolyte where a i , a j are the radii of ions i and j respectively with a ij = a i + a j .The special case of a i = a j is the restricted primitive model (RPM) electrolyte.
Working in the canonical ensemble, the mean electrostatic potential ψ({n};q) at the field point r q for ions fixed at r 1 ,. . ., r n is e k ν(k, q) + s e s ν(n + 1, q)ρ s ({n} ; n + 1) d(n + 1), (4) where the sum in s is over the ion species, n N and is the conditional number density of ions type s at r n+1 given ions at r 1 , . . ., r n and β = 1/k B T , k B Boltzmann's constant and T the absolute temperature.Interchanging q, n +1 in equation (4), e k ν(k, n + 1) + s e s ν(q, n + 1)ρ s ({n}; q) dq (6) and applying the Laplacian operator to r n+1 , gives the generalised Poisson equation where δ(r) is the three dimensional Dirac delta function.The first term on the right hand side of equation (7) gives the charge contribution from the n fixed ions while the second term is the mean charge at r n+1 from the mobile ions given n fixed ions.
For n = 1, equations ( 6), (7) reduce to where ρ s is the mean number density of ions type s and g is the pair distribution function for ions i and s.Equation ( 9), with formal solution equation ( 8), is the familiar Poisson equation of the DH theory.When n = 2, equations ( 6), (7) are The application of equations such as ( 8)- (11) depend upon knowing a practical relationship between the conditional densities and the mean electrostatic potentials.A convenient relationship is provided by the Kirkwood charging process for a single ion [12,13], although the BBGY formulism can be employed [14].We suppose that the ion i at r 1 has a charge λ 1 e i (0 λ 1 1) so that Then taking the natural logarithm of the number densities for {n} and {n+1} ions, differentiating each with respect to λ 1 , subtracting and using the definition of the mean electrostatic potential equation (6) gives for n 1.In particular for n = 1, equation ( 13) reduces to Writing ψ(1, 2; q) = ψ(1; q) + ψ(2; q) + φ(1, 2; q), where φ(1,2;q) is the fluctuation potential, then equation ( 14) simplifies to on using ρ s (1; 2) = ρ s g ij .The fluctuation potential describes the departure from linear superposition, for two fixed ions, of the mean electrostatic potential at r q from that of the two individual ions.Alternatively charging an ion j at r 2 , so that Poisson's equation ( 9) becomes where ζ i,s = g is (λ 2 = 0).

Poisson-Boltzmann theories
The Poisson-Boltzmann equation is derived from equation ( 18) by putting φ(1, 2; 2) = 0 and ζ i,s = 0 for r less than the distance of closest approach of i and s, and 1 otherwise.These two approximations were first delineated by Kirkwood [15] although his fluctuation term is not cast in a form suitable for further analysis.A drawback of using the classical Poisson-Boltzmann theory is its unsymmetrical g ij , which restricts any consistent use to ions of equal size and valence.This drawback stems from the unsymmetrical formulation of g ij in equations ( 16), (17).A symmetrical g ij can be derived by putting e i = 0 in equation ( 17) and e j = 0 in equation ( 16), substituting for ζ j,i , ζ i,j in equations ( 16) and ( 17) respectively, and then combining the results to give where ) and ψ 0 j = ψ j (2; 1 |e j = 0) .Neglecting the fluctuation potentials in equation ( 19) gives a symmetric mean field radial distribution function, and from equation ( 9) the corresponding set of symmetric PB (SPB) equations [16] ) for t running over the ion species.If all the ions have the same size then the discharged potentials ψ 0 t are zero.The SPB formulation enables a mean field treatment to be made of single, or mixtures, of electrolytes with unequal radii and valences.
The linear form of equation ( 20) can be solved for a single PM electrolyte, species p and q say, when ψ 0 p = ψ 0 q = 0 and g 0 pq = H(r − a pq ) where H(r) is the Heaviside unit step function [17,18].For r 2a qq , with a p a q , where κ = [(β/ε 0 ε r ) s e 2 s ρ s ] 1/2 is the DH constant and B, C are constants depending on the system parameters.When the ions have the same size C = 0.The first term on the right hand side is the classical screened potential with an effective charge while the second term is also screened but with a larger screening length.This second term arises because the different exclusion volumes of the ions leads to a coupling of the two ion species in the pair of differential equations.

Modified Poisson-Boltzmann theories
To improve upon the mean field formulation we need a procedure to calculate the fluctuation potential φ(1,2;3).Consider the triplet distribution function g ijk (1,2,3).Formally we can write where W ij , W ijk are the pair and triplet potentials of mean force respectively and w ijk is the departure from linear superposition of the pair potentials for three ions.
The unsymmetrical MPB closure is given by [12] w ijk = e k φ ij (1, 2; 3) (24) so comparing equations ( 8) and ( 23), the MPB closure is the replacement of the departure of linear superposition of the potential of mean force by the corresponding energy e k φ(1,2;3).This is in exactly the same spirit as the DH closure at the previous hierarchical level.The closures ( 24) and ( 25) are not symmetric, but symmetric closures can be obtained from the symmetric formulations The symmetric closures ( 26) and ( 27) do not imply fully symmetric g ij and g ijk respectively as the neglect of the discharged potentials means the distribution functions are inconsistent when an ion is discharged.Nevertheless, a symmetric g ij will be given by equation (19) even if an estimate of φ(1,2;3) is found using either (24) or (27).

Calculation of the fluctuation potential
For simplicity we now restrict ourselves to a RPM electrolyte with ions of diameter a. From equations ( 9) and (11) we have where equations ( 33), (34) are exact but equation (35) incorporates the MPB unsymmetrical closure (24).Subtracting equations ( 33) and ( 34) from (35) gives the differential equation for φ (= φ(1,2;3)), To make analytical progress we linearize equation (36) outside the exclusion volume of the two ions.Writing ω k as the exclusion volume r k3 < a and putting or using a more approximate equation outside the exclusion volumes where equations (39-41) are exact.Figure 1 illustrates the fluctuation potential problem for intersecting and non-intersecting exclusion volumes.The boundary conditions associated with the potential problem are (i) φ and the normal derivative ∂φ/∂n are continuous at the boundaries, (ii) φ → 0 as any of r 13 , r 23 , r 12 → ∞.
An approximate solution can now be found for large r (= r 12 ) when we can assume spherical symmetry around the individual exclusion volumes [14].Considering the exclusion volume centered at r 2 , then the general solution is where r 2 = r 23 and with θ 2 the angle between r and r 2 .So as φ(1, 2; 2) = lim we have from equation ( 17) that so that only B 0 is required.A similar analysis for the BBGY hierarchy requires only B 1 [14].Application of the boundary conditions (i) at r 2 = a gives where y = κa and c 0 = ∂c/∂r 2 .From equation ( 44) we can derive where u = rψ(1; 2).Similarly the corresponding B 0 for the exclusion volume centered at r 1 can be found, hence using the symmetric formulation (19) for g ij gives where The corresponding analysis for unequal ion sizes cannot be easily carried out as there is more than one exclusion volume region about an ion.For a single electrolyte with species p having the smaller radius, and ignoring the different exclusion regions about an ion so that equation (38) still holds [19], where with L 0 s = L s (e s = 0).The approximation outside the exclusion volumes becomes increasingly inaccurate as the size discrepancy increases and allows too great a weight to the approach of like ions.Putting L s = ψ s gives the SPB radial distribution function.
Combining the g ij RPM result (50) with Poisson's equation (9) gives which for a single electrolyte of species p, q is a pair of coupled equations for u p , u q .
If the valences are equal, then u p = −u q = u say and there is only one equation to solve.The PB equation is then derived from equation ( 54) by putting a = 0 in L(u) and g 0 st = H(r − a).
For low y there are two real roots with the lower corresponding to the DH solution.As y increases the two real roots converge, coalesce at a "Kirkwood value"of y C = 1.2412, then move off as complex conjugates.Eventually this complex conjugate pair merges and becomes purely imaginary at y I = 7.83.The other roots are complex conjugates, with larger real parts, and infinite in number.The linear equation thus predicts a damped solution for y < y C , a damped oscillatory solution for y C < y < y I and for y > y I , a purely oscillatory solution which is unjustifiable mathematically.The behaviour of the linear solutions for L(u) given by equations (57), (60), (61) are similar, with the corresponding values of y C being 1.032, 1.720, 1.642 respectively, and of y I being 2.79, 17.7, 9.17 respectively.Casting the mean spherical approximation (MSA) in terms of the mean electrostatic potential also leads to equation (56), with the MSA L(u) predicting y C = 1.229 and y I = ∞ [23].Another linear theory, the generalised DH of Lee and Fisher [24], predicts y C = 1.178 and y C = 6.652.All these linear theories predict the same value of y C and y I for both symmetric and asymmetric valences.Indeed, for a symmetric valence RPM, the value of y C for a particular linear theory is unaltered if ψ(1;2) is linear in charge or the direct correlation function c ij is given by where c 0 is short-ranged and both c 0 and f (r) are independent of e i and e j [25].
Investigations of the hypernetted chain (HNC) equation [8,[26][27][28] indicate that y C is dependent on the electrolyte valence.Also, from diagrammatic techniques [29] and dressed ion theory [30], the asymptotic behaviour of the screening length for small values of the coupling parameter Γ = βκe 2 /4πε 0 ε r , e elementary unit of charge, has a purely electrostatic contribution for a → 0. Unequal ion sizes further complicate the screening length behaviour [31].The unequal size L(u) of equation (53) predicts variations in y C and y I for different ion sizes and valences [19], but its inherent approximations limit its applications.
The combination of the neutrality condition with the second moment condition of Stillinger and Lovett [32] predicts "charge oscillations" for y > √ 6.This value of y = √ 6 is also the limiting value for which a single screening length is possible which allows these two conditions to hold [33].A simple derivation of this limitation is given [25] by considering ψ i = A 1 exp(−κ 1 r) and then finding κ 1 for which ψ i satisfies both the neutrality condition and the PM mean potential reformulation of the second moment condition [34] Numerical solutions of the non-linear SPB and MPB equations ( 20), ( 52), (54) have been found for various PM systems.Once g 0 st is specified, then for m ionic species there are only m coupled equations to solve for ψ s , s = 1, . . ., m, rather than the m(m + 1)/2 equations in a standard integral equation approach.Early work utilized a density expansion, but now the Percus-Yevick hard sphere distribution function for mixtures, with the Verlet-Weis correction [35], is used.Note that g 0 st is not the uncharged hard sphere radial distribution function but the pair distribution for two discharged ions in a sea of fully charged ions.In the BBGY potential formulation the analogous distribution function for an uncharged ion is given by the Kirkwood superposition approximation in the hierarchy.The numerical solution of the SPB and MPB equations have been based on a quasi-linearization technique [36], with a fixed point iteration process used for the MPB in the BBGY hierarchy [37].The quasi-linearization technique is very efficient and robust, with none of the usual difficulties associated with the Coulomb potential because of the damped nature of the mean electrostatic potential.
For a single electrolyte the MPB structural and thermodynamic properties are on a par with those of the benchmark HNC integral equation when compared with MC simulations, while in some other situations the MPB is more accurate and can be solved for a greater parameter range.Thus in contrast to the HNC equation, the MPB (and SPB) equations can predict a gas-liquid co-existence region for the PM electrolyte, and can be solved for a larger neutral species concentration range in an electrolyte + neutral species mixture [42].As a general principle the SPB theory is most suited to those physical situations involving univalent ions [9].However, in cases where neutral species are present, the SPB can usually provide at least qualitative, if not quantitative, predictions.

Main problems
The two major approximations in the present MPB theory for the PM electrolyte are: (a) the linearization of equation ( 36 Clearly the ideal result is a numerical solution of equation ( 36), using perhaps a finite element technique.For consistency a fully symmetric closure should be used although it is expected that either of the closures ( 24), ( 27) would be adequate.
The linear analysis provides physical insight into the MPB approach.Within the present linear analysis the most important progress would be to (i) treat the intersecting exclusion volume problem and (ii) adequately treat unequal ion sizes.The neglect of (i) leads to incorrect weighting for the short range interaction of two ions, and hence shortfalls when concepts such as ion pairing are important.An adequate treatment of (i) should give the required variation of y C with different valences, and further corrections could be given by including in the fluctuation potential problem the exclusion volume term (30) or the non-linear term for unequal valences [12].The point ion contribution in the screening length for small Γ should come from a solution of the fluctuation potential problem using the more accurate equation (37) rather than equation (38).This is expected as in the corresponding inhomogeneous MPB theory for a planar double layer [48,49], there is a point ion contribution in the singlet distribution function which can readily be shown to contribute to the screening length.The present flawed unequal size analysis compounds the error introduced in (i).Ionic criticality clearly exposes the limitation of (i) and (ii) where, for example, the RPM critical density is too low and for the PM the variation of the critical parameters is not in accord with simulation predictions [50,51].A satisfactory treatment of (i) and (ii) is also important in colloidal applications where large variations in size and valences can occur.
A further interesting problem is finding the MPB closure in terms of the direct correlation function c ij .Given the accuracy of the MPB theory, even with the approximations after the closure, such an insight could well prove useful for approaches based on the Ornstein-Zernike equation.An associated problem is whether or not the MPB closure implies that the second moment condition is satisfied.The present numerical solution indicates that it is nearly satisfied for typical electrolyte solutions, but the agreement diverges in the neighbourhood of the PM co-existence region.Of course, neutrality is automatically satisfied in the MPB formulation through the boundary conditions.A simple linear theory satisfying both neutrality and the second moment condition (65) can be constructed from the first two terms of equation (62) by the appropriate choice of A 1 and A 2 [25].
Emphasis has been placed here on PM electrolytes.Suffice it to say, it is extremely important to consider models which treat the solvent in a realistic fashion.The application of the MPB approach to treat such models has been hampered by some of the present PM limitations.Applications have only been made to the neutral solvent PM [42] and ion-dipole mixtures [52,53].There is clearly great scope for the application of mean electrostatic potential theories, based on the MPB closure, to improved electrolyte models.
) outside the exclusion volumes, (b) the restriction of the solution of the fluctuation problem to large separations of the ions.