11678
J. Phys. Chem. B 2001, 105, 11678-11682
Breakdown of the Capillarity Approximation in Binary Nucleation: A Density Functional Study† Ari Laaksonen*,‡ and Ismo Napari§ Department of Applied Physics, UniVersity of Kuopio, P.O.Box 1627, Kuopio, Finland 70211, and Department of Physics, UniVersity of Helsinki, P.O. Box 64, Helsinki, Finland 00014 ReceiVed: May 1, 2001; In Final Form: August 7, 2001
The binary classical nucleation theory is known to produce unphysical predictions of gas-liquid nucleation in systems in which the surface tension has a strong composition dependence. It has been shown earlier, using thermodynamic methods, that the unphysical predictions are related to the capillarity approximation, i.e., the assumption that the critical nucleus surface tension equals the surface tension of a flat interface. In this paper, we apply the density functional theory of nucleation to a surface active system and show by numerical calculations that the capillarity approximation breaks down severely for critical nuclei of intermediate composition. In these nuclei, the concentration profile of the surface active component shows a strong maximum in the nucleus surface layer.
I. Introduction The classical nucleation theory (CNT) is based on the assumption that the surface tension of a small liquid droplet equals the surface tension of a flat surface: the capillarity approximation. Although various attempts have been made in order to include the size dependence of surface tension1 into phenomenological nucleation theories, the resulting theories2-5 have not been much more successful than the CNT. In onecomponent systems, the CNT describes the gas-liquid transition often amazingly well, predicting nearly correct supersaturation dependences of the nucleation rate, although the predicted temperature dependences mostly disagree with experimental findings. However, in binary gas-liquid nucleation, the CNT breaks down completely, producing unphysical predictions of the properties of critical nuclei in systems where the surface tension has a strong composition dependence. The binary CNT was first derived by Neumann and Do¨ring,6 who applied it to the water-ethanol system and discovered that at certain vapor activities the theory predicts that the total number of water molecules in the critical nucleus becomes negative. The work of Neumann and Do¨ring was forgotten for several decades after World War II. The next development in binary CNT was the seminal 1950 paper by Howard Reiss,7 who considered saddle surfaces formed by binary clusters in the free energy space. Perhaps luckily, Reiss only studied systems which did not have compositional gradients of surface tension, and he thereby avoided tackling the unphysical predictions. After Wilemski8 had rederived the binary CNT equations of Neumann and Do¨ring, it was soon found9 that activity plots (plots showing how the gas phase activities of the two species can be varied keeping the nucleation rate constant) show strange humps instead of monotonic behavior in systems with strong surface tension gradients. These humps were subsequently10 shown to correspond to negative excess number of particles over the uniform vapor in critical nuclei. †
Part of the special issue “Howard Reiss Festschrift”. * To whom correspondence should be addressed. ‡ University of Kuopio. § University of Helsinki.
We have recently11 shown, using thermodynamic argumentation, that the unphysical predictions result directly from the failure of the capillarity approximation. In the present paper, we make density functional calculations on the properties of critical nuclei in a surface active binary system and corroborate the conclusion that the capillarity approximation fails badly when the surface tension depends strongly on droplet composition, leading to predictions of negative excess number of particles for the species with higher surface tension. II. Density Functional Theory The binary fluid system under consideration consists of spherical hard-core atoms (monomers) and molecules composed of two fused hard spheres (dimers). To study this system, we employ density functional theory (DFT) together with the interaction site model in the local density approximation. In this approach, the grand potential functional Ω of the inhomogeneous fluid is written as12
βΩ[F0(r), F1(r), F2(r)] )
∫ dr F0(r) ln F0(r) +
2
∫ dr Fj(r) ln fj(r) + ∫ dr Ψ[η(r)]Fs(r) ∑ j)1 ∫ ∫ dr dr′ s(2)(|r - r′|) f1(r) f2(r) + β
2
∑ ∫ ∫ dr dr′φLJjk (|r - r′|)Fj(r) Fk(r′) -
2 j,k)0
2
β
µj ∫ Fj(r) dr ∑ j)0
(1)
where β ) 1/kBT with T the absolute temperature and kB Boltzmann’s constant. Here F0(r) and Fi(r) are the number densities of monomers (index 0) and the two parts of the dimer molecule or sites (indices 1 and 2). The total density is Fs(r), 2 VFi(r), where and the packing fraction is given by η(r) ) ∑i)0 2 V ) π/6σ is the molecular volume for a hard-sphere, which is the same for all components, because we only consider hard
10.1021/jp0116454 CCC: $20.00 © 2001 American Chemical Society Published on Web 10/09/2001
Breakdown of the Capillarity Approximation
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11679
spheres of the same size σ. The chemical potential of the hard sphere component i is denoted by µi. The auxiliary functions fi(r) are related to the activities of the molecular sites. The first two terms on the right-hand side of eq 1 comprise the ideal-gas part, whereas the third and the fifth terms are contributions from the repulsive hard-sphere system and the attractive perturbation, respectively. The excess free energy per particle over the ideal gas Ψ is given by the formulas of Carnahan and Starling.13 The perturbative system is treated according to Weeks-Chandler-Andersen scheme14 with the attractive interactions obtained from the Lennard-Jones (LJ) potential. The bonding between the two particles in the dimer molecule cause a decrease in entropy which is included in the fourth term, where s(2) is an intramolecular correlation function. The dimers are assumed to have a rigid bond length σ; therefore, the correlation function is represented by s(2)(r) ) (4πσ2)-1δ(r - σ). Minimizing eq 1 with respect to functions Fi(r) and fi(r) gives a set of integral equations12 from which the density profiles of planar gas-liquid interfaces or spherical droplets can be solved in respective geometries.15 For homogeneous states these equations reduce to particularly simple forms. Noting that the densities of the dimer sites are now equal (F1 ) F2 ) F), we can write the chemical potentials as
µ0 ) kBT ln F0 + U0
(2)
µ ) µ1 + µ2 ) kBT ln F + U1 + U2
(3)
where functions Ui are defined as 2
βUi ) Ψ + Ψ′VFs - β
RijFj ∑ j)0
(4)
with Rij ) -∫ drφLJ ij (r). Pressure is obtained from
P ) kBT(F0 + F) + kBTΨ′VF2s -
1
2
∑ RijFiFj
2 i,j)0
(5)
The surface tension of a planar gas-liquid interface γ∞ is calculated from the equilibrium density profiles as
γ∞A ) Ω[{Fi(r)}] - Ωe
(6)
where A is the area of the interface and Ωe is the grand potential evaluated at the coexisting bulk densities which are obtained by equating the chemical potentials from eqs 2 and 3 and pressures from eq 5 for vapor and liquid states. In a similar fashion, the work of formation of a critical droplet ∆Ω is found from
∆Ω ) Ω[{Fi(r)}] - Ωv
(7)
where Ωv corresponds to the grand potential of a uniform vapor at given gas-phase activities ai ) exp(µi/kBT). The desired behavior of the fluid system is obtained by providing appropriate parameters for the LJ potential
φLJ ij (r) ) 4ij
[( ) ( ) ] σij r
12
-
σij r
6
(i and j ) 0, 1, and 2) (8)
We assume the standard arithmetic mean rule for the length parameters σij; however, the energy parameters ij for the cross interactions are modified as
ij ) (1 - kij)xiijj
(9)
where for the present case 11 ) 22 ) 0.6, k12 ) 0, k01 ) -0.3, and k02 ) 0.374 (the energy parameter of the monomer component 00 is denoted as ). This choice results in an enhanced surface activity of the dimer component. Furthermore, the dimer site 1 is strongly attracted to monomers (“monophilic” site), whereas the site 2 experiences much weaker attraction (“monophobic” site), thus inducing a preferred orientation of dimers with respect to monomers at gas-liquid interfaces; however, the bulk liquid phases of this system are homogeneous at all liquid compositions within the local density approximation. All our calculations are performed at temperature T/ ) 0.7 where the liquid is fully miscible but just on the brink of immiscibility. III. Thermodynamics We consider a two-component liquid droplet with a volume V in equilibrium with the vapor phase. The droplet contains g1 and g2 molecules of species 1 and 2, respectively, it is modeled as spherical, and it has a sharp boundary (Gibbs dividing surface) between the uniform liquid and vapor phases. The differences between the real nucleus with a smoothly varying density profile and the droplet model nucleus are expressed as surface excess quantities. The total numbers of molecules, g1 and g2, are then independent of the choice of dividing surface:
gi ) nli - nVi + nsi
(10)
Here nli ) VFli, and Fli denotes the density of species i in the droplet. Correspondingly, nVi ) VFVi, with FVi as the density of species i in the vapor phase; thus, gi represents the excess number of molecules in the droplet over the uniform vapor. The difference between the numbers of molecules obtained by integration over the actual and over the droplet model density profiles is given by nsi. The pressure difference between the vapor and the droplet interior is given by the Laplace equation
∆P ) Pl - PV )
2γ Rt
(11)
where γ is the surface tension of the droplet and the subscript t of the radius Rt denotes the surface of tension. The reversible work of formation of the droplet is
W ) -∆PVt + 4πR2t γ
(12)
We now proceed to derive equations for the droplet radius and work of formation that make use of chemical potential rather than pressure differences. The liquid-phase chemical potentials are related to pressure via
dµli ) Vli dPl
(13)
at constant temperature and liquid composition; Vli denotes the partial molecular volume of species i. Making the assumption of incompressible liquid, we can integrate eq 13 from PV to Pl to obtain
Vli∆P ) -∆µi
(14)
where ∆µi ) µli(Pl) - µli(PV) ) µVi(PV) - µli(PV), and the latter equality follows from the equilibrium between the droplet and the vapor. From eq 14, it follows that
11680 J. Phys. Chem. B, Vol. 105, No. 47, 2001
∆µ1/Vl1 ) ∆µ2/Vl2
Laaksonen and Napari
(15)
which can be used to numerically determine the nucleus composition. Substituting for ∆P in eqs 11 and 12, we obtain the Kelvin equation for the droplet radius and an equation for the work of formation as
Rt )
2γVli ∆µi
(16)
W ) (4π/3)Rt2γ
(17)
The usual form of the classical nucleation theory is now obtained by making the assumption that the surface tension of the droplet, γ, equals the surface tension of a flat interface, γ∞. We have recently shown11 that in binary systems the requirement that surface tension is independent of droplet curvature is equivalent to the requirement that the surface of tension coincides with the special equimolar surface denoted by radius Re at which the following condition is fulfilled:
ns1V1 + ns2V2 ) 0
Figure 1. Activity plots for two different constant nucleus free energies. See text for details.
(18)
When eq 18 holds at the surface of tension, we can write R ) Rt ) Re and
R)
2γ∞Vli ∆µi
(19)
W ) (4π/3)R2γ∞
(20)
which are the binary CNT equations for the radius and work of formation of the critical nucleus. IV. Results and Discussion The unphysical behavior produced by the binary CNT for our model system can be seen in Figure 1. On the curves of Figure 1, the free energy of the critical cluster rather than the nucleation rate is kept constant; however, the variation of the rate along a curve is modest. The density functional curves behave monotonically and have a negative slope throughout the activity space, whereas the upper classical curve exhibits a region of positive slope. The CNT prediction is counterintuitive: if the dimer activity is kept constant (e.g., at aD ) 1.5) and monomer vapor is added continuously into the system, the nucleation rate must first increase, then decrease, and then increase again in order for the activity plot to show a hump. Furthermore, in a binary system, the nucleation theorem16,17 leads to:18
( ) ∂µV1 ∂µV2
∆Ω
)-
g2 g1
(21)
indicating that a positive slope in an activity plot corresponds to a negative excess number of particles gi in a critical cluster, which of course is not only counterintuitive but also unphysical. In our recent paper,11 we made classical nucleation calculations on the water-ethanol system and traced the reason for the appearance of the negative excess number of particles gi back to the assumption that Rt ) Re (or, in other words, that γ ) γ∞). In Figures 2-4, we present direct numerical evidence from DFT calculations supporting this conclusion. The figures present surface tension ratios γ/γ∞ and differences between the two radii, (Re - Rt)/Rt. The surface tension of the droplet, γ,
Figure 2. Ratios of nucleus surface tension to flat interface surface tension as a function of gas-phase monomer activity. The markers correspond to nuclei shown in Figure 1.
and the radius to the surface of tension, Rt, were obtained from eqs 11 and 12. Note that the liquid pressure, Pl, applied in these calculations was obtained from the equation-of-state relations 2-5, giving the pressure of a liquid phase in equilibrium with the supersaturated vapor, and it does not necessarily have the same value as the pressure in the center of the actual density functional nucleus. Similarly, when calculating the radius to the equimolar surface, Re, the reference density and composition were taken to be those of a liquid in equilibrium with the supersaturated vapor, which generally are somewhat different from the density and the composition in the center of the nucleus. The equilibrium composition obtained from eqs 2-5 was also applied in calculation of γ∞. Figure 2 shows the surface tension ratios γ/γ∞ for the critical nuclei corresponding to the activity plot of Figure 1. With the ∆Ω/ ) 60 curve, the ratio decreases starting from am ) 0, and the decrease becomes quite steep at monomer gas-phase activities between 1.25 and 1.5. As can be seen from Figure 1, the slope of the activity plot becomes positive at am ) 1.25. All in all, the classical and density functional activity plots are quite far from each other in the monomer activity range of 1.253.5, where the surface tension ratio is below 0.95. The γ/γ∞
Breakdown of the Capillarity Approximation
Figure 3. Difference of the radii to the equimolar surface and to the surface of tension, normalized by nucleus size, as a function of monomer activity. The markers correspond to the nuclei shown in Figure 1.
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11681 can be expected based on Figure 2, the difference for ∆Ω/ ) 60 increases rapidly at the monomer activity range 1.25-1.5. The maximum difference between Re and Rt is roughly 15% of the nucleus radius, which in absolute terms equals 0.65σ. For ∆Ω/ ) 400, the maximum value of Re - Rt is 0.45σ, i.e., not very much smaller than that for ∆Ω/ ) 60. However, this time the difference is only 4% of the nucleus radius, indicating that in relative terms the capillarity approximation is reasonable, unlike with the smaller nuclei. The variation in radius as a function of aM does not have any qualitative effect on Figure 3 because Rt is a monotonically decreasing function; if we consider clusters of constant size, Re - Rt still shows a maximum. Figure 4 shows density profiles of the critical nuclei at four different monomer activities of Figure 3. The three curves shown in each panel correspond to monomer density (M) and to the densities of monophilic (D1) and monophobic (D2) ends of the dimer molecules. At dimer-rich vapors, the profiles indicate that the surface layers of the nuclei are depleted of monomers, the dimer densities are monotonically decreasing from nucleus center toward the vapor, and the dimers show very little orientational ordering. Above aM ) 1.0, a flip-over of the density profiles in the center of the nuclei takes place, and at aM ) 1.75, most of the dimers are already concentrated in the surface layers with monophobic ends oriented toward the vapor. The worst failure of the capillarity approximation takes place with clusters having a high concentration of dimers on the nucleus surface and a depletion of dimers in the interior. The orientation of the dimers may contribute to the large nonzero values of Re - Rt somewhat, but it cannot be the main source, because the unphysical behavior can also be seen in surface active systems in which both molecular species are spherical.19 V. Conclusion We have studied the surface tensions of binary critical nuclei in a system in which the planar surface tension exhibits a strong compositional gradient. The ratio of cluster surface tension to planar surface tension exhibits a strong minimum at intermediate nucleus compositions. The deviation of the nucleus surface tension from the planar value is related to the deviation of the surface of tension (Rt) from the equimolar surface (Re) specified by the condition ns1V1 + ns2V2 ) 0. Our calculations show that the minimum in the surface tension ratio is accompanied by a strong maximum in the ratio (Re - Rt)/Rt. The maximum occurs in nuclei in which the surface active component shows a density maximum in the nucleus surface layer. The breakdown of the capillarity approximation in this manner leads to the well-known unphysical predictions produced by the binary classical nucleation theory. Acknowledgment. A.L. acknowledges financial support by the Academy of Finland (Project No. 44278). References and Notes
Figure 4. Density profiles of monomers (M) and the monophilic (D1) and monophobic (D2) ends of dimer molecules in critical nuclei at ∆Ω/ ) 60.
curve for ∆Ω/ ) 400 shows also a minimum, but a less pronounced one; this time the surface tension ratio is mostly above 0.95. However, it can be noted from Figure 1 that the CNT and DFT activity plots start departing somewhat as the decrease of the surface tension ratio becomes steeper. Figure 3 shows the differences in the surface of tension and equimolar radii as functions of monomer gas-phase activity. As
(1) Tolman, R. C. J. Chem. Phys. 1949, 17, 333. (2) Dillmann, A.; Mayer, G. E. A. J. Chem. Phys. 1991, 94, 3872. (3) Laaksonen, A.; Ford, I. J.; Kulmala, M. Phys. ReV. E 1994, 49, 5517. (4) Kalikmanov, V. I.; van Dongen, M. E. H. J. Chem. Phys. 1995, 103, 4250. (5) Talanquer, V. J. Chem. Phys. 1997, 106, 9957. (6) Do¨ring, W.; Neumann, K. Z. Phys. Chem. A 1940, 186, 193. Neumann, K.; Do¨ring, W. Z. Phys. Chem. A 1940, 186, 203. (7) Reiss, H. J. Chem. Phys. 1950, 18, 840. (8) Wilemski, G. J. Chem. Phys. 1984, 80, 1370. (9) Garnier, J. P.; Mirabel, P.; Migault, B. Chem. Phys. Lett. 1985, 115, 101.
11682 J. Phys. Chem. B, Vol. 105, No. 47, 2001 (10) Laaksonen, A.; Kulmala, M.; Wagner, P. E. J. Chem. Phys. 1993, 99, 6832. (11) Laaksonen, A.; McGraw, R.; Vehkama¨ki, H. J. Chem. Phys. 1999, 111, 2019. (12) Napari, I.; Laaksonen, A.; Strey, R. J. Chem. Phys. 2000, 113, 4476. Napari, I.; Laaksonen, A.; Strey, R. J. Chem. Phys. 2000, 113, 4480. (13) Mansoori, G. A.; Carnahan, N. F.; Starling, K. E.; Leland, T. W. J. Chem. Phys. 1971, 54, 1523.
Laaksonen and Napari (14) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237. (15) Napari, I.; Laaksonen, A. J. Chem. Phys. 1999, 111, 5485. (16) Viisanen, Y.; Strey, R.; Reiss, H. J. Chem. Phys. 1993, 99, 4680. (17) Oxtoby, D. W.; Kashchiev, D. J. Chem. Phys. 1994, 100, 7665. (18) Oxtoby, D. W.; Laaksonen, A. J. Chem. Phys. 1995, 102, 6846. (19) Laaksonen, A.; Oxtoby, D. W. J. Chem. Phys. 1995, 102, 5803.
