A new family of bridge functions for electrolyte solutions

We present a new set of closures for restricted models of electrolyte solutions at the McMillan-Mayer level that improve upon the Hypernetted Chain prediction for the ion-ion pair correlation functions. The improvement is accomplished by proposing simple functional forms for the bridge functions and the specification of certain adjusting parameters according to several criteria. Under the new closures, and unlike the HNC case, the “sum” direct correlation function, which is crucial for determining the stability of the solution with respect to phase separation, remains finite at thermodynamic states along the spinodal and at the critical point.


Introduction
The accurate prediction of the structure and thermodynamic properties of electrolyte solutions over a wide range of thermodynamic states has been the focus of an enormous number of studies.In spite of these efforts, there are still a number of issues that require more development, especially those concerning the description of the structure and thermodynamic properties of an electrolyte solution close to criticality (i.e., the critical point associated with the unmixing of an ionic solution into two liquid phases of different electrolyte concentration).Another relevant issue is the need of an improved description of the structure and thermodynamic properties for aqueous solutions of 2-2 electrolytes (and also for electrolytes of higher charge) at low concentrations.
The vast majority of the efforts have focussed on the calculation of the structure and other properties of electrolyte solutions by statistical mechanical methods at the McMillan-Mayer level [1]: integral equations and computer simulation.In the first approach the ion-ion pair correlation functions g jl (r) are calculated by solving an appropriate set of Ornstein-Zernike (OZ) equations, that is supplemented with a set of additional relations between the indirect {h jl (r)} and the direct {c jl (r)} correlation functions [2,3].Because of their role in setting the mathematical problem determinate, the additional relations are referred to as closure relations or integral equation closures.
Among the closure relations proposed for ionic systems, the hypernetted chain (HNC) closure [2,3] has probably been the closure most thoroughly investigated.The HNC closure predicts the structural and thermodynamic properties of 1-1 electrolytes quite accurately up to concentrations greater that 1 mol dm −3 .For 2-2 electrolytes, however, this closure predicts a spurious peak in the like-sign (i.e., cationcation or anion-anion) correlation functions, which is absent in the correlation functions extracted from computer simulation calculations [4].At the same time, results obtained from molecular dynamics computer simulations indicate that the height (g max 12 ) MD of the peak in the cation-anion correlation function is considerably larger than the corresponding quantity (g max 12 ) HNC obtained with the HNC closure [5,6].A somewhat less known problem is that the direct correlation functions c jl (r) under the HNC closure become long-ranged as the state of the solution approaches the critical point (or, more generally, the stability limit of the solution with respect to phase separation) [7].
Recently we proposed an integral equation closure for general models of binary electrolytes that gives direct correlation functions c jl (r) that remain of finite range at every thermodynamic state [7].The new closure is formulated by decomposing the density fluctuation of each ionic species into contributions to two mutually orthogonal types of fluctuations: the charge density fluctuations and the neutral density fluctuations.The source in the ion-ion direct correlation functions responsible for the long-range behaviour close to the critical point is eliminated in the closure by appropriately choosing the bridge function.The particular form of the bridge function, which is reminiscent of the Percus-Yevick (PY) type bridge function [see equation (18) below], explains the name, HNC/PY, assigned to the closure.A shortcoming of the proposed HNC/PY closure is that it systematically gives (g max 12 ) HNC/PY smaller than (g max 12 ) HNC and (g max 12 ) MD .The purpose of this paper is to propose and investigate a new family of integral equation closures for electrolyte solutions that are constructed using the HNC/PY closure as the starting building block.The closures are obtained by simply extending the bridge function associated with the HNC/PY closure with supplementary terms of carefully selected form.In particular, the functional forms of the additional bridge function terms are extracted from the asymptotic behaviour of the direct correlation functions under the HNC/PY closure.Clearly this strategy guarantees the property (shared with the HNC/PY closure) that the new closures remain of finite range close to the critical point.
Although the HNC/PY closure was originally formulated for general solution models of binary electrolytes at the McMillan-Mayer level [7], here we confine the presentation to the simpler case of solutions of restricted or symmetric binary electrolyte models (see the next section for details).Furthermore, with respect to the performance of the new closures, in this work we consider only temperature and concentration states where the electrolyte solution is stable.
The outline of the rest of the paper is as follows.In section 2 we discuss some theoretical features of the HNC and the HNC/PY integral equation closures.Unlike in the original presentation [7], here we use the language of the sum and difference ion-ion correlation functions that is so convenient for problems involving symmetric binary electrolytes.In this section we also introduce the new integral equation closures HNC/λ 1 and b D /λ 1 by adding new terms to the bridge functions in the HNC/PY closure.In section 3 we report calculations that test the performance of the HNC/λ 1 closure for models of aqueous solutions of a 2-2 electrolyte at 298.15 K.The results are compared with the results obtained from calculations with the HNC and HNC/PY closures, as well as with the results obtained from molecular dynamics computer simulations of the same solution models by Smith et al. [5].Finally, in section 4 we briefly discuss the results obtained under the HNC/λ 1 closure, and suggest alternatives for further improvement (the b D /λ 1 closure and another extension, the HNC/λ 2 closure).We also compare our approach with previous work on the development of bridge functions appropriate for electrolyte solutions.

Theory
Our starting point is the well known McMillan-Mayer level expression [2,3] for the ion-ion correlation functions, where the indices j and l refer to the ionic species: 1 for the cations and 2 for the anions.As usual, β = (k B T ) −1 has the meaning of the inverse of the temperature expressed in energy units (k B is the Boltzmann constant).We recall that the indirect correlation function h jl (r) is related to the ion-ion radial distribution function g jl (r) [which is proportional to the probability that two ions, one of species j and the other of species l, are separated by a distance r in the ionic solution] according to the equation h jl (r) = g jl (r) − 1.The other functions in equation ( 1) are ūjl (r), the solvent-averaged potential between two ions of species j and l separated by a distance r at infinite dilution [1,2]; t jl (r) ≡ h jl (r) − c jl (r), in which c jl (r) is the ion-ion direct correlation function; and finally the ion-ion bridge function b jl (r) [2,3].The sets of functions {h jl (r)} and {c jl (r)} are further interrelated through (see below) the Ornstein-Zernike equations [2,3].
As mentioned in the introduction, in this study we focus on the case of a symmetric (or restricted) model of an electrolyte solution.In such a model the electrolyte dissociates into equal numbers of cations and anions.Furthermore, the cation and anion species are entirely equivalent, except for their charges z 1 e and z 2 e, which are of the same magnitude but of opposite signs (z 1 = −z 2 ; e is the charge of a proton).Quite generally, then, we write the McMillan-Mayer level ion-ion solvent averaged potentials ūjl (r) as the sum of a short-range contribution ū * (r), that is the same for every jl ionic pair-type, and a long-range Coulombic contribution whose sign depends on whether the labels j and l refer to the same or to different ionic species: It should be noted that in the McMillan-Mayer level description, the Coulomb interaction between the charges is attenuated by the dielectric constant ǫ 0 of the pure solvent.A particular but important case of equation ( 2) is the restricted primitive model, for which ū * (r) is a hard sphere potential.In this study, however, the discussion will not be limited to that model, but instead it will consider more general forms of ū * (r).
For the restricted electrolyte models represented by equation ( 2) it is advantageous to express any ion-ion function F jl [such as any of the functions that occur in equation ( 1)] in terms of the "sum" (subscript S) and difference (subscript D) functions: With this notation, equation (1) for jl =11, 12, 21, and 22 may be cast in the more succinct form where we have introduced the auxiliary function Notice that because ū * (r) is the same function for every jl pair, we have used in writing equations ( 4) and (5) that ūS (r) = ū * (r), the short-range part of the ion-ion potential energy of interaction.Correspondingly, we note that the nature of ūD (r) is purely Coulombic: ūD (r) = z 2 e 2 /ǫ 0 r .Equations ( 4) and (5) may also be rearranged to display the sum and the difference direct correlation functions explicitly: In terms of the Fourier transforms of the sum and difference functions, the Ornstein-Zernike equations at the McMillan-Mayer level are where ρ t = ρ 1 + ρ 2 is the total number density of ions in the solution (for the restricted electrolyte systems discussed here ρ 1 = ρ 2 applies).We recall that the sum/difference decoupling suggested by equations ( 9) and ( 10) is only apparent, as the set of functions {h S , c S } and {h D , c D } are still entangled by the relations (7) and (8).

HNC closure
In the set of equations ( 7)- (10), which provide the means for the calculation of the structure functions of the electrolyte solution, we still have to specify the sum and difference bridge functions b S (r) and b D (r).To do so is one of the goals of this study.In this section, however, we discuss the closure relations that result from the simplest possible approximation for these functions, namely when both bridge functions can be ignored (i.e., b S (r) = b D (r) = 0).This is the well-known Hypernetted Chain (HNC) approximation, for which equations ( 4) and ( 5) simplify to while the auxiliary function D(r) [cf.equation ( 6)] is given by It is interesting to analyze here a shortcoming of the HNC closure that is not widely appreciated.This may be seen most straightforwardly from equations (11) and (12).These equations can be formally solved for the function D 0 (r): in terms of the auxiliary function Equation ( 14) for D 0 (r) can then be substituted back into equation (11).Recalling that in the resulting equation t S (r) = h S (r) − c S (r) in the argument of the exponential, we can then solve for the sum combination c S (r) of the ion-ion direct correlation functions.We obtain, for values of r large enough that we can make the approximation exp(−β ū * (r)) ≃ 1, The particular form of the second term in this equation follows straightforwardly from the identity cosh arctanh The feature, revealed by equation ( 16), that c S (r) under the HNC closure is expressed by a sum of two separate contributions, one that depends only on h S (r) while the other depends only on χ(r), is the source of difficulties when the closure is used to examine the material stability limit of the electrolyte solution.We recall [8] that the stability of the solution against phase separation (unmixing) is determined exclusively by the sum combination c S (r) of the direct correlation functions.(Briefly, the solution is stable or metastable against phase separation as long as 1−ρ t c S (0) 0, where c S (0) denotes the Fourier transform of c S (r) specialized at k = 0).On the other hand, the signature of the states on the spinodal line (and, of course, at the critical point) is that at those thermodynamic states the sum combination h S (r) becomes infinitely long ranged (or, equivalently, h S (k) diverges as k → 0 at these states).In contrast, the range of χ(r), which is proportional to the difference combination h D (r), remains finite at the stability limit of the electrolyte solution [9].
At large values of the interionic distance r, where both h S (r) and χ(r) are small, we can approximate c S (r) in equation ( 16) by the first few terms of its power series expansion with respect to h S (r) and χ(r): When the thermodynamic state approaches the stability limit and, correspondingly, h S (r) becomes more and more long-ranged, equation (17) reveals that, under the HNC closure, the sum combination c S (r) (due to the presence of terms with "pure" [h S (r)] n powers) also becomes infinitely long-ranged .

HNC/PY closure
As a first step in designing an integral equation relation that improves over the HNC closure, it is natural to introduce into the closure a bridge function b S (r) that eliminates from c S (r) the terms that contain pure [h S (r)] n powers of the sum combination function h S (r).This can be accomplished most simply by choosing b S (r) of the form [7] b We still insist on the approximation b D (r) = 0.When these approximations are introduced into the general equations ( 4) and ( 5) they become with D 0 (r) given by equation ( 13).
In the next section we show that the choice of equation ( 18) for the sum combination b S (r) bridge function makes c S (r) under this closure short-ranged (which is the correct behaviour) at the critical point or any other state along the spinodal line of the electrolyte solution.
Equation ( 18) is the well-known form that the bridge function for a pure substance takes under the Percus-Yevick (PY) integral equation closure.In fact, the PY flavour of this closure becomes immediately obvious in equation ( 19) if we ignore the charge-fluctuation factor cosh D 0 (r).On the other hand, through the choice b D (r) = 0 for the difference bridge function combination, this closure retains all the features of the HNC closure in what concerns the difference combination direct correlation function c D (r).For these reasons we use the acronym HNC/PY to refer to this closure.A version of this closure for general models (i.e., not restricted) of electrolyte solutions has been reported [7].
As the calculations reported in section 3 show, the performance of the HNC/PY closure for dilute solutions of 2-2 restricted electrolytes is somewhat disappointing, with the HNC closure giving (broadly speaking) better results.The goal of this paper is, using the HNC/PY closure as the foundation, to propose a new family of improved integral equation closures for solutions of symmetric electrolytes.

Beyond the PY form for the sum bridge function b S (r)
To obtain improved results but, at the same time to preserve the feature that c S (r) is of finite range at every thermodynamic state, we slightly modify the PY form of the sum combination b S (r) of the bridge functions [cf.equation ( 18)].We do this by adding a corrector function λ(r) to the argument of the logarithmic term: For generality we do not assume at this point that b D (r) = 0, but instead we leave this function momentarily unspecified.The general equations ( 4) and ( 5) take the form where D(r) was introduced in equation ( 6).
We follow now the discussion in section 2.1.We begin by solving equations ( 22) and (23) for the auxiliary function D(r): where χ(r) is given in equation (15).At interionic distances such that exp(−β ū * (r)) ≃ 1, equations ( 22) and (24) may be solved for the sum combination c S (r) of direct correlation functions, with the result c S (r) = λ(r) To derive this equation we have used the identity Equation (25) shows that the contribution λ(r) coming from the sum bridge function b S (r) is separated from the remaining terms.Except for the multiplicative factor (1 + h S (r)), these terms depend only on the function χ D (r), which is shortranged at every thermodynamic state.Furthermore, the equation does not explicitly involve the bridge function b D (r): this expression of c S (r) in terms of λ(r), h S (r) and χ(r) is not affected by whichever approximation we choose for b D (r).
At large values of r, where χ(r) is small, we can approximate equation (25) for c S (r) by in which h S (r) does not appear explicitly Leaving aside for a moment the function λ(r), equation ( 26) shows that at large r the sum combination c S (r) involves terms from which pure powers [h S (r)] n are absent; such powers of h S (r) (that originate from the expansion of the denominator of χ(r)) are always multiplied by powers [h D (r)] m of the difference combination function h D (r).The latter function is of finite range at every thermodynamic state [9].When λ(r) = 0, equation (26) coincides with the large-r behaviour of c S (r) under the HNC/PY closure (cf.section 2.2).This demonstrates the statement made in that section that c S (r) under the HNC/PY closure is of finite-range at the stability limit of the ionic solution.
In view of the extreme complexity involved in the direct calculation of λ(r), here we simply consider this function to be a functional of the functions χ(r) and h D (r) (or, equivalently, given equation ( 15), a functional of h S (r) and h D (r)).Furthermore, we propose that this functional has a form that closely resembles the first terms in the asymptotic expansion of c S (r) under the HNC/PY closure [which is given by equation (26) after making λ(r) = 0].With these assumptions it follows that the simplest candidate for λ(r) has the form where the value of the bridge parameter B S needs to be specified.Replacing this form of λ(r) back into equation ( 26) gives the asymptotic behaviour of the sum direct correlation function c S (r) at large r with the parameter A S is related to B S by It is clear that this approach for constructing λ(r), that borrows from the asymptotic expansion of c S (r) under HNC/PY, guarantees that c S (r) under the new closure [that incorporates λ(r)] remains of finite range at the stability limit of the electrolyte solution.It is also worth noting that the strategy just outlined for the construction of λ(r) is completely independent on which approximation we choose for the difference bridge function b D (r).

Approximation for the difference bridge function b D (r)
As in the case of λ(r), we propose to construct an approximation for the difference bridge function b D (r) by considering the first few terms of the asymptotic form of c D (r) under the HNC/PY approximation.
From equations ( 6) and ( 24) we find for c * D (r) ≡ c D (r) + β ūD (r) (the short-range part of the difference direct correlation function) The last term is, of course, equivalent to arctanh { χ(r) }.This expression is explicitly linear in h D (r).Expanding the logarithmic term in powers of χ(r) we obtain The second and third terms can be combined together, taking into account equation (15), so that the large-r asymptotic expansion of c * D (r) becomes It is important to note that this result applies irrespective of any approximation made to b S (r).By taking b D (r) = 0 in equation ( 31) we obtain the asymptotic expansion of c D (r) under any closure that approximates the difference bridge function to zero, like the HNC or HNC/PY closures.We may generate an approximate expression for b D (r) by proceeding as we did for constructing λ(r) in the last section; i. e., we assume that b D (r) is proportional to the first term of the asymptotic expansion of c * D (r) under the HNC/PY closure.We then have b in terms of the bridge parameter B D .With this approximation c * D (r) [equation ( 31)] behaves at large r as where A D = 1+B D .Clearly, this expansion does not contain terms with pure powers [h S (r)] n ; such powers appear to be always multiplied by powers [h D (r)] m that are of finite range for every thermodynamic state.Finally, it is straightforward to show from equations ( 28) and ( 33), together with the expected small-k behaviour of the Fourier transforms of h S (r) and h D (r) that the direct correlation functions c S (r) and c D (r) under the closures proposed in this section have the form that is required for the Stillinger-Lovett second moment condition [10] to be satisfied.

Closure relations
The development in section 2 motivated specific forms for the sum b S (r) and difference b D (r) bridge functions [equations ( 21) and ( 27) for b S (r); equation (32) for b D (r)].Having this information, we may think of several approximate implementations of equations ( 7) and ( 8): HNC closure.In this closure we have b S (r) = b D (r) = 0 at every r.The sum and difference direct correlation functions c S (r) and c D (r) are given by relations similar to equations ( 7) and ( 8), but with b S (r) absent from the argument of the exponential, and with D 0 (r) in place of D(r) [cf.equations ( 13) and ( 6)].
HNC/PY closure.In this closure λ(r) = b D (r) = 0 at every r.Hence c S (r) and c D (r) are given by expressions similar to equations ( 7) and ( 8), but with b S (r) given by equation ( 18) and, again, with D 0 (r) in place of D(r).
HNC/λ 1 closure.In this closure b S (r) in the argument of the exponentials in equations ( 7) and ( 8) is given by equation ( 21 b D /λ 1 closure.In this closure c S (r) and c D (r) take the form of equations ( 7) and ( 8), with b S (r) given by equation ( 21), λ(r) given by equation ( 27), and b D (r) given by equation (32).
It is obvious that the implementation of closures HNC/λ 1 and b D /λ 1 require the specification of the bridge parameters B S and B D .We address this important issue in the case of the HNC/λ 1 closure in the next subsection.
In what remains of section 3 we report numerical results for the structure and thermodynamic properties of dilute solutions of a model 2-2 electrolyte as calculated with the HNC, HNC/PY, and HNC/λ 1 closures.We briefly comment on closure b D /λ 1 in section 4.

Numerical algorithm, electrolyte model, and calculations
To solve the Ornstein-Zernike integral equations [cf.equations ( 9) and ( 10)] under any of the above closure relations we have used the familiar Picard iterative approach.The technical difficulties associated with the long-range nature of the difference direct correlation function c D (r) are avoided by means of Ng's method [11,12].
As a representative of a restricted model of an electrolyte solution, we consider the model of a 2-2 electrolyte studied by Duh and Haymet [6].In this model the ions interact through a McMillan-Mayer level potential of the form of equation ( 2), with ū * (r) given by ū * where B = 5377(ze) 2 Å K and σ = 2.8428 Å.In the Coulombic term ūD (r) we specify z = 2 and ǫ 0 = 78.358, the value of the dielectric constant of water at 298.15 K.The calculations reported here for this model electrolyte parallel those in the study by Duh and Haymet [6].We have solved each of the integral equation closures summarized in the previous section for electrolyte solutions at several concentrations in the range 0.001 mol dm −3 c s 0.5625 mol dm −3 , where c s is the molar concentration of the electrolyte.
In addition to the structure functions [the cation-cation g 11 (r) and the cationanion g 12 radial distribution functions, obtained from g S (r) and h D (r) by inversion of equations ( 3)], we have calculated the following thermodynamic properties: (a) the reduced excess internal energy per ion u ≡ βU ex /N t : (b) the osmotic coefficient φ ≡ βπ/ρ t of the solution: (c) the generalized compressibility: In the above equations π is the osmotic pressure of the solution.It may be shown that the generalized compressibility equals where µ s and γ ± are, respectively, the chemical potential and the mean activity coefficient of the electrolyte, while ρ s is the number density of the salt (ρ s = ρ t /ν, where ν = ν 1 + ν 2 is the sum of the cation and anion stoichiometric coefficients).
It is important to remember that equations (37)-(39) correspond to different thermodynamic routes [energy route for u, virial route for φ, and compressibility route for β(∂π/∂ρ t ) T ], and when implemented using correlation functions g S (r) and h D (r) calculated from approximate integral equation closures, their results are not necessarily thermodinamically consistent.

Results
As a representative of the results obtained at each concentration examined, we first discuss in some detail the results obtained under the HNC, HNC/PY, and HNC/λ 1 closures for the model electrolyte solution at concentration c s = 0.02 mol dm −3 at temperature 298.15 K.In figure 1 we report the pair correlation functions g 11 (r) and g 12 (r) at this concentration as calculated under the HNC and the HNC/PY integral equation closures.The first thing to notice is that (g max 12 ) MD = 48.0, the peak height for the cation-anion correlation function g 12 (r) calculated from molecular dynamics computer simulation [5], is substantially larger than the corresponding HNC and HNC/PY results, (g max 12 ) HNC = 37.6 and (g max 12 ) HNC/PY = 33.7.From this perspective, and in spite of the fact that for this closure b S (r) = 0, the HNC/PY closure does not improve over the HNC predictions.We notice, however, that there is no evidence of the spurious peak at r ≃ 8.5 Å in the like-sign ion-ion correlation function g 11 (r) calculated under the HNC/PY closure, contrary to the case for g 11 (r) under HNC closure.The simulation study [5] indicates that such a peak is an artifact of the approximate HNC closure.
Thus, except for the case of the osmotic coefficient and the absence of the peak in g 11 (r), at this concentration of the 2-2 electrolyte the structure and thermodynamic properties obtained with the HNC/PY closure differ from the simulation results by a larger extent than the predictions of the HNC closure.This pattern in the performance comparison between HNC and HNC/PY is also observed, in different latter route involves numerical computation of the indicated partial derivative, with the function βπ at each density calculated with the relation βπ = ρ t φ, in which the osmotic coefficient is obtained with the help of the (virial route)-based equation (38).
The other thermodynamic consistency test considered follows from the Maxwell relation between the reduced internal energy and the reduced chemical potential of the electrolyte [8]: Notice that the left hand side of the relation requires the full reduced internal energy per ion; the term 3/2 is just the value of the non-interacting contribution βU id /N t to the reduced internal energy.By taking the derivative of each side of equation ( 41) with respect to the total ionic density ρ t we derive, after taking into account the first equality of equation (40), It should be noted that the partial derivatives on the left hand side of this equation are carried at constant temperature, while those on the right hand side are carried at constant density of the ions.To use this equation as a of the thermodynamic consistency of closure HNC/λ 1 , we calculate each side of the equation following a different thermodynamic route.Specifically, for the left hand side we calculate u at several values of ρ t using the energy route, equation (37).For the right hand side we calculate β(∂π/∂ρ t ) T at several temperatures using the compressibility route, equation (39).The results are shown in figure 3b.
Figure 3a shows that at B S = 0, which corresponds to the HNC/PY closure, the difference between the values of the generalized compressibility β(∂π/∂ρ t ) T calculated under the compressibility and virial routes is the largest, with the value calculated by the compressibility route larger than the one obtained through the virial route.The thermodynamic inconsistency decreases, however, when B S is increased, and the figure shows that at B S ≃ 0.395 the two routes in fact give the same value for β(∂π/∂ρ t ) T .Although this is satisfying, this thermodynamic consistency is achieved at a value of B S where g max 12 is too large and, furthermore, for which the shape of g 11 (r) differs considerably from the like-sign correlation function extracted from molecular dynamics simulation (see below).
Figure 3b illustrates the performance of the HNC/λ 1 closure with respect to the second thermodynamic consistency test, equation (42).It is evident from the figure that under this closure the (energy route)-based left hand side of the relation is different at every B S from the (compressibility route)-based right hand side.Furthermore, the magnitude of the thermodynamic inconsistency is practically independent of the value of the bridge parameter B S .
Finally, in figure 4 we report the ion-ion correlation functions g 11 (r) and g 12 (r) for the 2-2 electrolyte at concentration c s = 0.02 mol dm −3 .In addition to the HNC and HNC/PY results already shown in figure 1, we also display the correlation functions bridge parameter B S in the HNC/λ 1 closure.An appealing alternative for establishing the value of B S that gives under HNC/λ 1 improved results (relative to HNC and HNC/PY) is suggested by the correlation functions g jl (r) reported in figure 4 using B S = 0.083.For this value of the bridge parameter (at c s = 0.02 mol dm −3 ) u HNC/λ 1 ≃ u HNC = −1.280.Inspection of figure 4a reveals that at this value of B S there is only a minimal trace of the spurious peak in g 11 (r), with the overall shape of the function in better agreement with the simulation result [5,6].At the same time figure 4b shows that g 12 (r), although not as good as the result at B S = 0.14, is still an improvement over the predictions of HNC and HNC/PY [(g max 12 ) HNC/λ 1 = 41.75 versus (g max 12 ) HNC = 37.6, (g max 12 ) HNC/PY = 33.7,and (g max 12 ) MD = 48.0].For the osmotic coefficient we find φ HNC/λ 1 = φ HNC/PY = 0.709, which is in better agreement with φ MD = 0.706 than the HNC result φ HNC = 0.713.
We have tested this strategy for fixing the B S parameter of closure HNC/λ 1 for solutions of the 2-2 electrolyte in the concentration range 0.001 − 0.5625 mol dm −3 .The results obtained are summarized in table 1.The second column in the table reports the value of B S for which u HNC/λ 1 = u HNC at each concentration.It is interesting that this "optimun" value of the bridge parameter does not change monotonically with the electrolyte concentration.
By following this procedure, the HNC/λ 1 results for the reduced excess energy u will be obviously equivalent to those of HNC; they are systematically smaller in magnitude than the corresponding molecular dynamics results.For the osmotic coefficient, HNC and HNC/PY are quite comparable, and the table shows that these two closures perform marginally better against simulation than the HNC/λ 1 closure.On the other hand, the last part of the table reveals that HNC/λ 1 performs considerably better than HNC and HNC/PY with respect to the height g max 12 of the peak in the cation-anion radial distribution function.Especially encouraging are the results at the smallest concentrations c s = 0.001 and 0.005 mol dm −3 .To illustrate this point, in figure 5 we compare the g jl (r) at c s = 0.001 mol dm −3 calculated under closures HNC, HNC/PY, and HNC/λ 1 with B S = 0.047.From table 1 we see that (g max 12 ) HNC/λ 1 = 231 compares very well with (g max 12 ) MD = 234.Inspection of figure 5b reveals that the peak in g 12 (r) predicted by HNC/λ 1 is slightly narrower than the peak calculated under HNC (and it is also narrower than the peak extracted from the molecular dynamics data [6]).With respect to the like-sign correlation function g 11 (r), the HNC/λ 1 closure gives, compared against the simulation results, better results than HNC/PY.Unfortunately, for B S = 0.047, at which u HNC/λ 1 = u HNC , g 11 (r) under closure HNC/λ 1 displays the spurious peak closeby to r ≃ 8 Å.The size of the peak, however, is much smaller than the peak in g 11 (r) under the HNC closure.

Discussion
Starting from the previously reported HNC/PY integral equation closure for solutions of binary electrolytes [7], we have proposed a very simple strategy for improving the closure performance concerning the structure and thermodynamic properties of the solution.Basically, the improvement is accomplished by adding new terms to the ion-ion bridge functions b jl (r) of the HNC/PY closure.The general idea is to borrow from the first few terms in the asymptotic expansion of the HNC/PY c jl (r) the specific functional form of the new terms in the bridge functions b jl (r).We note, however, that although the original HNC/PY work was concerned with general McMillan-Mayer level models of binary electrolytes, the present extension of the HNC/PY closure has been implemented in the present work only for restricted (or symmetric) models of electrolytes.In fact, by exploiting the symmetry of such electrolyte models, the discussion was actually cast in the language of the sum and difference functions introduced in equation ( 3).
The calculations presented in section 3.3 tested the predictions of the simplest extension, the HNC/λ 1 closure, for aqueous solutions of a restricted 2-2 electrolyte.For this closure, only the sum combination b S (r) of ion-ion bridge functions is modified relative to the corresponding HNC/PY bridge function; the difference bridge function b D (r) for this closure, like in the HNC and HNC/PY closures, is assumed to vanish at every r.The explicit of expression for b S (r) is given by equations ( 21) and ( 27), and contains an unspecified parameter B S .Initially we had hoped that the optimal value of B S would be fixed by demanding thermodynamic consistency for β(∂π/∂ρ t ) T when calculated through the compressibility and virial routes.While at every concentration studied we succeeded finding a B S for which the thermody-namic consistency is realized, the calculated correlation functions g jl (r) at those B S deviate substantially from the structure functions g jl (r) extracted from computer simulations.We were therefore forced to adopt the less elegant (but quite practical) criterion of choosing B S such that u HNC/λ 1 = u HNC is satisfied.
The performance of the HNC/λ 1 closure based on this criterion for fixing B S was examined in section 3.3, where we concluded that, except for the osmotic coefficient, this new integral equation provides a substantial improvement over the HNC/PY and the HNC closures.Here we notice that, while still insisting on the approximation b D (r) = 0, there is room for improving the closure relation by further extending b S (r) with a 2-parameter λ(r) function: The functional form of each term is borrowed from the respective term in the asymptotic expansion of c S (r) under the HNC/PY closure [the second term in equation (26)].We refer to this closure by the acronym HNC/λ 2 , as its implementation requires the specification of the two bridge parameters B S and B ′ S .When considered from the viewpoint of the ion-ion bridge functions b jl (r), the assumption b D (r) = 0 in the HNC/λ n family of closures implies that under these closure relations b 11 (r) = b 12 (r) applies.Ichiye and Haymet [13] and Duh and Haymet [6] have extracted the bridge functions for 2-2 electrolytes at several concentrations from computer simulations, and their results indicate that b 11 (r) and b 12 (r) have different sign, with b 12 (r) > 0. (See also Kjellander and Mitchell [14]).The bridge functions in the Ionic-Percus Yevick (IPY) integral equation closures developed by Ichiye and Haymet [6,13,15] specifically account for this feature, and are able to account for the simulation results to a better degree than the HNC/λ 1 closure described in section 3.3.It should be mentioned that the IPY closures have no free parameter.To our knowledge, the behaviour of these closures in the context of thermodynamic consistency between the various routes has not been examined.(The exponential exp[b 12 (r)] in the IPY2 closure [15] contains a term with the factor exp[−τ 12 (r)], where τ 12 (r) is an optimally designed t 12 (r) function.The presence of this exp[−τ 12 (r)] factor suggests that the range of the corresponding direct correlation function would become infinite as the electrolyte solution approaches the limit of material stability).
Only when we allow b D (r) = 0 we can have b 11 (r) and b 12 (r) of different sign.This leads us to the b D /λ 1 integral equation closure introduced in section 3.1.This closure requires the specification of the two bridge strength parameters B S and B D .The criterions for fixing these parameters and the performance of this new integral equation are the subject of ongoing research in our group.
We conclude with the mention that, for the dilute systems studied here, the results of the HNC/PY closure are practically identical to the results calculated with the Livshits-Martinov integral equation relation [16].
), with λ(r) represented by (27).The auxiliary function D(r) again takes the form of D 0 (r) [i.e., in this closure we again assume that b D (r) = 0 at every r].

Figure 1 .
Figure 1.Ion-ion radial distribution functions for an aqueous solution of a restricted 2-2 electrolyte at c s = 0.02 mol dm −3 and 298.15 K.The correlation function are calculated under the integral equation closures HNC and HNC/PY.(a): Like-sign (cation-cation or anion-anion) correlation function g 11 (r); (b): cationanion correlation function g 12 (r).r is in units of Å.

Figure 5 .
Figure 5. Ion-ion radial distribution functions for an aqueous solution of a restricted 2-2 electrolyte at c s = 0.001 mol dm −3 and 298.15 K.The correlation function are calculated under the integral equation closures HNC, HNC/PY, and HNC/λ 1 with the bridge parameter B S = 0.047, for which u HNC/λ 1 = u HNC .(a): Like-sign (cation-cation or anion-anion) correlation function g 11 (r); (b): cationanion correlation function g 12 (r).r is in units of Å.