Thermodynamics of primitive model electrolytes in the symmetric and modified Poisson-Boltzmann theories. A comparative study with Monte Carlo simulations

Osmotic coefficients, individual and mean activity coefficients of primitive model electrolyte solutions are computed at different molar concentrations using the symmetric Poisson-Boltzmann and modified Poisson-Boltzmann theories. The theoretical results are compared with an extensive series of Monte Carlo simulation data obtained by Abbas et al. [Fluid Phase Equilib., 2007, 260, 233; J. Phys. Chem. B, 2009, 113, 5905]. The agreement between modified Poisson-Boltzmann predictions with the"exact"simulation results is almost quantitative for monovalent salts, while being semi-quantitative or better for higher and multivalent salts. The symmetric Poisson-Boltzmann results, on the other hand, are very good for monovalent systems but tend to deviate at higher concentrations and/or for multi-valent systems. Some recent experimental values for activity coefficients of HCl solution (individual and mean activities) and NaCl solution (mean activity only) have also been compared with the symmetric and modified Poisson-Boltzmann theories, and with the Monte Carlo simulations.


Introduction
Designing many of the important industrial chemical processes involving aqueous electrolytes requires an understanding of the thermodynamics of these systems. Ion-ion or ion-surface charge interactions of electrolyte solutions in water are fundamental processes also in areas such as biology and physical chemistry (see for example, references [1][2][3][4][5][6][7][8][9][10]). These interactions relate to the activity of water and the individual activities of the ions in solution. Two of the more important properties of such Coulomb fluids, which describe the thermodynamics, are the osmotic coefficient φ and the activity coefficient γ. These are experimentally measurable quantities [2] that are widely used in physical chemistry and electrochemistry to quantify the deviation of experimental measurements from theoretically predicted ideal cases.
Theoretically, the classical approach to estimating the osmotic and activity coefficients is based on the Debye-Hückel (DH) theory [11]. However, due to its mean-field nature and, more importantly, due to its linear character, the theory is somewhat restricted in its applications. The Debye-Hückel limiting law (DHLL), however, provides useful limiting expressions for many physical quantities of interest [12]. A lasting legacy of the DH theory is the underlying physical model that it portrays, that is, the point ions moving in a dielectric continuum. By imparting a size to the ions we arrive at the primitive model (PM) of electrolytes, viz., charged hard spheres in a dielectric continuum. If the sizes of the spheres are the same, then we have a restricted primitive model (RPM). The RPM and the PM have been extensively used over the years in formal statistical mechanical theories to describe the structure and thermodynamics of charged fluids in the bulk and near charged interfaces. Mention can be made of liquid structure integral equations such as the hypernetted chain [13] and the mean spherical approximation [13][14][15], and potential based approaches, viz., the symmetric Poisson-Boltzmann (SPB) [16][17][18][19] and modified Poisson-Boltzmann (MPB) [19][20][21][22] theories. Parallel numerical simulations, for example, using Monte Carlo (MC) and Molecular Dynamics simulations [23] have been valuable in theoretical developments.
A few years ago, Abbas et al. [24,25] made extensive MC simulations for a series of symmetric (in ion size and valence) and asymmetric PM electrolytes and reported their osmotic coefficient, single ion and mean activity coefficient results. In total, they treated 104 electrolyte systems covering a wide range of concentration and different valence combinations. The object of their work was to see how valid the PM was in representing the experimental results of salt solutions. Ionic activity and osmotic coefficients were calculated for 1 + :1 − ,2 + :1 − , and 3 + :1 − electrolytes and best fitting ionic radii determined for describing the experimental results. In this paper, we focus on obtaining the corresponding SPB and MPB results for these systems [24,25], and compare the theoretical predictions against the exact MC data.
Experimental measurements of osmotic and mean activity coefficients of electrolytes also abound in the literature [2][3][4]. These two quantities are related by the famous Gibbs-Duhem equation of physical chemistry [1], which has now been extended to multicomponent fluids [26]. The advent of sophisticated experimental techniques over the past two decades has made determination of activity of individual ionic species in electrolytes increasingly feasible [27]. Here, we will also make comparisons of the SPB and MPB predicted activity coefficients with the experimental results for HCl due to Sakaida and Kakuichi [28] and for NaCl taken from the literature [2].

Model and Methods
As suggested in the Introduction we have treated here the PM of the electrolyte, that is, the constituent ions are mimicked by rigid spheres of arbitrary radii with an embedded charge of arbitrary valence at the centre of each sphere. The charged hard spheres move in a continuum solvent characterized by a single parameter, viz., a relative permittivity ǫ r . The model is the one used in the MC studies of Abbas et al. [24,25].
The pair interaction potential in the Hamiltonian is given by where r s and Z s are the radius and valence of ion species s, e is the proton charge, r is the separation between the centers of ions i and j, and ǫ 0 the vacuum permittivity. In a special case when the ions are of the same size, that is, r i = r j , we have the RPM. The above PM (or the RPM) was treated by the symmetric Poisson-Boltzmann and the modified Poisson-Boltzmann theories. The development and the formulation of the SPB and the MPB have been detailed elsewhere in the literature (see references [16][17][18][19] for the SPB and references [19][20][21] for the MPB) and will not be repeated here.

Results
The SPB and MPB equations were solved numerically using a quasi-linearization iteration scheme [30]. This is a well tested, robust method, which has been successfully utilized in earlier works involving the SPB and MPB [16,20].
The results presented in this paper correspond to 1:1, 2:1, and 3:1 valence electrolytes at temperature T = 298 K in a water-like solvent with a dielectric constant ǫ r = 78.5, and for a wide range of concentration. These values are consistent with some experiments we have compared our results with [2,28] and the MC simulations of Abbas et al. [24,25]. The ion sizes used are also from their simulations, and although ions are generally not spherical, the ion sizes represent optimized values against experimental osmotic coefficients and mean activity coefficients.

Comparison with experiments
We begin this discussion by comparing the SPB and MPB activity coefficient results against some experimental data for HCl and NaCl. There is some debate in the literature over the experimental determination of individual ion activities [27,29]. Wilczek-Vera and Vera [27] have compared some experimental results for different salt systems and have found, in many situations, a reasonable agreement between the results for the salts using different techniques.

HCl
Sakaida and Kakiuchi have measured the individual activity coefficients of HCl (and hence the mean activity coefficient) at room temperature using an ionic liquid salt bridge [28]. We have calculated the corresponding SPB and MPB results for this system using the H + , Cl − radii taken from (a) the Abbas et al. [24,25] MC simulation data,viz., r H + = 2.20 × 10 −10 m, r Cl − = 1.81 × 10 −10 m, and (b) from Fraenkel's [31] theory, r H + = 0.58 × 10 −10 m, r Cl − = 1.81 × 10 −10 m. The experimental data are in "molal" scale, whereas the MC and theoretical data are in "molar" scale. For consistency, we have converted the experimental data into "molar" units. However, the difference between the two sets of data are negligible since the range of concentration used in the experiment is small. The experimental and theoretical results for the individual and mean activity coefficients for this system at various concentrations are shown in figures 1 and 2, respectively. Following literature custom we have plotted the logarithms of the activity coefficients as function of the square root of the salt concentration. In each figure, the top panel shows the individual activity coefficients, while the bottom panel gives the mean activity coefficients.
In figure 1 (top panel) the experimental, theoretical, and simulation individual activity coefficient data are close together up to about c ∼ 0.07 mol/dm 3 . Indeed, in the limit of low concentration, they tend  [24,25]. The experimental data are from [28]. to the same limiting value consistent with the DHLL. However, as the concentration increases beyond c ∼ 0.07 mol/dm 3 , the simulation and the theoretical activities deviate from the experimental results, especially for the cationic activity coefficient. It is noted that for the range of concentration plotted, the SPB and MPB results follow the simulations rather well. Although the MPB difference is almost quantitative with the MC, the SPB shows deviations at higher concentrations, which can be attributed to the neglect of ionic correlations in the mean-field theory. Such correlations become more significant for dense solutions. In contrast to the situation in the top panel, the experimental, the MC, and the theoretical mean activity coefficients in the bottom panel remain relatively close to each other throughout, which suggests some error cancellations for the MC and the theories.
The MC data are conventionally taken to be "exact" for a given physical model. These results suggest that the PM probably does not lead to an adequate physical description of HCl at high concentrations. Fraenkel [31,32] took this idea further and developed the Smaller ion Shell (SiS) theory, which is based on an improved DH type theory with different ion sizes, and gives a better representation of the HCl experimental individual activity coefficients [28]. The SPB and MPB activities using ionic radii from Fraenkel's theory [31] are given in figure 2. Although MC data are not available at these ionic radii, from the consistency of the simulation data with the theories for all the 104 cases for which the MC data are available, it is a fair conjecture that the SPB and MPB curves in this figure will be very close to the "exact" results for the model. The relative behaviour of the theoretical curves relative to the experimental curves are similar to that in figure 1. At these radii too, the shortcomings of the PM are clearly seen.

NaCl
We have plotted the SPB and MPB results for the mean activity coefficients at three sets of ionic parameters: (i) r Na + = 0.97 × 10 −10 m, r Cl − = 1.81 × 10 −10 m [31], (ii) r Na + = 1.68 × 10 −10 m, r Cl − = 1.81 × 10 −10 m [24,25], and (iii) r Na + = 1.58 × 10 −10 m, r Cl − = 1.91 × 10 −10 m [24,25]. The experimental values of the mean activity coefficients are from Robinson and Stokes [2], and are given in terms of γ ± as function of concentration, hence the MC, SPB, and MPB plots of γ ± for consistency. The results are shown in figure 3. At the Fraenkel ion-size parameters (lower panel), the theoretical mean activities are qualitative with the experiment, with no MC data being available at these parameters. However, in the upper panel of the figure, for both sets of ionic radii, the MPB mean activities are in excellent agreement with the MC data with the SPB showing deviations (from the MC and MPB) as concentration increases. The theories and the simulations are again qualitative and relatively closer to the experimental trends than with the Fraenkel parameters..

Comparison with simulations
Abbas et al. [24,25] simulated osmotic coefficients, individual and mean activity coefficients of a total of 104 salt solutions. We have utilized their MC parameters to obtain the SPB and MPB results for the same systems, which encompass a wide range of concentration, ionic sizes, and valence combinations of 1:1, 2:1, and 3:1 [33].
We will present here, a sampling of the results for 1:1, 2:1, and 3:1 valence systems, which will involve both the RPM and PM. For many of the electrolyte systems, Abbas et al. [24,25] obtained a single common ionic radius optimized through comparing the simulated and experimental osmotic coefficient. We have chosen such systems because this affords a comparison of the RPM and PM for the same experimental system. For the 1:1 case, the results for LiOH salt are shown for the RPM (figure 4) and PM ( figure 5). LiOH is a widely used salt in rechargeable batteries and hence the importance of its thermodynamic properties. It is noted that for the symmetric valence RPM situation in figure 4, we have ln(γ + ) = ln(γ − ) = ln(γ ± ). In both RPM and PM cases, the MPB φ, ln(γ i ), and ln(γ ± ) are almost quantitative with the simulations, while not unexpectedly, the SPB values are qualitative and deviating at higher concentrations.
For the 2:1 case, we have chosen the MgCl 2 salt system, with the RPM and PM results being displayed in figure 6 and figure 7, respectively. Again, the MPB remains nearly quantitative, and the SPB is qualitative but close to the simulations.
For a still higher asymmetric case of a 3:1 system, we have the AlCl 3 electrolyte. The RPM results  The MC data are from references [24,25]. for this system are in figure 8 and the PM results are in figure 9. In these figures, the behaviour of the curves seen for the 1:1 and 2:1 cases continues, although for this high valence asymmetry situation, even the MPB reveals some deviations from the MC. The deviations shown by the SPB are bigger although the theory is still qualitative. The behaviour pattern seen in these figures for the SPB and MPB curves vis-a-vis the simulations is repeated for the remaining 95(104-9) cases not shown here, the most conspicuous general feature being the consistent overall agreement between the theories and the MC results for almost all the simulations, this being true for both the individual and mean activity coefficients.

Conclusions
We have employed the symmetric Poisson-Boltzmann and the modified Poisson-Boltzmann theories of statistical mechanics to characterize the thermodynamics of electrolyte solutions. The osmotic coefficient, the individual activity coefficients and hence the mean activity coefficients of 104 primitive model electrolyte systems with arbitrary ionic sizes and ionic valences were calculated and the results were compared with the corresponding Monte Carlo simulations. In addition, the theoretical predictions were also contrasted with the experimentally measured activity coefficients of HCl and NaCl solutions.  We have found that overall the SPB and MPB theories reproduce the MC simulation results to a remarkable accuracy for 1:1 electrolytes, and to a slightly lesser extent for the asymmetric 1:2 and 1:3 valence systems. At concentrations higher than approximately 2 mol/dm 3 , discrepancies between SPB and MPB tend to become relatively more prominent even for 1:1 valence cases. This can be clearly appreciated in every graph of MC vs MPB and SPB and is rooted in the neglect of interionic correlations in the classical mean field approximation [34]. Although at dilute solution concentrations, the influence of correlations is minimal, they can be substantial at higher concentrations. These findings lead to the main conclusion of this work: The agreement of the MC data to the SPB, MPB results is excellent at small concentrations, whereas the SPB and MPB theories tend to deviate from each other at higher concentrations. Overall, the MPB reproduces the MC results semi-quantitatively or better.
The comparison of the experimental individual activity coefficients of HCl solutions with the corresponding simulations and the theoretical results reveals some shortcomings of the primitive model itself. Indeed, the discrepancy between the MC and the experiments must be due to the inadequacy of the model. However, the consistency of the mean activity coefficients implies some cancellation of errors.
Perhaps the most serious approximation in the PM is the continuum approximation of the solvent, this limitation having long been known [2]. Besides steric effects, a dissimilar polarization of the water molecules in different ion hydration shells would be expected to have a critical bearing on the individual activity coefficients. The classical Born solvation model uses the continuum solvent, so it cannot give the necessary molecular description. Born models are still useful though, a fairly recent example being the use of simulated PM results with a Born description of the ion-water interaction [35].
The 104 electrolyte solutions, for which MC data are available, cover a wide range of concentration, ion sizes and valence combinations. The agreement of the MPB predictions, and to a lesser extent the agreement of the SPB predictions [33] with the relevant simulation data, suggest the viability of these theoretical approaches in describing the thermodynamics of charged fluids in the electrolyte solution regime within the framework of the primitive model.