Thermodynamically Consistent Adaptation of ... - ACS Publications

Scaled particle theory (SPT) is adapted to an arbitrary hard-sphere equation of ... particular, the Carnahan-Starling EOS, the resulting work of cavit...
1 downloads 0 Views 139KB Size
Download PDF
Ind. Eng. Chem. Res. 2006, 45, 5489-5500

5489

Thermodynamically Consistent Adaptation of Scaled Particle Theory to an Arbitrary Hard-Sphere Equation of State Daniel W. Siderius and David S. Corti* School of Chemical Engineering, Purdue UniVersity, West Lafayette, Indiana 47907-2100

Scaled particle theory (SPT) is adapted to an arbitrary hard-sphere equation of state (EOS) while utilizing an SPT interpolation function that satisfies all the exact conditions of SPT and, as such, is thermodynamically consistent. The new modification of SPT is tailored to predict with high accuracy the reversible work of forming a cavity in a hard-sphere fluid, but it also estimates the surface tension and excess surface adsorption. The application of additional SPT conditions also provides estimates of the first two derivatives of the hardsphere radial distribution function (RDF) at contact. When various accurate equations of state are used, in particular, the Carnahan-Starling EOS, the resulting work of cavity formation is nearly identical to accurate Monte Carlo simulation data and the predictions of the recent SPT6 [Heying and Corti, J. Phys. Chem. B 2004, 108, 19756]. The estimates of the surface tension and surface adsorption at moderate to high packing fractions deviate somewhat from both previous theoretical techniques and the simulation. Meanwhile, EOSbased predictions of the initial slope of the RDF are very close to simulation estimates. Predictions of the initial curvature, like those of SPT6, deviate from simulation, but not drastically so. The overall accuracy of our version of SPT demonstrates that a thermodynamically consistent modification of SPT provides correct predictions of the work of cavity formation and is, therefore, suitable for use as the ideally repulsive reference state in various solvation theories of liquids. 1. Introduction In many solvation theories of liquids, the free energy difference upon addition of a solute molecule is typically obtained via a two-step process.1-4 The first step is the creation of a molecular-sized cavity (a region devoid of solvent particle centers) in the liquid, and the second step is the activation of the intermolecular potential between the solute and solvent molecules. It is well-known that the energy cost of forming the molecular-sized cavity within the liquid is a nontrivial quantity. For example, in some cases, cavity formation in model liquids is developed utilizing a framework modeled on the perturbation theory of Weeks, Chandler, and Andersen,2,4,5 in which the work of cavity formation is divided into a repulsive reference contribution and an attractive perturbation. The repulsive term is often replaced by the hard-sphere potential as an ideal repulsive reference.6,7 Rigorous theoretical or semiempirical methods that provide values for the work of cavity formation are of constant interest for use in the prediction of solvation free energies. In particular, generating accurate expressions for the work of cavity formation in hard-sphere fluids is an important area of research. Accurate predictions of the work of forming cavities in hardparticle fluids is also relevant to the study of depletion interactions in colloidal dispersions.8-10 In colloidal fluids, interactions between colloidal particles or between a single colloid and a surface can sometimes be approximated by the hard-sphere potential.11 As we discuss later, the equivalency of hard-sphere colloids and cavities allows for the substitution of a colloid particle with an appropriately sized cavity. Thus, the work of forming the cavity is equal to the excess chemical potential of the equivalent hard-sphere colloid. Because the depletion potential of a colloid particle at a certain position is * To whom correspondence should be addressed. Tel.: (765) 4966064. Fax: (765) 494-0805. E-mail: [email protected].

equal to the difference in the free energy between that position and a reference position far away from the other particle or surface, any calculation of the depletion potential requires an accurate estimation of the work of cavity formation. Accuracy is of particular importance when cavities are allowed to intersect a hard wall, where knowledge of the hard-sphere surface tension and line tension are needed.9,10 Many theories for predicting the work of cavity formation stem from the scaled particle theory (SPT) of Reiss, Frisch, and Lebowitz (RFL).12 SPT predicts the work of cavity formation in the hard-sphere fluid exactly for cavities smaller than a hardsphere solvent particle and develops an interpolation function that spans from small cavity radii to the limit of macroscopic size. The interpolation is based on a number of exact conditions imposed by geometry and thermodynamics and yields a hardsphere equation of state (EOS). SPT is quite accurate in predicting the work of cavity formation, even up to large cavity sizes, when many of the exact conditions are used. The most recent version of SPT uses six conditions,13 providing estimates of the work of cavity formation that are much more accurate than the original version of SPT, which used only three conditions.12 This six-condition version of SPT also produces an EOS that is comparable in accuracy to the semiempirical Carnahan-Starling EOS.14 The usefulness of the more advanced versions of SPT, however, is limited by the reasonably intensive computations required to obtain the SPT interpolation. For example, both the five-15 and six-condition13 versions of SPT require the solution of an integral equation. Because of the numerical challenges inherent in highly accurate versions of SPT, some effort has been directed toward developing approximations of the work of cavity formation that are similar in structure to SPT (i.e., invoke some of the SPT conditions) but are based on existing equations of state (rather than on a self-consistent set of equations that are then used to generate an EOS16,17). Basing an approximation for the work

10.1021/ie051038t CCC: $33.50 © 2006 American Chemical Society Published on Web 01/06/2006

5490

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

of cavity formation on an existing EOS is attractive, not only because many empirical equations of state are highly accurate, leading to good estimates of the work of cavity formation, but also because the resulting cavity function is very simple (no additional numerical solution is required). As discussed in the next section, previous approaches that have incorporated a given EOS did not, however, use all the possible information describing cavities in hard-sphere fluids or were not constructed via a thermodynamically consistent set of equations. Admittedly, the numerical inaccuracy brought about by using a set of thermodynamically inconsistent set of equations is small. Yet one could argue that the cavity function, required to be as accurate as possible, should satisfy all the fundamental conditions of SPT as well as of thermodynamics (e.g., the GibbsDuhem equation). Given the inconsistencies in some previous approaches, we therefore present a method of adapting SPT to an arbitrary hardsphere EOS that incorporates all known exact SPT conditions and results in a thermodynamically consistent expression for the work of cavity formation. Our modification of SPT to an existing EOS is specifically targeted to emulate the accuracy of the sixcondition form of SPT, but in a simpler implementation. All of the geometric and thermodynamic conditions presently available from SPT are used. Like all versions of SPT, our thermodynamically consistent form of SPT also has uses beyond the prediction of the work of cavity formation. SPT also provides information about the surface thermodynamics of cavities and the local structure of the hard-sphere fluid.12,15 By being incorporated into the SPT framework, the basis EOS is now able to describe properties of the hard-sphere fluid beyond just the pressure or chemical potential. SPT indeed proves to be a powerful tool for studying the thermodynamics of the hard-sphere fluid. This paper is organized as follows. Section 2 offers a review of the basic concepts of SPT and provides the foundation on which to develop our modification of SPT. This section also describes the theoretical underpinnings of previous modifications of SPT and discusses the assumptions and interpolation choices each one made to obtain the work of cavity formation. In Section 3, our adaptation of SPT to a given EOS is described. Section 4 compares the thermodynamic predictions of our modification of SPT to molecular simulation results and previous theoretical approaches. Conclusions are provided in Section 5. 2. Review of Scaled Particle Theory To provide a proper foundation for adapting SPT to an arbitrary hard-sphere EOS, we introduce the terminology of SPT and list some of its exact expressions. The number of exact expressions that can be derived within SPT is, in principle, infinite in number. Only a few, however, yield conditions that are immediately useful, not requiring further approximations when applied. We begin by considering a system of hard spheres of diameter σ at a temperature T and a number density F ) N/V, where N is the number of hard-sphere particles and V is the system volume. The fluid is assumed to be in the thermodynamic limit so that it has no knowledge of its container or, equivalently, that spatial correlations between particles are zero for large particle separations. A cavity of radius λ, or a spherical region devoid of hard-sphere centers, is now introduced into the fluid. As demonstrated in the original SPT paper,12 this cavity is equivalent to a hard-sphere solute of diameter σs, where λ and σs are related according to

λ)

σs + σ 2

(1)

Figure 1. Illustration of the relationship between the hard-sphere solvent and solute diameters σ and σs, respectively, and the cavity radius λ. The centers of the solvent particles are excluded from the spherical regions denoted by dashed lines, which correspond to the equivalent cavities of the hard particles. The local density of solvent particles at the cavity surface is defined as FG(λ,F), where F is the bulk number density of the solvent.

Equation 1 arises because the closest approach between two rigid spheres is the sum of their radii, as illustrated in Figure 1. From eq 1, a cavity of radius λ < σ/2 corresponds to the seemingly nonphysical case of σs < 0. These “particles” of negative diameter do not interact with the solvent hard sphere via a hard potential and, consequently, can penetrate the hard core of a solvent sphere. The local density of solvent hard spheres at the cavity surface is defined as FG(λ,F).12 G is the central function of SPT and is related to a number of hard-sphere thermodynamic properties. As we discuss later, G(λ,F) is known explicitly for λ e σ/2 but is nonanalytic for all other values. One can argue, however, that despite its nonanalytic nature, G(λ,F) can be wellrepresented by a smoothly varying function for λ > σ/2.12 In terms of G(λ,F), the reversible work, W, of inserting or growing a cavity of radius λ is given by12

W(λ) ) 4πFkT

∫0λG(r,F)r2 dr

(2)

where k is the Boltzmann constant. W is simply obtained by integrating the product of the kinetic pressure on the cavity surface, FG(λ,F)kT, and the differential volume 4πλ2 dλ. G(λ,F) can also be written in terms of the probability of randomly observing, or inserting, a cavity of radius of at least λ, denoted by P0(λ). Using fluctuation theory,18 one can write

P0(λ) ) exp[-βW(λ)]

(3)

where β ) 1/kT. By introducing eq 2 into the above expression and taking a derivative with respect to λ, one arrives at12

G(λ,F) )

-1 ∂ ln P0(λ) ∂λ 4πFλ2

(4)

With eq 4, obtaining a form of G depends on identifying P0(λ), followed by straightforward manipulation. This is a more difficult task than it might seem, however, because P0(λ) is known exactly for only a small range of λ. For example, cavities of λ e σ/2 may contain at most one hard-sphere center because the hard cores of two spheres may not overlap. In this case, the probability that the cavity contains a single particle center is (4/3)πFλ3. Therefore, P0(λ) is equal to 1 - (4/3)πFλ3, providing the following exact expression.12

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5491

G(λ,F) )

1 4 1 - πFλ3 3

(λ e σ2)

The description of G for a cavity whose radius exceeds σ/2 requires the introduction of arguments relating P0 to multiparticle correlation functions, because at least two particle centers may occupy the cavity.12 Using these expressions for P0, RFL rigorously showed that P0 and its first two derivatives are continuous at λ ) 1/2.12 Since the ith derivative of G depends on the (i + 1)th derivative of P0 (see eq 4), both G and ∂G/∂λ are continuous at λ ) σ/2, with their values being equal to12

1 1 G ,η ) 2 1-η ∂G 6η ) ∂λ λ)1/2 (1 - η)2

( )

|

(6) (7)

where the dimensionless packing fraction η ) πFσ3/6 has been introduced for notational simplicity and λ has been scaled by σ. (These dimensionless variables will be used throughout the remaining analyses.) Further manipulation of P0 and its derivatives indicates that the third derivative of P0 is discontinuous at λ ) 1/2.12 Entering the discontinuity of ∂3P0/∂λ3 into the second derivative of eq 4 yields the third exact condition on G,12

∂2G ∂λ2

|

λ)1/2+

)

24η + 48η2 48η g(1+) 3 1 η (1 - η)

(8)

where the derivative is noted to be taken at a value infinitesimally larger than 1/2. The first term on the right side is equal to ∂2G/∂λ2, as evaluated from below 1/2, and the second term is the discontinuity across 1/2. We note that the discontinuity is proportional to g(1+), the value of the hard-sphere radial distribution function at contact. Two other exact conditions on G may be obtained through geometric and thermodynamic arguments. Both result from the equivalence of cavities and hard-sphere solute particles as noted in eq 1. For λ ) 1, the cavity is equivalent to a solvent particle and so G(1,η) is equal to g(1+). (This equality allows for the replacement of g(1+) with G(1,η) in eq 8 and elsewhere.) By combining this condition with the virial equation of state (EOS) for hard spheres, the fourth exact condition of SPT is12

p ) G(∞,η) ) 1 + 4ηG(1,η) FkT

µex )

(5)

(9)

where the equality G(∞,η) ) p/FkT has also been introduced. As λ f ∞, the cavity becomes indistinguishable from a hard wall so that the pressure exerted by the collisions of hard spheres with the cavity surface becomes identical to that exerted by the hard-sphere fluid on a flat surface. Hence, FG(λ,η)kT, the pressure on the cavity surface, will equal the fluid pressure. Another condition results from the equivalence of solvent particles and cavities of radius λ ) 1. As W(1) is the work required to “grow” a cavity of radius λ ) 1, it is also the free energy required to scale up a hard particle from zero size to its usual diameter of σ (hence the name “scaled” particle theory). Thus, W(1) contains all nonideal contributions to the free energy of adding another hard particle to the fluid and is equal to µex, the excess chemical potential of the hard-sphere solvent. Thus, by introducing a form of the Gibbs-Duhem relation,

( ) ex

∫0FF′1 ∂p∂F′

dF′

(10)

and combining the above with both eq 2 and eq 9, one can write the fifth condition on G as12

-ln(1 - η) + 24η

∫1/21G(r,η)r2 dr ) G(∞, x) - 1 dx + G(∞,η) - 1 ∫0η x

(11)

Equation 11 is the Gibbs-Duhem integral equation condition on G and ensures thermodynamic consistency between G(∞,η) (essentially the EOS) and the integral of G up to λ ) 1. Despite being a nonanalytic function for λ > 1/2, the use of a series representation for G that interpolates between the available conditions has been quite successful.12,13,15,19 This series is based on the thermodynamic expression for the formation of a macroscopic cavity, where the work of forming a cavity is given by

2δ 4 + ‚‚‚ W(λ) ) πλ3p + 4πλ2γ∞ 1 3 λ

(

)

(12)

In the above equation, the first term is the pressure-volume work contribution and the second term is the surface tension contribution. γ∞ is the fluid surface tension (more properly the boundary tension) of a cavity of zero curvature, and δ is the Tolman length, providing a first order correction for surface curvature.20 Equation 12 also notes that higher-order corrections to the surface tension may be included. When eq 12 is compared to eq 2, one is motivated to express G as a Laurent series, with the generic form being n

G(λ,η) )

∑ i)0

Ri(η) λi

(13)

The fitting coefficients Ri are functions of the fluid density and can be obtained by demanding that eq 13 satisfy the exact conditions described previously. It is straightforward to show that R0(η) ) p/FkT from the relationship between G(∞,η) and the fluid pressure. The number and identities of the other coefficients depend on the number of exact conditions used to solve for Ri. For n g 3, one is required to choose R3 ) 0 to prevent the appearance of nonphysical logarithmic terms which would arise from the integration of eq 13 in eq 2.19,21 It should be noted that this constraint on R3 is not a true condition on G, but it is required because of the selection of eq 13 as an interpolating function. In the original SPT paper, the series coefficients of G were obtained for n ) 2 along with three exact conditions, those being the continuity of G and its first derivative at λ ) 1/2 and the virial condition shown in eq 9.12 For this set of matching conditions, the Laurent series coefficients are

R0(η) ) R1(η) )

1 + η + η2 (1 - η)3 -3η(1 + η)

R2(η) )

2(1 - η)3 3η2 4(1 - η)3

(14)

(15)

(16)

The resulting EOS (given by the value of R0(η) in eq 14) is

5492

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

curiously identical to the Percus-Yevick EOS for hard spheres, which was derived later from a fundamentally different approach.22 Also, by equating the integral of eq 13 (in eq 2) to eq 12, the surface tension and Tolman length are given by12

γ∞σ2 -9η2(1 + η) ) kT 2π(1 - η)3

(17)

δ η ) σ 2(1 + η)

(18)

Tully-Smith and Reiss later revisited the three-condition form of SPT (now known as SPT3) and showed that it actually included fiVe conditions on G, in that it assumed R3 ) 0 and satisfied the Gibbs-Duhem integral equation (eq 11) without including the additional fitting coefficients.19 SPT3 did not adequately represent simulation data at high packing fractions, motivating Mandell and Reiss to revisit SPT and employ all five exact conditions.15 The resulting EOS, known as SPT5, was an improvement over SPT3, though it still deviated from simulation data at high packing fractions.15 Two additional relations for G, originally exact though requiring approximate, yet very accurate, closure conditions, were recently developed by Heying and Corti and incorporated in the SPT6 and SPT7 versions of scaled particle theory.13 These conditions will be described in Section 4.3. The SPT6 version proved to be the most accurate SPT to date, providing an equation of state comparable with the highly accurate semiempirical hard-sphere equations of state of Carnahan and Starling14 and Kolafa et al.23 Prior to SPT6 and SPT7, the predictions of W(λ) from SPT differed from simulation results at moderate to high densities by statistically significant amounts. Because of the importance of hard-sphere cavity functions as reference states for various liquid state and solvation theories, interest in developing accurate cavity functions has not diminished over the years. Since SPT remains a convenient framework for describing cavities in hardsphere fluids, as well as providing a large number of exact conditions, its ideas have been incorporated into several approaches to generate more accurate cavity functions. In particular, W(λ) or G(λ,η) has been modified to be consistent with the pressure of a more accurate EOS, such as the Carnahan-Starling EOS (henceforth, CS-EOS).16,17 For example, Matyushov and Ladanyi developed a function which fit, using the usual SPT notation, W(λ) by a cubic series in λ.16 The series coefficients in their method were chosen such that G(1,η) and G(∞,η) were consistent with the CS-EOS and the exact value of G(λ,η) at λ ) 1/2. The result is similar to SPT in that the cavity work function contained no terms higher in order than λ3. All conditions on G were not, however, included, and so the method is not thermodynamically consistent. In particular, their chosen interpolation does not match the chemical potential of the hard-sphere solvent according to the CS-EOS and, consequently, overpredicts the work of cavity formation for mid-to-high packing fractions. The deviation, though, is quite small. In a more recent paper, Sun adapted the CS-EOS to an SPT-like interpolation series in a more traditional manner.17 This interpolation series is identical to eq 13 for n ) 3. The pressure is fixed according to the given EOS, and the values of G(1/2,η) and G′(1/2,η) are set according to the usual SPT eqs 6 and 7. A fourth condition used is the consistency of G(1,η) and the EOS (see eq 9). Like the method of Matyushov and Ladanyi, this interpolation does not correctly predict the chemical potential and, as such, is not thermodynamically consistent. In addition, this approach includes an O(λ-3) term in the interpolation series (i.e., R3 * 0). As previously mentioned, this term is nonphysical and should not

be included in the Laurent series representation of G.19,21 The underprediction in the chemical potential arising from this term is again small, especially compared to the predictions of SPT3, but this underprediction persists for larger cavities. And so, a need still arises for a method that utilizes both a given EOS and all the exact conditions on G and that enforces thermodynamic consistency between the cavity function and the chosen EOS. In the following section, we present such a method, i.e., a full, thermodynamically consistent adaptation of SPT to any existing hard-sphere equation of state. 3. Adaptation of SPT to an Equation of State Our adaptation of SPT to an arbitrary hard-sphere EOS begins, like previous modifications of SPT, with the decision to impose an EOS on the SPT function G rather than solve for the EOS. In addition, we enforce consistency between the series representation of G and all possible thermodynamic conditions implicit in the chosen EOS. These consistency conditions are then combined with the statistical geometric conditions of SPT to provide a set of equations that, when solved, yield expressions for the various Ri. From the relationship between G and the hard-sphere fluid pressure, the first condition of our SPT modification is, as for others,

G(∞,η) )

p )Z FkT

(19)

where Z is the compressibility factor and is now considered to be known, that is, set equal to some chosen EOS. Although eq 19 is an exact condition, specifying the form of Z a priori makes this relation “pseudoexact”. The second condition follows from eq 9 in which

G(1,η) )

G(∞,η) - 1 Z - 1 ) 4η 4η

(20)

Like the use of the virial equation in the solution of a selfconsistent SPT, eq 20 ensures thermodynamic consistency between the values of G(1,η) and G(∞,η). The last thermodynamic consistency condition provided directly by the EOS is the relationship between G and the solvent’s excess chemical potential. Since the chemical potential is known explicitly from the EOS, eq 11 is transformed to

βµex(η) ) -ln(1 - η) + 24η

∫1/21G(r,η)r2 dr

(21)

Eq 21 serves to constrain the form of G(λ,η) between λ ) 1/2 and λ ) 1. µex is sometimes available explicitly as a function of η (as in the case of the CS-EOS), but it is otherwise calculated by numerically integrating the EOS as shown by the right-hand side of either eq 10 or eq 11. Equations 19-21 are the three necessary thermodynamic consistency conditions of our current SPT modification. The other conditions on G that should be used are the three exact conditions at λ ) 1/2. As a review, these conditions are

1 1 G ,η ) 2 1-η ∂G 6η ) | ∂λ λ)1/2 (1 - η)2

( )

|

∂2G ∂λ2

) λ)1/2+

24η + 48η2 48η G(1,η) 1-η (1 - η)3

(22) (23) (24)

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5493

where G(1,η) has been substituted for g(1+) in eq 24. No other exact conditions can be written for G without supplying additional conditions on the radial distribution function g(r).13 Consequently, we have a total of six conditions on G (eqs 1924) that can be used to evaluate the coefficients in a chosen interpolation series, like that shown in eq 13 for n ) 6. Because the Laurent series is chosen as the representation of G beyond λ ) 1/2, we must again have R3(η) ) 0. Hence, our thermodynamically consistent SPT is obtained by simultaneously solving eqs 19-24 with the Laurent series of eq 13 reduced to

G(λ,η) ) R0(η) +

R1(η) R2(η) R4(η) R5(η) R6(η) + 2 + 4 + 5 + 6 (25) λ λ λ λ λ

Solving for all Ri(η) requires fairly straightforward, though tedious, calculations (again we note that R0(η) ) Z). The resulting expressions for Ri(η) are quite large and are provided in the Appendix. In general, the expressions for each Ri(η) are given in terms of Z and µex, and they will not yield expressions explicit in η if µex is obtained by numerical integration. For equations of state with an explicit relation for µex, such as the CS-EOS, the Ri(η) are fully known in terms of η (see the Appendix). In the following section, we discuss the use of the above interpolation scheme in predicting various thermodynamic quantities via the CS-EOS,14 the Kolafa EOS (K-EOS),23 and the 4,3-Pade´ approximation of Sanchez (S-EOS; based on the known virial coefficients up to the eighth),24 three equations of state considered to be among the most accurate expressions available to date. Again, our interpolation scheme is not limited to these equations of state but can be applied to any existing hard-sphere EOS.

4. Predictions of Thermophysical Properties from SPT In addition to the work of cavity formation that it is initially tailored to calculate, SPT (more general than just our method) can provide a number of other useful thermodynamic properties. For example, by revisiting the surface thermodynamics of cavities in hard-sphere fluids,25,26 one can obtain predictions of the surface tension and the excess surface adsorption of the hardsphere fluid. SPT can also be used to generate an estimate of the initial slope and curvature of the radial distribution function, g(r), via the newer exact SPT conditions derived by Heying and Corti13 as well as g(r) over a limited range. Such SPT predictions are presented in successive subsections and compared to both simulation results and calculations from previous versions of SPT and other sources. 4.1. G(λ,η) and the Work of Cavity Formation. Since G is most often used to calculate the free energy cost of cavity formation, an initial test of any new expression for G is its ability to accurately predict W(λ). Since our method enforces consistency between W(1) and the solvent excess chemical potential, as specified by some chosen EOS, the accuracy of both G(λ,η), at least up to λ ) 1, and its corresponding prediction of W(1), relative to simulation, depends solely on the EOS. (An examination of various simulation results in the literature points to the 4,3-Pade´24 as the best EOS for predicting µex,13 though not necessarily the hard-sphere pressure.) For larger cavities where the particular value of G(∞,η) becomes important, we expect to see deviations between the various SPT-based predictions of G and W and the corresponding values generated by simulation.

Figure 2. SPT function G(λ) for η ) 0.35. The solid line denotes the prediction of our method based on the CS-EOS. Predictions from SPT312 are shown by the dotted line and those from SPT613 are denoted by a dashed-dot line that is nearly indistinguishable from the CS-EOS-based results. The long-dashed line is the prediction of M&L.16 The result of Sun17 is given by the short dashed line, which is masked by the CS-EOS prediction. For this packing fraction, G(∞,η) as calculated from simulation is ∼5.206, which is nearly identical to G(∞,η) predicted by the CS-EOS and SPT6, which are 5.206 and 5.193, respectively. SPT3 overpredicts G(∞,η) with a value of 5.362.

For example, we plot G(λ,η ) 0.35) in Figure 2, obtained from our thermodynamically consistent method using the CSEOS. Also included are the self-consistent versions of SPT3 and SPT6,13 the asymptotically correct method of Matyushov and Ladanyi (M&L),16 and the CS-EOS fit of Sun.17 (We have not included predictions from our adaptation of SPT to the K-EOS or S-EOS, because their predictions are indistinguishable from the CS-EOS on the scale of Figure 2.) At λ ) 3, for example, the values of G from the CS-EOS, SPT6, SPT3, M&L, and Sun are 4.364, 4.357, 4.539, 4.390, and 4.363, respectively. Our method, using the CS-EOS, yields a G that is similar to what is predicted by SPT6 and is much more accurate than that by SPT3 (as should be expected since the SPT3 pressure is not as accurate as the CS-EOS). The CS-EOS-based predictions differ slightly from those of both M&L and Sun, though all methods, other than SPT3, fall close to the G generated by SPT6. The differences between the various modifications to SPT can be attributed to the interpolation scheme used in each method, since our results and the methods of M&L and Sun all use the CS-EOS. For example, the asymptotically correct expansion of M&L yields slightly higher values of G, because its interpolation function overpredicts µex at this particular packing fraction. Sun’s SPT modification yields slightly lower G values, a property that can be attributed to the chosen interpolation scheme that overpredicts µex for all but the smallest packing fractions. (We also observed differences in G when the K-EOS and S-EOS were used as the basis EOS, with the differences being negative and positive, respectively, as compared to the SPT6 results.) Overall, the various methods produce SPT functions clustered close to the accurate predictions of SPT6. The small deviations in G between the various versions of SPT are further reflected in the predictions of W(λ) (since W is an integral of G, see eq 2). Figure 3 displays the predictions of W(2) and W(3) as a function of η for all the methods shown in Figure 2. (Again, we have not included predictions from our adaptation of SPT to the K-EOS or S-EOS, because their predictions are nearly indistinguishable from the CS-EOS-based

5494

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

constant needed to match W to its exact value at λ ) 1/2. By rearranging the above equation, the surface tension of a planar interface is identified as

(

)

4 W(λ) - πλ3p 3 γ∞ ) lim λf∞ 4πλ2

(27)

Although γ∞ should be correctly referred to as the planar, or hard-wall, surface (boundary) tension, it is usually just called the surface tension. When the interpolation form for G (eq 13) is entered into eq 27, the surface tension is found to be simply proportional to R1(η) and is given by

γ∞σ2 3η ) R1(η) kT π Figure 3. Reversible work of cavity formation, W(λ)/kT, as a function of packing fraction, η, for cavity radii λ ) 2 and 3. The identities of the lines are the same as those in Figure 2. The solid circles denote results from Monte Carlo simulations.

results.) Also shown in the figure are the results of Monte Carlo (MC) simulations, which were generated with a Metropolistype MC algorithm using biased umbrella sampling to reach large values of W(λ). This simulation technique has been described in previous papers.9,13,27 Our method, using the CSEOS, yields works of cavity formation that are nearly indistinguishable from SPT6 and much more accurate than SPT3 (again as should be expected). The CS-EOS-based predictions differ slightly from those of both M&L and Sun, but all methods, other than SPT3, fall very close to simulation values. At λ ) 3 and η ) 0.35, the values of W/kT from the CS-EOS-based SPT modification, SPT6, SPT3, M&L, and Sun are 303.8, 303.2, 316.1, 305.5, and 303.6, respectively. With MC simulations yielding a value of 308.1, all predicted values differ from simulation by no more than 1.6% (certainly within simulation error margins) for any method other than SPT3. These results indicate that our adaptation of SPT to the CS-EOS accurately predicts the work of cavity formation up to very high packing fractions. Similar accuracy is found when SPT is adapted to the K-EOS23 or the S-EOS,24 where W(3σ) for η ) 0.35 differ from the CS-EOS predictions by -0.07kT and 0.44kT, respectively. The SPT predictions of W are highly accurate whenever the EOS itself is accurate. Ultimately, a more stringent test of each G would be to compare each W with simulation results for very large cavities. Because either the equations of state show small differences or the chosen interpolation schemes are not the same, plots of W should show greater deviations at larger λ. These predictions would be difficult to validate, however, since determining W for large cavities is computationally prohibitive. 4.2. Surface Properties from G(λ,η). Since the chosen interpolation of G is based on an extension of traditional surface thermodynamics8,12,28 to cavities, SPT also provides expressions for the surface (boundary) tension, Tolman length, and excess surface adsorption. To obtain these properties, we begin with the reversible work of cavity formation given in eq 12, rewritten as

W(λ) ) 2δ δ2 δ3 δ4 δ5 4 3 πλ p + 4πλ2γ∞ 1 + 2 + 3 + 4 + 5 + K (26) 3 λ λ λ λ λ

(

)

In eq 26, δ is again the Tolman length and each δi represents another higher-order correction to the surface tension. K is a

(28)

The general form of R1(η), fit to any chosen EOS and corresponding µex, is given in eq 44 of the Appendix. Since γ∞ e 0 for the hard-sphere fluid, we note that R1(η) e 0 as well. To provide an initial comparison of the surface tension as calculated by various SPT-based techniques, we include numerical values for γ∞ in Table 1 for two self-consistent SPT methods, Sun’s SPT modification, Henderson and Plischke’s semiempirical equation,29 and three thermodynamically consistent adaptations of SPT. In general, we see excellent agreement between SPT6 and our three EOS-based adaptations of SPT, confirming that similar equations of state produce similar predictions for γ∞. When compared to the SPT3 prediction of γ∞ and the semiempirical equation of Henderson and Plischke,29 however, our SPT modification predicts a larger magnitude (more negative) for the surface tension. Sun’s SPT modification predicts a surface tension between SPT3 and the more accurate equations of state. Because of the lack of thermodynamic consistency and the inclusion of nonphysical terms in Sun’s interpolation, however, this result does not necessarily indicate a better prediction of the surface tension (given that the use of the R3 term should most likely make these surface thermodynamic predictions unreliable). The difference between the EOSbased predictions and SPT3 is not as easy to argue away because previous simulation predictions of γ∞30-32 largely agree with the SPT3 surface tension when listed or plotted as a function of the packing fraction. The difference between the surface tension obtained from SPT3 and more accurate equations of state can be interpreted in a slightly different manner than is usually done. Since the pressure is known to vary significantly between the various forms of SPT and other equations of state, we note that R1(η) and any surface tension derived from some version of SPT depend on the chosen interpolation scheme (and on which conditions on G were used). Consequently, any SPT-based prediction of the surface tension should be considered in light of the EOS arising from the SPT interpolation. To examine the relationship between the predicted surface tension and the basis EOS, we plot γ∞ versus the pressure for several SPT and simulation sources in Figure 4. When plotted in this manner, the large amount of variation in γ∞ between the various forms of SPT is somewhat surprising in light of Table 1. The similarity of our adaptation of SPT to the CS-EOS and SPT613 follows from the close correspondence of the pressure and chemical potential of these two methods. Most interesting is the large difference between our CS-EOS-based method and that of Sun,17 considering their fairly equal predictions of W(λ). As mentioned before, this difference can be attributed to the inclusion of R3(η) in Sun’s interpolation for G. Despite matching

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5495 Table 1. Hard-Sphere Fluid Surface Tension, γ∞σ2/kT, Calculated from Various Forms of SPTa η

SPT312

SPT613

Sun17

H-P29

CS-EOS

K-EOS

S-EOS

0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.49

-0.0216 -0.0604 -0.1343 -0.2653 -0.4886 -0.8626 -1.4854 -2.5279 -3.8631

-0.0218 -0.0617 -0.1380 -0.2738 -0.5070 -0.9034 -1.5825 -2.7796 -4.4365

-0.0227 -0.0640 -0.1432 -0.2829 -0.5187 -0.9073 -1.5420 -2.5802 -3.8800

-0.0219 -0.0614 -0.1365 -0.2683 -0.4905 -0.8575 -1.4588 -2.4472 -3.6917

-0.0218 -0.0611 -0.1369 -0.2729 -0.5091 -0.9142 -1.6099 -2.8206 -4.4378

-0.0222 -0.0627 -0.1413 -0.2828 -0.5282 -0.9472 -1.6613 -2.8904 -4.5127

-0.0222 -0.0628 -0.1415 -0.2831 -0.5290 -0.9510 -1.6794 -2.9677 -4.7460

a The column marked H-P lists the semiempirical predictions of Henderson and Plischke.29 The last three columns are our adaptations of SPT to the CS-EOS,14 K-EOS,23 and S-EOS,24 respectively.

series. Expanding the γ∞ from our adaptation of SPT to the CS-EOS (see eq 50) provides the series

γ∞ σ2 ) - 1.4324η2 - 5.8357η3 + ‚‚‚ kT

(29)

which compares quite well with Bellemans’ “exact” result from cluster theory33

γ∞ σ2 ) - 1.4324η2 - 6.0979η3 + ‚‚‚ kT

(30)

The SPT3 surface tension yields the series approximation

γ∞ σ2 ) - 1.4324η2 - 5.7296η3 + ‚‚‚ kT

Figure 4. Reduced surface tension, γ∞σ2/kT, of the hard-sphere fluid versus the reduced pressure, pσ3/kT. The solid line represents our prediction of SPT adapted to the CS-EOS. The prediction of SPT312 is shown by the dotted line, and the prediction of SPT613 is given by the dashed-dot line. Sun’s SPT prediction17 is given by the short dashed line. Simulation results of Attard and Moule,31 Henderson and van Swol,30 and Heni and Lo¨wen32 are given by the closed circles, open diamonds, and closed triangles, respectively. Error bars associated with simulation data points represent one standard deviation.

the pressure and approximating the chemical potential satisfactorily, Sun’s series includes nonphysical terms that lead to the difference in γ∞. We have also included in Figure 4 simulation results for γ∞ to provide another test of our SPT modification. In general, good agreement with simulation results is observed for our adaptation of SPT to the CS-EOS, but at high pressure our method overestimates the magnitude of γ∞. We note, however, that these predictions fall within the error bars of Henderson and van Swol’s simulation results.30 On the other hand, the simulation measurements of Attard and Moule31 and of Heni and Lo¨wen32 are larger than both our predictions and those of SPT6.13 Overall, the MC simulation results agree better with SPT3 than the more accurate equations of state, a somewhat puzzling result. Both SPT6 and our adaptation either use or generate a highly accurate EOS yet produce surface tension values quite different from existing literature data. Nevertheless, one can argue that the predictions of SPT613 and our SPT adaptation may in fact represent more accurate estimates of γ∞ than those of SPT3, mainly because the former methods yield significantly more accurate predictions of p and µex (as well as include a large number of exact conditions on G). This conclusion is at least supported by considering the expansion of γ∞ in a virial-type

(31)

where the O(η3) coefficient is smaller than that in the other expansions. Now, eqs 29-31 open (cautiously) the possibility that SPT6 and the accurate EOS-based versions of SPT generate better predictions of γ∞, while SPT3 and the current simulation results for γ∞ are less accurate than previously thought. The simulation data points (from all sources in Figure 4) show a significant amount of scatter when plotted versus pressure, making the correct trend difficult to visualize. Furthermore, both SPT6 and EOS-based adaptations of SPT incorporate all known exact conditions on G in contrast to SPT3 and Sun’s SPT. With SPT6 and our EOS-based SPT yielding highly accurate predictions of W (along with very good estimates of the initial slope and curvature of g(r), see the following section), one is initially hesitant to accept that the estimates of the surface tension are significantly in error. If one were to evaluate the accuracy of various versions of SPT using only the surface tension, one would have to conclude that Sun’s SPT modification is superior to SPT6 or the current SPT modification. This conclusion would be incorrect, though, since Sun’s interpolation function for G is not thermodynamically consistent and contains the nonphysical R3 term. Given these new insights provided by SPT6 and our EOS-based modification of SPT, we suggest that the question as to what is the true value of the surface tension of the hard-sphere fluid remains open and the discrepancy between present simulation results and EOS-based theoretical predictions of γ∞ warrants further investigation. Another surface property that can be calculated via SPT is the excess surface adsorption. From the Gibbs adsorption isotherm, the excess surface adsorption is given by28

Γ(λ) ) -

( ) ∂γ(λ) ∂µ

T

(32)

where we make note that, in the above equation, γ and, thereby, Γ are functions of the cavity radius λ. By replacing ∂µ via the

5496

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Figure 5. Reduced excess surface adsorption, Γ∞σ2, versus the reduced pressure, pσ3/kT. The data sources and labels are the same as those in Figure 4.

Gibbs-Duhem relation, letting λ f ∞, and then introducing a chain rule expansion, the excess adsorption on a planar interface is found to be15

Γ∞ )

-6η (∂γ∞/∂η) πσ3 (∂p/∂η)

(33)

Unfortunately, eq 33 does not lend itself to equations of state that do not provide expressions for γ∞ in terms of η. Estimates of Γ∞ from SPT6 and equations of state that do not provide explicit µex will therefore require numerical derivatives of γ∞. Because γ∞ is directly connected to the EOS, we again base our comparisons with existing theory and simulation on the relationship between Γ∞ and p. In Figure 5, we again plot the prediction from our SPT adaptation of the CS-EOS along with other estimations and simulation results. Like the results for γ∞, our predictions of Γ∞ agree well with those of SPT6. The differences of Γ∞ between the various forms of SPT are similar to what is observed for γ∞, pointing to the relationship between Γ∞ and γ∞. Unfortunately, we again observe a statistically significant difference between our theory and previously published simulations results. In light of eqs 29-31, one would also expect the Γ∞ predictions from SPT6 or the CS-EOS to be fairly accurate. By combining Bellemans' series representation of γ∞ with the virial series for Z24 in eq 33, one obtains the following series for Γ∞

Γ∞σ ) 2.8648η - 4.6246η + ‚‚‚ 2

2

3

(34)

Likewise, our CS-EOS adaptation yields,

Γ∞σ2 ) 2.8648η2 - 5.4113η3 + ‚‚‚

|

∂3G ∂λ3

)

[

]

λ)1/2+

48η 1 + 16η + 10η2 2 - 26η + g(1 ) - 2g′(1+) (37) + 3 1-η 1 η (1 - η) In eq 37, the first term in brackets is identical to the third derivative of G at λ ) 1/2- (using eq 5) and the two other terms represent the discontinuity across λ ) 1/2 (much like the discontinuity of ∂2G/∂λ2 at λ ) 1/2+ that is proportional to g(1+)12). Rearrangement of eq 37 and substitution of G(1,η) for g(1+) provides the following relationship between the initial slope of g and the SPT function G

g′(1+) )

1 + 16η + 10η2 1 - 13η G(1,η) + 1-η 2(1 - η)3 1 - η ∂ 3G 96η ∂λ3

( | )

(36)

Comparing eq 34 to eq 35 suggest that the CS-EOS should provide a better match to simulation data than SPT3. It is interesting to note that Henderson and van Swol combined the SPT3 expression for γ∞ and the CS-EOS in eq 33 to produce a “more accurate” expression for the surface excess adsorption,30 potentially suggesting that Γ∞ from SPT3 is not as accurate as once thought. We argue that γ∞ and p should be obtained from

(38)

λ)1/2+

When combined with some version of SPT, eq 38 generates an approximation for g′(1+). The second derivative at contact is connected to the fourth derivative of G at λ ) 1/2+ in the following manner13

|

∂4G ∂λ4

[

192η 30η + 102η2 + 30η3 1-η λ)1/2+ (1 - η)4 6 + 27η + 147η2 + 15η g′(1+) - g′′(1+) (39) g(1 ) 2 1 η (1 - η) )

]

(35)

with SPT3 giving

Γ∞σ2 ) 2.8648η2 - 5.7296η3 + ‚‚‚

the same EOS when calculating Γ∞ and that Henderson and van Swol’s mixed result is more convenient than correct. Many of the arguments we discussed regarding the true value of the surface tension apply here as well. Considering the amount of exact conditions imposed on G by SPT6 and our EOS-based modification of SPT, we are initially reluctant to disregard their surface adsorption predictions. Again, further investigation of simulation results and theoretical predictions of Γ∞ appears to be required to determine which expression for Γ∞ is truly correct. 4.3. Limiting Derivatives of g(r). Another group of properties that may be extracted from SPT is the derivatives of the radial distribution function, g(r), at contact. The equations that relate G and these derivatives were derived by Heying and Corti13 in their development of SPT6 and SPT7 (one of these equations was given, though not derived, in the original SPT paper12). As we explain later in more detail, our adaptation of SPT to an EOS utilizes these equations in a manner opposite to Heying and Corti; while they used approximations of the initial derivatives of g(r) to constrain G(λ,η), we use G(λ,η) from the adaptation of SPT to an EOS to estimate the derivatives of g(r). By taking derivatives of P0 beyond those used to connect the second derivative of G to g(1+), one finds that13

As in the second and third derivatives of G, the first term in brackets is the fourth derivative of G at λ ) 1/2- (see eq 5), and the terms dependent on g and its derivatives are proportional to the discontinuity across λ ) 1/2. Rearrangement of eq 39 and substitution of eq 38 provides the following expression for the initial curvature of g

g′′(1+) )

45η - 36η2 - 90η3 6 + 48η + g(1 ) + 1-η 2(1 - η)4 5 ∂3G 1 - η ∂ 4G (40) 3 32 ∂λ λ)1/2+ 192η ∂λ4 λ)1/2+

( | )

( | )

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5497 Table 2. Limiting Slope of the Radial Distribution Function, g′(1+), as a Function of the Packing Fraction as Calculated via Simulation and by Various Theoretical Approachesa η

MC113

MC240

SPT613

L40,41

Sun17

CS-EOS

K-EOS

S-EOS

0.157 0.209 0.262 0.314 0.340 0.367 0.393 0.419 0.445 0.471

-1.59(3) -2.85(2) -4.91(2) -8.29(4) -10.77(3) -13.98(2) -18.21(2) -23.81(3) -31.27(4) -41.35(4)

-1.55(2) -2.80(3) -4.90(7) -8.17(10)

-1.58 -2.85 -4.90 -8.26 -10.71 -14.02 -18.12 -23.74 -31.32 -41.75

-1.58 -2.85 -4.88 -8.20

-0.71 -1.17 -1.84 -2.80 -3.43 -4.20 -5.15 -6.30 -7.73 -9.51

-1.47 -2.82 -5.05 -8.69 -11.31 -14.67 -19.01 -24.62 -31.94 -41.57

-1.98 -3.52 -5.93 -9.74 -12.43 -15.84 -20.19 -25.77 -32.99 -42.41

-2.00 -3.54 -5.97 -9.88 -12.68 -16.31 -21.03 -27.25 -35.56 -46.85

-13.86(13) -23.64(27) -41.26(1.34)

-13.67 -22.93 -39.04

a

The final three columns are our predictions from the labeled equations of state. Numbers in parentheses indicate the error in the final significant figure.

Like the recent paper by Robles and Lo´pez de Haro34 that generates the derivatives of g based on the rational function approximation (RFA)35,36 and a basis EOS, eqs 38 and 40 yield estimates of the initial derivatives of g given an interpolation scheme for G along with, in the case of our adaptation of SPT, a chosen basis EOS. We note that eqs 38 and 40 are fairly straightforward, whereas the RFA expressions for g′(1+) and g′′(1+) are quite complex.34 In their development of SPT, Heying and Corti developed physically intuitive approximations of g′(1+) and g′′(1+) using an SPT integral equation,37-39 to obtain the six- and sevencondition forms of SPT.13 We, however, already have a sixcondition “form” of SPT. Instead of introducing additional approximations, we use our current adaptation to solve eqs 38 and 40 directly for these derivatives at contact. Once the interpolation coefficients Ri are obtained by the interpolation scheme described in Section 3, the third and fourth derivatives of G at λ ) 1/2+ are easily calculated, leading to the initial derivatives of g. In Table 2, we list the numerical values of g′(1+) from several sources. The first two columns are the simulation results of Heying and Corti (MC1)13 and Smith et al. (MC2).40 The next columns are the predictions from SPT6 and the integral equation of Labı´k and co-workers (L),41,40 followed by Sun’s modification of SPT to the CS-EOS.17 The last three columns are the predictions from the current adaptation of SPT to several equations of state (CS-EOS, K-EOS, and S-EOS). In general, we observe that our modification of SPT to several equations of state produces values of g′(1+) similar to previous theoretical estimates and simulation results. The percentage difference between SPT6 and our EOS-based results decreases with increasing packing fraction, confirming that the third derivative of each source’s G is predicted with reasonable accuracy. Like SPT6, the CS-EOS-based estimations of g′(1+) appear to be more accurate than the RFA-based prediction (which is not shown in the table). When compared to the integral equation of Labı´k and co-workers and the simulation values, we see similar agreement at lower packing fractions, but the difference at higher packing fractions exceeds 10% only for the S-EOS. The predictions obtained with Sun’s method are significantly different from any of the simulation results or theoretical predictions, pointing to errors in the third derivative of G. These errors may be attributed to the lack of thermodynamic consistency in Sun’s method and the inclusion of the nonphysical R3 term in his interpolation series. The better agreement between our predictions (regardless of the EOS) and the simulation or previous theory reinforces the requirement that a hard-sphere cavity function be thermodynamically consistent.34

Table 3. Limiting Curvature of the Radial Distribution Function, g′′(1+), as a Function of the Packing Fraction as Calculated via Simulation and by Various Theoretical Approachesa η 0.157 0.209 0.262 0.314 0.367 0.393 0.419 0.445 0.471

MC113 MC342 SPT613 3.0(5) 7.3(5) 16.2(4) 35.8(9) 78.7(6) 117(1) 176(1) 264(1) 400(1)

1.7 29.1 64.8 90 143 197 293

11.2 22.3 42.5 79.5 153 204 283 394 557

Sun17 -6.3 -6.90 -12.9 -23.4 -41.8 -55.8 -74.5 -99.6 -134

CS-EOS K-EOS S-EOS 8.08 21.2 46.0 90.5 169 229 308 414 557

21.3 39.7 70.0 120 203 263 343 447 585

21.8 40.2 71.1 123 215 285 381 514 702

a The final three columns are our predictions from the labeled equations of state.

Table 3 lists the values of g′′(1+) for several of the same sources used in Table 2, though an additional column of simulation results, labeled MC3 and generated by Groot et al.,42 is included. As noted by Heying and Corti, there are significant discrepancies between their simulation data13 and that of Groot et al.42 They point out, however, that the accuracy in their values for g′(1+) as compared to theory suggests that their simulation estimates are more accurate than those of Groot et al.42 Nevertheless, our EOS-based predictions of g′′(1+) vary considerably from the simulation data. Our predictions are similar in magnitude to the simulation results but show some variation from EOS to EOS. As compared to the SPT6, the CS-EOS seems to best approximate g′′(1+), while the K-EOS and S-EOS predict larger values of the curvature. This is similar to the RFA-based approximations (again not shown here) for g′′(1+), where the CS-EOS provides better estimates than the S-EOS. The similarity between values of g′′(1+) from SPT6 and CS-EOS reflects the similarity between their basis equations of state and chemical potentials. The simplicity with which we calculate the initial curvature of g, however, makes our method useful in providing an estimate of g′′(1+). In contrast, Sun’s method yields an initial curvature of g that has the wrong sign, corresponding to an incorrect shape of g(r) in the vicinity of r ) 1. Once again, this result stresses the importance of requiring that any modification of SPT be thermodynamically consistent and not include nonphysical terms in the interpolation series. The calculation of g′(1+) is less sensitive to the nonphysical R3 term in Sun’s method, but the errors brought about by its inclusion are amplified when calculating g′′(1+). Finally, the radial distribution function can be approximated over the interval 1/2 e λ < 1/x3 using an equation that appears in the original SPT paper. Given the function

5498

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

F(λ) )

)

P0(λ) - 1 + 8ηλ3 24η2λ3 exp[-24η

∫0λG(r,η)r2 dr] - 1 + 8ηλ3 24η2λ3

(41)

where F(λ) is related to the two-particle contribution to P0, RFL showed that12

g(2λ) )

(

1 ∂ 3F 1 ∂2F 8 ∂F + 7 + 96 ∂λ3 96λ ∂λ2 λ ∂λ

) (

1 1 for e λ < 2 x3

)

(42) Equation 42 provides an exact connection between g(r) and G from particle contact at r ) 1 up to r ) 2/x3. Since G is not known exactly beyond λ ) 1/2, eq 42 yields an approximation of g(r) that depends on the chosen representation of G. In Figure 6, we plot the resulting approximation of g(r) obtained from eq 42 for various versions of SPT and MC simulation at η ) 0.393. Over a very short range, the accurate SPT6 and our adaptation of SPT to three accurate equations of state provide a reasonably good estimation of g(r) differing slightly from simulation between r ) 1.05 and 1.15. The g(r) obtained from Sun’s SPT modification is quite different in that it does not capture the shape of g(r) for any r, even predicting an incorrect g(1+). This was foreshadowed by the incorrect sign of g′′(1+) that Sun’s method predicted earlier (see Table 3), but is also indicative of the conditions used to obtain Sun’s interpolation for G. Because the discontinuity in ∂2G/∂λ2 at λ ) 1/2+ (which is proportional to g(1+)) is not used as a matching condition (see eq 8) in Sun’s approach, eq 42, which is consistent with this discontinuity at λ ) 1/2, will predict a value of g(1+) that differs from the value of g(1+) that is obtained only by using the condition of G(1,η) ) g(1+). In other words, the equality of g(1+) and G(1,η) predicts one value and the discontinuity of ∂2G/∂λ2 across λ ) 1/2 predicts another, unless both of these conditions are enforced with SPT. Combining an incorrect contact value of g with an initial curvature whose sign is incorrect clearly generates an inaccurate prediction of g(r). The more accurate versions of SPT again highlight the importance of a cavity function being thermodynamically consistent and not containing nonphysical terms. Reasonably good predictions of the initial shape of g(r) are obtained, however, when SPT is adapted to an arbitrary hardsphere EOS in a thermodynamically consistent manner. 5. Conclusions We have presented a method of adapting SPT to an arbitrary hard-sphere EOS that utilizes all known exact conditions on G and ensures thermodynamic consistency between the pressure and the chemical potential obtained from the integral of G. This approach produces very accurate values of the work of cavity formation at all packing fractions, comparing favorably with those of SPT6 and other EOS-based modifications of SPT.16,17 Estimates of the initial slope and curvature of the hard-sphere radial distribution function were also developed within the framework of SPT. Our EOS-based predictions of these initial derivatives compare favorably to the predictions of SPT6, though both yield initial curvatures that differ somewhat from the simulation data. The most surprising result of this work is the difference between both the CS-EOS-based and SPT6 predictions of the surface tension and the traditionally accepted estimates of SPT3,

Figure 6. Radial distribution function, g(r), of the hard-sphere fluid over a limited range for η ) 0.393. The closed circles represent MC simulation data points. Predictions based on SPT6 are shown by the dashed-dot line. The short dashed line is based on Sun’s modification of SPT. Our predictions for the CS-EOS, K-EOS, and S-EOS are shown by the solid line, longdashed line, and the dotted line, respectively.

which largely agrees with the simulation data currently found in the literature.30-32 We argued that SPT6 and the CS-EOSbased results may still, in fact, represent the more correct values of the surface tension. In contrast to SPT3 (and even the thermodynamically inconsistent method of Sun), both SPT6 and our EOS-based SPT approach incorporate all known exact conditions on G, yielding superior estimates of W. The accuracy of the predictions of W, along with the very good estimates of the initial slope and shape of g(r), indicates that the coefficients used in the Laurent series provide a faithful representation of G over a large range of λ. And so, we are reluctant to immediately discredit the EOS-based predictions of γ∞ (and Γ∞). We may, of course, be proven wrong. In the end, further investigation is clearly needed to better understand why the surface thermodynamic predictions of our SPT modification differ from the simulation results. Improvements to the current work may follow from the introduction of additional matching conditions into the SPT interpolation function (like SPT6 and SPT7) or by proposing a new series representation of G. For example, if the surface tension estimate from SPT3 is proven to be accurate, this value of γ∞ could be used as another matching condition (just like the imposition of a given EOS). Alternatively, the SPT6 approximation of the initial slope of the radial distribution function (which is quite accurate) could provide another matching condition. New interpolation functions for G, suggested by new insights into the thermodynamics of hard spheres or by extending the current Laurent series approximation, can be proposed. For example, one can use two connected, but separate, Laurent series to represent the SPT function,43 though the chosen cavity radius at which the two series are linked is arbitrary. Finally, all of the ideas presented above can be straightforwardly extended to the two-dimensional version of SPT,43,44 though the divergence of ∂2G/∂λ2 at λ ) 1/2+ presents a larger challenge than the simple discontinuity in ∂2G/∂λ2 at λ ) 1/2 presents in the three-dimensional SPT. Acknowledgment This work was supported by the National Science Foundation under Grant No. 0133780.

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5499

Appendix: Analytic Expressions for ri(η) In Section 3, the interpolation scheme used for adapting SPT to an arbitrary equation of state is described. The coefficients for the Laurent series representation of G are expressed most generally in terms of some given Z ) p/FkT and βµex and are obtained from simultaneously solving eqs 19-24. The results are as follows.

R0(η) ) Z

R1(η) )

R2(η) )

R4(η) )

(43)

-48 + 251η - 295η2 + 110η3 + Z(1 - η)2(48 - 67η + 16η2) 6η(1 - η)3

-

288 - 1577η + 1858η2 - 695η3 - Z(1 - η)2(288 - 361η + 52η2) 24η(1 - η)3

R5(η) )

108 - 656η + 763η2 - 287η3 - 2Z(1 - η)2(54 - 65η + 5η2)

R6(η) )

24η(1 - η)3 -24 + 149η - 172η2 + 65η3 + Z(1 - η)2(24 - 29η + 2η2) 32η(1 - η)3

+

-

6η(1 - η)3

+

-768 + 4541η - 5320η2 + 1997η3 + Z(1 - η)2(768 - 925η + 82η2) 96η(1 - η)3

40(1 - η)3[βµex + ln(1 - η)]

256(1 - η)3[βµex + ln(1 - η)]

-

24η(1 - η)3 736(1 - η)3[βµex + ln(1 - η)] 96η(1 - η)3

106(1 - η)3[βµex + ln(1 - η)] 24η(1 - η)3 24(1 - η)3[βµex + ln(1 - η)] 32η(1 - η)3

(44)

(45)

(46)

(47)

(48)

We again note that the use of the Laurent series (eq 13) requires that R3 ) 0. With this set of equations, one need only insert the pressure and excess chemical potential to predict any of the thermodynamic properties discussed in Section 4. Explicit expressions in terms of η may be obtained when the equation of state and chemical potential are both explicitly known functions of η. For example, the Carnahan-Starling equation of state yields an equation for βµex in terms of η. When these CS-EOS expressions for the pressure and chemical potential are entered into eqs 43-48, the following coefficients are obtained.

R0(η) )

R1(η) )

R2(η) )

R4(η) )

1 + η + η2 - η3 (1 - η)3

-40η + 131η2 - 174η3 + 93η4 - 16η5 - 40(1 - η)4 ln(1 - η) 6η(1 - η)4

128η - 448η2 + 558η3 - 243η4 + 26η5 + 128(1 - η)4 ln(1 - η) 12η(1 - η)4

-368η + 1288η2 - 1563η3 + 609η4 - 41η5 - 368(1 - η)4 ln(1 - η)

R5(η) )

48η(1 - η)4

106η - 371η2 + 450η3 - 171η4 + 10η5 + 106(1 - η)4 ln(1 - η)

R6(η) )

24η(1 - η)4 -12η + 42η2 - 51η3 + 19η4 - η5 - 12(1 - η)4 ln(1 - η) 16η(1 - η)4

(49)

(50)

(51)

(52)

(53)

(54)

5500

Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006

Note Added after ASAP Publication. In the version of this paper that was published on the Web 1/6/2006, there was an error in the numerator of the first term in eq 40. The correct version was published on the Web 1/24/2006. Literature Cited (1) Pratt, L. R.; Pohorille, A. Theory of hydrophobicity: Transient cavities in molecular liquids. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 2295. (2) Lum, K.; Chandler, D.; Weeks, J. D. Hydrophobicity at small and large length scales. J. Phys. Chem. B 1999, 103, 4570. (3) Henderson, J. R. Solvation of a solvophobic sphere. J. Chem. Phys. 2002, 116, 5039. (4) Ho¨finger, S.; Zerbetto, F. The costly process of creating a cavity in n-octanol. Theor. Chem. Acc. 2004, 112, 240. (5) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54, 5237. (6) Matyushov, D. V.; Voth, G. A. A perturbation theory for solvation thermodynamics: Dipolar-quadrupolar liquids. J. Chem. Phys. 1999, 111, 3630. (7) Punnathanam, S.; Corti, D. S. Work of cavity formation inside a fluid using free-energy perturbation theory. Phys. ReV. E 2004, 69, 036105. (8) Corti, D. S.; Reiss, H. Depletion force between a colloid particle and a wall: Simple determination by means of scaled particle theory. Mol. Phys. 1998, 95, 269. (9) Siderius, D. W.; Corti, D. S. Extension of scaled particle theory to inhomogeneous hard particle fluids. I. Cavity growth at a hard wall. Phys. ReV. E 2005, 71, 036141. (10) Siderius, D. W.; Corti, D. S. Extension of scaled particle theory to inhomogeneous hard particle fluids. II. Theory and Simulation of fluid structure surrounding a cavity that intersects a hard wall. Phys. ReV. E 2005, 71, 036142. (11) Dinsmore, A. D.; Yodh, A. G.; Pine, D. J. Phase diagrams of nearly hard-sphere binary colloids. Phys. ReV. E 1995, 52, 4045. (12) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. Statistical mechanics of rigid spheres. J. Chem. Phys. 1959, 31, 369. (13) Heying, M. D.; Corti, D. S. Scaled Particle Theory Revisited: New conditions and improved predictions of the properties of the hard sphere fluid. J. Phys. Chem. B 2004, 108, 19756. (14) Carnahan, N. F.; Starling, K. E. Equation of state for nonattracting rigid spheres. J. Chem. Phys. 1969, 51, 635. (15) Mandell, M. J.; Reiss, H. Scaled particle theory: Solution to the complete set of scaled particle theory conditions: Applications to surface structure and dilute mixtures. J. Stat. Phys. 1975, 13, 113. (16) Matyushov, D. V.; Ladanyi, B. M. Cavity formation energy in hard sphere fluids: An asymptotically correct expression. J. Chem. Phys. 1997, 107, 5815. (17) Sun, J. Scaled Particle Theory Improved by Using the CarnahanStarling Equation. Chin. J. Chem. Phys. 1999, 12, 427. (18) Tolman, R. C. Principles of Statistical Mechanics; Oxford University Press: London, 1938. (19) Tully-Smith, D. M.; Reiss, H. Further development of scaled particle theory of rigid sphere fluids J. Chem. Phys. 1970, 53, 4015. (20) Tolman, R. C. The effect of droplet size on surface tension. J. Chem. Phys. 1949, 17, 333. (21) Stillinger, F. H.; Cotter, M. A. Free Energy in the Presence of Constraint Surfaces. J. Chem. Phys. 1971, 55, 3449. (22) Wertheim, M. S. Exact solution of the Percus-Yevick integral equation for hard spheres. Phys. ReV. Lett. 1963, 10, 321.

(23) Kolafa, J.; Labık, S.; Malijevsky, A. Accurate equation of state of the hard sphere fluid in stable and metastable regions. Phys. Chem. Chem. Phys. 2004, 6, 2335. (24) Sanchez, I. C. Virial coefficients and close-packing of hard spheres and disks. J. Chem. Phys. 1994, 101, 7003. (25) Vieceli, J. J.; Reiss, H. Thermodynamics of curved boundary layers. J. Chem. Phys. 1972, 57, 3754. (26) Mandell, M. J.; Reiss, H. Thermodynamics of curved boundary layers. J. Stat. Phys. 1975, 13, 107. (27) Punnathanam, S.; Corti, D. S. Homogeneous bubble nucleation in stretched fluids: Cavity formation in the superheated Leonard-Jones liquid. Ind. Eng. Chem. Res. 2002, 41, 1113. (28) Gibbs, J. W. The Collected Works of J. Willard Gibbs; Longmans, Green and Co.: New York, 1928; Vol. I. (29) Henderson, D.; Plischke, M. Sum rules for the pair-correlation functions of inhomogeneous fluids: Results for the hard-sphere-hard-wall system. Proc. R. Soc. London, Ser. A 1987, 410, 409. (30) Henderson, J. R.; van Swol, F. On the interface between a fluid and a planar wall. Theory and simulations of a hard sphere fluid at a hard wall. Mol. Phys. 1984, 51, 991. (31) Attard, P.; Moule, G. A. A force-balance Monte Carlo simulation for surface tension of a hard-sphere fluid. Mol. Phys. 1993, 78, 943. (32) Heni, M.; Lo¨wen, H. Interfacial free energy of hard-sphere fluids and solids near a hard wall. Phys. ReV. E 1999, 60, 7057. (33) Bellemans, A. Statistical mechanics of surface phenomena. 1. A cluster expansion for surface tension. Physica 1962, 28, 493. (34) Robles, M.; Lo´pez de Haro, M. On the contact values of the derivatives of the hard-sphere radial distribution function. J. Chem. Phys. 1997, 107, 4648. (35) Bravo Yuste, S.; Santos, A. Radial distribution function for hard spheres. Phys. ReV. A 1991, 43, 5418. (36) Bravo Yuste, S.; Lo´pez de Haro, M.; Santos, A. Structure of hardsphere metastable fluids. Phys. ReV. E 1996, 53, 4820. (37) Reiss, H.; Casberg, R. V. Radial distribution function for hard spheres from scaled particle theory, and an improved equation of state. J. Chem. Phys. 1974, 61, 1107. (38) Reiss, H.; Ellerby, H. M.; Manzanares, J. A. Ornstein-Zernikelike equations in statistical geometry: Stable and metastable systems. J. Phys. Chem. 1996, 100, 5970. (39) Chua, L. P. New Integral Equation for the Structure of the Hard Sphere Fluid Based on Scaled Particle Theory. M.S. Thesis, Purdue University, West Lafayette, IN, 2002. (40) Smith, W. R.; Labı´k, S.; Malijevsky, A.; Sedlbauer, J. A new geometrically based integral-equation hierarchy for hard-sphere systems. 2. Pair and triplet background correlation-functions. Mol. Phys. 1994, 83, 1223. (41) Labı´k, S.; Malijevsky, A.; Smith, W. R. A new geometrically based integral equation hierarchy for hard-sphere systems. 1. Basic theory and thermodynamic results Mol. Phys. 1994, 83, 983. (42) Groot, R. D.; van der E erden, J. P.; Faber, N. M. The direct correlation function in hard sphere fluids. J. Chem. Phys. 1987, 87, 2263. (43) Cotter, M. A.; Stillinger, F. H. Extension of scaled particle theory for rigid disks. J. Chem. Phys. 1972, 57, 3356. (44) Helfand, E.; Frisch, H. L.; Lebowitz, J. L. Theory of the two- and one-dimensional rigid sphere fluids. J. Chem. Phys. 1961, 34, 1037.

ReceiVed for reView September 15, 2005 ReVised manuscript receiVed November 11, 2005 Accepted November 21, 2005 IE051038T