J. Phys. Chem. 1996, 100, 1415-1419
1415
Phase Separation of Ionic Fluids: An Extended Ebeling-Grigo Approach Showguei Yeh Department of Chemistry, State UniVersity of New York at Stony Brook, Stony Brook, New York 11794-3400
Yaoqi Zhou Department of Chemical Engineering, North Carolina State UniVersity, Raleigh, North Carolina 27695-7905
George Stell* Department of Chemistry, State UniVersity of New York at Stony Brook, Stony Brook, New York 11794-3400 ReceiVed: August 21, 1995; In Final Form: October 16, 1995X
The Ebeling-Grigo (EG) improvement of the mean-spherical approximation (MSA) via the Bjerrum theory of ionic association in the restricted primitive model (RPM) is reviewed and extended. We find that, in terms of reduced density and temperature, the critical point (F* c, T* c) of the EG approximation (EGA) is located at (0.026, 0.0789), rather than the (0.018, 0.0837) quoted by EG. The corrected result is closer to the most recent simulation estimates of the RPM critical point. A further improvement is obtained by incorporating the activity coefficient of the associated ionic pair to yield a dipole-ion (DI) contribution. The new approximation, EGA+DI, has a critical point at (0.027, 0.0716), which is in still better agreement with the simulation estimates. A critical discussion comparing these results with recent simulation results and other recent theoretical work is given.
I. Introduction The mean-spherical approximation (MSA) has been widely used in ionic systems because it provides simple analytical solutions and correctly treats the hard-core condition.1 However, the MSA is a linearized approximation which does not take into account ionic association. Various methods of improving it have been suggested.2-4 One of the methods, due to Ebeling and Grigo4 is to explicitly include the ionic pairing contribution via the second ionic virial coefficient5 and the mass-action law, using an elegant extension of the Bjerrum (BJ) theory6 of ionic association. In an earlier paper,7 we improved the MSA via a new approach, similar to the approach of Ebeling and Grigo4 but using our reference cavity theory of association8 adopted to ionic association instead. In the simplest version of our theory,7 which we call the pairing MSA 1 (PMSA1), we neglect the activity coefficient of the fully associated ion pairs, which are regarded as a separate dipolar species, and obtain the critical point (F* c, T*c) at (0.025, 0.0748). Here, F* ) 2F0σ3 and T* ) kBTσ/(ze)2 with σ, the ionic diameter, , the solvent dielectric constant, z, the ionic charge number, e, the electronic charge, T, the temperature, 2F0, the total number density, and kB, the Boltzmann constant. In a second PMSA (or PMSA2), we include the activity coefficient of these dipolar particles at the MSA level. The new critical point is located at (0.023, 0.0733). In a third PMSA (or PMSA3), we further include the effect of the presence of the dipolar-particle cores. The final critical point is located at (0.025, 0.0745). These critical points are considerably closer than the MSA result3,9 (0.014, 0.079) to the most recent Monte Carlo estimates of F* c from 0.025 to 0.04 and T* c from 0.053 to somewhat over 0.057.10 Conceptually our PMSA1 is on much the same level of approximation as the Ebeling-Grigo approximation (EGA),4a which also neglects the activity coefficient of the fully associated ionic pairs except that the EG theory uses a modified Bjerrum theory of ionic association. However, somewhat surprisingly, not only is the critical point value (0.018, 0.0837) that EG give X
Abstract published in AdVance ACS Abstracts, December 15, 1995.
0022-3654/96/20100-1415$12.00/0
in considerably worse agreement with simulation than the PMSA1 critical point,7 but their critical temperature is in even worse agreement than the MSA’s critical temperature (0.079).9 In order to understand this difference, we have recalculated the critical point using their approximation. We find that the critical point of the EGA is in fact located at (0.026, 0.0789). The critical point is thus not far from the PMSA1 result (0.025, 0.0748), which is what one would anticipate on conceptual grounds. (A similar reassessment of EG theory has also been made independently by Guissani and Guillot.11) We have gone on to improve the EGA by incorporating the activity coefficient of associated ionic pairs at the MSA level to obtain a dipoleion (DI) contribution in much the same spirit as that of Fisher and Levin in their Debye-Hu¨ckle level of theory.12 The new approximation, EGA+DI, further improves the agreement with simulations. The critical point of EGA+DI is located at (0.027, 0.0716). In Section II, the EGA is reviewed and the new approximation (EGA+DI) is developed. In Section III, a critical discussion is given in which our recent results are compared with other approximations and simulation results. II. Theory A. Ebeling-Grigo Approximation. Consider the restricted primitive model (RPM), which is a classical system of N charged hard spheres, N/2 cations and N/2 anions with the same hardcore diameter σ and opposite charge numbers (+z, -z). In the EGA, the pressure (P) and Helmholtz free energy (A) consist of three parts: ideal (ID), hard core (HC), and electrostatic (EL), calculated by applying the MSA to the nonassociated ions, where the degree of association is determined via a law of mass action. In this case, the pressure and the free energy can be written as
βP ) β(PID + PHC + PEL)
(2.1)
βA/N ) β[AID + AHC + AEL]/N
(2.2)
where β ) 1/kBT. © 1996 American Chemical Society
1416 J. Phys. Chem., Vol. 100, No. 4, 1996
Yeh et al.
The ideal part after association is
βP ) 2F0 - RF0 ID
(2.3)
βAID/N ) βA°/N - (1 - R){1 - ln[(1 - R)F*]} (R/2){1 - ln[(RF*)(2σ3/K0)]} (2.4) where R is the degree of association which satisfies the law of mass action:
γ+γR ) K0 2 γM+F0(1 - R)
(2.5)
Here γ+ and γ- are the activity coefficients of free (unassociated) cations and anions, respectively, γM+- is the activity coefficient for associated ionic pairs, F0 is total number density of cations or anions, including those constituting the associated pairs, and K0 is the ideal association constant. The EGA evaluates γ(, the mean activity coefficient of free ions, in the MSA and approximates the γM+- as 1. EGA EL γ( ) γ(
(2.6A)
EGA )1 γM +-
(2.6B)
As a result, the degree of association, R0, in the EGA can be obtained from
R0 F0(1 - R0)
EL 2 ) K0[γ( ] 2
(2.7)
The EGA hard-core part is obtained from the Carnahan and Starling14 equation, which is known to be very accurate for the pressure and free energy of hard spheres.
βPHC ) 2F0 βAHC/N )
2η(2 - η) (1 - η)
∑
b2m
m)2(2m)!(2m
(2.9)
(1 - η)2
]
1/2 z2λ 1 + x - (1 + 2x) σ x
(2.10)
βA/N ) β[AID + AHC + AEL + ADI]/N
x ) κσ, k2 ) 4πβFi(ze)2/
(2.11)
and Fi is the number density of free cations and anions. The MSA also gives1
βAEL/N ) -
((1 + 2x)1/2 - 1)3 24πσ3
3x2 + 6x + 2 - 2(1 + 2x)3/2 24πF0σ3
(2.12)
βADI/N )
Fd ln γM+2F0
(1 + η + η2 - η3) (1 - η)3
-
((1 + 2x)1/2 - 1)3 24πσ3
(2.17)
with Fd the number density of ion pairs and γM+- the ion-pair activity coefficient at the infinite dilute limit of the pair. Here, we shall use an accurate expression for the dipolar-dumbbell activity coefficient γM+- obtained in a previous paper7 of ours
ln γM+- )
[
]
1/2 2 z2λ (2x + 2)(1 + 2x) - x - 4x - 2 σ x2
(2.18)
where x is defined in eq 2.11. Equation 2.17 provides an additional contribution to the mean activity coefficients of free ions DI ln γ( )
[
]
∂(βADI)/V ∂Fi
(2.19) Fd
[
]
∂(-βADI) (2.20) ∂V T As a result, the equation of state for our new approximation, EGA+DI, is (1 + η + η2 - η3) ((1 + 2x)1/2 - 1)3 βP ) 2F0 + (1 - η)3 24πσ3 βPDI - R1F0 (2.21) with the degree of association, R1, given by EL 2 DI 2 (γ( ) (γ( ) ) K 0 2 γM+F0(1 - R1)
(2.22)
where the ideal association constant K0 is the same as in the EL DI , γ( , and γM+- are obtained from eqs EGA (eq 2.15) and γ( 2.10, 2.19, and 2.18, respectively. Here, we shall continue to assume that the hard-core part is independent of the degree of ionic association. III. Quantitative Results and Discussion
(2.13)
The final equation for pressure of the EGA is
βP ) 2F0
(2.16)
where the ion-pair/ion contribution is approximately obtained from
R1
where
βPEL ) -
(2.15)
- 3)
where b ) β(ze)2/(σ) is the Bjerrum length. B. Improvement over the Ebeling-Grigo Approximation. One natural way to improve the EGA is to include the effect of the ion-pair/ion interaction by taking into account the ion-pair (dipole) activity coefficient.7,12 In this case, the free energy can be rewritten as
βPDI )
η(4 - 3η)
[
∞
K0 ) 8πσ3
and also the contribution to the total pressure
(2.8)
3
where the packing fraction η ) πF0σ3/3. The hard-core part is assumed to be independent of the degree of association. From the MSA we have1 EL ln γ( )-
eq 2.7 and the ideal association constant K0 in eq 2.7 is chosen so that eq 2.14 gives an exact second ionic virial coefficient,13
- R0F0 (2.14)
where R0, the degree of association of EGA, is calculated from
A. Quantitative Results. To investigate the critical phenomena, it will be easier to deal with the reduced parameters. The reduced pressure is defined as
P* )
σ4 P ≡ φF*T* (ze)2
(3.1)
with T*, the reduced temperature ()1/b), F*, the reduced density ()2F0σ3), and φ, the osmotic coefficient [)βP/(2F0)].
Phase Separation of Ionic Fluids
J. Phys. Chem., Vol. 100, No. 4, 1996 1417
TABLE 1: Osmotic Coefficient for 2-2 Electrolytes Calculated in Different Theories (σ ) 4.2 Å, E ) 78.358, T ) 298.16 K) EGA C (mol/L)
R0
ref 4a
this work
MSA
HNC15
0.0625 0.25 0.5625 1.0 2.0
0.3826 0.3915 0.3687 0.3463 0.3193
0.6741 0.5907 0.5965 0.6654 0.9136
0.6120 0.5740 0.6024 0.6810 0.9342
0.6303 0.5543 0.5732 0.6520 0.9112
0.6396
a
a
0.5966 0.6586 0.8733
Our values of R0 agree with those of EG.4a Our R0 is their 1 - R.
Nonlinear algebraic equations of R0 and R1 and the derivatives in the EGA+DI approximation have been evaluated numerically by us using IMSL routines. When the association term in EGA is turned off, i.e., R0 ) 0, we recover the simple MSA results as expected.1 However, we are unable to reproduce the pressure and the critical point of the EGA obtained in ref 4a. In Table 1, we compare our EGA result for osmotic coefficients with the results obtained in the original EG paper;4a our values of the degree of association agree with theirs. The agreement between the HNC results15 and our EGA results is not as good as the agreement between the HNC and original EGA results. We also find that the critical point of EGA is located at (0.026, 0.0789) instead of (0.018, 0.0837).4a This new critical point is close to the PMSA1 result (0.025, 0.0748), as one would expect on conceptual grounds. A similar re-evaluation of the EG result was also made independently by Guissani and Guillot.11 In Table 2, the critical points in the EGA and EGA+DI are compared with some of earlier theories and most recent computer simulations. The critical point of EGA+DI is located at (0.027, 0.0716), which is in somewhat better agreement with the most recent Monte Carlo estimates than the PMSAs.7 The phase diagrams can be calculated through the constructions involving the intersection of the Gibbs free energy and pressure at each temperature. The coexistence curves for EGA and EGA+DI are compared with the most recent simulation and PMSA3 in Figure 1. Both EGA and EGA+DI have coexistence curves of reasonable shape, and EGA+DI is in the best agreement with simulation. In Figure 2, the reduced chemical potential of EGA+DI is compared with that of MSA and of PMSA3 and computer simulation data. The reduced chemical potential of EGA+DI, βµ′, is calculated from the equation
Figure 1. Coexistence curves of the RPM in the EGA (-‚-); EGA+DI (s); PMSA3 (- -). See text for further details. The simulation data10d is indicated by the square symbol (0). Critical points in various approximations are indicated by dots (b).
(3.2)
Figure 2. Reduced chemical potential, βµ′, as a function of the logarithm of the reduced density (log F*) at various temperatures. The EGA+DI (s), PMSA3 (- -), and the MSA (- - - -) results are compared with simulation data [0, T* ) 0.147, ref 10a; O, T* ) 0.07 and 0.045, ref 10d]. From top to bottom, T* ) 0.147, 0.07, and 0.045, respectively.
β∆A βA βA° ) - ln(2F0σ3) (3.3) N N N As Figure 2 shows, the EGA+DI chemical potential is in
excellent agreement with the simulation data, in particular, in the low-temperature and low-density region. It is in slightly worse agreement in the high-temperature and high-density region. B. Significance of Results. A question naturally arises as to the reason for the significant disparity between the EGA+DI
βµ′ ) ln(F0σ3) +
β∆A βP -1 + N 2F0
with
TABLE 2: Comparison of Different Theoretical Results with Simulations
a
theory
ref
T* c
F*c
P* c
Rc (%)
MSA EGA EGA EGA+DI PMSA3 FL(DH+Bj+DI) FL(DH+Bj+DI+HC) MC, 1994 MC, 1994 MC, 1992
9 4a a a 7 12 12 10d 10c 10b
79 × 10-3 83.7 × 10-3 78.9 × 10-3 71.6 × 10-3 74.5 × 10-3 57 × 10-3 (55-58) × 10-3 53 × 10-3 g57 × 10-3 56 × 10-3
0.014 0.018 0.026 0.027 0.025 0.028 0.025-0.026 0.025
9.7 × 10-5 35.4 × 10-5 60.0 × 10-5 26.5 × 10-5 18.4 × 10-5
41.8 60.8 56.2
This work.
0.04
1418 J. Phys. Chem., Vol. 100, No. 4, 1996 coexistence curve and the coexistence curve as determined by the most recent simulations.10b-d (We consider only these, rather than the earlier, results of ref 10a because of the consensus that has emerged from the analysis given in refs 1c, 10b-d, 12, and 16 that the later simulation work has yielded more reliable coexistence-curve data.) What is left out of the theory that may be responsible for the lack of agreement with the RPM simulation results? One feature that is left out is ionic association of the charged spheres into clusters that are larger than dimerssespecially neutral tetramers, hexamers, etc. Results of two recent simulation studies16,17 show that the thermodynamic contributions of such neutral clusters are large at low densities and low temperatures. Neither study attempted to take RPM cluster statistics very close to the critical point, but an extrapolation of the results of ref 16 into the critical region suggests that the larger clusters make a significant contribution there too. Although the relative contribution from such larger clusters hinges on the precise definition of ionic cluster that is being used, it is clear that their presence appreciably lowers the unassociated ion population that one estimates from the EGA+DI, although the effect that the neglect of these larger clusters has on the T*c and F* c of the EGA+DI is harder to estimate. The neglected large-cluster contribution is an intracluster contribution, in contrast to the ion-ion and dipole-ion contributions that are accounted for in the EGA+DI. All intercluster contributions beyond these ion-ion and dipole-ion terms are also left out of the EGA+DI. Although the results of refs 16 and 17 suggest that in the critical region the neglected intracluster terms dominate the neglected intercluster terms, there is some compelling theoretical evidence in ref 11 that the neglected intercluster dipole-dipole interaction, at least, is also non-negligible, primarily as a result of its effect on the dielectric constant of the solvent. Guissani and Guillot add to a EG-type description both a dipole-ion (DI) term (via a generalized Born free-energy term computed in the MSA) and a dipole-dipole (DD) contribution (again computed in the MSA, as a dipolarsphere contribution). The resulting coexistence curve obtained in ref 11 can be described as our EGA+DI curve upon which is superimposed a rather broad bump at F* ) 0.064, which causes the coexistence curve to have its maximum of T* ) 0.0766 there, rather than at F* ≈ 0.027, at which T* ≈ 0.072. This new maximum becomes the new critical point, at which there is very little association (with Rc ) 0.025). The GuissaniGuillot result is an EG-type theory to which the intercluster effects of both a DI and a DD term have been added. However, it does not represent a quantitative improvement over the EGA+DI result. It seems plausible on the basis of the results of refs 11, 16, and 17 that the intracluster contribution from larger clusters is the primary contribution missing from the EGA+DI result and that it acts in the opposite direction from that of the missing intercluster contributions such as the DD contribution with respect to values of both Rc and T* c. Hence, although it seems most profitable to consider adding the missing intracluster contribution as the next step in improving the EGA+DI, in light of the results in ref 11 one must be prepared to find that one must add both the missing intercluster and at least the dipole-dipole intracluster terms in order to obtain substantial improvement. It turns out that the RPM, despite the analytic simplicity of its Hamiltonian, is by no means a simple fluid, as a result of the high degree of ionic association that takes place in the vicinity of the coexistence curve. Recently, Fisher and Levin12 have had remarkable success in reproducing the simulated coexistence curve over a consider-
Yeh et al. able range of density about F* c by adding a dipole-ion term to the linearized Debye-Hu¨ckel theory augmented by Bjerrum’s theory of ion pairing. We shall refer to the result as the DH+Bj+DI theory, for which T*c ) 0.057 and F* c ) 0.028. The question arises as to why this theory predicts T*c so much more accurately than EGA+DI, which yields a T* c about 26% higher than DH+Bj+DI. The latter yields a T*c that may itself be a bit on the high side but is consistent with the most recent simulation results within the latter’s range of considerable uncertainty (see Table 2). Conceptually, the linearized DebyeHu¨ckel theory upon which the Fisher-Levin result builds is far less satisfactory than the MSA. In particular, the linearized Debye-Hu¨ckel theory is a theory of charge-charge correlation but not density-density correlation. Hence, in its thermodynamic prediction it includes no hard-core term that persists when the ions are discharged. Moreover, its electrostatic treatment of the core term in describing the charge-charge correlation gives rise to a violation of the Stillinger-Lovett second-moment condition. The Fisher-Levin treatment of the DI term also appears to be somewhat more crude than the DI treatment used in the EGA+DI theory.18 (On the other hand, the quantitative difference between the Bjerrum treatment used by Fisher and Levin and the EG treatment of association is quantitatively negligible in the critical region, as Fisher and Levin have stressed.) As we have already discussed in ref 7, the addition of an accurate hard-core contribution to the DH+Bj+DI free energy changes the location of the critical point only slightlysit lowers T*c a bitsbut it narrows the predicted coexistence curve considerably, causing it to fall out of good agreement with the simulated curve on the high-density side. Fisher and Levin already had considered the addition of a hard-core term themselves, but it was chosen to give rise to appropriate thermodynamic behavior at a close-packing density many times higher than T* c, and it makes too small a contribution in the critical region. The DH+Bj+DI+HC listed in our Table 2 refers to the two slightly different versions of this theory mentioned in ref 12. In ref 7 we concluded that the accuracy of the DH+Bj+DI may be largely fortuitous, resulting in a substantial cancellation of the contributions from terms lacking in the Fisher-Levin treatment. We believe this conclusion is buttressed by new work reported in ref 11, which adds to DH+Bj+DI a DD term (as well as an accurate HC term) appropriately tailored to the Fisher-Levin approach. Guissani and Guillot find that the added terms give rise to huge effects on the coexistence curve that throw it into very much worse agreement with the simulation results. It becomes much narrower, with a T* c increasing to 0.0784 and F* c decreasing to 0.0179. The Rc is 0.08. On the basis of such observations it seems to us that the DH+Bj+DI does not enjoy its outstanding quantitative accuracy as a result of having faithfully incorporated the most important contributions to the thermodynamic behavior but rather because all the contributions it leaves out tend to cancel one another. We hasten to add that this does not diminish its value as a means of accurately describing thermodynamic behavior of the RPM over a rather large range of densities in the vicinity of the coexistence curve. For this purpose neither the EGA+DI nor any other approximation procedure of which we are aware can rival it. Acknowledgment. Helpful discussions and correspondence with Dr. Guillot are gratefully acknowledged. We are also indebted to Prof. J. C. Rasaiah for calling to our attention the most accurate HNC evaluations available for our Table 1. S.Y. and Y.Z. acknowledge the support of the National Science Foundation, and G.S. acknowledges the support of the Division
Phase Separation of Ionic Fluids of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, U. S. Department of Energy. References and Notes (1) Waisman, E.; Lebowitz, J. L. J. Chem. Phys. 1972, 56, 3086; 1972, 56, 3093. (2) Recent reviews: (a) Pitzer, K. S. Acc. Chem. Res. 1990, 23, 333. (b) Levelt-Sengers, J. M. J.; Given, J. A. Mol. Phys. 1993, 80, 899. (c) Fisher, M. E. J. Stat. Phys. 1995, 75, 1. (d) Stell, G. J. Stat. Phys. 1995, 78, 197. (e) For a recent work, also see: Bhuiyan, J. L.; Outhwaite, C. W.; Molero, M.; Gonza´lez-Tovar, E. J. Chem. Phys. 1994, 100, 8301. (3) (a) Stell, G.; Wu, K. C.; Larsen, B. Phys. ReV. Lett. 1976, 37, 1369. (b) Stell, G.; Larsen, B. J. Chem. Phys. 1979, 70, 361. (c) Larsen, B.; Stell, G.; Wu, K. C. J. Chem. Phys. 1977, 67, 530. (d) Hafskjold, B.; Stell, G. In Studies in Statistical Mechanics; Montroll, E. W.; Lebowitz, J. L., Eds.; North-Holland: Amsterdam, 1982; Vol. VIII. (4) (a) Ebeling, W.; Grigo, M. Ann. Phys. 1980, 37, 21. (b) Ebeling, W.; Grigo, M. J. Solution Chem. 1982, 11, 151. (5) For example: Friedman, H. L. Ionic Solution Theory; Interscience: New York, 1972. (6) Bjerrum, N. K. Dan. Vidensk. Selsk. Mat.-Fys. Medd. 1926, 7, 1.
J. Phys. Chem., Vol. 100, No. 4, 1996 1419 (7) Zhou, Y.; Yeh, S.; Stell, G. J. Chem. Phys. 1995, 102, 5785. (8) Zhou, Y.; Stell, G. J. Chem. Phys. 1992, 96, 1504; 1992, 96, 1507. Stell, G.; Zhou, Y. Fluid Phase Equilib. 1992, 79, 1. (9) Medina-Noyola, M.; McQuarrie, D. A. J. Stat. Phys. 1978, 18, 445. (10) (a) Valleau, J. P. J. Chem. Phys. 1991, 95, 584. Graham, I. S.; Valleau, J. P. J. Phys. Chem. 1990, 94, 7894. (b) Panagiotopoulos, A. Z. Fluid Phase Equilib. 1992, 76, 97. (c) Caillol, L. M. J. Chem. Phys. 1994, 100, 2161. (d) Orkoulas, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 1994, 101, 1452. (11) Guissani, Y.; Guillot, B. Mol. Phys., in press. (12) Fisher, M. E.; Levin, Y. Phys. ReV. Lett. 1993, 71, 3826. (13) Ebeling, W. Z. Phys. Chem. 1968, 238, 400. (14) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969, 51, 635. (15) Rasaiah, J. C. In The Liquid State and Its Electrical Properties; Kunhardt, E. E., Chrostophorou, L. G., Luessen, L. H., Eds.; NATO ASI Series Vol. 193; Plenum Press: New York, 1989. (16) Caillol, J. M.; Weis, J. J. J. Chem. Phys. 1995, 102, 7610. (17) Bresma, F.; Lomba, E.; Weis, J. J.; Abascal, J. L. F. Phys. ReV. E 1995, 51, 289. (18) The DI part of ESA + DI is already discussed in ref 7, where a detailed comparison with the treatment of ref 7 is given.
JP952412D
