Comments on the linear modified Poisson-Boltzmann equation in electrolyte solution theory
C. W. Outhwaite L. B. Bhuiyan


CaDIS: Cataract Dataset for Image Segmentation
    Maria Grammatikopoulou^∗, Evangello Flouty^∗, Abdolrahim Kadkhodamohammadi, Gwenolé Quellec, Andre Chow, Jean Nehme, Imanol Luengo and Danail Stoyanov

This work was supported by the Wellcome/EPSRC Centre for Interventional and Surgical Sciences (WEISS) at UCL (203145Z/16/Z), EPSRC (EP/P012841/1, EP/P027938/1, EP/R004080/1) and the H2020 FET (GA 863146). Danail Stoyanov is supported by a Royal Academy of Engineering Chair in Emerging Technologies (CiET1819/ 2/ 36) and an EPSRC Early Career Research Fellowship (EP/P012841/1). 

M. Grammatikopoulou, E. Flouty, A. Kadkhodamohammadi, A. Chow, J. Nehme, I. Luengo and D. Stoyanov are with Digital Surgery LTD, 230 City Road, EC1V 2QY, London, UK (emails: {maria.grammatikopoulou, evangello.flouty, rahim.mohammadi, andre, jean, imanol.luengo, danail.stoyanov}@touchsurgery.com). G. Quellec is with Inserm, 29200, Brest - France (email: gwenole.quellec@inserm.fr). Corresponding author: Danail Stoyanov.

^∗ M. Grammatikopoulou and E. Flouty contributed equally to this work.

    December 30, 2023
=======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





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ü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.

oncentrated electrolytes, restricted primitive model, screening length,
linear modified Poisson-Boltzmann theory

82.45.Fk, 61.20.Qg, 82.45.Gj




§ 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 <cit.>. By contrast, the non-linear MPB
approach along  with the hypernetted chain (HNC) <cit.>
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  <cit.>. It can be shown that
the MSA is the linearized  version of the HNC theory (see for example,
the references <cit.>). 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 κ <cit.>
and the analogous MSA parameter Γ=[√((1+2κ a))-1]/(2a)
(a is the ionic diameter) <cit.>. 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
<cit.>. 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 <cit.>, and that
of Debye and Hückel <cit.> 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. <cit.>.
The experimental screening length decreases at the higher electrolyte
concentrations in contrast to that predicted by the DH theory. Smith et al.
<cit.>, Lee et al. <cit.> have recently investigated the  experimental
screening length of concentrated electrolytes. A scaling analysis <cit.>
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_ij 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

    ∇ ^2ψ _i(1;2)=-1/ε _0ε_r∑ _se_sn_sg_is ,

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
<cit.>. Taking the ion j to have a charge λ _2e_j  (0 ⩽λ _2⩽  1),
then the charging process gives <cit.>

    g_ij = g_ij (λ_2 = 0)exp{-β e_j[ψ (1;2) + ∫_0^1ϕ (1,2;2)λ _2]},

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 (<ref>) for g_ij is not symmetric, but a symmetric pair distribution function is
obtained by charging up the ion i and combining the result with equation (<ref>) <cit.>.

    The MPB theory rests upon the closure introduced in <cit.>.  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:

    W(1,2;3) = W(1,3) + W(2,3) + w(1,2,3),

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),

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 <cit.>.

	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 <cit.> has derived, for
a RPM,

    g_ij = g_ij^0exp{-β/2[e_j L(u_i)+e_i L(u_j )]},

where g_ij^0 = g_ij (e_i = e_j = 0),

    L(u)= 1/2r(1+y)[u(r+a)+u(r-a)+ κ∫_r-a^r+au(R) R].

Here, u = r ψ(1;2), r = r_ij, y = κ a with
κ = [β/(ε_0ε_r)∑_s e_s^2 n_s]^1/2
the DH parameter. Substituting for g_ij in equation (<ref>) 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:

    ^2u(r)/ r^2 = g_ij^0κ ^2 r L(u),      r ⩾ a,

with

    u(r) = r/au(a)-e_i (r-a)/4ε_0ε_r a ,    0 ⩽ r ⩽ a.

Substituting the linear solution equation (<ref>) into equation (<ref>) means that

    L(u)     =  κ ^2/4(1+y)[2u(r+a)+(r-2a)(2+2y-κ r) ( u/ r)_r = a
    +2(1+2y-κ r)u(a) + 2κ∫_a^r+au(R) R], 
         a ⩽ r ⩽ 2a.

Taking the Laplace transform of equation (<ref>), with g_ij^0 = 1, gives the
general solution to be of the form

    u(r) = ∑_n=1^∞ A_nexp(-z_n r/a),

where the sum is over the roots z_n = α _n + β _n
of the transcendental equation

    z cosh(z) + y sinh(z)= z^3 (1+y)/y^2.


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 <cit.>, with

    ^2u(r)/ r^2 = κ ^2 r L_M(u),     r ⩾ a,

where now

    L_M = 1/2r{(1-b)^2 [u(r+a)+u(r-a)] + 2b/a∫ _r-a^r+a(1- b/a |r-R|)u(R) R},

with b = [y+1-√((1+2y))]/y.  The general solution is the form of
equation (<ref>), where now z satisfies the transcendental equation

    z^2 - 2aΓ z + 2(aΓ )^2 [1-exp(z)] = 0,

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 <cit.>. 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 <cit.>. Figure <ref> displays the behaviour of the MPB and MSA
roots, while a table of some of the roots is given in <cit.>. A critical value of
y_c = 1.032 was first predicted by Kirkwood <cit.> in his analysis of the
potentials of mean force. His transcendental equation, corresponding to those of
equations (<ref>) and (<ref>), 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 <cit.>.



We now consider u(r) to be given by the two roots with the smallest real parts so that

    u(r) = A_1exp(-z_1 r/a) + A_2exp(-z_2 r/a),    r ⩾ 2a,

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 <cit.>.
So, taking the complex conjugate pair z_1, z_2 to be α±β, equation (<ref>)
can then be written as follows:

    u(r) = A exp(-α r/a)cos(β r/a-B),      r ⩾ 2a,

where A_1 = A̅_2 = X +  Y,   A = 2 √(X^2 + Y^2),
tanB = Y/X. Equation (<ref>) 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 <cit.>.
Equation (<ref>) predicts, for both theories, that the screening α /a decreases and
the frequency parameter β /a increases for increasing y.

The approximate solution (<ref>) involves two unknown constants. These cannot be
determined in a unique fashion as equation (<ref>) is simply the two leading terms of
the general solution (<ref>). 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 ( u/ r)_r = a. In this general case, the boundary
conditions at r = a give the neutrality condition

    ∑ _n=1^∞A_n (1 + z_n) exp(-z_n) = e_i/4ε _0ε_r.


Another general condition is that of Stillinger-Lovett <cit.>, which is
true away from the critical region of the RPM electrolyte. In terms of the mean electrostatic
potential, the condition is <cit.>

    β∑_se_s n_s∫ψ _s V = 1.


    Assuming that the solution (<ref>) 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 <cit.>, is for the
solution to satisfy neutrality and equation (<ref>) at r = a. A third possibility,
amongst others, is to evaluate u(r) in a ⩽ r ⩽ 2a
 by substituting the solution (<ref>)
for r ⩾ 2a into equation (<ref>) and integrating twice. In this case, the boundary
conditions are u(r), u(r)/ r 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 <ref> 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 <cit.>. 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.15 K, a = 4.25× 10^-10 m, relative
permittivity ε _r = 78.381 with 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 <ref> and <ref> for a uni-univalent
electrolyte. Figure <ref> illustrates the behaviour of ψ ^* as y increases
through  y_c from y = 1.243 to y = 1.978, and figure <ref> 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 (<ref>). 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 (<ref>). 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 <ref>.
Lee et al. <cit.> 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 <ref>
of <cit.>). 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 <ref>, 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 <cit.>. 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 (<ref>). 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 α±β 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 (<ref>), 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.

10
blum Blum L., In: Theoretical Chemistry, Advances and Perspectives, Vol. 5, Eyring H., Henderson H.
        (Eds.), Academic Press, New York, 1980, 1.
friedman Friedman H.L., A Course in Statistical Mechanics, Prentice Hall, New Jersey, 1985.
vlachy Vlachy V., Annu. Rev. Phys. Chem., 1999, 50, 145, 10.1146/annurev.physchem.50.1.145.
bhuiyan Bhuiyan L.B., Vlachy V., Outhwaite C.W., Int. Rev. Phys. Chem.,
        2002, 21, 1, 10.1080/01442350110078842.
adrielQuiñones A.O., Bhuiyan L.B., Outhwaite C.W., Condens.
        Matter Phys., 2018, 21, 23802,
 10.5488/CMP.21.23802.
henderson1 Henderson D., Blum L., J. Electroanal. Chem. Interfacial Electrochem., 1978, 93, 151,
 10.1016/S0022-0728(78)80228-9.
henderson2 Henderson D., Blum L., J. Chem. Phys., 1979, 70, 3149, 10.1063/1.437813.
debye Debye P., Hückel E., Z. Phys., 1923, 24, 185.
gouy Gouy G., J. Phys. Theor. Appl., 1910, 9, 457, 10.1051/jphystap:019100090045700.
chapman Chapman D.L., Philos. Mag., 1913, 25, 475, 10.1080/14786440408634187.
rotenberg Rotenberg B., Bernard O., Hansen J.-P., J. Phys.: Condens. Matter, 2018, 30, 054005, 
 10.1088/1361-648X/aaa3ac.
smith Smith A.M., Lee A.A., Perkin S., J. Phys. Chem. Lett., 2016, 7,
        2157, 10.1021/acs.jpclett.6b00867.
lee Lee A.A., Perez-Martinez C.S., Smith A.M., Perkin S., Phys. Rev. Lett., 2017, 119, 026002,
 10.1103/PhysRevLett.119.026002.
kirkwood1 Kirkwood J.G., Chem. Rev., 1936, 19, 275, 10.1021/cr60064a007.
outh1 Outhwaite C.W., In: Statistical Mechanics, Vol. 2, Singer K. (Ed.),
    The Chemical Society, London, 1975,  188–255, 10.1039/9781847556936-00188.
outh2 Outhwaite C.W., Condens. Matter Phys., 2004, 7, 719, 10.5488/CMP.7.4.719.
outh3 Outhwaite C.W., Molero M., Bhuiyan L.B., J. Chem. Soc., Faraday Trans., 1991, 87, 3227,

    10.1039/FT9918703227.
outh4 Outhwaite C.W., J. Chem. Phys., 1969, 50, 2277, 10.1063/1.1671378.
kirkwood2 Kirkwood J.G., J. Chem. Phys., 1934, 2, 767, 10.1063/1.1749393.
outh5 Outhwaite C.W., Hutson V.C.L., Mol. Phys., 1975, 29, 1521, 10.1080/00268977500101331.
hafskjold Hafskjold B., Stell G., In: The Liquid State of Matter: Fluids, Simple and Complex,
        Montroll W.L., Lebowitz J.L. (Eds.), North-Holland Publishing Company, Amsterdam, New York, Oxford, 1982, 175.
martynov1 Martynov G.A., Sov. Electrochem., 1965, 1, 285.
martynov2 Martynov G.A., Sov. Electrochem., 1965, 1, 484.
rasaiah Rasaiah J., Chem. Phys. Lett., 1970, 7, 260, 10.1016/0009-2614(70)80303-7.
kjellander Kjellander R., J. Phys. Chem., 1995, 99, 10392, 10.1021/j100025a048.  
lee2 Lee B.P., Fisher M., Europhys. Lett., 1997, 39, 611, 10.1209/epl/i1997-00402-x.
leote Leote de Carvalho R.J.F., Evans R., Mol. Phys., 1994, 83, 619, 10.1080/00268979400101491.
stillinger Stillinger F.H., Lovett R., J. Chem. Phys., 1968, 48, 3858, 10.1063/1.1669709.
outh6 Outhwaite C.W., Chem. Phys. Lett., 1974, 24, 73, 10.1016/0009-2614(74)80216-2.
henderson Henderson D., Smith W.R., J. Stat. Phys., 1978, 19, 191, 10.1007/BF01012511.
outh7 Outhwaite C.W., Mol. Phys., 1971, 20, 705, 10.1080/00268977100100671.





§ ANALYTIC EXPRESSIONS FOR THE THREE LMPB THEORIES





 §.§ 0  ⩽ y ⩽ 1.2412



Roots of equation (<ref>), z_1 = α _1, z_2 = α _2.

    ψ _i = [A_1exp(-α _1r/a)+A_2exp(-α _2r/a)]/r,     r ⩾ a.


(a) LMPB1

    A_1   =  e_i/4ε _0ε _rexp(α _1)[α _1^2G_2-ω (1+α _2)]/D,    
    
    A_2   = e_i/4ε _0ε _rexp(α _2)[-α _2^2G_1+ω (1+α _1)]/D,     
    
    D     =  α _1^2(1+α _1)G_2-α _2^2(1+α _2)G_1 , ω = 6α _2^2/y^2,  
    
    G_j   =  α _j^3+3α _j^2+6α _j+6,   j = 1,2.


(b) LMPB2, see also reference <cit.>.

    A_1   =  e_i/4ε _0ε _rb_2exp(α _1)/D,   
    
    A_2   =  -e_i/4ε _0ε _rb_1exp(α _2)/D ,  
    
    D     =  b_2(1+α _1) - b_1(1+α _2),   
    
    b_j    =  2α _j^2/λ-(2+y)α _j-2(1+y)-2y/α_j-2(1-y/α _j)exp(-α _j),   j = 1,2,               
    λ   =  y^2/2(1+y).


(c) LMPB3. Not given as results very detailed, see for example subsection <ref> below.



 §.§ 1.2412 < y < 7.83



Roots of equation (<ref>), z_1 = α + β,  z_2 = α -β.

    ψ _i   = (e_i/4ε _0ε _rr)Aexp[-α (r/a-1)]cos[β(r/a-1)-B],    r ⩾ a,
    
    A    = √(X^2+Y^2)/D,      B = tan ^-1(Y/X),
    
    R    = α ^2 - β ^2,      S = α ^2 + β ^2.


(a) LMPB1

    G     = G_1 +  G_2 ,   
    
    G_1    = α ^3 -3αβ ^2+3α ^2-3β ^2 +6α +6,  
    
    G_2    = 3α ^2β -β ^3 +6αβ +6β,     
    
    H     = 6/y^2S^2,  
    
    X     = β H-RG_2+2αβ G_1 ,


    Y     = (1+α )H-RG_1+2αβ G_2 ,    
    
    D     = (1+α )X-Yβ.


(b) LMPB2

From <ref> (b), b_1=Y+ X,  b_2=Y- X, giving

    Y     = 2R/λ-(2+y)α -2(1+y)-2yα/S-2exp(-α )(pcosβ +qsinβ ),   
    
    X     = 4αβ/λ-(2+y)β +2q-2exp(-α )(qcosβ  -psinβ  ),    
    
    D     = (1+α )X-Yβ ,    
    λ    = y^2/2(1+y) ,  p=1-yα/S ,   q=qy/S.


(c) LMPB3

    ψ _i    =  λ e_i/4ε _0ε _r[DX+EY+μ (r)+c_8η (r)],  a ⩽ r  ⩽ 2a,                   
    ψ _i    =  e_i/4ε _0ε _r A exp(-α r/a)cos(β r/a -B),  r ⩾ 2a,            
    
    A     =  2√(X^2+Y^2),   B = tan ^-1(Y/X),   
    λ    =  y^2/2(1+y) ,      p = 1-yα/S ,     q = yβ/S ,  
    
    R     =  α ^2 - β ^2,      S = α ^2 + β ^2,     
    
    x     =  r-2a/a ,     
             
    c_1    = 6y^3/12+12y-y^3 ,              
    
    c_2    = 120y^2/240+240y+y^2(60+15c_1+40y+4yc_1) ,              
    
    c_3    = 5c_1/24 ,   
    
    c_4    = 1/24(4c_3+5)+y/120(5c_3+6)+d_2c_8 ,  
    
    c_5    = 1-y^2(3+2c_8)/6-y^4c_4/2(1+y) ,    
    
    c_6    = 19/120+29y/720+c_3/8+yc_3/30 ,   
    
    c_7    = 1-y^2[3/2+2y/3+c_1(4+y)/16+c_8c_9],   
    
    c_8    = -c_224+4c_1+12y+yc_1/24 ,     
    
    c_9    = 200+110y+3c_1(15+4y)/240 ,    
    
    d_1    = p(1+α/S)+q^2/y+3y/2 ,    
    
    d_2    = βy(S+2α )-S/S^2 ,         
    
    D         = K + D_1+ Mx +[c_1V/2+yexp(-2α )B_4/S]x^2 + M_3η (r),


    E         = L+E_1+ Nx +[c_1W/2+yexp(-2α )B_5/S]x^2 + N_3η (r),    
    
    D_1    = 2exp[-(r/a+1)α](Rs_1+2αβ s_2)/S^2 ,      
    
    E_1    = 2exp[-(r/a+1)α](2αβ s_1-Rs_2)/S^2 ,      
    
    s_1    = pcos[β(r/a+1)]+qsin[β(r/a+1)],   
       s_2    = qcos[β(r/a+1)]-psin[β(r/a+1)],   
    η (r)     = 1/240[c_1(15+4y)+120+40(1-y)x-10x^2]x^2,   
    μ (r)     = 1/48[c_1(4+y)+24]x^2-y/6x^3,   
    
    K     = 2/λexp(-2α )cos2β-2/s^2(RB_2+2αβ B_3)exp(-3α ), 
         L     = 2/λexp(-2α )sin2β-2/s^2(2αβ B_2-RB_3)exp(-3α ), 
    
    B_2    = pcos3β +qsin3β,  B_22 = pcos2β +qsin2β,   
    
    B_3    = qcos3β -psin3β,  B_32 = qcos2β -psin2β,   
    
    B_4    = αcos2β - βsin2β ,   
    
    B_5    = βcos2β + αsin2β ,   
    
    B_6    = (1+α/S)cos2β -β/Ssin2β ,   
    
    B_7    = (1+α/S)sin2β +β/Scos2β ,   
    
    B_8    = (2+α/S)cos3β -β/Ssin3β ,   
    
    B_9    = (2+α/S)sin3β +β/Scos3β ,   
    δ _1    = α ^3-3αβ ^2,   δ _2=β ^3-3βα ^2,  
    
    M     = 2/S(α B_2+β B_3)exp(-3α )-2/λ(αcos2β +βsin2β ) exp(-2α ),   
    
    N     = 2/S(β B_2-α B_3)exp(-3α )-2/λ(αsin2β -βcos2β ) exp(-2α ),   
    
    P     = {y/3SB_4+2/S^3[δ _1B_22-δ _2B_32-(δ _1B_2-δ _2B_3)
    exp(-α)]}exp(-2α),   
    
    Q     = {y/3SB_5+2/S^3[-δ _2B_22-δ _1B_32+(δ _2B_2+δ _1B_3)
    exp(-α)]}exp(-2α),   
    
    V     = K -M/2+P,   
    
    W     = L -N/2+Q,   
    
    V_1    = 2/Sexp(-2α)(d_1B_4+d_2B_5),   
    
    W_1    = 2/Sexp(-2α)(d_1B_5-d_2B_4),   
    
    M_2    = 2/S{β [-qcos2β +(y+p)sin2β ]-α [(y+p)cos2β +qsin2β ] }exp(-2α ),


    N_2    = 2/S{α [qcos2β -(y+p)sin2β ]-β [(y+p)cos2β +qsin2β ]}exp(-2α ),   
    
    M_3    = c_2(M-c_1V+M_2),     N_3 = c_2(N-c_1W+N_2),   
    
    K_2    = {5y/12SB_4+2/S^3[δ _1(pB_6+qB_7)-δ _2(qB_6-pB_7)] 
        - exp(-α)[δ _1(pB_8+qB_9)+δ _2(pB_9-qB_8)]}exp(-2α),   
    
    L_2    = {5y/12SB_5+2/S^3[-δ _2(pB_6+qB_7)+δ _1(qB_7-pB_6)]
    
        + exp(-α)[δ _2(pB_8+qB_9)-δ _1(pB_9-qB_8)]}exp(-2α),   
    
    K_3    = 3K/2-2M/3+5c_3V/24+K_2 ,   
    
    L_3    = 3L/2-2N/3+5c_3W/24+L_2 ,   
    
    M_4    = λ(K-2M+3c_1V/2+V_1+c_9M_3) ,  
    
    N_4    = λ(L-2N+3c_1W/2+W_1+c_9N_3),   
    
    F     = 2/Sy^2exp(-2α )[2B_4+1/s(Rcos2β -2αβsin2β )],  
    
    G     = 2/Sy^2exp(-2α )[2B_5+1/s(Rsin2β +2αβcos2β )],  
    
    F_1    = y^2M_3/3+λ y^2(K_3+c_6M_3)+F,   
    
    G_1    = y^2N_3/3+λ y^2(L_3+c_6N_3)+G,   
    
    X         = N_4c_5-G_1c_7/D_0 ,    
    
    Y         = F_1c_7-M_4c_5/D_0 ,    
    
    D_0    = N_4F_1-G_1M_4.












	label1  , i i, i, i
	label2i  i, i  i,  i, 





       
	  -   ' .
	ֳ       .
	˳     
	      -  
     
	  .     
	,      .
	       
	      NaCl.
	

	
	   ,  ,    -






