Comments on the linear modified Poisson-Boltzmann equation in electrolyte solution theory

Three analytic results are proposed for a linear form of the modified Poisson-Boltzmann equation in the theory of bulk electrolytes. Comparison is also made with the mean spherical approximation results. The linear theories predict a transition of the mean electrostatic potential from a Debye-H\"{u}ckel type damped exponential to a damped oscillatory behaviour as the electrolyte concentration increases beyond a critical value. The screening length decreases with increasing concentration when the mean electrostatic potential is damped oscillatory. A comparison is made with one set of recent experimental screening results for aqueous NaCl electrolytes.


Introduction
Little attention has been paid to the linear modified Poisson-Boltzmann (MPB) equation in the electrolyte solution theory, mainly due to the ready availability of the mean spherical approximation (MSA) analytical results [1].By contrast, the non-linear MPB approach along with the hypernetted chain (HNC) [2] theory have proved two of the more successful theories in predicting the structure and thermodynamic properties of the primitive model (PM) (arbitrary sized charged hard spheres moving in a continuum dielectric) electrolyte solution [3][4][5].It can be shown that the MSA is the linearized version of the HNC theory (see for example, the references [6,7]).However, unlike the MSA, the linear MPB theory remains relatively less well explored.
In a broad sense, the usefulness of a linear theory is twofold.
Firstly, from a theoretical point of view, linear theories can and often provide valuable insights into the physics and chemistry of a situation, which might otherwise remain obscure when analyzed through their non-linear counterparts only.A case in point is the similarity (or even equivalence) between the Debye-Hückel (DH) parameter κ [8] and the analogous MSA parameter Γ = [ (1 + 2κa) − 1]/(2a) (a is the ionic diameter) [1].Although the HNC gives a more accurate solution generally, it does not lead to a physically significant quantity such as Γ.Similarly, in our case, the linear MPB (LMPB) predicts the asymptotic form of the non-linear equation and gives useful information regarding the screening length.Being based on the mean electrostatic potential approach, the linear theory is a logical extension of the Debye-Hückel (DH) theory [8].In particular, a natural transition is given of the mean electrostatic potential behaviour in going from a DH damped exponential to a damped oscillatory as the concentration increases beyond some critical value.
Secondly, from a practical perspective, linear solutions are normally easier to obtain and are often analytic.The latter feature makes a linear solution particularly convenient to use as an initial guess in an iterative algorithm to obtain the corresponding non-linear solution.
Screening has played an important role in analysing charged systems since the electrical double layer work of Gouy and Chapman [9,10], and that of Debye and Hückel [8] in the electrolyte bulk.The phenomenon is significant in shaping the structure and thermodynamics in these systems.A recent theoretical work on the subject is due to Rotenberg et al. [11].The experimental screening length decreases at the higher electrolyte concentrations in contrast to that predicted by the DH theory.Smith et al. [12], Lee et al. [13] have recently investigated the experimental screening length of concentrated electrolytes.A scaling analysis [13] suggests that the decay length increases linearly with the Bjerrum length at higher concentrations.When applied to the restricted primitive model (RPM), the MSA and LMPB theories both predict a decrease in the screening length above a critical y c (y = κa).Neither theory can predict the scaling result, a probable factor being the neglect of solvent effects.
In this paper we will present the mean electrostatic potential results due to three different versions of the LMPB, and compare them with the corresponding results from the MSA, the MPB, and the DH theories.The comparative behaviour of the screening properties predicted by the MPB and MSA theories will also be examined.

MPB theory
Improvements to the classical PB theory rest upon a more accurate pair distribution function g i j for two ions i and j.In the mean electrostatic potential approach, this implies an improvement in the mean electrostatic potential ψ(1; 2) at the field point r 2 for an ion i at r 1 through the solution of the Poisson equation where e s is the charge on an ion s at r 2 , n s is the mean number density of ions of type s, ε r is the solvent relative permittivity and ε 0 is the vacuum permittivity.One technique to express the pair distribution function in terms of mean electrostatic potentials is to use Kirkwood's charging process [14].Taking the ion j to have a charge λ 2 e j (0 λ 2 1), then the charging process gives [15,16] where φ(1, 2; 3) is the fluctuation potential defined by ψ(1, 2; 3) = ψ(1; 3) + ψ(2; 3) + φ(1, 2; 3), with ψ(1, 2; 3) the mean electrostatic potential at r 3 for ions i and j fixed at r 1 and r 2 , respectively.Thus, the fluctuation potential φ(1, 2; 3) describes the departure from linear superposition of the singlet potentials.Equation (2) for g i j is not symmetric, but a symmetric pair distribution function is obtained by charging up the ion i and combining the result with equation ( 2) [17].
The MPB theory rests upon the closure introduced in [18].Consider the conditional potential of mean force W(1, 2; 3), which measures the work done in bringing an ion s to r 3 , given the fixed ions i and j at r 1 and r 2 , respectively.In an analogous fashion to that for the mean electrostatic potentials, it can be written as follows: where w(1, 2, 3) is the corresponding departure from linear superposition of the singlet potentials of mean force.The MPB closure is w(1, 2, 3) = e s φ(1, 2; 3), (5) which is analogous, but at the next hierarchical level, to the DH closure W(1; 3) = e s ψ(1; 3).Attempts to improve the DH theory by using the closure W(1, 2; 3) = e s ψ(1, 2; 3) fails as this closure neglects terms of the order being calculated [15,16,18,19].
The governing set of differential equations for the fluctuation potential can be derived by using the Poisson equation for one or two fixed ions.This gives for the PM a two sphere potential problem which enables an approximate solution to be found for φ (1, 2; 3).Given an appropriate fluctuation potential, an improved pair distribution function can now be found.Previous analysis [15][16][17] has derived, for a RPM, where g 0 i j = g i j (e i = e j = 0), Here, u = rψ(1; 2), r = r i j , y = κa with κ = [β/(ε 0 ε r ) s e 2 s n s ] 1/2 the DH parameter.Substituting for g i j in equation ( 1) gives the RPM non-linear MPB equation.

The linear MPB and MSA equations
The linear form of the MPB equation, for the special case of a single electrolyte with equal valences, is as follows: with Substituting the linear solution equation ( 9) into equation ( 8) means that Taking the Laplace transform of equation ( 8), with g 0 i j = 1, gives the general solution to be of the form where the sum is over the roots z n = α n + iβ n of the transcendental equation The roots with the smallest real part correspond to the physical situation.For small y, there are only two real roots with the smallest corresponding to the DH solution.As y increases, the two real roots coalesce at y c = 1.2412.Above this value, the two roots move off as a complex conjugate pair until they become purely imaginary at y I = 7.83.
The MSA theory is one of the most successful and widely used theories of charged fluid systems.Although of a linear nature, its appeal lies in its analytical results, which have been applied to numerous electrolyte models and theories.The mean electrostatic potential has been studied in [20], with where now with b = [y + 1 − (1 + 2y)]/y.The general solution is the form of equation (11), where now z satisfies the transcendental equation with aΓ = y(1 − b)/2.There are two real roots for small y with the smallest corresponding to the DH solution, analogous to the MPB theory.Again, as y increases, the two real roots coalesce, but at the slightly smaller value of y c = 1.229 [20,21].Above this critical value of y, the two roots become a complex conjugate pair, but unlike the MPB situation, they do not become purely imaginary for large y [1]. Figure 1 displays the behaviour of the MPB and MSA roots, while a table of some of the roots is given in [15].A critical value of y c = 1.032 was first predicted by Kirkwood [14] in his analysis of the potentials of mean force.His transcendental equation, corresponding to those of equations ( 12) and ( 15), has a behaviour pattern similar to the MPB theory.Unfortunately, his value of y I = 2.79 means that his theory is restricted to much lower concentrations than the MPB and MSA theories.Critical values of y c have been predicted by other theories [22][23][24][25][26].
We now consider u(r) to be given by the two roots with the smallest real parts so that with the first term corresponding to the DH value for small y.Both the MPB and MSA predict, above their respective y c , that the solution is damped oscillatory [15,16].So, taking the complex conjugate pair z 1 , z 2 to be α ± iβ, equation ( 16) can then be written as follows: where 17) holds for the MPB when 1.2412 < y < 7.83, while for the MSA when y > 1.229.Analysis of the Ornstein-Zernike equation indicates that the asymptotic form of the correlation function takes on a similar form [27]. Equation (17) predicts, for both theories, that the screening α/a decreases and the frequency parameter β/a increases for increasing y.
The approximate solution (16) involves two unknown constants.These cannot be determined in a unique fashion as equation ( 16) is simply the two leading terms of the general solution (11).A consistent analysis requires that the general solution satisfies the integro-difference differential equation in a r 2a and the continuity of u(a) and (du/dr) r=a .In this general case, the boundary conditions at r = a give the neutrality condition Another general condition is that of Stillinger-Lovett [28], which is true away from the critical region of the RPM electrolyte.In terms of the mean electrostatic potential, the condition is [29] Assuming that the solution ( 16) holds for r a, then neutrality and the Stillinger-Lovett condition allows A 1 and A 2 , or equivalently A and B, to be calculated.Another approach, following Kirkwood [14], is for the solution to satisfy neutrality and equation ( 16) at r = a.A third possibility, amongst others, is to evaluate u(r) in a r 2a by substituting the solution ( 16) for r 2a into equation ( 8) and integrating twice.In this case, the boundary conditions are u(r), du(r)/dr continuous at r = a and 2a, with neutrality and the Stillinger-Lovett condition satisfied.The three approximations are called LMPB1, LMPB2, and LMPB3, respectively.See appendix A for details.The above approximate analysis to determine the LMPB mean electrostatic potential is unnecessary for the MSA as an alternative approach gives an analytical result [20,30].However, this alternate formulation does not give immediately a direct identification of the MSA screening and frequency parameters.

Results and discussion
The physical parameters used in all the calculations were as follows: the temperature T = 298.15K, a = 4.25 × 10 −10 m, relative permittivity ε r = 78.381with the concentration varying from c = 0.79 M to 18.4 M, so the range of y is 1.243 to 6. Plots of the dimensionless mean electrostatic potential ψ * (r/a) [= β|e|ψ(r/a)] for the 6 theories MPB, LMPB1, LMPB2, LMPB3, MSA and DH are shown in figures 2 and 3 for a uni-univalent electrolyte.Figure 2 illustrates the behaviour of ψ * as y increases through y c from y = 1.243 to y = 1.978, and figure 3 shows the behaviour at the higher values of y = 3, 4, 6.The non-linear MPB has been shown to accurately predict both the structural and thermodynamic properties of the RPM for 1:1 electrolytes, and thus is a good metric to assess the accuracy of the linear theories.For y < y c , except for y approaching y c , all the theories are qualitatively very similar to the DH result, and hence, are not shown here.Above y c , the theories demonstrate a damped oscillatory behaviour which cannot be predicted by the DH theory.Overall LMPB2 is in closest agreement with MPB, while surprisingly LMPB3 does relatively poorly.It would have been expected that the LMPB3 would be the most accurate because of the treatment of the region a r 2a via equation (10).Comparison with the MPB, when the exclusion volume term is unity, brings the LMPB3 and MPB into almost quantitative agreement for c 18.4 M.This indicates that the improved mean electrostatic potential treatment of the LMPB3 analysis needs to be balanced by incorporation of the uncharged hard sphere term in equation ( 8).An interesting feature of the LMPB3 is that besides damped oscillation terms, there is a (r − 2a) quadratic contribution in the region a r 2a.
Experimental results for a wide class of salts indicate that the decay length λ s of concentrated electrolytes increases with concentration, in contrast to that of the DH prediction.The LMPB and MSA also predict that their decay length a/α, which corresponds to the experimental λ s , increases at higher concentrations, but at too low a rate -see figure 4. Lee et al. [13] have analysed the electrolyte screening length behaviour at high concentration using a scaling analysis.They find that the decay length increases linearly with the Bjerrum length, or alternatively λ s /λ D ∼ (a/λ D ) 3 where λ D = 1/κ, in general agreement with experimental results for a range of 1:1 valency ionic liquids.We note that in the present paper λ s /λ D = y/α and a/λ D = y so that the experimental trend thus provides a test of theories through a study of y/α ∼ y 3 .Unfortunately, no realistic comparisons can be made with the MPB and MSA predictions due to the decay length being too small.For example, at y = 3, the LMPB y/α = 1.95, which is an order of magnitude too small (see figure 1 of [13]).Only for y close to y I is the LMPB decay length of the correct order of magnitude.The most probable reason for the shortfall is the lack of any treatment of the solvent.At low to medium solution concentration, the solvent contribution can be partially treated by using an hydrated ion diameter.In figure 4, the experimental ratio λ s /λ D for varying a/λ D is given for an aqueous solution of NaCl for a = 2.94 × 10 −10 m and a = 5.2 × 10 −10 m.The use of the hydrated diameter a = 5.2 × 10 −10 m brings the experimental points closer to the MPB results.The high concentration of the solvent relative to the solute ions means that the structure around the ions is driven by the solvent.A step in this direction is to consider the solvent primitive model where the solvent is modelled by uncharged hard spheres moving in a constant permittivity.Recently, the underscreening in this solvent primitive model has been studied in the MSA [11].The ratio λ s /λ D increases at higher concentrations as required by experiment, but at too slow a rate.

Summary
The potential approach of the MPB theory provides a natural extension to the DH theory.At lower electrolyte concentrations, the predictions mimic the DH theory but qualitative differences occur at higher concentrations.In the region y c < y < y I , the MPB predicts a damped oscillatory potential which is in accordance with simulation and other theoretical work.For this interval of y, the linear potential is of the form u(r) = A exp (−αr/a) cos (βr/a − B), r 2a, where A, B are constants and α, β given by the solution of equation (12).There is no unique way of determining the constants, so three methods are considered, in each case the neutrality condition being satisfied.The Stillinger-Lovett condition is used in both the LMPB1 and LMPB3 approaches with the solution holding at r = a and r = 2a, respectively, while the LMPB2 is derived by assuming that the equation is satisfied at r = a.Comparison of the three theories with the accurate non-linear MPB for 1:1 salts indicates that overall the simple LMPB2 provides the best approximation.The relative poor behaviour of LMPB3 is deceptive.The LMPB3 potential is nearly identical to that of the non-linear MPB when its exclusion volume term is unity.This means that realistic improvements to the linear MPB lie through treating the exclusion volume term in the LMPB3 approach.
Experimental results have shown that the screening parameter, as interpreted via a DH type potential, initially increases as the electrolyte concentration rises, but then decreases at higher concentrations.In the LMPB theory at lower concentrations, there are two pure damped exponential terms with screening parameters α 1 /a, α 2 /a with the smaller α 1 /a corresponding to the DH screening.The parameters α 1 , α 2 coalesce at y c as the concentration increases and become complex conjugate roots α ± iβ for y c < y < y I .These complex conjugate roots lead to a damped oscillatory behaviour for the electrostatic potential with screening parameter α/a.In this region of y, the screening parameter reduces with concentration, as seen in experiments over a similar range of y.The experimental results also suggest that the decay length increases linearly with the Bjerrum length.Such a behaviour is not seen in the LMPB as the rate of a decrease of the screening length is too small, only at values of y close to y I is the screening length of the correct order of magnitude.The MSA mean electrostatic potential has features in common with the LMPB.It becomes damped oscillatory at a value of y c slightly less than that of the MPB, but in contrast there is no finite value of y I .This means that its screening factor behaves in an analogous fashion to that of the MPB as y increases through y c .Above y c , the MSA screening parameter is greater than that of the MPB and slowly reduces to zero as y tends to infinity.Close agreement with experiment at the higher electrolyte concentrations is unlikely with the present PM model.A critical feature of the PM is the lack of any adequate treatment of the solvent.Models such as the solvent primitive model have provided insight as to the importance of solvent steric effects, but both the solvent and solute molecules need to be adequately modelled and treated on an equal footing.
The linear MPB equations studied here complement the full, non-linear MPB approach to the electrolyte solution theory.From a theoretical perspective, these solutions offer a valuable physical insight by painting a clearer picture of the onset of oscillations in the mean electrostatic profile and screening in electrolytes.The LMPB solutions transcend the simplistic DH view, compare favourably with the more well known MSA while retaining the simplicity and convenience of their analytic nature.Thus, the LMPB theories might well prove useful as a basis for interpreting the thermodynamic properties of concentrated electrolytes.Again, the theories could be applied to studying the liquid-gas transition in simple electrolytes in an analogous fashion to the MSA.Although the LMPB1 and LMPB3 satisfy the Stillinger-Lovett condition, this does not prohibit the LMPB1 and LMPB3 from having a solution in the liquid-gas region, as seen with the MSA.The solution of the LMPB theories is basically dependent on y through the soluton of equation ( 12), so the theories can be solved in the simulated liquid-gas coexistence region.The pattern of the LMPB results for the range of concentration treated suggests these to be potentially very useful as initial input in the numerical solution of the full MPB equation and alternative non-linear theories.

Figure 1 .
Figure 1.(Colour online) Roots of the MPB transcendental equation (12).For y < y c there are two real roots α 1 , α 2 with the black (lower) curve corresponding to the DH κ (= α 1 /a) for low y.For y c < y < y I , the black curve represents the screening parameter α, and the blue curve represents the screening frequency | β|.The MSA roots of equation (15) are nearly identical to those of the MPB for y < y c , where now y c = 1.229, and are hence not shown.For y > y c , the MSA screening parameter α is given by the upper green curve.The MSA does not have a finite y I .The two vertical dashed lines at y c = 1.2412 and y I = 7.831 indicate the two MPB critical values.

Figure 4 .
Figure 4. (Colour online) The experimental ratio λ s /λ D for an aqueous solution of NaCl as functions of a/λ D compared with the MPB y/α 1 , y/α 2 for y < y c , y/α for y c y y I .The filled squares and the filled circles are the experimental results for the unhydrated (a = 2.94 × 10 −10 m) and hydrated (a = 5.2 × 10 −10 m) NaCl, respectively.The experimental results are taken from the Supporting Information of Smith A.M., Lee A.A., Perkin S.; electronic link given in reference [13].The upper curve for y < y c is that for α 1 which would be unity if it was the DH value.The MSA screening length is given by the lower red curve and is indistinguishable from the MPB value on the graphical scale for y 4. The MSA does not have a finite y I .The two vertical dashed lines at y c = 1.2412 and y I = 7.831 indicate the two MPB critical values.