J. Phys. Chem. 1995,99, 2857-2864
2857
Phase Separation of Asymmetric Binary Hard-Sphere Fluids: Self-consistent Density Functional Theory? Yaakov Rosenfeld Nuclear Research Centre-Negev, P.O. Box 9001, Beer-Sheva 84190, Israel Received: June IO, 1994; In Final Form: September 27, I994@
It was recently shown (Phys. Rev. Lett. 1994, 72, 383 1) that the “test-particle-self-consistent” fundamentalmeasure free energy density functional for general inhomogeneous hard-sphere fluids (J. Chem. Phys. 1993, 98, 8126) predicts phase separation in bulk binary mixtures with large size ratios, R ~ R >I -4, when the packing fractions for the two species are comparable, in qualitative agreement with recent experiments on “nearly hard-sphere” colloidal particles. These preliminary calculations are focused on the spinodals and not on the equilibrium coexistence curve, which is much harder to calculate. To prepare the necessary tools for undertaking the complete thermodynamic calculation, the present paper provides a full description of this self-consistent density functional for mixtures and clears many fine points in the theory. It also gives a more detailed presentation of the results relevant for phase separation. These results should be useful for studying external field effects (e.g., sedimentation and confining walls) on the entropically driven demixing transition, as encountered in real experiments on colloids.
Phase separation in a fluid binary mixture of hard-spheres is a long-standing fundamental problem in statistical mechanics,’ which is the object of renewed theoretical and experimental In particular, it was recently shown9 that the “test-particle-self-consistent” fundamental-measure free energy density functional for general inhomogeneous hard-sphere predicts phase separation in bulk binary mixtures with large size ratios, Rz/Rl > -4, when the packing fractions for the two species are comparable, in qualitative agreement with recent experiments on “nearly hard-sphere” colloidal particles. These preliminary calculations are focused on the spinodals and not on the equilibrium coexistence curve, which is much harder to calculate. To prepare the necessary tools for undertaking the complete thermodynamic calculation, the present paper provides a full description of this self-consistent density functional for mixtures and clears fine points in the theory. It also gives a more detailed presentation of the results relevant for phase separation. Fluid phase separation in simple (“atomic”) mixtures is usually induced by enthalpic driving forces: steric effects due to nonadditive effective diameters, and van der Waals attractions. The hard-sphere demixing, however, can be driven only by entropic effects and is part of the general question of how phase separation can be induced by purely repulsive interactions.’* The physical origin of hard-sphere fluid phase separation is the osmotic depletion effect:I3 when the surfaces of two large spheres approach within a distance less then the diameter of a small sphere (the small spheres are expelled from between the large ones), the pressure exerted by the small spheres on the outer surfaces is no longer compensated, leading to an effective attraction between the large spheres. Ever since Lebowitz and Rowlinson’ showed that within the Percus-Yevick (PY) closure of the Ornstein-Zernike equations, hard-sphere fluid mixtures ~~~
~~
~~
It is my great pleasure to contribute to this issue in honor of Stuart Rice. The research described below has originated in part from work which I have done at the University of Chicago during a Sabbatical year with Stuart’s group. I have benefitted, like many others, from conversations with Stuart and from his wise guidance. That earlier work could not have been pursued without Stuart’s generosity, and his patience with the viruses which I left on his personal computer. @Abstractpublished in Advance ACS Abstracts, February 1, 1995.
are completely miscible for all concentrations and size ratios, it has been generally believed that hard-sphere fluids never phase separate. This was contested recently by Biben and Hansen,2 who found that certain more accurate closure^'^ predict phase separation for dense binary hard-sphere fluid mixtures when the size ratio differs considerably from 1, and the partial packing fractions of the two species are comparable. Motivated by this result, the existence of an entropically driven demixing transition was demonstrated in simple, exactly solvable, two-dimensional hard-core lattice model^.^,^ Direct tests of the results from these approximate theories, using numerical simulations of the relevant systems with highly asymmetric size ratios, are precluded, however, by severe ergodicity p r o b l e m ~ even ’ ~ for a ratio R2/ Rl = 3. The hard-sphere demixing is supported by recent experiments on asymmetric “nearly hard-sphere” sterically stabilized6 and charge-~tabilized’.~ colloidal suspensions. Bulk phase separation was observed into two disordered phases, and a new ordered phase was found to be located on the cell walk8 The analysis of these experimental results requires a theoretical tool which is able to take into account the effects of external fields, in particular gravity and confining walls. This paper offers a new and general theoretical framework for addressing the question of hard-sphere phase separation in extemal fields, by employing a comprehensive free energy functionallo.’ for inhomogeneous hard-sphere fluid mixtures. Since a description of the density functional theory can be found in excellent review articles (e.g. ref 16), section I1 focuses on the general new method in density functional theory,” namely, the “test-particle” optimization, by which any free energy model can be optimized at second order by imposing consistency in the test-particle limit, and provides the appropriate extensions of ref 11. In section I11 this general method is applied to bulk hard-sphere fluid mixtures by using, specifically,the fundamentalmeasure free energy After testing the high accuracy and self-consistency of this model in describing pair correlations for bulk hard-sphere mixtures, I confirm that bulk fluid mixtures of hard-spheres are unstable against phase separation for sufficiently asymmetric size ratios. The results are compared with the experimental data and discussed in section IV. To
0022-365419512099-2857$09.00/0 0 1995 American Chemical Society
2858 J. Phys. Chem., Vol. 99, No. 9, 1995
further supplement the previous9 short account (without most of the details) of the phase-separation work and in order to make the paper self contained, explicit expressions for the model free energy functional and for the model bridge functional are given in the Appendix. The present calculations are in qualitative agreement with the experimental results on colloidal suspensions, correctly reproducing their main characteristics. These results open the way (by using the same free energy functional) for studying the entropically driven demixing also for nonuniform hardsphere fluids and, in particular, for investigating the effects of gravitational settling and confining walls, as encountered in experiments with colloidal particles.
Rosenfeld
-
where gi(7) = ei(F)/ei,o,cf3”[{em,o};(17 - r’l)] is the bulk direct correlation function as obtained from the second functional derivative of the excess free energy functional, eq 5, in the uniform fluid limit, and
11. Density Functional Theory and ‘Test-Particle” Self-consistency
A general powerful method for both uniform and nonuniform fluids is the density functional theory.I6 The density profiles {e“} for the fluid subject to external potentials {umO} which couple to the particles of type {m, m = 1,2, ...,M } are obtained by solving the Euler-Lagrange equations:
which correspond to the minimization of the grand potential Q[{emO>l:
where pi are the chemical potentials. The ideal-gas free energy is given by the exact relation
where Ai = (h2/27”kBT)”2 are the de Broglie wavelengths. The central quantity in the density functional theory for nonuniform fluids is the excess (over “ideal gas” contributions) free energy, Fex[(em(F)}], which originates in interparticle interactions, and from which many equilibrium properties of the fluid (e.g., tensions at interfaces, solvation forces for confined fluids, phase transitions for different inhomogeneities) can be derived. It is a unique functional of the spatially varying one particle densities, {em@)},which is in general unknown. A hierarchy of direct is given by functional derivatives correlation functions c(”sFD) (FD) of the excess free energy functional. In particular the one particle ( n = 1, the excess chemical potential) and pair ( n = 2) direct correlation functions are
defines the “Bridge” functional for the hard sphere mixture. pi,exe em,^}] is the uniform, bulk fluid limit of the excess chemical potentials pi,ex[{pm@)};71. The bridge functional Bi[{em(F));a contains all the terms beyond the first order in the functional Taylor expansion, around the uniform fluid limit, of the excess chemical potential functional pi,ex[{@mO};a/kBT,is thus related to the sum of all terms beyond second order in the corresponding functional Taylor expansion of the excess free energy. The density profile equations (6, 7) have the form of the modified hypernetted-chain equations for the density profile.” By truncating the expansion of the excess free energy after second order, the bridge functional vanishes, and the density profile equations (6, 7) then have the hypemetted chain approximation form. There is an underlying connection, originally due to Percus,I6 between the structure of the nonuniform fluid and that of the corresponding bulk, uniform fluid. This connection is made in the “testparticle limit” of the density profile equations: When the external potential is obtained by fixing a test particle of type tat the origin, ui(r) = &i(r), where &i(r) is the corresponding pair potential between particles of types t and i in the fluid, then the density profiles normalized to unity at large r, correspond to the pair distribution functions in the bulk uniform fluid g,i(r) = ei(r)/ei,o,where {em,o}are the average densities of the bulk fluid. The test particle limit of the density profile equations takes the form“
where hi(r) = gtj(r) - 1, and cf,”(r) are the uniform fluid, bulk limit of the direct correlation functions as obtained from the second functional derivative (FD) of the excess free energy functional. The symmetrized bridge function
is obtained as the appropriate weighted bulk average of the bridge functions bg(r), which are derived from the bridge functional Bi[(Qm(F)};aby using ei(r) = gi,ogfi(r): For a fluid in contact with a reservoir bulk fluid, of average densities (em,o}, the density profile equations (1) can be written in the following form:
W )= Bi[{em,ogrm(r>>;r1 The bridge functional” is obtained from the sum of all the terms
J. Phys. Chem., Vol. 99, No. 9, 1995 2859
Asymmetric Binary Hard-Sphere Fluids beyond the second order in the functional Taylor expansion, around the uniform fluid limit, of the excess free energy functional. The exact free energy functional must obey the “testparticle self-consistency”: the exact gij(r)’s as obtained from the solution of the coupled density profile equations (8, 9) are identical to those obtained from the Omstein-Zemike relations:
When we consider an approximate free energy model, its accuracy can be inferred from comparison of the pair correlations, { g y ( r ) } , as obtained from the solution of the “testparticle” density profile (DP) equations (8, 9), with those obtained by the Omstein-Zemike (OZ) equations (lo), { g y F D ( r ) }using , the direct correlation functions {cFsFD)(r)}. This consistency property will be referred to as “the test particle self-consistency”. It is somewhat analogous to the “pressurecompressibility consistency” of integral equation theories for bulk correlations. The test-particle consistency can be exactly satisfied only if we use the exact free energy functional. Given a model free energy, this approximate free energy functional can be optimized at second order by keeping the bridge functional, while imposing the test-particle self-consistency (SC). The resulting coupled equations (8-lo), for both {gij(r)} and {cv(r)},in which the bridge functions are obtained from given fixed bridge functionals, eq 7, represent a new self-consistent method for calculating pair correlations for bulk fluid mixtures, (gic(r)}, which also provide the required input (equation of state ~ ~ ’ [ { Q ~ {, , ~L ~} ~] , ,~ ’ [ { Q ~ and , ~ } ]pair } , direct correlations, { c ~ c ’ [ { e m , o } ; r } )for the optimized model excexss free energy of the nonuniform fluid. Starting with a model excess-freeenergy functional F e x [ { ~ m min } ]eq 2, then after the test-particle optimization it takes the form
~ ’ [ { Q ~ -F ~Ie xI[ {Ie m ( 7 > 1 1
= AFex[{em,O)l
kBT
kBT
+
i
The optimized free energy functional is obtained by replacing the first three terms in the expansion of the nonoptimized functional with the self-consistent terms. There is no attempt to impose any specific structural-thermodynamic consistency relations, everything is predetermined by the quality of the given free energy functional. In particular, the fundamental-measure free energy model, which we apply, responds well to this optimization.
111. Optimized Fundamental-Measure Free-Energy Model-Functional The present approach employs a free energy functional for the inhomogeneous hard-sphere fluid mixture which is based on the fundamental geometric measures of the As a special case, this functional derives in a unified way the Percus-Yevick (PY) direct correlation functions’ and the scaledparticleIg equation of state for the uniform hard-sphere fluid mixture. This free energy model has already been tested very successfully by comparison with computer simulations of density profiles for a large variety of situations where size or packing effects play an important role.’O~ll~ls In view of ref 1 this free energy without test-particle optimization does not predict phase separation for bulk hard-sphere mixtures. Using the optimization by the test-particle method,” as described above, this free energy is presently applied to bulk binary mixtures of hardspheres, of radii R1, R2 and partial densities ei,o ( i = 1, 2), with special emphasis on large size asymmetry ( R I > q,only the small sphere test-particle results are relevant. The results are summarized in Figure 1 where A for the size ratios s = R2/R1 = 10 is plotted as function of 3 2 for different values of 7’1,showing the same general behavior as found in ref 2. A similar picture (Figure 2 ) is obtained for other values of s. As 72 is increased for 3’1 fixed below some critical value, the ratio A may drop below the ideal mixture value of 1 but then it increases again. As 7 2 is increased for a 7’1 fixed above
-
J. Phys. Chem., Vol. 99, No. 9, 1995 2861
Asymmetric Binary Hard-Sphere Fluids
-
0.5/-,
0.0
0.1
0.2
0.3
0.4
?“large Figure 3. Spinodal lines in the (rl,r2)= (rsmdrqarge) plane for different values of the size ratio s = Rz/R,. The open diamonds represent the s = 6 experimental results for silica particles (covering triangles and squares in Figure 1 of ref 6). The ants and open triangles represent the s GZ 7 experimental bulk-phase separation results for polystyrene spheres (same symbols as between the middle and upper lines in Figure 3 of ref 8). The reference dotted line, = V I 72 = 0.54, is provided as a visual aide to mark the s = 1 melting, and to show that the spinodals obey 7 = V I 72 = constant.
r
diagnostic A toward its instability value of A 0. Yet, the main change from the Percus-Yevick results, namely, the enhanced depletion attraction is achieved already by the starting value gij,o(r),with further iterations producing only relatively small refinements. We observed that the bridge functional is not sensitive to details of the pair functions, b,(r) = b,,o(r), and noted the slow variation of the cavity distribution function ) b ~ ( rE ) along the iterations process, Ho(r) hO(r) - c ~ ( r h,PY( I ) - c V ( r ) - b~,o(r).These observations of the generally slow variation of the bridge function and of the cavity distribution functions throughout the iteration process, provide the basis for the present simplified analytic calculation of the spinodal. From these we obtain the following first-order correction to the Percus-Yevick direct corrections:
q j ( 4 = cY(r>+ (gij,o(r>- gY(r>> Inserting this approximation into the diagnostic equation
D,(O) = 0
+
+
some critical value, concentration fluctuations become very strong and S,,(O) tends to diverge causing severe numerical problems, and we could not make the numerical iterative procedure for obtaining gLJ,lto converge beyond some maximal value of 72. But any reasonable extrapolation of the data indicates that S,,(O) will go to zero. Similar trends are observed for the calculations at constant pressure. Extrapolating the results for A down to A = 0, we obtained the spinodal lines (Figure 3). For any given value of the size ratio s = R21R1 the critical pressure P,, below which no phase separation occurs, corresponds to the point with largest value of 72. For s = 10, the case considered in detail in ref 2, we predict a much higher P, ‘ =vs 0.01 < P*, < O.l), corresponding (P*, = ~ P , R I ~ I ~ B T0.3 to a much lower critical concentration of the large spheres (XI 2 0.002 vs x2 E 0.02), Le., a narrow range, 72 < 0.2, in comparison with the wide range, 72 < 0.45, predicted in ref 2. The “self-consistent” spinodals obey approximately 7 = 81 -t 72 constant. For s 5 4 the spinodal region for ourfluid model is in the region where the system is expected to be solid, but a phase-separation instability of the fluid is predicted by our calculations for s 2 4. Like Biben and Hansen, we could not find the second spinodal at higher concentrations of the large spheres.
IV. First-Order Density Functional Correction to the Percus-Yevick Theory The depletion attraction of the large spheres, including its enhancements close to a wall, is built into the PY closure and into the scaled-particle theory which, nevertheless, predict no phase separation. The PY depletion attraction is enhanced by the thermodynamically consistent closures used here in ref 2, enough to produce a spinodal instability. This mechanism becomes clearer through the following simplified density functional calculation. The PY functions are the obvious starting point (recall cf,”(r) = c y ( r ) )for the iterative solution of the self-consistency equations. Inserting q,(r) = g r ( r ) into the right-hand side of eqs 8 and 9, we obtain the first-order results b,,o(r) and gI,,o(r). The iterative process through the solution of the modified hypernetted-chain equations (eqs 8-10) is required in order to gradually build the numerical instability which is associated with the corresponding glide of the
(21)
(22)
we calculate the corresponding approximation for the spinodal. This approximation can be rewritten as Ac&r) = Ag&r)
(23)
and reinterpreted as (e.g.) a mean-spherical approximation:
A ~ ~= -( A~Q ~ ) ’(~>/IC~T
(24)
where the additional depletion interactions on top of the reference PY result is evaluated as
The calculation of the bridge functions using our model bridge functional involves only integrations of the density profiles with simple weight functions. Thus, the whole calculation based on eqs 21 and 22 using eqs 8 involves analytic input functions and several Fourier transforms which are easily performed numerically. The results are presented in Figure 4. The unstable regions, encircled by the spinodals, decrease with decreasing size ratio, s = R21R1 = 10, 6, 4 and eventually disappear for s 3.8. Note that the spinodals from the selfconsistent theory seem to merge with the first-order results at small packing fractions of the large spheres. The s = 10 results from another semianalytic model5 (the solution to their eq 27, compare with their Figure 2) partially agree with our first-order density-functional correction to the PY theory (see Figure 4). Yet, when evaluated also for other values of the size ratio, the spinodals from eq 27 in ref 5 are nearly s independent and with an opposite trend as function of s when compared to the present model. In all these models and in ref 2, the spinodals are driven by what amounts to small corrections (recall that XI
