Structure and thermodynamics in the linear modified Poisson-Boltzmann theories in restricted primitive model electrolytes

Structure and thermodynamics in restricted primitive model electrolytes are examined using three recently developed versions of a linear form of the modified Poisson-Boltzmann equation. Analytical expressions for the osmotic coefficient and the electrical part of the mean activity coefficient are obtained and the results for the osmotic and the mean activity coefficients are compared with that from the more established mean spherical approximation, symmetric Poisson-Boltzmann, modified Poisson-Boltzmann theories, and available Monte Carlo simulation results. The linear theories predict the thermodynamics to a remarkable degree of accuracy relative to the simulations and are consistent with the mean spherical approximation and modified Poisson-Boltzmann results. The predicted structure in the form of the radial distribution functions and the mean electrostatic potential also compare well with the corresponding results from the formal theories. The excess internal energy and the electrical part of the mean activity coefficient are shown to be identical analytically for the mean spherical approximation and the linear modified Poisson-Boltzmann theories.


Introduction
One of the more enduring theories in the physics and chemistry of Coulomb fluids over the past (nearly) hundred years has been the theory of Debye and Hückel (DH) [1], which is the linearized form of the classical Poisson-Boltzmann (PB) theory. The intuitive simplicity of the DH concept together with the ease of its implementation have been the theory's main attractions. For instance, almost all variables required for a structural and thermodynamic description of an electrolyte solution occur in closed forms in the DH and the Debye-Hückel Limiting Law (DHLL) theory [2]. Formal statistical mechanical analysis (see for example, [3]) and subsequent machine simulations [4][5][6][7][8][9] over the years have brought out the deficiencies of the DH, the principal ones being the neglect of the ionic exclusion volume and the ionic correlation terms. Some of the more recent, salient references, and reviews are given by [10][11][12][13][14][15].
The potential approach to the theory with its origins in the DH mechanism has evolved over the decades through the pioneering work of Kirkwood in the 1930s [3] and later through the works of other authors [16][17][18][19][20][21][22] to the modified Poisson-Boltzmann (MPB) equations of today (see for example, reference [22]). A popular alternate route is based on the liquid structure integral equations such as the hypernetted chain (HNC) [23,24] and the mean spherical approximation (MSA) [25][26][27]. The density functional theory (DFT) has also been explored [28,29].
A widely used physical model used in conjunction with the above statistical mechanical theories in studies of electrolytes has been the primitive model (PM), viz., arbitrary sized charged rigid spheres moving in a dielectric continuum [10]. The solvent is thus structureless being characterized by a dielectric constant or relative permittivity . If the ion sizes are equal, then we have the restricted primitive model (RPM). The RPM is also the underlying model of the DH theory, since a RPM with the vanishing ion radius, except perhaps for the size of the central ion, would lead to the latter. The HNC and the MPB have been two of the most successful theories of PM or RPM electrolytes in the electrolyte solution regime having been applied to a wide variety of situations under different physical conditions. The MPB formalism set in the PM (or the RPM) yields a highly non-linear differential equation whose solution requires involved numerical techniques [18] and thus may not be readily available. Fortunately, as Outhwaite [30] showed a linear version of the MPB equation is tractable analytically leading to a closed form expression for the mean electrostatic potential .
Since the analysis of Outhwaite stated above, little has been reported in the literature on linear theories based on the MPB, although a lot of work has been done with the MSA, viz., the works by Blum (see for example, references [25,27]), and by Outhwaite and Hutson [26], apart from the obvious DH. In addition to being easier to use, linear theories can offer valuable insights and understanding of the properties of systems being examined albeit at the cost of a little accuracy. Linear solutions can also be valuable in iterative numerical solution of corresponding non-linear equations. In a recent paper, Outhwaite and Bhuiyan [31] studied the linear MPB (LMPB) in some detail and formulated three versions of the equation, LMPB , the index i (i = 1,2,3) referring to special characteristics of a particular equation. Significantly, these equations yielded an analytical solution for the , which, for 1:1 valency RPM electrolytes, compared well with that from the MSA and MPB for a range of concentrations, and symmetric Poisson-Boltzmann (SPB) theory [32] at low concentrations. Linearization retained the aspects of the MPB fluctuation potential terms since the linear 's also showed damped oscillations at higher concentrations beyond the critical (= ) = 1.2412 ( being the Debye-Hückel constant and the common ionic diameter). In simple terms, fluctuation potential is the potential formulation of the inter-ionic correlations and such oscillations are signatures for inter-ionic correlations, which are not seen with the mean-field DH or the SPB.
In view of the comparative behaviour of LMPB with that from the SPB, MPB, and MSA theories for a broad range of concentrations of RPM electrolytes seen in [31], we thought it of interest to apply the LMPB approach to an analysis of the structure and thermodynamics in these systems. As we will see later, the LMPB expressions for thermodynamic quantities such as the osmotic coefficient and the electrical contribution to the mean activity coefficient (el) ± develop into closed analytical forms, which make their numerical evaluation straightforward. Addition of the hard-core component (HS) ± to the (el) ± leads to the (full) ± . It is worth mentioning here that the knowledge of and ± has practical significance in many chemical processes involving electrolyte solutions in industry and bio-sciences (see for example, [33,34]). Relevant to this work, we note also that recently Quiñones et al. [35] made extensive comparisons of the SPB and MPB and ± with the corresponding RPM or PM Monte Carlo (MC) data of Abbas et al. [36,37] for a wide range of solution concentrations. The MPB, and to a lesser extent, the SPB showed a very good agreement with the simulations. It would be interesting to see how well the LMPB predictions compare with these results. Being experimentally measurable [38], the measured and ± can also provide a good assessment of theories.
The organization of this paper is as follows. In the following section we outline the principal equations of the LMPB theories pertinent to the calculation of structure and thermodynamics of electrolytes. In section 3 we present and discuss the results of this work, while in section 4 some general conclusions are drawn.

Model
The model electrolyte system employed in this work is an aqueous RPM electrolyte at around room temperature. This is consistent with one of the models used by Abbas et al. [36,37] in their MC simulations, the other being the PM.  In the theoretical calculations, the common ionic radius is Na + = F − = 1.435 ×10 −10 m, taken from the MC simulation data of references [36,37]. The MC data are from the same references.
The various particle-particle interaction potentials are Here, and are the valency and radius of ion species , while is the separation between two ions of types and , respectively. Parameter 0 is the vacuum permittivity and | | is the magnitude of the electronic charge. For the RPM we have = /2 for all .

Methods
The linearization of the MPB equation and the subsequent development of the LMPB equations for the RPM electrolyte have been discussed in details in [31] and will not be repeated here. We restrict ourselves to outlining the salient equations in these theories.
The Poisson equation for the (1; 2) at the field point r 2 in presence of an ion at r 1 is Here, (1, 2) is the radial distribution function for the ion pair and separated by (= ) = |r 1 − r 2 |, with being the mean number density of ions of type . The operator ∇ operates on the coordinates of 23801-3 the field point 2. It is convenient to use the transformation = (1; 2), whence the above equation transforms to In the MPB approximation, the has been developed for the RPM as (see for example, reference [22]) The quantity 0 above is the exclusion volume term, which is the radial distribution function for two discharged ions in a sea of fully charged ions, viz., 0 = ( = = 0), = 1/( B ), with B the Boltzmann's constant and the temperature. The operator ( ) is given by (3)] above yields the following

Linearization of the non-linear MPB equation [equation (4) substituted in
Within the charge free space 0 < < , the solution of the Laplace equation d 2 ( ) = 0 is as follows: Equation (6) is valid for both symmetric and asymmetric valency systems. For instance, for the linear theory with equal ion sizes we have = , so that in using equation (4) in equation (3) with the expression (5) for ( ), we have the terms such as ( ), which can be written as ( ). Hence, the equation follows.
The general solution of the linear MPB equation is governed by the roots of a transcendental equation (cf. reference [31]) Taking the first two roots with the smallest real part, as these give the physical solution, we can see that for small there are two real roots that coalesce at the critical = 1.2412 for symmetric valencies. For > the roots form a complex conjugate pair becoming imaginary at a second critical point = 7.83. Apart from the continuity of and d /d at boundaries such as = , there exists the exact condition, viz., the local electroneutrality Another useful relation is the Stillinger-Lovett (SL) second moment condition [39], although it does not hold near the critical point of the electrolyte. The SL condition can be written as [40] | | In [31], the solutions to the LMPB equation were classified depending on the boundary conditions and the exact conditions a solution satisfies. In the theoretical calculations, the common ionic radius is Na + = Cl − = 1.745 ×10 −10 m, taken from the MC simulation data of references [36,37]. The MC data are from the same references.

LMPB1, LMPB2 and LMPB3 equations
The LMPB1 equation satisfies the electroneutrality and the SL conditions, while the LMPB2 satisfies the neutrality and the continuity of ( ) at = . In addition to the electroneutrality, and the SL conditions, the LMPB3 satisfies the continuity of ( ) and d /d at both = and = 2 . This occurs since a more accurate solution for ( ) can be derived in the region 2 , for example, by using the linear solution (7) in equation (5) to obtain in 2 . We refer the reader to [31] for further details. For , the LMPB1 and LMPB2 solutions can be written as where 1 , 2 are the two real roots of the transcendental equation (8). The constants A 1 , A 2 take on different forms for LMPB1 and LMPB2. For instance, for LPMB1 we have while for LMPB2 we have The LMPB3 solution for this range of has not been given in [31]. For the range < < , again the LMPB1 and LMPB2 solutions have the common form where and are now the real and imaginary parts of the complex conjugate pair of roots of equation (8). Furthermore, =   [36,37]. The MC data are from the same references.
For the LMPB2 we have from equation (12) For the LMPB3 and for the range < < , the solution is where = 2 √︁ ( 2 + 2 ) and = tan −1 ( / ). The quantities , , , , ( ) and ( ) are very involved and we refer the reader to [31] for details.

Structure and thermodynamics
For the LMPB formulations, the pair distribution can be constructed as which gives a non-linear . Another possibility would be to use a linearized version of the above, viz., The expression (18) is analogous to the DHX and EXP theories [2,4,18,41,42], while the expression (19) is analogous to the MSA [26]. To avoid confusion, in the rest of this paper we refer to this linear form as LMPBi (linear) ( = 1, 2, 3). A consistency check on the 's may be carried out through the integral We note that differentiating this twice with respect to immediately yields the requisite Poisson's equation for¯. Hence, the extent of agreement between and¯for an LMPB theory would be a measure of the consistency of the particular .
The osmotic coefficient can be calculated from the relation (see for example, reference [18]) where = and (ex) is the excess internal energy, viz., Calculation of the activity coefficient is most conveniently achieved through the Günteberg charging process where the ion at the origin is charged up from zero to its full charge in a sea of charged ions. The individual ionic activity has been derived as [18,20] ln = ln (HS) + ln (el) .  The hard sphere part ln (HS) was analyzed by Ebeling and Scherwinski [43], while the charging process gives for the electrical part [20] ln Using now equations (7), (24) with the LMPB solutions (11), (14) and (17) it is straightforward to calculate the electrical part of the individual activity coefficient for these theories. For example, for the LMPB1 and LMPB2 we have where Γ = (| | 2 )/(4π 0 ), with the rest of the constants being given in equations (12), (15) for LMPB1 and equations (13), (16) for LMPB2.

Equivalence of ln
(el) In the course of these calculations we have observed an equivalence of the excess internal energy with the natural logarithm of the electrical component of the mean activity coefficient for the LMPB and the MSA theories, viz., where = . Here, we outline an argument why this is expected for a linear theory. Such an equivalence has recently been observed by Kjellander with regard to his Multiple-Decay Extended Debye-Hückel MDE-DH theory of electrolytes [44].
From equation (7) we have Using the Gauss's law for the electric field −(d /d ) = , and assuming that ( ) is linear in (= | | ), equation (29) can be written as where is independent of . Combining now equations (22) and (31) we get for a single electrolyte. In writing the above we have invoked the global neutrality = 0. Now, from equations (24) and (31), we immediately have Hence, using equation (27) the mean activity can be written as Equation (28) now follows upon eliminating between the equations (32) and (34).

Results and discussion
In presenting the results we include data from the SPB, MPB, and the MSA theories for comparison purposes. While the MSA results were obtained from the relevant analytic expressions (see, for example, references [27,45,46]), the non-linear SPB and MPB equations were solved numerically using a quasilinearization technique [47] used successfully in earlier works [19][20][21]32]. The hard sphere part 0 was approximated by the Percus-Yevick uncharged pair distributions [48,49] and their corrections due to Verlet and Weis [50]. In evaluating the 0 we used an efficient numerical technique developed by Perram [51]. Similarly, the hard sphere individual activity coefficient ln (HS) was determined using the formulations of Ebeling and Scherwinski [43].
Except for one case of asymmetric 2:1 valency system, all calculations reported here are for 1:1 symmetric valency systems at temperature = 298 K and relative permittivity = 78.38, which corresponds to a water-like solvent. These values are in line with that used in the MC simulations of Abbas et al. [36,37]. The other physical parameters like the common ionic diameter and the concentration were variable and were fitted to the MC system being compared to.

Thermodynamics
We begin this discussion by considering the results for and ln ± in four 1:1 RPM salt solutions, viz., NaF, NaCl, HCl, and LiI shown in figures 1, 2, 3, and 4, respectively. Abbas et al. [36,37] actually simulated over 100 PM and RPM salts with different valencies 1:1, 2:1, and 3:1 and covering a wide range of ionic sizes and solution concentrations. We have chosen these four since these encompass a fair range of ionic size starting from = 1.435 × 10 −10 m for NaF (figure 1) to = 2.325 × 10 −10 m for LiI (figure 4). We have actually carried out calculations for more 1:1 salts at their MC parameters, where the results show similar characteristics and are hence not shown here for brevity. The four sets of results being displayed also constitute a good representative sample.
A striking feature in figures 1-4 is the remarkable consistency of the LMPB curves both among themselves and with the MC data. Indeed, the results of all the theories including that from the MSA and the MPB are in very close agreement with each other with only the SPB ln ± showing some deviation at the two lower diameters (figures 1 and 2). This behaviour pattern carries over to the 2:1 valency (MgCl 2 ) situation in figure 5 where again the LMPB theories are at par with the formal theories in reproducing the MC data.

Structure
The structural results are presented in figures 6-10. For these calculations, we have used a fixed value of = 4.25 ×10 −10 m. Figure 6 shows the LMPB1, LMPB2, SPB, MPB, and MSA results for the reduced mean electrostatic potential * (= | | = | | ) / ) (upper panel) and the (lower panel) at the electrolyte concentration = 0.01 mol/dm 3 ( = 0.13985). For the * , the various curves are indistinguishable from each other, which is also the case with the 's except for the MSA curves. This occurs owing to the linear nature of the MSA. The contact values of the MSA 's are underestimated and generally the MSA curves lie below the others. The unphysical negative values of the MSA coion around the contact are noticeable. This is of course a known shortcoming of the MSA at low concentrations. In the lower panel, we have also plotted the corresponding LMPB1(linear) and LMPB2(linear) 's. It is interesting, although perhaps not surprising, to note that these linear 's follow the MSA very closely leading to the coion 's also becoming negative at and near the contact.
The results at the higher concentrations of = 1 mol/dm 3 ( = 1.3985) and = 4 mol/dm 3 ( = 2.797) are displayed in figures 7 and 8, respectively. In both situations > , and we see the beginnings of oscillations in the LMPB, MPB, and MSA profiles in figure 7, which, expectedly, become more pronounced in figure 8. Although not quite as quantitative as they are at the lower = 0.01 mol/dm 3 in figure 6, the LMPB predictions continue to show good qualitative agreement overall with that from the MPB at these enhanced concentrations. Nonetheless, the differences in the structure can produce discrepancies in various properties such as thermodynamics via different routes and transport properties, especially at higher concentrations and valencies.
The oscillations in the LMPB curves are due to the fact that linearization of the MPB equation retains the aspects of the fluctuation potential terms. By contrast, no such terms occur in the classical SPB and hence no oscillations.
A consistency check on the LMPB 's and a comparison of these 's with their linear version are shown in figures 9 and 10 at = 1 mol/dm 3 ( = 1.3985) and = 4 mol/dm 3 ( = 1.3985), respectively. For illustrative purposes we only show the results for LMPB2 since the other LMPB1 and LMPB3 results are very similar. In the upper panels, the LMPB2 * from equations (14) and (20) are plotted. The consistency of these * at both of these concentrations is noteworthy, and points, in turn, to the consistency of the approximation made in equation (18). A similar effect is observed with regard to the non-linear LMPB2 and the linear LMPB2(linear) 's, viz., equations (18) and (19), in the lower panels with the linear 's showing some discrepancy only near the contact.

Conclusion
This work represents a continuation of our earlier study, [31], on the applicability of linear modified Poisson-Boltzmann theory to the electrolyte solution theory. The main achievement of this paper is the characterization of thermodynamics of RPM electrolytes using the set of three linear modified Poisson-Boltzmann theories proposed in [31]. The osmotic and mean activity coefficients predicted by the LMPB1, LMPB2, and LMPB3 theories for model electrolytes mimicking NaF, NaCl, HCl, LiI, and MgCl 2 solutions are consistent among themselves, and show a very good agreement with those from the MPB and the MSA theories. The linear results also reproduce the MC data for these systems to a comparable degree of accuracy.
We have also studied structural aspects of RPM electrolytes at different concentrations through the mean electrostatic potential and the radial distribution functions as revealed by the LMPB theories. Again, the results show an overall qualitative or better level of correspondence with the MPB and MSA data. A notable feature of the LMPB and * curves is that they show oscillations at higher solution concentrations. These oscillations are manifestations of inter-ionic correlations and occur since linearization of the MPB equation retains the aspects of the (MPB) fluctuation terms in the LMPB equations as seen in [31]. In contrast, the classical mean-field theories such as the DH and SPB do not incorporate ionic correlations and as such fail to capture such oscillations.
An interesting finding of the present work is the equivalence of the excess internal energy and the electrical contribution to the mean activity coefficient for linear theories. In the course of our calculations we have found this to be true of the LMPB1, LMPB2, LMPB3, and the MSA -all linear theories, and under all physical conditions. This has also been observed by Kjellander [44] with his MDE-DH theory of electrolytes. We have outlined here a general argument as to why this should necessarily be so for all linear theories. This exact result is likely to be useful when working with such theories. For the non-linear SPB and MPB cases, however, we have found the numerical values of these two quantities to be close, but not identical.
The results of this study give a practical relevance to the LMPB approach. In some sense, the LMPB theories are intermediate between the DH/DHLL theories and the MPB or other formal statistical mechanical theories. There are parallels to the DH/DHLL in that many of the thermodynamic and structural quantities of interest in the LMPB theories are analytical. However, unlike the former theories, the latter incorporate the ionic exclusion volume and correlation terms and hence improve the accuracy of their results relative to the DH/DHLL. These features of the LMPB theories can be useful in the routine, everyday analysis of experimental data. Similarly, in the analysis of numerically intensive theoretical problems (see for example, reference [52]) the use of LMPB analytical expressions as initial input in iterative processes can prove to be useful. Furthermore, the success of the LMPB approach in analyzing the thermodynamics of 1:1 valency systems makes this a potentially attractive method that can be used to explore more complex situations such as higher valency systems [44] and/or systems with a variable dielectric constant [53]. Studies of higher valency and mixed electrolyte systems would also be useful in order to see the limitations of the LMPB theories.