Charge reversal and surface charge amplification in asymmetric valence restricted primitive model planar electric double layers in the modified Poisson-Boltzmann theory

The modified Poisson-Boltzmann theory of the restricted primitive model double layer is revisited and recast in a fresh, slightly broader perspective. Derivation of relevant equations follow the techniques utilized in the earlier MPB4 and MPB5 formulations and clarifies the relationship between these. The MPB4, MPB5, and a new formulation of the theory are employed in an analysis of the structure and charge reversal phenomenon in asymmetric 2:1/1:2 valence electrolytes. Furthermore, polarization induced surface charge amplification is studied in 3:1/1:3 systems. The results are compared to the corresponding Monte Carlo simulations. The theories are seen to predict the"exact"simulation data to varying degrees of accuracy ranging from qualitative to almost quantitative. The results from a new version of the theory are found to be of comparable accuracy as the MPB5 results in many situations. However, in some cases involving low electrolyte concentrations, theoretical artifacts in the form of un-physical"shoulders"in the singlet ionic distribution functions are observed.


Introduction
Electric double layers play a key role in a wide variety of phenomena, in the fields diverse as colloids, physical and biological systems, and electrochemistry [1][2][3][4][5][6][7][8]. Recent years have seen an increasing relevance of double layers in DNA sequencing [9,10], drug delivery [11,12], and green technology, viz., designing novel energy storage devices as in super-capacitors (see for example, reference [13]). An electric double layer appears when a system of charged particles in a fluid is in the neighborhood of a charged surface. Various models and theories have been developed depending upon the physical system being treated. A simple model of an electrolyte next to a plane charged surface is a basic model in many areas. Gouy and Chapman [14,15] analysed this model, based on the Poisson-Boltzmann (PB) equation, where the electrolyte was treated as a system of point ions moving in a medium of constant permittivity. Later on, Stern [16] proposed to lower the permittivity in the region immediately next to the planar electrode (Stern layer) in order to bring the computed double layer capacitance in line with the experiments. At not too high surface charge and electrolyte concentration for univalent ions, the Gouy-Chapman-Stern (GCS) theory describes reasonably well the diffuse part of the electric double layer. However, double layer research over the last four decades involving formal statistical mechanical

Model and theory
The electric double layer is modelled by a restricted primitive model electrolyte in the neighborhood of a uniformly charged plane electrode. The solute molecules are hard spheres of diameter a with point charges at their center moving in a medium of constant relative permittivity r , while the electrode has a constant relative permittivity of w . The analysis presented here builds upon the earlier formulations of MPB4 [38,39] and MPB5 [40] versions of the theory and clarifies the relation between them. The present approach is a bit more general since the MPB5 situation will be seen to follow as a special limiting case. The connections to the MPB4 are also explored.
At a perpendicular distance x from the electrode into the solution, the mean electrostatic potential ψ(x) satisfies Poisson's equation where the sum is over the ion species, 0 is the vacuum permittivity, e s is the charge of an ion of species s, n s and g s (x) are the mean number density and singlet distribution function for an ion of species s, respectively. Kirkwood's charging process [41] for an ion i at x 1 gives where is the fluctuation potential at r 2 for an ion at r 1 , r = r 12 , and β = 1/(k B T) (k B is the Boltzmann's constant and T is the absolute temperature).
The fluctuation potential can be shown to satisfy the following set of differential equations [40] where x = x 2 and κ(x) is the local Debye-Hückel parameter defined by

33801-2
Equations (4) and (5) are exact, while equation (6) involves two approximations. These are Loeb's closure [42], which gives a closed system of equations for φ i , and the linearization of the equation involving Loeb's closure. The boundary conditions associated with equations (4)-(6) are φ i and ∂φ i /∂n continuous (n denotes normal) at x = a/2 and r = a, while at the electrode φ i is continuous but Approximate solutions of equations (4)- (6) have been considered for both point ions [43] and finite size ions [40]. The point ion solution suggests an appropriate solution in the regions x < 0, 0 < x < a/2. From the analysis of reference [43], the limiting case of having no distance of closest approach and k → 0 gives where f = ( r − w )/( r + w ) and r * is the distance of the image of e i in the electrode from the field point. These two results suggest that an appropriate approximation for the ion size result for the equation (4) is which were used in reference [40]. Previous work [38,40,[44][45][46][47] has shown that a good approximation to the solution of equation (6) is where κ = κ(x 1 ). The point ion result, with a distance of closest approach, indicates a complex situation such that B, B * , C can contain the imaging factor f . The use of the approximate solutions (10)- (12) mean that the boundary conditions cannot be satisfied. An approximate fitting of the boundary conditions was considered in [40], giving the MPB5 theory. Here, we adopt the approach of reference [38]. The application of Green's second identity to φ * i and 1/r in the exclusion sphere V, where V is a truncated sphere if a/2 < x < 3a/2, gives, on using equation (5) 4π lim Here, S W and S I are respectively the wall and double layer surface portions of the truncated sphere. The integrals over S W vanish for x > 3a/2, and S I is then the complete spherical surface. For reference, the values of the integrals appearing in equation (13), using the approximations for φ i , are given in the appendix.
To determine the relations between B, B * and C we use Gauss's theorem for ion i Thus, combining the expressions resulting from the evaluation of equations (13), (14) with equation (2) 33801-3 gives where for a/2 < x < 3a/2 and for x > 3a/2 where F 0 = lim x→∞ F, y = κa and s = √ a 2 + 2xa. At x = 3a/2, we have that F is continuous but ∂F/∂ x is discontinuous in general. The expression for F depends upon the ratios C/B, B * /B for a/2 < x < 3a/2 and on the single ratio B * /B for x > 3a/2. Putting C = B exp(−y), B * = B f exp(−y + κs) in a/2 < x < 3a/2 and B * = B f exp y in x > 3a/2 gives the MPB5 theory [40]. Early calculations of F did not consider the exclusion region 0 < x < a/2 for φ i so that F is defined in the two regions a/2 < x < a and x > a, for example the MPB4 theory [38,39]. In these earlier works, the value of F differs from that of the MPB5 when x > 3a/2. On putting B * = f B in equation (19), the F of MPB4 for x > a is derived. This suggests an alternative F, analogous to that for MPB4, which is given by C = B in equations (17), (18) and B * = f B in equations (17)- (19). There are many possibilities of choosing the ratios C/B, B * /B, depending upon the chosen criterion. However, a better approach may be to improve upon the analysis of equations (4), (5) and the nonlinear version of equation (6).
Various expressions have been used for evaluating the exclusion volume term ζ i (x) in equation (2) [38,39,48]. The most accurate exclusion volume terms are those based on the BBGY expansion [39] or charging up the wall and introducing the non-uniform direct correlation function [48]. Here, we will use the BBGY approach, where for x > a/2 with ζ i (x) = 0 for x < a/2. The contact value ζ i j (a) between ions i and j is that proposed by Fischer [49] while Φ(y, X) = φ(1; 2/e s = 0/r 12 = a).

Monte Carlo simulations
The Monte Carlo simulations were performed in the canonical (N, V, T) ensemble using the standard and widely used Metropolis algorithm [50,51]. The techniques were the same as those used in some of our earlier works (see for example, references [52,53]). The MC cell was a rectangular parallelepiped of dimensions L x , L y , and L z , with L y = L z . One of the y-z faces was a charged wall, viz., the electrode 33801-4 with a uniformly distributed surface charge at x = 0, while the other at x = L x was a neutral wall. The surface charge density on the electrode was determined from where N + and N − are the number of cations and anions, respectively. It is to be noted that with the above definition of σ, the local electro-neutrality at the wall is built in. In order to simulate the semi-infinite system, conventional minimum image techniques along with periodic boundary conditions in the y and z directions were implemented [51]. The effect of the long range nature of the Coulomb interactions was treated by the parallel charge sheets method pioneered by Torrie and Valleau [17]. This procedure was later examined and improved upon by Boda et al. [54], who also confirmed the validity of the method for charge hard spheres. The situation when the permittivity of the medium of the charged wall differs from that of the electrolyte solution has also been treated by Torrie at al. [55] and their procedure adopted here. These authors regarded the polarizability of the electrode as a classic electrostatic image problem, viz., the total surface charge density σ has a contribution from the fictitious image charge − f σ due to the polarization of the electrode. In passing we would like to note that the Torrie and co-workers did their simulations in the grand canonical (µ, V, T) ensemble unlike the canonical ensemble in the present case.
In the canonical ensemble, the bulk target concentration can only be achieved by a trial-and-error method of adjusting the cell length L x and/or the number of particles. We used both of these methods where necessary and imposed a tolerance of less than 2% error in reproducing the bulk concentration. Typically, around 400 particles were simulated and approximately 10 8 configurations sampled out of which the first 10 7 (10% of the total number of configurations) were used for system equilibration before statistics were taken.

Results
The properties of the diffuse double layer are found by solving the Poisson equation, equations (1) and (15), for the mean electrostatic potential and hence from equation (15) the singlet distribution function. Three different values of F were treated in detail, these corresponding to those of MPB4, MPB5, and the new F given by C = B, B * = f B in equations (16)- (19), respectively. A quasi-linearisation technique was used [38,56] which has been successfully employed in past investigations. The parameters chosen were those used by Wang [37] who simulated charge asymmetric 1:3 and 3:1 RPM electrolytes with imaging. For the RPM electrolyte, we treat 1:2, 2:1 and some 1:3, 3:1 valences with a = 4 · 10 −10 m, r = 80, T = 298 K while the wall is characterized by w = 2, 80 or ∞. When w r , imaging occurs which is repulsive for w < r and attractive for w > r . The value w = 2, which is relevant for a medium such as silica, while w = ∞ represents a metallic electrode.
We start by looking in detail at the predictions of the MPB5 theory for 2:1 electrolytes with and without imaging. Figures 1 and 2 treat the repulsive ( f = 0.9513) and attractive ( f = −1) cases, respectively, at electrolyte concentration c = 0.24 mol/dm 3 and surface charge density σ = −0.02 C/m 2 . Displayed are the singlet distribution functions g s (x/a), the dimensionless mean electrostatic potential ψ * (x/a) = |e| βψ(x/a) and the integrated charge density q(x/a) defined by The integrated charge gives a measure of the total charge within a distance x of the electrolyte and enables a succinct illustration of the two phenomena CR and SCA. For no imaging, the simulated singlet distribution functions apparently display a monotonous behaviour of the co and counter ions, with a corresponding monotonous behaviour of the mean electrostatic potential and the integrated charge. The MPB5 theory has a small shoulder in the counterion g s (x) at x ∼ 3a/2 arising from the discontinuity of ∂F/∂ x at x = 3a/2. This is an unfortunate feature at low 33801-5 The symbols represent MC data, while the lines represent the MPB results. The notations "co" and "cntr" stand for "coion" and "counterion", respectively. Legend as given in the figure. concentrations which has been observed earlier [40]. The influence of this feature is not observed in ψ * (x) since ψ * (x) and its derivative are continuous everywhere. The principal effect of imaging in the distribution functions for the repulsive situation is that the divalent counterion is attracted for large x but reduces and becomes small near the wall, while a slight reduction is seen in the monovalent coion. Conversely for the attractive situation, the coion profile increases to a value greater than 1 near the electrode, and the counterion profile rapidly increases in the neighborhood of the electrode. These imaging effects are reflected in the reduction (repulsive case) and increase (attractive case) over the no imaging situation, of the mean electrostatic potential and integrated charge. A very small CR is seen when f = −1. The MPB5 accurately predicts ψ * (x) and q(x), and apart from the shoulder feature, the g s (x) are quantitatively or qualitatively correct.
At the higher concentration c = 1 mol/dm 3    the bulk Debye-Hückel constant. A damped oscillatory behaviour in ψ * (x) is predicted in the linear MPB theory for y 0 > 1.2412 [57] for symmetric valences and for y 0 > 1.15 [45] for 2:1/1:2 valences. However, the latter is probably too high since in simulation and the HNC, the g s (x) value for 1:2 homogeneous electrolytes implies y 0 ≈ 0.76 [58]. Here, y 0 = 2.26, while the previous concentration of 0.24 mol/dm 3 gave y 0 = 1.11, so that a qualitative change in the behaviour between the two concentrations is predicted by the MPB theories. With no imaging, the initial maximum in ψ * (x) is at x ∼ a with CR starting at about the same value of x and having a maximum at x ∼ 5a/4. Both the repulsive and attractive cases are qualitatively similar to that for no imaging, showing small but different deviations from f = 0. For example, with the mean potential, the first maximum is shifted to a larger x for the repulsive case, but closer to the electrode with a larger maximum for the attractive case. At this higher concentration, there is no shoulder in the MPB5 g s (x) at x ∼ 3a/2, and the MPB5 gives an accurate picture of all three functions. Clearly, the effect of discontinuity in ∂F/∂ x at x = 3a/2 gets masked at high concentrations. The MPB4 theory seems a possible alternative to the MPB5 theory as it does not have this discontinuity at x = 3a/2. Figures 5 and 6 show the MPB4 for the 2:1 attractive cases. Clearly, the theory is inadequate, especially for the lower concentration. The main inaccuracy is the unphysical upturn in the coion distribution function close to the wall. A further possible alternative theory is that given by C = B, B * = f B in 33801-7  equations (16)- (19). Looking again at the 2:1 attractive cases, figures 7 and 8, this theory incorrectly predicts (i) an unphysical rise in the coion g s (x) near the wall, (ii) a shoulder in the counterion g s (x) for no imaging at the lower concentration. Apart from the unphysical behaviour (i), the theory based on the new F gives a good representation of the mean potential and integrated charge at c = 1 mol/dm 3 . Unfortunately, neither the MPB4 nor that based on the new F is an improvement on the MPB5. Figures 9 and 10 consider the prediction of the MPB5 theory for the 1:2 electrolyte where now the coion is divalent. Counterion attraction with the wall plays a dominant role in shaping the structural properties of the electric double layer and hence many of the 1:2 double layer properties are similar to those of the 1:1 case, while the 2:1 properties resemble that for the 2:2 situation (see for example, reference [59]). Thus, it comes as no surprise that the predicted structure for the 1:2 cases are relatively closer to the simulations than that seen earlier for the 2:1 cases. Only attractive imaging is treated in figures 9, 10 since any deviation from no imaging tends to be greater than those for the repulsive case. At a low concentration, figure 9, the coion is attracted near the wall, and since it is divalent, it has a greater contact value than in the monovalent case in figure 2. The difference in the valence of the counterion means that the deviations from no imaging for the mean potential and an integrated charge are less and oppositely directed. At a higher concentration, figure 10, the imaging has little effect. There is still a

33801-8
Charge reversal and surface charge amplification using the modified Poisson-Boltzmann theory  damped oscillatory feature in the functions, with CR occurring, but it is not as pronounced as for the divalent counterion. Overall, the MPB5 theory accurately predicts the 1:2 double layer, apart from the behaviour of the coion in the neighborhood of x = 3a/2. The alternative MPB4 and new F theories perform well in the 1:2 case.
Lastly, we briefly look at the MPB theory for a 1:3 and 3:1 electrolyte with the imaging at c = 0.24 mol/dm 3 and σ = −0.02 C/m 2 . Replacing the divalent coion by a trivalent ion leads to CR ( f = 0.9513, f = 0) and also to SCA when there is an attractive imaging ( f = −1), viz., figures 11, 12. Alternatively, when the trivalent ion is the counterion there is again CR and SCA, but now SCA occurs when there is a repulsive imaging and CR for an attractive imaging, viz., figure 13. MC simulations [37] indicate that the MPB5 prediction of SCA for no imaging is incorrect. The CR is expected as y 0 = 1.5636, while the counterintuitive SCA arises solely through polarization forces. In the 1:3 situation, the attraction of the coion image overcomes the repulsion of the surface charge, while for the 3:1 electrolyte, a much larger repulsive force acting on the counterion is sufficient to neutralise the screening in the neighborhood of the wall. SCA does not occur for either 1:2 or 2:1 electrolytes at these parameters. Higher valences provide a stringent test for theories. The discrepancies in g s (x) mean that the MPB5 and the new F theory predict both CR and SCA with varying degrees of success.

Conclusions
In this paper, we have looked again at the formulation of a modified Poisson-Boltzmann equation for the planar electric double layer formed by a restricted primitive model electrolyte. The derivation of the MPB equation follows the pattern and techniques utilized in some of our earlier works [38,39] but is a bit more general. Our focus was on obtaining a general expression for the quantity F, a key ingredient in the theory, which incorporates important aspects of inter-ionic correlations. We have shown that an earlier version of the theory, viz., MPB5 [40], follows as a special case of the generalized F. Although the MPB4 version does not follow directly, the relationship between these two versions can be clearly assessed. With such features built in, the theory is more flexible and it is likely that the MPB equation will become more adaptable to particular physical situations.
The MPB theory does not satisfy exactly the contact condition [60][61][62] which holds for no imaging. Previous work [63][64][65] has shown that the MPB5 theory fairly accurately satisfies the contact condition for 1:1 electrolytes and gives a good agreement for 1:2/2:1 and 2:2 electrolytes. A similar conclusion holds in general when the new F is used, but the MPB4 theory can be poor. Adapting the MPB approach to satisfy the contact condition would be useful for inner consistency. However, for the present univalent and divalent RPM electrolytes, the MPB5 theory has shown to be in good overall agreement with both structural and thermodynamic simulation properties [59].   We have studied the effect of imaging on 1:2, 2:1 valence electrolytes with the emphasis on the phenomenon of charge reversal. Also briefly treated were 1:3, 3:1 electrolytes where, besides the charge inversion, polarization can produce the counterintuitive surface charge amplification. The conclusions are mainly twofold: (a) structural and charge reversal behaviour are generally well described by the various MPB theories 33801-11 compared to MC simulations. For the 1:3, 3:1 systems, the theories are poor relative to the 1:2, 2:1 cases, which is not unexpected at the higher multi-valence situations. Furthermore, the theoretical predictions tend to be closer to the simulation data for monovalent counterions rather than for multivalent counterions, a feature that is well known in the double layer literature.
The new formulation of F allows for different choices to be made for F. One of the new F's analyzed gives reasonable results that in some situations are at a par with MPB5. However, there are issues, the main one being the appearance of an nonphysical shoulder in the singlet distribution functions at lower concentrations, similar to that seen in MPB5.
The shoulders are artifacts of the theory, which come from discontinuities in the slope of the fluctuation potential φ at x = 3a/2. In the analysis outlined, although φ is continuous across the planar surface of the truncated exclusion sphere (cf. section 2), the normal derivative is not continuous. The continuity of the normal derivative may well be achieved by an appropriate choice of B, B * , C in equations 16,19. It is to be seen if such choices in the present approach will give a simple means of eliminating the shoulder in the distribution function. Further analytic work along the lines of MPB5 may be difficult, so that eventually the best approach may well be a numerical evaluation of the fluctuation potential.