J. Phys. Chem. B 2001, 105, 11893-11900
11893
Comparison of Cluster Definitions in Homogeneous Nucleation Rate Calculations† Ranjit Bahadur and Richard B. McClurg* Department of Chemical Engineering and Materials Science, UniVersity of Minnesota, Minneapolis, Minnesota 55455 ReceiVed: June 11, 2001; In Final Form: August 10, 2001
We compare cluster definitions used in the framework of classical nucleation theory to estimate nucleation rates. Cluster definitions influence the magnitude of the partition functions used to calculate the equilibrium constants that relate aggregation and fission rate constants. Therefore, cluster definitions affect calculated nucleation rates. Geometric and energetic definitions are simple to apply but are prone to systematic errors. We recommend the dynamic cluster definition as the most likely to lead to calculated rates agreeing with measurements. The definitions are illustrated using a simple xenon model.
I. Introduction When a vapor is cooled and/or its partial pressure increased just beyond the saturation curve, it becomes metastable. Then the thermodynamic equilibrium state is a combination of both liquid and vapor phases. Formation of the liquid phase does not occur instantaneously, however. Instead, density fluctuations larger than a condition-dependent critical size are required to nucleate the liquid phase. The rate at which droplets of the condensed phase form per unit volume and time is called the nucleation rate, (J). If one could continuously observe a metastable vapor at Angstrom spatial resolution and femtosecond time resolution, the majority of the vapor would appear very much like any vapor phase with frequent molecular collisions between individual molecules. Only rarely would one find a group of molecules that are in close proximity and remain so for more than a few rotation/vibration periods. We refer to such groups as clusters. Clusters collide with other molecules on a much longer time scale. Many of these collisions serve only to redistribute energy among the cluster’s translational, rotational, and vibrational degrees of freedom. Occasionally the redistribution of energy is such that a molecule from the cluster leaves the cluster in a fission event. Alternatively, a condensable molecule can undergo an inelastic collision resulting in the growth of the cluster in an aggregation event. Clusters smaller than the critical size tend to undergo fission more readily than aggregation, so only rarely do they grow large enough to become macroscopic liquid droplets. The combination of the rarity of clusters and their propensity toward fission rather than aggregation makes direct simulation of nucleation difficult under realistic conditions. An enormous number of molecules must be simulated for a very long time to witness a rare event like nucleation under reasonable conditions. Rather than integrating the equations of motion for a metastable phase, Langer developed a method to calculate the nucleation rate by focusing on the rate-limiting step for nucleation.1,2 His method is based on an expansion of the system Hamiltonian near the critical fluctuation. Such an approach is most fruitful when the system is described by a relatively small number of coordinates and conjugate momenta, such as in Ising models of magnetism or order-disorder phenomena in alloys.2 A suitable Hamiltonian for vapor condensation includes all of †
Part of the special issue “Howard Reiss Festschrift”.
the coordinates of the cluster and the molecules that collide with it. This enormous number of coordinates makes Langer’s method unwieldly for analysis of nucleation from supersaturated vapor. In our description of metastable vapors, we described three time scales. The vibration/rotation time scale is by far the shortest. Collisions with the background gas occur on an intermediate time scale. Finally, aggregation/fission events occur on the longest time scale. Instead of considering all of the system coordinates and momenta, which evolve on all three time scales, it is preferable to develop a model that focuses on clusters in the aggregation/fission time scale. In the so-called Classical Theory,3-5 nucleation is modeled as a reversible aggregation one monomer at a time:
C1 + Ci h Ci +1
(1)
It is implicitly assumed that the state of the system at any time following the initial conditions (0 e t < ∞) is fully specified by the concentration of clusters of all sizes (Ci), independent of the history of the concentration distribution. This is appropriate if the clusters are rapidly thermalized through collisions with a background gas. This process should be rapid relative to the aggregation/fission time scale. Langer refers to this limit as overdamping.2 Also, the rates of i-mer aggregation and fission processes are assumed to be proportional to the i-mer concentration. Since the state and transition probabilities of the system are determined by the cluster distribution, the commonly accepted model conforms to the definition of a Markov process.6 Although the precise path through cluster concentration space is random, the asymptotic dynamics for a large system evolving via the set of reactions in eq 1 are fully determined by initial concentrations and by rate constants for aggregation and fission. In practice, the aggregation rate constants are determined from the kinetic theory of gases while the ratio of aggregation and fission rate constants for each reaction is fixed using detailed balancing at saturation conditions. Thus the model of the nonequilibrium process is reduced to a kinetic part and an equilibrium thermodynamic part. Details of the derivation of the nucleation rate expression are given elsewhere.7 The resulting expression is ∞
J ) 1/
a sat i +1 -1 (k1,i C sat ) ∑ 1 Ci S i)1
10.1021/jp012201a CCC: $20.00 © 2001 American Chemical Society Published on Web 11/01/2001
(2)
11894 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Bahadur and McClurg
where J is the nucleation rate expressed as number of droplets formed per unit volume and time, ka1,i is the aggregation rate is the saturation constant for the reaction in eq 1, C sat i concentration of i-mers, and S is the monomer activity relative to saturation. In this paper, we compare a variety of methods for determining the cluster concentrations at saturation (C sat i ) and suggest which of the methods is consistent with the assumed Markov process. We contend that it is this consistency requirement that dictates the proper cluster definition.
Evidently, λ ) λ1 and the partition functions Q1 and Q100‚‚‚ are the same quantity: the partition function of a system containing one monomer. Likewise, the partition functions from eqs 3 and 7 are related by
QN )
Q1 ) Q100‚‚‚
Cluster concentrations at saturation are calculated either using a droplet model based on bulk properties or using atomistic simulations based on intermolecular potentials. Droplet models have been discussed elsewhere.7 It is widely recognized that droplet models are inaccurate when applied to small clusters rather than macroscopic droplets. Atomistic models are potentially more accurate, but they require statistical mechanics to relate intermolecular potentials and cluster concentrations via partition functions. In this paper, we follow the atomistic approach and begin by summarizing several relevant results relating cluster size distributions to the partition functions of individual clusters.8 This serves to define nomenclature used elsewhere and to explain the necessity of a cluster definition. The grand partition function (Ξ) of the vapor is
QN λN ∑ Ng0
(3)
where QN is the canonical partition function for a system containing N monomers and λ is the absolute activity. This formulation implicitly considers clusters since QN is a weighted average over all possible configurations including those in which subsets of the N monomers form clusters. Now, consider an alternative approach in which clusters appear explicitly. Let λi be the absolute activity of a cluster of size i and Ni the total number of such clusters. Then eq 3 becomes
Ξ)
(∏λNi )QN ∑ Ng0 i
(4)
i
(10)
Similarly, the partition function for a two-monomer system (Q2) is divided into contributions from separate monomers (Q200‚‚‚) and one dimer (Q010‚‚‚),
Q2 ) Q200‚‚‚ + Q010‚‚‚
(11)
and the classical partition function for the three-monomer system (Q3) is broken into three terms,
Q3 ) Q300‚‚‚ + Q110‚‚‚ + Q0010‚‚‚
(12)
The purpose of a cluster definition is to unambiguously determine the division of phase space into these contributions. Whereas thermodynamic properties can be calculated using either eq 3 or 7 and are therefore independent of a cluster definition, kinetics are sensitive to such definitions. This is easily shown by computing the equilibrium concentration of i-mers from the grand potential (eq 7)
C sat i )
N hi V
sat
)
( )
λi ∂ ln Ξ V sat ∂λi
(13)
Expressions for individual cluster concentrations can be obtained by using the definition of Ξ from eq 4 and differentiating with respect to the cluster activities as indicated in eq 6. At equilibrium, the activities of the larger clusters can be related to the activity of the monomer through eq 6. These steps are carried out in greater detail by Hill.8 The resulting average concentrations of monomers and dimers are
C sat 1 )
where
(9)
where the sum is over all cluster sets N conforming to eq 5. Therefore, the partition function for a one-monomer system is
II. Statistical Mechanics Background
Ξ)
QN ∑ N
1 2 [λQ100‚‚‚ + λ2(2Q200‚‚‚ - Q100‚‚‚ ) + ‚ ‚ ‚] (14) V sat
N
N)
iNi ∑ i)1
(5)
and QN is the partition function for a set of N monomers in a particular aggregation state N ) {N1, N2, . . . , Ni, . . .}. At equilibrium, the cluster and monomer activities are related by
λi ) λi1
(6)
We now have
Ξ)
λ1ΣiN QN ∑ Ng0 i
(7)
Equations 3 and 7 give the grand potential of the same system at the same conditions and therefore must be identical. The correspondence can be established by taking the limit as λ and λi approach zero
Ξ ) 1 + λQ1 + O(λ2) ) 1 + λ1Q100‚‚‚ + O(λ21, λ2) (8)
C sat 2 )
1 2 [λ Q010‚‚‚ + λ3(Q110‚‚‚ - Q100‚‚‚Q010‚‚‚) + ‚ ‚ ‚] (15) sat V
Equations 14 and 15 demonstrate that cluster concentrations at saturation depend on the division of the partition function for an N-monomer system (QN) into partition functions for various aggregation states (Qijk‚‚‚). The nucleation rate (J) is, in turn, a function of the cluster concentrations as shown in eq 2. Therefore, the calculated nucleation rate depends on the cluster definition used. At this point we should stress that equilibrium thermodynamic properties of the vapor do not depend on cluster definitions as long as each point in phase space is assigned to one of the partition functions. This is demonstrated by Hill for calculation of the second virial coefficient,9 but it is true of any equilibrium property. The nucleation rate is nonzero only for metastable states and is therefore not an equilibrium property. It is tempting to regard the set of cluster concentrations at saturation as an equilibrium property, but equilibrium thermodynamics does not provide a unique way to identify clusters. Therefore the cluster
Cluster Definitions in Nucleation Rate Calculations
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11895
concentrations at saturation in eq 2 cannot be determined from thermodynamics alone, but must be based on a cluster definition that is consistent with the assumed Markov process. There are several clusters definitions in use which we compare in section III. Before comparing various cluster definitions, we summarize the partition functions that we use repeatedly to illustrate the various definitions. The classical partition function for an N-monomer system is
QN )
∑Qijk... ) σh3N ∫ ‚ ‚ ‚ ∫ 1
exp
dq ‚ ‚ ‚ dq ( -H kT ) 1
N
dp1 ‚ ‚ ‚ dpN (16)
where H is the Hamiltonian for the system with generalized positions qi and momenta pi. The integral is performed over all of phase space and therefore includes contributions from all possible combinations of bound states. Using relative coordinates and integrating with respect to the center-of-mass and its conjugate momentum gives
QN )
1 σh
3N-3
( )∫ ∫ ( ) V Λ3M
exp
‚‚‚
-HNint dq1 ‚ ‚ ‚ dqN-1 dp1 ‚ ‚ ‚ dpN-1 (17) kT
where V is the system volume, ΛM is the deBroglie wavelength based on the total mass (M), the internal Hamiltonian is int
H
|pc|2 )H2M
(18)
Partition functions for larger clusters are written similarly in mobile coordinates. III. Cluster Definitions In the following subsections we describe a variety of cluster definitions, provide representative references in which the definitions have been used, and illustrate the application of each for the calculation of partition functions for bound dimers (Q010‚‚‚) and trimers (Q0010‚‚‚). For simplicity in the illustrations, we assume that monomers behave classically, are spherical, and lack internal structure. Generalization to monomers with internal structure (molecular rather than atomic monomers) is straightforward but needlessly complicates the nomenclature. A. Geometric Definitions. The simplest and most widely applied cluster definitions depend only on the system coordinates (q) and not the momenta (p). Stillinger considered a cluster to consist of a group of monomers, each of which lies within a distance rc of at least one other monomer in the group.10 Rao et al. used this definition to study nucleation with molecular dynamics and Monte Carlo techniques.11 Early work by Reiss et al. defined a cluster to be a group of monomers contained within a sphere of radius rc about the center of mass.12 (More recently, he has introduced another cluster definition that is discussed below.) Lee et al. followed this definition in their Monte Carlo calculations.13 We refer to these definitions collectively as geometric cluster definitions. In the geometric cluster definitions, the conjugate momenta of the monomers do not determine if they comprise a bound cluster or not (the only constraint being placed on coordinates), so that the corresponding integrals in eq 17 may be carried out,
Q0...010... )
1 σ
(∏ ) V
j
and pc is the center-of-mass momentum. For the two-monomer system,
Q2 )
( )∫
1 V σh3 Λ3M
‚‚‚
H2int )
∫
( )
-H2int exp dqdp kT
|p|2 + U(r) 2µ12
µ12 )
m1m2 m1+m2
(19)
(20)
( )∫
1 V Q3 ) 6 3 σh ΛM
‚‚‚
∫
( )
-H3int exp dq1dq2dp1dp2 (22) kT
|p1|2 |p2|2 + + U(q1,q2) H3int ) 2µ12 2µ123 µ123 )
µ12m3 µ12 + m3
(23)
(24)
Here, σ is one if the monomers are different, two if any two are indistinguishable, and six if all three are indistinguishable.
Λj
[ ]
-U(q) dq1 ‚ ‚ ‚ dqN-1 kT (25)
where Λj is the deBroglie wavelength based upon the mass of monomer j. Primes on the integrals above indicate that only configurations conforming to the distance cutoffs are to be included. The potential energy depends only upon the separation of the monomers in our illustration, so the orientation integrals in eq 25 can also be carried out. The bound dimer criterion then reduces to |r| < rc, and the partition function (eq 25) becomes
(21)
where r is the distance between monomers and σ is the symmetry number for the cluster, which equals 1 if the monomers are different and 2 if they are indistinguishable. For the three-monomer system (using mobile coordinates),
∫′ ‚ ‚ ‚ ∫′ exp 3
Q010‚‚‚ )
4πV σΛ31 Λ32
∫0r r2 exp c
[ ]
-U(r) dr kT
(26)
The final integral may be evaluated numerically for a given interaction potential. Similarly, the partition function for the bound trimer is
Q0010‚‚‚ )
σΛ31
V Λ32 Λ33
∫′ ‚ ‚ ‚ ∫′ exp
[
]
-U(q) dq1dq2 (27) kT
with primes indicating that the integration is to be carried out over all configurations with at least two of the separations less than the cutoff value. The remaining integrations may be evaluated using Monte Carlo integration.14 Extension to larger clusters is straightforward. B. Energetic Definition. Hill8 defined a cluster as a collection of monomers that are each energetically bound to at least one other monomer. He considered a pair of monomers to be energetically bound if the sum of their interaction energy and their relative kinetic energy is negative. Barrett15 used this
11896 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Bahadur and McClurg
definition in calculating equilibrium cluster populations. Alternatively, one could construct a cluster definition in which the Hamiltonian of a cluster is negative.16 These definitions coincide for the dimer but differ for larger clusters. We refer to these definitions as energetic cluster definitions. We illustrate energetic cluster definitions using Hill’s definition and spherical monomers. Therefore, we begin by focusing on a dimer. As long as the potential energy is negative [U(q) < 0], there is a limited range of momenta such that the internal kinetic energy is less than the negative of the potential energy [K(p,q) < -U(q)] and the internal Hamiltonian is also negative
H
int
) K(p,q) + U(q) < 0
(28)
The dimer partition function may be calculated by integrating over all configurations with negative interaction potentials and limiting the range of momentum integration such that the Hamiltonian is also negative.
Q010‚‚‚ )
4πV σΛ31 Λ32
∫U(r)e0r2F(r) exp
[ ]
-U(r) dr kT
Γ[3/2] - Γ[3/2, -U(r)/kT] Γ[3/2]
int
) K(p,q) + U eff(q,L)
int
) K(p,q) + U eff(q,L) < H
int max(L)
(32)
Therefore, the partition function may be calculated by integrating over all allowed values of the net angular momentum and limiting the configurations and momenta such that the Hamiltonian is less than its maximum value,
(30)
Q0...010... )
and Γ[n, x] is the incomplete Γ-function. Details of the integrations leading up to eqs 29 and 30 are given elsewhere.8 The remaining integral in eq 29 may be evaluated numerically for a given interaction potential. For systems larger than the dimer, Monte Carlo integration may be used to evaluate the partition function.14 Positions and momenta are chosen for each monomer and are accepted if all of the monomers are bonded to at least one other monomer according to eq 28 and rejected otherwise. Dahler has shown that calculation of the trimer partition function under the alternate energetic criteria, H int < 0, reduces to a twodimensional quadrature.16 C. Dynamic Definition. In the introduction, we described clusters as a group of molecules that are in close proximity and remain so for more than a few vibrational periods. This intuitive description contains both proximity and stability components. Although geometric definitions require proximity, they do not ensure stability. Energetic definitions are an attempt to ensure stability, but they are overly restrictive. There are classically bound states with sufficient energy to dissociate but cannot due to conservation of angular momentum. These are known as resonant states. We call the cluster definition that includes all classically bound states the dynamic cluster definition because it includes all states whose internal dynamics are bounded. Lushnikov and Kulmala have used a dynamic cluster definition in their analysis of nucleation when dimerization is the ratelimiting step.17 To include the contribution of resonant states in the cluster partition function, we express the internal Hamiltonian as a function of the relative coordinates (q), the net angular momentum (L), and the remaining momenta (p). Those terms that do not depend on nonconserved momenta can be considered as part of an effective potential (U eff), and the internal Hamiltonian becomes
H
H
(29)
where
F(r) )
a dissociated state. Systems can tunnel through the barrier, but at a rate that decreases exponentially with the square root of the mass.18 Tunneling is probably negligible except for clusters of quantum gases (e.g., H2 or He), and we therefore neglect tunneling in the subsequent discussion. To evaluate partition functions including resonant states, we note that for sufficiently large angular momentum (L), there is no local minimum on the effective potential. This places an upper limit on the magnitude of the cluster angular momentum (Lmax). As long as the angular momentum is less than the maximum value (L < Lmax), there is an upper limit to the energy for the cluster to remain bound [H int max(L)]. If the effective potential is less than this value, there is a limited range of momenta such that the internal Hamiltonian is less than its maximum value
(31)
If there is a local minimum in the effective potential, a cluster can remain classically bound even if it has the same energy as
1 σ
(∏ ) V 3 j Λj
∫0L ∫′ ‚ ‚ ‚ ∫′F(q,L)
exp
max
[
]
-U eff(q,L) dq1 ‚ ‚ ‚ dqN-1 dL (33) kT
where primes indicate that the integration is over configurations with U eff(q,L) < H int max and
F(q,L) )
1 h3N-3
( )
∏ j Λ3j ∫′ ‚ ‚ ‚ ∫′ Λ3M
exp
[ ]
-K(p,q) dp1 ‚ ‚ ‚ dpN-1 (34) kT
where primes indicate that the integration is over momenta such eff that K(p,q) < [H int max - U (q,L)] and the magnitude of the angular momentum is constrained to L. Again, coefficients in eq 34 are chosen such that F(q,L) is nondimensional. As with the previous definitions, we illustrate the dynamic definition for calculating the partition functions of bound dimers and trimers of spherical monomers. The effective potential and Hamiltonian for two bodies interacting with central forces is discussed in standard texts.19 The effective potential is
U 2eff ) U(r) +
L2 2µr2
(35)
For sufficiently large angular momentum (L), the effective potential does not have a local minimum and there are no bound states. Below Lmax, there is a local minimum and three classes of states. First are high-energy states that dissociate in a single vibrational period. Second are low-energy states with large separations that are outside of the potential well and have insufficient energy to enter. Finally, there are low energy states that are inside of the potential well. The third class is the only one that remains a cluster on the aggregation/fission time scale and therefore should contribute to the cluster partition function. Thus the corresponding partition function is
Cluster Definitions in Nucleation Rate Calculations
Q010‚‚‚ )
( )∫
16π2 V σh3 Λ3M
∫r
Lmax
0
rmax(L)
min(L)
∫
pmax r (r, L) L 0
exp
(
)
-H int dpr dr dL (36) kT
where the limits on the separation (r) and the radial momentum (pr) result from the upper limit on the Hamiltonian to remain in the potential well. For clusters larger than the dimer, there are saddles on the energy surface for L < Lmax. In the range of energy between the lowest saddle and the highest turning point, the cluster’s dynamic stability is uncertain. This range corresponds to clusters with sufficient energy to dissociate given their angular momentum, but it is still plausible that a given cluster remains in a bounded trajectory never passing through the saddle(s). Examples of such bounded trajectories are the KAM tori.20,21 We have chosen not to include such states for two reasons, one pragmatic and one physical. First, there is no simple criterion to determine if a point in phase space belongs to a bounded trajectory in this range. Apart from the conservation laws used here, we know of no way to determine if a given point in phase space belongs to a bounded trajectory short of integrating the equations of motion. The dynamics of three or more bodies is known to exhibit chaotic behavior so that within this range of energies, even nearby points in phase space may have very different behavior. Fortunately we expect that the bounded trajectories in this range of energies are a very special subset that happen to rebound off of the “hills” on the potential surface without ever passing over a saddle. Recall that the classical nucleation theory implicitly assumes that clusters are rapidly thermalized. For vapor condensation, this occurs through frequent collisions with a background gas. We expect that clusters with sufficient energy to dissociate given their angular momentum but residing in one of these special bound trajectories are likely to experience a collision that shifts the dynamics into a nearby unbound state. Therefore, our second reason for excluding points in this energy range is that they are unlikely to remain bounded on the aggregation/fission time scale. Neglecting the bound trajectories in the ambiguous energy range results in a rigorous lower bound to the partition function, but we expect it to be a close lower bound. Application of the dynamic cluster definition to clusters larger than the dimer is not trivial, even when the ambiguous energy range is neglected. The difficulty lies in locating all of the energy surface saddles and determining if a given point in phase space is inside or outside of the energy well. The required steps are presented for the trimer in the Appendix. The Hamiltonian is given in eq 42 and the effective potential is given in eq 47. The saddle points on the effective potential surface were located by numerically following the local minima of the effective potential as a function of angular momentum until the local minima vanish at saddle points. The lowest saddle point was identified by comparing the effective potential values at each of the saddles. For our model potential, the lowest saddle corresponded to a configuration in which one of the monomers is on the verge of leaving the cluster along a line perpendicular to the remaining dimer and the angular momentum vector is parallel to the dimer. Obviously there are equivalent saddles with different choices for the departing monomer. The partition function [eq 33] is then integrated over 0 e L e Llsp subject to the limitations in eqs 46 and 48. Also, the momenta are limited such that the internal Hamiltonian is less than the value at the saddle. These integrations were carried out using Monte Carlo
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11897 sampling,14 rejecting points that either had too much energy to remain bound or are outside of the potential well. We anticipate that extension to larger clusters is progressively more challenging because there are increasingly more saddles to identify and because identifying which phase points are inside of the well is more difficult in a higher dimensional space. D. Modified Liquid Drop. The modified liquid drop model proposed by Reiss and co-workers seems to overcome the limitation of arbitrarily selecting a cutoff distance (rc) in the geometric definition by using two identifying parameters for a cluster. In addition to the cluster size (i), the cluster volume (V) is considered.22 They suggest that the work required to create a cluster containing i monomers in a volume V has a saddle point as a function of its two variables. This saddle point corresponds to the cluster in their model. Since the work of formation depends on momenta as well as positions, this is not a purely geometric definition. The saddle point can be located by a series of simulations of various system sizes as a function of volume. Note that this saddle point is not the same as the saddle point discussed in the context of the dynamic cluster definition. In the dynamic definition, the saddle point is on the effective potential surface for a single cluster while this saddle is on the size-volume surface for a series of clusters sizes. Thus the specification of a bound cluster depends on the free energy of clusters of adjacent sizes in this definition. The dependence of the definition of one cluster on the thermodynamics of other clusters makes if difficult to compare this model with other cluster definitions. IV. Sample Calculations We illustrate the discrepancies between the partition functions described in section III using a simple model for xenon. We chose xenon because it has an orientation independent potential and the heavy nuclei minimize quantum effects, justifying the use of classical partition functions. The xenon interaction potential was modeled using the three-parameter Kihara potential23
U(r) ) 4
[(σr -- 2a2a) - (σr -- 2a2a) ] 12
6
(37)
with ) 2.003 kJ/mol, σ ) 3.963 Å, and a ) 0.190 Å. The parameters for the potential were chosen to match the bulk crystal cohesive energy and lattice constant24 and the dimer vibration frequency.25 Figure 1 presents partition functions for monomers, dimers, and trimers using the geometric, energetic, and dynamic cluster definitions. Quantum partition functions, in the rigid rotor/ harmonic oscillator approximation, are also shown for comparison.9,26 The temperature range covers the triple point to the critical point of xenon.27 A cutoff distance of three times the location of the potential minimum was used for the geometric definition, as suggested by Stillinger.10 Numerical integrations over complicated domains were performed using a simple Monte Carlo method, which is not to be confused with Metropolis sampling Monte Carlo integration.14 The integrations were converged to at least (30%. The various dimer partition functions span an order of magnitude while the trimer partition functions span several orders of magnitude. The effect of the partition function discrepancies on saturation cluster concentrations is illustrated in Figure 2. Since the small but nonzero concentrations of large clusters at saturation influence the nucleation rate in eq 2, the wide discrepancies between the trimer concentration estimates are significant. We
11898 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Bahadur and McClurg
Figure 3. Phase space included by different cluster definitions for a typical dimer interaction potential. The solid line is the effective potential with a moderate total angular momentum, while the dashed line is the bare interaction potential. Geometric definitions include energy states up to an arbitrary cutoff (g). The energetic definition includes only negative energy states below e-e′. The dynamic definition includes energy states inside the effective potential below d-d′. For high-dimensional energy surfaces, d corresponds to the lowest saddle.
Figure 1. Partition functions for bound xenon clusters for the geometric (dashed), energetic (dot-dashed), and dynamic (solid) cluster definitions. The rigid rotor/harmonic oscillator approximation is given for comparison (dotted). The partition functions coincide for one monomer, differ by an order of magnitude for bound dimers, and differ by several orders of magnitude for bound trimers.
Figure 2. Saturation cluster concentrations for xenon. The key is as in Figure 1. The cluster size labels are each above and to the left of the corresponding curves. Discrepancies between the definitions mirror those in the partition functions.
expect that the model discrepancies continue to grow with increasing cluster size. V. Discussion Figure 3 illustrates the domain of integration for the dimer in the various cluster definitions. The dotted line is the interaction potential, the solid line is the effective potential for a
moderate net angular momentum, and the distance between them is the rotational kinetic energy. The distance above the solid line indicates kinetic energy associated with stretching. In the darkly shaded region, below the line e′-e, the Hamiltonian is negative. This is the region included in the energetic definition. In the moderately shaded region, below d′-d, the Hamiltonian is positive but insufficient to escape the effective potential well. Clusters in this region are classically bound and are properly included in the dynamic definition. Finally, the lightly shaded region out to the arbitrary cutoff distance g is included in the geometric definition. The corresponding diagram for larger clusters is similar except that the domains are high-dimensional, making them impossible to draw and d corresponds to the lowest saddle point on the effective potential surface rather than a local maximum. As discussed in section I, classical nucleation rate calculations are based on a step-growth mechanism and a Markov process analysis. This analysis focuses on the aggregation/fission timescale and assumes that the internal cluster dynamics are thermalized on a much shorter time scale. An important assumption in the analysis is that the transition probability from one cluster concentration distribution to another is independent of the system history. This assumption dictates the appropriate cluster definition. Consider a high-energy state resulting from a collision. Such unbound configurations must not be included as clusters since they will certainly dissociate on the rotation/vibration time scale and thus the probability of growth is negligible. Such configurations are merely in the process of fission. On the other hand, all configurations with a significant probability of growth must be counted as clusters. The dynamic definition best satisfies this criterion since all classically bound clusters are included but those with sufficient energy to dissociate are excluded. Geometric definitions suffer from arbitrary cutoff distances. Since the proper cutoff distance changes with net angular momentum and the relative importance of high angular momentum states changes with temperature, a single cutoff distance does not adequately define clusters. As a consequence of this effect, molecular simulations of low angular momentum clusters with vibration energy approaching the dissociation limit would alternately be counted as dimers or as a pair of monomers as they vibrate. Also, geometric definitions do not limit the momenta. Thus transient states associated with collisions are (inappropriately) included as clusters in this definition. The net effect of these effects is that geometric definitions underestimate partition functions of clusters at low temperatures and overes-
Cluster Definitions in Nucleation Rate Calculations timate them at high temperatures. For the xenon model described above, the geometric model came closest to the dynamic model. This indicates that the cutoff chosen by Stillinger is a good tradeoff between the two effects, at least for the chosen system, cluster sizes, and conditions. Note, however, that the temperature dependence of the partition function is incorrect. Energetic definitions undercount clusters by excluding resonant states. Unless tunneling of monomers through the effective potential barrier is significant, resonant states should be counted as bound clusters. Therefore, partition functions calculated using the energetic definition based on the cluster Hamiltonian (H < 0) provide a rigorous lower bound to those using the dynamic definition. Hill’s energetic definition is also prone to underestimating the partition function through neglect of resonant states, but by using a pairwise criterion it does include some states that are not bound. Thus Hill’s definition may not provide a rigorous lower bound. Quantum partition functions in the rigid rotor/harmonic oscillator approximation are included in Figure 1 for comparison with the cluster definitions. Since the harmonic oscillator approximation is based on a quadratic expansion to the potential surface, clusters cannot undergo fission in this model. Therefore, the rigid rotor/harmonic oscillator approximation is not a cluster definition but is used in lieu of a definition. Instead of providing criteria to distinguish a cluster, this approximation leads to a partition function estimate for a collection of monomers assumed to comprise a cluster. The rigid rotor approximation is quite accurate except for clusters of quantum gases. The more significant approximation is the quadratic expansion to the potential surface, which limits the accuracy at elevated temperatures. Therefore, we recommend the rigid rotor/harmonic oscillator approximation at temperatures well below the critical point. Nucleation calculations using eq 2 require cluster concentrations at saturation for clusters much larger than trimers. Typically, clusters containing one to several dozen monomers are required at reasonable conditions. However, by comparison of the partition functions for monomers, dimers, and trimers, we have illustrated that there are significant discrepancies between the various classes of cluster definitions. The cluster definition most likely to yield nucleation rates consistent with experiment is the one that is most consistent with the assumptions of the assumed Markov process. In particular, the transition probabilities between clusters are assumed to be related to the cluster concentrations. Definitions that either include states that are in the process of fission or exclude states that have the potential to aggregate are not consistent with this assumption. The definition adhering to the assumptions inherent in a Markov process is the dynamic definition. Unfortunately, this definition is the most difficult to apply. The geometric definition is the simplest to apply, but its accuracy depends on a fortuitous balance between tendencies to overcount and undercount clusters. The energetic definition is also relatively simple and can provide a rigorous lower bound to the partition function using the dynamic definition. Therefore, we recommend the rigid rotor/harmonic oscillator approximation at low temperature and the dynamic definition at temperatures approaching the critical temperature. Further method development for application of the dynamic definition to larger clusters is needed. Acknowledgment. Our thanks go to Professor John Dahler who provided a manuscript in preparation (ref 16). We gratefully acknowledge support for this research from the Shell Oil Foundation and from the Graduate School of the University of Minnesota.
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11899 Appendix There are several steps in applying the dynamic definition to a cluster of a given size. First the Hamiltonian must be expressed with the net angular momentum given explicitly. The form of the internal Hamiltonian used for the trimer in the dynamic cluster definition (section IIIC) was derived starting with the reduced form in section 158 of ref 28. To eliminate most of the cross terms in the momenta, we applied two additional contact transformations. The first, defined by
∂W ∂W , pi ) ∂Pi ∂qi
Qi )
(38)
W ) P1(q1 - q3) + P2(q2 - q4) + P3(m1q1 + m2q3) P4(m1q2 + m2q4) + (39) m1 + m2 m1 + m2 transforms to Cartesian mobile coordinates in the invariant plane. The second contact transformation, defined by
qi )
∂W ∂W , Pi ) ∂pi ∂Qi
(40)
W ) p1Q1 sin(Q3) + p2Q1 cos(Q3) + p3Q2 sin(Q4) + p4Q2 cos(Q4) (41) further transforms to radial mobile coordinates. The resulting reduced Hamiltonian is
H
int
p12 p22 p32 p42 ) + + + + 2µ12 2µ123 2µ q 2 2µ q 2 12 1 123 2
[
[L2 - (p3 + p4)2] cos(q4)2 sin(q3 - q4)2
2µ12q12
+
cos(q3)2 2µ123q22
]
+ U(r12, r13, r23) (42)
where we have used the reduced masses defined in eqs 21 and 24. The separation distances between monomers are
r13 )
r23 )
x x
m22q12
(m1 + m2)
2
m12q12
(m1 + m2)
2
r12 ) q1
(43)
+
2m2cos(q3 - q4)q1q2 + q22 m 1 + m2
(44)
-
2m1 cos(q3 - q4)q1q2 + q22 m 1 + m2
(45)
The definitions of the generalized momenta and positions in eq 42 require that they are limited to the following ranges:
|p3 + p4| e L {q1, q2} g 0
(46)
0 e {q3, q4} < 2π In eq 42, p1 and p2 are stretching momenta, p3 - p4 is the bending momentum, L is the magnitude of the total angular momentum, L⊥ ) p3 + p4 is the angular momentum perpendicular to the plane of the monomers, and L| ) L - (p3 + p4) is the angular momentum parallel to that plane. Note that L is
11900 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Bahadur and McClurg
conserved but L⊥ and L| are not. Terms that do not depend on the stretching and bending momenta in eq 42 belong to the effective potential. If L⊥ ) 0, which is the case for the lowest saddle point using our potential, the effective potential is
U eff 3 )
[
]
cos(q4)2 cos(q3)2 L2 + + U(r12, r13, r23) sin(q3 - q4)2 2µ12q12 2µ123q22 (47)
There is a local minimum on this effective potential surface for total angular momenta up to a maximum value (Llsp) corresponding to the lowest saddle point. The effective potential at the saddle sets an upper limit on the Hamiltonian, ensuring that the cluster has insufficient energy to escape from the potential well. In addition to the energetic limit, geometric criteria are needed to ensure that the cluster is inside of the potential well. The border between bound and unbound states is fixed by the lowest saddle on the effective potential surface for a given angular momentum (L). At this saddle point, the largest separation between monomers is Rmax(L) and the internal energy is H int max(L). Then the following conditions
H
