An analysis of the fluctuation potential in the modified Poisson-Boltzmann theory for restricted primitive model electrolytes

An approximate analytical solution to the fluctuation potential problem in the modified Poisson-Boltzmann theory of electrolyte solutions in the restricted primitive model is presented. The solution is valid for all inter-ionic distances, including contact values. The fluctuation potential solution is implemented in the theory to describe the structure of the electrolyte in terms of the radial distribution functions, and to calculate some aspects of thermodynamics, viz., configurational reduced energies, and osmotic coefficients. The calculations have been made for symmetric valence 1:1 systems at the physical parameters of ionic diameter $4.25 \times 10^{-10}$ m, relative permittivity 78.5, absolute temperature 298 K, and molar concentrations 0.1038, 0.425, 1.00, and 1.968. Radial distribution functions are compared with the corresponding results from the symmetric Poisson-Boltzmann, and the conventional and modified Poisson-Boltzmann theories. Comparisons have also been done for the contact values of the radial distributions, reduced configurational energies, and osmotic coefficients as functions of electrolyte concentration. Some Monte Carlo simulation data from the literature are also included in the assessment of the thermodynamic predictions. Results show a very good agreement with the Monte Carlo results and some improvement for osmotic coefficients and radial distribution functions contact values relative to these theories. The reduced energy curve shows excellent agreement with Monte Carlo data for molarities up to 1 mol/dm$^{3}$.


Introduction
One of the more consistently active areas of research in the statistical mechanics of fluids over the years has been in the field of Coulomb fluids. These encompass among others, electrolytes, ionic liquids, molten salts, colloids, and polyelectrolytes, the practical relevance of which extend from biological systems to industrial chemical processes. The literature on this is vast and theoretical progress was limited until the application of liquid state theory [1][2][3][4][5] based on classical statistical mechanics. We would like to cite here a few of the recent reviews on the subject [6][7][8].
A widely used model used in the development of formal statistical mechanical theories of ionic solutions treats the solvent as a structureless, continuous dielectric medium with a relative permittivity ǫ r , and the solute particles as charged hard spheres of arbitrary diameters d i and charges Z s e with Z s being the valence of species s. This is the so-called primitive model (PM) of ionic solutions. When the ions are of the same size, it is called the restricted primitive model (RPM). Computer simulations of the RPM and PM over the years (see for example, references [9][10][11][12][13][14]) have shown the usefulness of these models in interpreting experimentally determined structures and thermodynamics of charged fluid systems. Furthermore, the simulation data have proved to be invaluable in theoretical development.
The statistical mechanics of primitive models in liquid state physics has followed two broad paths: In the first, the focus is on computing the pair correlation function or the radial distribution function g i j (r i , r j ) from the inter-molecular pair potential u i j (r i , r j ) starting from the Ursell-Mayer cluster expansion [1][2][3], or the distribution function method [3,5]. Two main routes are used,viz., the Kirkwood, Bogolubov, Born, Green, Yvon (KBBGY) hierarchies (see for example, reference [5]) and the Ornstein-Zernike (OZ) equation [2,3,5]. The KBBGY hierarchies relate correlation functions for n and n + 1 fixed particles, the molecular potential, and a charge parameter ξ. To evaluate the pair correlation function, for example, a closure relation between the pair correlation function g i j (r i , r j ) and the next higher order correlation function, that is, the triplet correlation g i jk (r i , r j , r k ) must be provided to break the hierarchy. One such relation is the superposition approximation [3]. In the OZ approach, the total correlation between two ions is considered to consist of two parts: the direct correlation function c i j (r i , r j ) between the two particles, and the indirect correlation h i j (r i , r j ), which takes into account the presence of a third particle. This is clearly shown by the OZ equation (see for example, reference [3]), which is often regarded as a definition of the direct correlation function. To solve the OZ equation, a closure relation between the direct and the total correlation functions is required. Among the more well known closures are: the Percus-Yevick (PY) [15], the Hyper-netted chain (HNC) [16], and the mean spherical approximation (MSA) [17].
In the second method, which is our interest in the present work, the focus is on obtaining the same g i j (r i , r j ), but through a potential approach to the theory based on the Poisson's equation. The classical theoretical analysis of electrolyte solutions in this regard is that of Debye and Hückel (DH) [18], which is a linearized version of the corresponding non-linear Poisson-Boltzmann (PB) equation. A key theoretical paper on an assessment of the inherent approximations in the Poisson-Boltzmann (PB) equation, and hence in the linearized DH equation is due to Kirkwood [19]. Kirkwood showed through a statistical mechanical analysis that the main approximations in the classical theories are the omission of (i) ionic exclusion volume effects, and (ii) the fluctuation potential term, which involves the inter-ionic correlations. There have been many attempts since Kirkwood to improve upon the PB/DH theory notable among which has been the extensive work done by Outhwaite and co-workers (see for example, references [20][21][22][23][24][25][26][27][28][29]), who within the framework of the PM, have analyzed Kirkwood's methods and obtained estimates for the fluctuation term. The resulting modified Poisson-Boltzmann (MPB) approach to ionic solutions is thus based on extending the classical mean electrostatic potential approach of DH theory by expressing the distribution functions in the Kirkwood, Bogolubov, Born, Green, Yvon (KBBGY) hierarchies in terms of mean electrostatic potentials. Essentially, the MPB improves upon the classical PB theory by incorporating (i) ionic exclusion volume effects, and (ii) inter-ionic correlation effects. This potential procedure solves for the mean electrostatic potential ψ(r) as opposed to the integral equations that attempt to solve directly for the radial distribution function g i j (r i , r j ). Outhwaite and co-workers [22][23][24][25][26][27][28] have further symmetrized the classical PB theory and the MPB theory so that the Onsager relation, g i j (r) = g ji (r) is satisfied for a homogeneous fluid. They have also coupled an exclusion volume term to the symmetrized PB theory, and call it the symmetric Poisson-Boltzmann (SPB) theory [25][26][27].
In the MPB theory, the mean electrostatic potential is expressed in terms of the fluctuation potential φ(1, 2; 3) (see for example, reference [28]) (3 is the field point, while there are fixed ions at 1 and 2), which measures deviations from the superposition principle of Kirkwood [19], and, therefore, contains information on the interionic correlations in the theory. Expressed in terms of the mean potentials, the fluctuation potential is given by [25,28] This equation is a statement that the mean potential at field point 3 is the sum of the direct potentials of particles fixed at 1 and 2, and the correlated potential contribution at the field point from the simultaneous presence of particles at 1 and 2. As we will see in the next section, the fluctuation potential can be written in terms of distributions functions as where e s is the charge and ρ s ({n}; q) is the number density of the s-th species of ions at r q with n fixed particles at r i (i = 1, . . . , n) with the sum being over all species, ǫ 0 is the vacuum permittivity, and ǫ r the relative permittivity (dielectric constant) of the solvent.

43801-2
In the simplest language, the fluctuation potential is the inter-ionic correlations expressed in potential form. The fluctuation potential φ(1, 2; 3) obeys a system of partial, non-linear, differential equations, and for the RPM case, the linearized version of the equations is given in reference (see for example, reference [28,29]). An approximate solution, valid for large inter-ionic separation, under the assumption of spherical symmetry, was found by Outhwaite [21]. One of the main problems in present MPB theory is the restriction of the fluctuation potential for large inter-ionic separations, where approximate spherical symmetry is valid. In the present work, an approximate analytical solution to the fluctuation potential problem is found, that is valid for the whole range of interionic distances. This solution has an advantage of simplicity that can provide insight into the eventual fully numerical methods for solving this kind of problems. The approximate analytical solution for φ(1, 2; 3) can serve as a guide to solving the problem numerically without using the approximations of this research.
The organization of this paper is as follows. In the following section (section 2) we start by giving details of the interaction potentials of the model, a brief introduction to the PB equation and the MPB theory approach. We then proceed to the main theoretical development of this work based on the primitive models. In this part, the set of differential equations for the fluctuation potential in dimensionless form is developed and an approximate solution is found using ordinary electrostatics.
In section 3 we utilize solution of the fluctuation potential to present structural and thermodynamic results for a 1:1 valence RPM electrolyte. We start by showing three-dimensional plots of the fluctuation potential solution. The plots show the fluctuation potential at a planar slice passing through the center of the ions for two ionic separations and for the like and unlike ion cases. A physical interpretation of the results in terms of ionic correlation energy is presented. To further test the solution's validity, configurational energies, and osmotic coefficients are calculated and compared to the Monte Carlo (MC) simulation data of Card and Valleau [9], and Rasiah, Card, and Valleau [10].
In section 4 we present some conclusions out of this work and stress the importance of the approach for future work that may involve a full iterative process using the solution presented here but without the approximations made.

Molecular model
As indicated in the introduction, the model electrolyte system used in this study consists of a binary, symmetric valence RPM at room temperature.
The ion-ion interaction potential in the Hamiltonian is thus where Z s is the valence of ion species s, e is the magnitude of the fundamental charge, r is the distance between the centres of two ions of types i and j, respectively, and d is the common ionic diameter. The relative permittivity ǫ r is assumed to be uniform throughout the entire system.

Theory
The formulation of the SPB and the (traditional) MPB have already appeared elsewhere in the literature (see for example, references [22,[25][26][27]), and will not be repeated here. We will restrict ourselves to outlining the main steps leading to the equations governing the fluctuation potential and their solution.
We begin by formulating the fluctuation potential problem in the restricted primitive model for a symmetric valence electrolyte, viz., |Z + | = |Z − |, consisting of N ions and satisfying global electroneutrality s Z s ρ s = 0. We will closely follow the notations used in reference [28]. In the defining relation for the fluctuation potential in equation (1), the mean electrostatic potentials ψ(1; 3), ψ(2; 3), and ψ(1, 2; 3) can 43801-3 be formally written as and, where e 1 , e 2 are the charges of the fixed ions at 1 and 2, respectively, and the sum runs over all the ionic species. Figure 1 shows the mean electrostatic potential at a field point q due to 1 and 2 fixed ions, respectively. Subtracting the equations (4) and (5) from equation (6) leads to the earlier equation (2). The Poisson equations follow and, Here, the operator ∇ is understood to operate on the coordinates of the field point. These equations can also be expressed in terms of the distribution functions using for example, g 1α (1, q) = ρ α (1, q)/ρ α , and so on and so forth, with ρ α being the mean number density of ion species α. The distributions can, in turn, be defined in terms of the potentials of mean force W, viz., the doublet or the triplet

43801-4
where W i j , W i jk are the pair and triplet potentials of mean force, respectively. Also, β = 1/(k B T) with k B the Boltzmann constant and T the absolute temperature. Hence, the conditional distribution, The term w i jk is the potential of mean force associated with the departure from linear superposition of the pair potentials. A hierarchy of such equations can be constructed for higher order correlations. At the lowest order, the classical PB theory follows upon neglecting w i jk (1, 2; 3), and to improve upon the PB, we need a procedure to estimate this term.
In the MPB formulation, the hierarchy is broken at the triplet level by a closure condition that relates the w i jk with the fluctuation potential φ i j [28] It is of interest to contrast this MPB closure with the Debye-Hückel closure For the RPM system with a finite ion diameter d, the Poisson equations (7)- (9) can be expressed in terms of the potentials of mean force as where the MPB closure (13) has been used in equation (17). The equations (15) and (16) are exact, for one fixed ion in position 1 and 2, but equation (17) incorporates the deviation from the superposition principle in the form of the fluctuation potential term. To obtain an equation for the fluctuation potential [equation (1)], we subtract equations (7) and (8) from (9), Equation (18) is the base nonlinear equation in the fluctuation potential problem. The equation also suggests that the charge density source for fluctuation potential is associated with the charged atmospheres of the triplet and doublet densities.
To illustrate the geometry of the fluctuation potential problem, one can expand the summation over species as where a number with a superscript notation with a positive or negative sign represents the presence of the corresponding ion at the referred position in space. Figure 2, represents the geometry of the fluctuation potential system of equations with Ω being the total volume of the ionic solution, ω 1 and ω 2 represent the exclusion volumes of ion 1 and 2, respectively, ω * is the overlap volume, and 3 is the field point. Region I [Ω − (ω 1 + ω 2 )] is the bulk volume defined as the total volume minus the exclusion volumes of ions 1 and 2. Region II (ω 1 − ω * ) and III (ω 2 − ω * ) are the interior of the exclusion volumes of ion1 and 2 minus the overlap volume. Region IV is the overlap 43801-5 volume. The nonlinear system of equations governing the fluctuation potential are then given by the following expressions II : III : At this point it is convenient to work in terms of reduced (dimensionless) quantities. Here, the relevant ones are the reduced mean electrostatic potential Ψ = eβψ, the reduced fluctuation potential Φ = eβφ, and y 0 = √ 24Z + Z − ηΓ. Also, η = (π/6) s ρ s d 3 is the volume or packing fraction and Γ = Z + Z − e 2 /(4πε 0 ε r k B T d) is the plasma coupling parameter. After expressing the Laplacian in ionic diameter scale, and imposing global electro-neutrality, we have a set of dimensionless fluctuation potential equations for the size symmetric case II : III : The boundary conditions are that the fluctuation potential and its normal derivative are continuous across the boundaries. Denoting the right-hand sides of these equations by P, we can write them in a general form with a formal solution [30,31] Φ(1,

43801-6
Specifically, we have in the various regions II : III : IV : In order to make analytical progress, we approximate the radial distribution functions g(1, 3) and g(2, 3), in the various P's appearing in the above equations by their DH values where the subscript in Z represents the sign of the charge state of the ion at the field point 3. Inserting the radial distribution functions (34) in the integrals in equation (29) will render the contribution to the fluctuation potential in regions II and III as ordinary integrals in space.
To obtain an approximation for P in the bulk region I, outside ions 1 and 2, we use the properties of the radial distribution functions in the various regions, and expand the exponents up to linear terms, leading to For a small fluctuation potential, we neglect the right-hand side of equation (35). For example, for a symmetric valence 1:1 RPM electrolyte, the theme of this work, we have noted that the DH radial distributions in equation (34) are of the order unity for the physical parameters and the range of concentrations used. If the fluctuation potential is of the order 10 −2 or less, then the right-hand side of equation (35) will be of a similar order and can be neglected as a first approximation for such a system. Under these approximations, the fluctuation potential is given by, where and The integral in equation (36) needs to be calculated numerically. This was done by discretization of space, and will be discussed in the next section.

43801-7
A useful way of testing the fluctuation potential solution is through subsequent evaluation of the structure and thermodynamics of the electrolyte solution. We have utilized the MPB formulation in reference [28] to calculate the pair correlation functions, where the DH functions (34) are used for Ψ(1; 2) and an analytic expression for the Percus-Yevick (PY) radial distribution functions for hard spheres [4] have been used for the excluded volume term, ζ 12 . The integral implies charging up of the ion at r 2 .
For the calculation of osmotic coefficients φ and the reduced configurational energy U/(N k B T), we use equation (12) from reference [32], written in dimensionless reduced variables as, and where r ′ = r/d with g A and g B corresponding to like and unlike ions, respectively, and the argument 1 of g A and g B in (40) refers to the contact value.

Results
All calculations in this work pertain to (1:1) symmetric valence RPM electrolyte for ions of common diameter d = 4.25 × 10 −10 m, in a continuum dielectric medium of relative permittivity ǫ r = 78.5, and at temperature T = 298 K, which is akin to a water-like solvent at room temperature. We have utilized electrolyte concentrations of 0.1038, 0.425, 1.00, and 1.968 mol/dm 3 . One reason for using these physical parameters is that these have been used earlier in the literature (see for example, reference [29] and for which MC simulation data exist [9,10]. The SPB and the conventional MPB equations were solved numerically using a quasi-linearization iteration scheme [33]. The procedure has been used with much success in earlier works [24][25][26][27] and we refer the reader to these references for further details. The fluctuation potential was obtained numerically by solving the integral in equation (36). The radial distribution functions g i j (r) were then calculated using the fluctuation potential solution in equation (39), while the osmotic coefficient φ, and the reduced configurational energy −U/(N k B T) have been determined through equations (40) and (41), respectively. In what follows we will briefly describe the numerical procedure involved before taking up the discussion of the results.

Numerical solution
The calculation of the fluctuation potential Φ(1, 2; 3) was achieved by creating a Cartesian grid in space with scaled distance of 10% of the ionic diameter, which in our dimensionless units is 1, so that in the present context, the grid spacing is 0.1. This grid was created to represent the physical regions involved in the fluctuation potential problem as shown in figure 2. Those regions consist of the spherical regions ω 1 and ω 2 , which correspond to the boundaries of the ions 1 and 2 , and the rest of the solution region, which is denoted by Ω − (ω 1 + ω 2 ). The region (ω 1 + ω 2 ) is denoted as the ionic excluded volume. The quantities F 1 and F 2 [equations (37) and (38), respectively] represent the charge densities associated with the regions ω 1 and ω 2 in the integral of equation (36). The boundary of the rectangular Cartesian grid representing figure 2 was defined by a parameter Λ, which represents the distance from the boundary of the ions to the edge of the grid. This parameter was chosen in such a way that the fluctuation potential solutions tend to zero at the exterior boundary of the grid. Usually this parameter was between 3 and 5 ionic diameters for the highest concentration but was found to a lot larger than at 43801-8 the lower concentrations. The fluctuation potential solution is an integral over the regions ω 1 and ω 2 . The summation used to numerically calculate the integral included approximately eight thousand terms for a point inside regions ω 1 and ω 2 . To produce the figures 3-5, the fluctuation potential was calculated at each point in a planar slice passing through the centers of ω 1 and ω 2 . For contact distance between the regions ω 1 and ω 2 , and Λ = 5, this planar slice contains approximately ten thousand points. The simplicity of equation (36) and the approximation of the g i j in equation (34) in terms of the corresponding DH functions are what makes the calculations fairly tenable.
The evaluation of the pair correlation functions was performed in a similar grid as the one used for the three-dimensional figures but now the fluctuation potential was only required to be calculated at the center of region ω 2 (figure 2), and the solution used in equation (39), where the Kirkwood charge integral over the fluctuation potential is calculated. The calculation of osmotic coefficient and the reduced configurational energy was achieved using the formulae (40) and (41), respectively.

Fluctuation potential
We begin this discussion with the analysis of the three-dimensional representations of the fluctuation potential Φ(1, 2; 3) shown in figures 3-5. To our best knowledge, such representation of the fluctuation potential does not presently exist in the literature. The plots show the fluctuation potential Φ(1, 2; 3) obtained from equation (36) with the various g's approximated through equations (34). The behaviour pattern of the fluctuation potential in these figures can be understood in terms of the charge density associated with the quantities F 1 and F 2 , inside the regions ω 1 and ω 2 . Figure 3 shows the fluctuation potential for a planar slice passing through the centers of two positive ions of valence +1 each. The charge density contributed by the spherical region ω 1 due to the positive ion in this region is calculated using F 1 [equation (37)], which is a function of g(2, 3), where the point 3 is inside region ω 1 . The positive sign in the fluctuation potential in region ω 1 is given by the sign of −g(2, 3 + ) + g(2, 3 − ). Since the charge at position 2 is positive, the second term associated with unlike charges is greater in magnitude than the  first term in F 1 causing an overall positive fluctuation potential in region ω 1 . The positive sign in region ω 2 has similar origins and thus analogous interpretations. Figure 4 shows the fluctuation potential for a positive ion (valence +1) in region ω 1 and a negative ion (valence −1) in region ω 2 . In contrast to the situation in figure 3, in this case the functions g(1, 3) and g(2, 3) in F 1 and F 2 lead to the sign of the fluctuation potential in regions ω 1 and ω 2 to be opposite to the signs of the ions 1 and 2, respectively. To see this, we first look at the fluctuation potential in region ω 1 calculated through F 1 with the charge density given by −g(2, 3 + ) + g(2, 3 − ). As the ion in region ω 2 is negative, the first term associated with this unlike charge dominates giving an overall negative sign to the fluctuation potential in region ω 1 where the positive ion is located. On the other hand, the fluctuation potential in region ω 2 is calculated using F 2 where the charge density is given by −g(1, 3 + ) + g(1, 3 − ). The second (positive) term here is the larger one in magnitude again being linked to the unlike charge, and hence the positive sign of the fluctuation potential in region ω 2 . So, it can generally be stated that the fluctuation potential for like ions near the vicinity of these ions is of the same sign as that of the physical ions and is of the opposite sign for unlike ions. This peculiar behavior is a consequence of the fluctuation potential in ω 1 being related to the g(2, 3) centred at 2, and that the fluctuation potential in region ω 2 being related to the g(1, 3) centred at the opposite region ω 1 . This combined with the relative magnitudes of the g's in functions F 1 and F 2 explain the behavior of the polarities in Φ (1, 2; 3).

43801-9
The magnitude of the Φ(1, 2; 3) that we have noted in the course of the present calculations, is generally small, especially for large inter-ionic separations. The reasons for this can again be traced to the dominant charge density appearing in equation (36). For instance, the charge density in region ω 1 is a function of g(2, 3) where the field point 3 is in region ω 1 and the point 2 is at the center of region ω 2 , and similarly the charge density in region ω 2 is a function of g (1, 3) where the field point 3 is in region ω 2 and point 1 is at the center of region ω 1 . As the inter-ionic separation is increased, the dominant functions in F 1 and F 2 associated with the unlike ions decrease, while the g's associated with the like charges tend to 1. It is clear from equations (37) and (38) that both F 1 and F 2 tend to zero at large distances but increase at contact distances, as evident in figure 5. Significantly, the fluctuation potential for similar charges is 43801-10 seen to become quite large compared with that in figure 3. This suggests that for small separation of the ions, the fluctuation potential term becomes important in evaluating g i j . Figures 3-5 indeed show that the fluctuation potential is the largest in the immediate vicinity of ions 1 and 2.
Another point regarding the fluctuation potential worthy of note is the relationship between the fluctuation potential and the electrostatic energy of the ions. In figure 3 we have ions of the same sign, and clearly the fluctuation potential manifests as an increase in electrostatic energy of the ions since the fluctuation potential is of the same sign as the ions. For ions of opposite sign as in figure 4, the sign of the fluctuation potential is opposite to that of the ion in the vicinity. This leads to a decrease in electrostatic potential energy leading to attractive inter-ionic correlation in this case. This implies, consistent with what has been known in the literature, that the sign of the fluctuation potential in the vicinity of ion 1 is mostly due to the cloud of counter ion (from ion 2) and vice versa. It can be seen further from figures 3 and 4, that the fluctuation potential increases as the separation of the ions decreases, establishing the importance of having a solution that is valid at short distances. Our results also show that the fluctuation potential increases with electrolyte concentration.

Structure and thermodynamics
In figure 6, we present the radial distribution functions obtained in this work along with the corresponding curves for the SPB and MPB theories at 1 mol/dm 3 concentration. It is clear that the curves are very similar for distances larger than 2 ionic diameters. Importantly, the present results and the MPB results are almost identical. The contact values for the radial distribution functions for like ions, from the present theory, are slightly closer to the MC result [9] than that from the SPB and MPB. This is probably due to a better treatment of the fluctuation potential in this work. Table 1 shows contact values g i j (1) and for comparison purposes, the corresponding results from the SPB, the MPB, and the MC [9,10] data are also included. The contact values from the present theory are consistent with the other 43801-11  Table 1. Contact values of the radial distribution functions g i j (1) from different theories. The common diameter of the ions is d = 4.25 × 10 −10 m, the temperature T = 298 K, and the dielectric constant of the electrolyte ǫ r = 78.5. The MC values are from reference [9]. theories and show a very good agreement with the MC simulation data. Tables 2 and 3 show reduced configurational energies, and osmotic coefficients from the Debye-Hückel, SPB, MPB, and MC [9,10], and this work. These values are also presented in a graphic form as in figures 7 and 8, respectively. The reduced configurational energy curves ( figure 7) show an excellent agreement between the MPB and this work with the MC curve up to 1 mol/dm 3 concentration. At the highest 1.968 mol/dm 3 concentration, the MPB is a little closer to the MC. Figure 8 shows osmotic coefficients for the theories and the relevant MC data [9,10]. These curves show a generally very good agreement between the MC results and the theories.

Conclusions
In this study we have made an analysis of the fluctuation potential in the modified Poisson-Boltzmann theory of bulk electrolyte solutions. An approximate analytical solution of the fluctuation potential equation was obtained for symmetric valence 1:1 electrolytes in the RPM. This solution was later utilized to obtain structural and thermodynamic descriptions of the electrolyte in terms of ion-ion radial distribution functions, reduced excess energy, and the osmotic coefficients, respectively.

43801-13
The fluctuation potential is a central ingredient in a potential approach to the theory (of charged fluids) such as the modified Poisson-Boltzmann theory. The fluctuation potential solution developed in this work, albeit with approximations to make analytical progress and for symmetric 1:1 valence systems, is a preliminary attempt to assess the implications of such a solution. In such cases, due to the linearization of the fluctuation potential in the bulk region I [equation (35)] and the small magnitude of Φ(1, 2; 3), the P function in bulk region I can be taken to be zero, thus neglecting charge density for that region. A less approximate and nearly full treatment could be achieved by solving for the fluctuation potential in region I using equation (29) with P I being given by equation (35) in conjunction with equation (34). An intermediate procedure (between the above two situations) to obtain a better, viable, and still feasible approximation for Φ(1, 2; 3) in region I would be to solve equation (30) [with P I given by equation (35) where the quantity C contains the valencies Z + , Z − , and has spatial dependence through g(1, 3) and g(2, 3). Thus, although C is not a constant per se, it can be assumed to be approximately constant for the purposes of solution to equation (42). An approximate analytic form of Φ(1, 2; 3) in region I, whose value is not necessarily zero, would then be available. Equation (42) has some parallels to a similar equation for the fluctuation potential in the MPB formalism in the planar electric double layer [34]. Such a procedure will be useful for higher and multivalent electrolytes when the magnitude of the fluctuation potential in region I is likely to be significant and hence P I can no longer be neglected. This will be a focus of our future work. The MPB description of the electric double layer phenomenon is an area where the present techniques might have some significance since the fluctuation potential plays an equally important role in the theoretical framework for the inhomogeneous fluid at the interface. In the MPB approach to the double layer theory in planar [34,35], cylindrical [36][37][38], and spherical [39,40] symmetries, the form of the corresponding fluctuation potential used is rather approximate and generally suffers from similar defects as those vis-a-vis the traditional MPB theory for the bulk. The statistical mechanical methods used in this paper are quite general and can be extended and adapted to interfacial double layer geometry where an analogous fluctuation potential analysis might prove useful.
Another area of possible relevance for this study is in the theoretical analysis of charged fluid systems with a variable dielectric constant (relative permittivity). The topic has attracted a lot of recent research attention (see for example, references [41][42][43]) and has been shown to be relevant for important technological systems, viz., super-capacitors [44,45]. In the electric double layer, the MPB has been found to be capable of dealing with systems having an inhomogeneous dielectric constant [34,35]. Very recently, the MPB was applied to a double layer system with three different dielectric constants [46], although the associated fluctuation potential problem could only be solved for point ions. Thus, again a fluctuation analysis in such situations along the lines of the present work could be valuable.
The three-dimensional plots of the fluctuation potential give a valuable insight into the correlations between ions. Furthermore, the present structural and thermodynamic results point in the right direction and are indicative of the potential usefulness of a full solution of the fluctuation potential. The radial distribution functions, especially at contact distances between the ions, the reduced excess energy, and osmotic coefficients show an expected improvement over that from the PB (or SPB), and overall, tend to be in a very good agreement with the predictions from the traditional MPB theory and Monte Carlo simulations.
The fluctuation potential problem is a challenging one. A complete solution of the fluctuation potential equation, valid for a general case and for asymmetry in ionic size and/or valence will involve a numerical solution comprising an iterative algorithm. Our solution here might prove useful in such an involved procedure. Such a project is contemplated in the near future.