Kinetics of Multicomponent Nanosize Clusters on Solid Surfaces

7130

Langmuir 2003, 19, 7130-7140

Kinetics of Multicomponent Nanosize Clusters on Solid Surfaces David N. Brunelli and Rex T. Skodje* Department of Chemistry and Biochemistry, University of Colorado, Boulder, Colorado 80309 Received March 13, 2003. In Final Form: May 6, 2003 We discuss the kinetics of multicomponent liquidlike clusters adsorbed onto solid surfaces. In this work, the rate expressions governing the growth and decay of surface clusters coupled to a two-dimensional gas phase are derived based on previous treatments of unary systems. General expressions are presented for the equilibrium partial pressures, the monomer evaporation rates, and the monomer capture rates. The composition dependence of the rate expressions is accomplished using the theory of real solutions to represent the activity coefficients and the surface tension. Numerical results from dynamical Monte Carlo simulations are presented for a nontrivial model of two-dimensional binary clusters. A generalized Gibbs-Thomson expression for the partial pressures is found to be in excellent agreement with the numerical results even for very small cluster sizes. The species selective evaporation rates are found to obey a simple rate law that predicts strong deviations from a naı¨ve geometric scaling behavior. The growth and decay rates of binary clusters in the simulations were accurately modeled using the analytical rate expressions.

I. Introductions The kinetics of small atomic clusters adsorbed on solid substrates has attracted considerable attention1,2 because of its relation to important phenomena such as thin film growth. Experiment3-9 has demonstrated that clusters consisting of hundreds to thousands of monomers can exhibit collective diffusion, growth/decay through capture/ emission of atoms into a two-dimensional gas phase, and coalescence into larger clusters upon contact. The use of scanning tunneling microscopy (STM), atomic force microscopy (AFM), and other direct imaging methods has permitted the direct observation of the distribution function of cluster size and position as functions of real time in a variety of systems. The analysis of such detailed experimental data requires theoretical models capable of describing the effects of the underlying atomic level processes that govern the evolution of the adlayer.10-29 (1) Zinke-Allmang, M. Scanning Microsc. 1990, 4, 523. (2) Jenson, P. Rev. Mod. Phys. 1999, 71, 1695. (3) Komura, S.; Osamura, K.; Fujii, H.; Tanaka, T. Phys. Rev. B 1985, 31, 1278. (4) Kruger, D. W.; Savage, D. E.; Lagally, M. G. Phys. Rev. Lett. 1989, 63, 402. (5) Wen, J.-M.; Evans, J. W.; Bartelt, M. C.; Burnett, J. W.; Thiel, P. A. Phys. Rev. Lett. 1996, 76, 652. (6) Morgenstern, K.; Rosenfeld, G.; Comsa, G. Phys. Rev. Lett. 1996, 76, 2115. (7) Pai, W. W.; Swan, A. K.; Zhang, Z.; Wendelken, J. F. Phys. Rev. Lett. 1997, 79, 3210. (8) Morgenstern, K.; Rosenfeld, G.; Poelsema, B.; Comsa, G. Phys. Rev. Lett. 1995, 74, 2058. (9) Heinonen, J.; Koponen, I.; Merikoski, J.; Ala-Nissila, T. Phys. Rev. Lett. 1999, 82, 2733. (10) Chakraverty, B. K. J. Phys. Chem. Solids 1967, 28, 2401. (11) Villain, J. Europhys. Lett. 1986, 2, 531. (12) Voter, A. F. Modeling of Optical Thin Films: 20-21 August 1987, San Diego, California; Proceedings of SPIE, Vol. 821; International Society for Optical Engineering: Bellingham, WA, 1988; p 214. (13) Villarba, M.; Jonsson, H. Surf. Sci. 1994, 317, 15. (14) Sholl, D. S.; Skodje, R. T. Phys. Rev. Lett. 1995, 75, 3158. (15) Schulze Icking-Konert, G.; Giesen, M.; Ibach, H. Surf. Sci. 1998, 398, 37. (16) McLean, J.; Krishnamachari, B.; Peale, D.; Chason, E.; Sentha, J.; Cooper, B. Phys. Rev. B 1997, 55, 1811. (17) DeW. Van Siclen, C. Phys. Rev. Lett. 1995, 75, 1574. (18) Krishnamachari, B.; McLean, J.; Cooper, B.; Sentha, J. Phys. Rev. B 1996, 54, 8899. (19) Khare, S. V.; Bartelt, N. C.; Einstein, T. L. Phys. Rev. Lett. 1995, 75, 2148.

Given a reasonable potential energy surface, a brute force simulation of the atomic configuration of a representative sample of clusters can be carried out using standard molecular dynamics (MD) or Monte Carlo (MC) techniques. While much can be learned from such simulations, real experiments often take place on time scales far too long to maintain an atomistic description of the system. On the other hand, a purely thermodynamic analysis,30 such as provided by classical nucleation theory, is capable of reproducing certain characteristics of large clusters as they approach the bulk limit at extremely long times. Unfortunately, thermodynamics does not provide fundamental rates necessary for a kinetic model nor does it provide a particularly accurate description of the smaller clusters in the distribution. Therefore, intermediate level models are sometimes employed in the study of the longtime-scale evolution of adlayers. We note in particular the use of a quasichemical approximation sometimes referred to as the “rate equation description”. In the rate equation approach,31,32 the atomic level description of the adlayer is abandoned in favor of a set of concentrations that represent the number of clusters of various sizes. The cluster concentrations change according to standard chemical rate equations with the values of rate coefficients set by the atomic level kinetic processes such as evaporation and capture of monomers. The use of such a nonatomistic description of the adlayer kinetics can greatly extend the range of time scales accessible to simulation but requires the input values of the rate constants over a large range of cluster sizes and (20) Soler, J. M. Phys. Rev. Lett. 1997, 79, 4238. (21) Bartelt, N. C.; Theis, W.; Tromp, R. M. Phys. Rev. B 1996, 54, 11741. (22) Metiu, H.; Rosenfeld, G. Surf. Sci. Lett. 1997, 373, L357. (23) Shao, H.; Weakliem, P.; Metiu, H. Phys. Rev. B 1996, 53, 16041. (24) Evans, J. W.; Bartelt, M. C. Phys. Rev. B 2001, 63, 235408. (25) Sholl, D. S.; Skodje, R. T. Physica A 1996, 231, 631. (26) Sholl, D. S.; Fichthorn, K. A.; Skodje, R. T. J. Vac. Sci. Technol. 1997, 14, 1275. (27) Lo, A.; Skodje, R. T. J. Chem. Phys. 1999, 111, 2726. (28) Lo, A.; Skodje, R. T. J. Chem. Phys. 2000, 112, 1966. (29) Thiel, P. A.; Evans, J. W. J. Phys. Chem. B 2000, 104, 1663. (30) Jain, S. C.; Hughes, A. E. J. Mater. Sci. 1978, 13, 1611. (31) Venables, A. J. Philos. Mag. 1973, 27, 697. (32) Biham, O.; Barkema, G. T.; Breeman, M. Surf. Sci. 1995, 324, 47.

10.1021/la034435m CCC: $25.00 © 2003 American Chemical Society Published on Web 07/01/2003

Kinetics of Multicomponent Nanosize Clusters

Figure 1. A typical liquidlike binary cluster in a dilute twodimensional gas environment. The rate of evaporation, Re, and the rate of capture, Rc, are defined through the change in connectivity of monomers with the main part of the cluster.

conditions. If the intracluster configuration undergoes rapid relaxation, the rate expressions can be accurately parametrized using just the cluster size, N, giving the resulting rate equations a particularly simple form for single-component systems. When a potential energy expression is available, a sensible way to determine the rate constants is from MC simulation on a shorter time scale.28 Otherwise, it is quite useful to have analytical expressions for the rate constants that can be adapted to a specific system through experimental data or otherwise. The objective of the present work is to develop rate expressions that govern the growth of clusters consisting of more than one chemically distinct species. Most work on the kinetics of surface clusters has focused on the behavior of single-component adlayers. However, the characteristics of multicomponent clusters are quite important for the study of alloyed surfaces, doped semiconductors, and composite nanomaterials as well as for the analysis of surface reactions. In this report, we shall present a treatment of binary atomic clusters which focuses on two of the primary rates that control the cluster evolution: the evaporation of monomers from the cluster surface and the reverse process of the capture of monomers from the two-dimensional gas phase. We shall treat the case of two-dimensional surface clusters, see Figure 1, where the evaporation/capture processes can be defined through the connectivity of a monomer with the cluster. While the kinetics of multicomponent surface clusters has received relatively little attention, there exists a considerable body of work devoted to the study of mixed clusters in the wider context of classical nucleation theory (CNT).33-37 The behavior of atmospheric aerosol particles has provided a particularly important application of binary nucleation theory.38,39 On the basis of the classical Kevin (33) Reiss, H. J. Chem. Phys. 1950, 18, 840. (34) Doyle, G. J. J. Chem. Phys. 1961, 35, 795. (35) Wilemski, G. J. Chem. Phys. 1984, 80, 1370. (36) Laaksonen, A.; Talanquer, V.; Oxtoby, D. W. Annu. Rev. Phys. Chem. 1995, 46, 489. (37) Laaksonen, A.; Kulmala, M.; Wagner, P. E. J. Chem. Phys. 1993, 99, 6832.

Langmuir, Vol. 19, No. 17, 2003 7131

equation, the cluster free energy of CNT depends on the curvature of the interface through a term of the form ∼κΓ, where κ is the curvature and Γ is a surface tension. The competition between the surface term and the bulk term leads to the familiar free energy maximum at the critical size for nucleation. A self-consistent formulation of the nucleation problem has been accomplished within a purely thermodynamic framework, although it has been noted that there can be unphysical predictions such as negative surface concentrations. This has led to considerable discussion of the validity of the Kelvin equation for small clusters and of the proper interpretation of the surface tension function for multicomponent systems. In the numerical simulations presented below, we find that the Kelvin equation (as manifested through the GibbsThomson formula for vapor pressure) is remarkably accurate even for clusters smaller than 100 monomers. Instead, it is found that the most serious errors of bulk predictions occur in the kinetic evaporation rates that do not become proportional to surface area until the cluster is vastly larger than the critical nucleus. The physical characteristics of nanosized clusters depend critically on their temperature. A cluster can be solidlike or liquidlike as indicated by some reasonable measure of atomic mobility within the interior of the cluster. Unlike bulk systems, the phase transition between these states does not occur at a precise temperature because of the finite size,40 but clear cases can be easily identified. In this work, we shall focus on model systems that lie unambiguously in the liquidlike regime. For binary clusters, there is the additional issue of the miscibility of the two species. As with the liquid-solid phase transition, there can be size-dependent effects in the transition to miscibility. For simplicity, we shall consider only the temperature regime where the two liquids are miscible at all relevant cluster sizes. We shall not, however, assume that the resulting solution can be represented as ideal or regular. For nonideal solutions, the surface concentrations of species will generally be different from that of the bulk. This may lead to unusual behavior of the evaporation rates and other kinetic parameters as functions of size and composition. To provide a general scheme to treat such systems, we introduce a procedure to model the rates that characterize cluster growth using a general expansion of the activity coefficients in the variables representing mole fractions. The resulting expressions can then be fit to experimental data or simulation data to obtain the cluster kinetics for all sizes and compositions. The remainder of this work is organized as follows. In section II, the standard treatment of the equilibrium vapor pressure of a unary cluster is reviewed for two-dimensional surface clusters. The Gibbs-Thomson formula is generalized to the case of binary clusters. An expansion procedure based on the theory of nonideal solutions is employed to parametrize the resulting expressions in the mole fraction variable that describes the composition of the cluster. In section III, rate expressions for the capture and emission of monomers into the two-dimensional gas phase are presented for binary clusters. In section IV, a numerical simulation is presented for the vapor pressure of binary clusters. The model is based on a dynamical Monte Carlo simulation of two-dimensional liquidlike binary clusters with an interaction potential provided by a Lennard-Jones potential. The analytic expressions derived in section II are tested against the numerical simulations and found (38) Raes, F.; Saltelli, A.; van Dingenen, R. J. Aerosol Sci. 1992, 23, 337. (39) Kathmann, S. M.; Hale, B. N. J. Phys. Chem. B 2001, 105, 11719. (40) Berry, R. S. Nature 1998, 393, 212.

7132

Langmuir, Vol. 19, No. 17, 2003

Brunelli and Skodje

to be quite accurate. In section V, the evaporation and capture rates of binary clusters are simulated using the same model presented in section IV. The analytical rate expressions are tested and found to be reliable. A brief conclusion is presented in section VI. II. Vapor Pressure of Clusters: The Gibbs-Thomson Equation The equilibrium vapor pressure surrounding a cluster will increase as the cluster size decreases. The physical origin of this behavior can be traced to the increase in surface free energy of a curved interface that is quantified using the Kelvin equation. To review the standard thermodynamic treatment,41 the free energy of formation of a cluster consisting of N monomers at constant temperature is given by

∆G )

∫0Nµcl dN + ∫N0µg dN

(2.1)

where µcl and µg are the chemical potentials of the cluster and gas phases, respectively. If we specialize to the case of two-dimensional circular clusters, µcl can be written using the Kelvin equation as

µcl ) µ∞cl +

ΓΩ r

(2.2)

where µ∞cl is the bulk potential for a flat interface, Γ is the line (surface) tension (free energy per unit length), Ω is the surface area of one monomer, and r is the cluster radius which satisfies

Λ t πr2 ) NΩ

(2.3)

defining the area of the cluster Λ. We consider the model with constant line tension. If the two-dimensional gas surrounding the cluster is assumed to behave ideally, that is, pA ) ngkT, then µg can be written as

µg )

µ0g

+ kT ln p

(2.4)

Figure 2. A schematic diagram of the free energy of formation, ∆G, for the formation of a cluster of size N from a twodimensional ideal gas. The upper panel shows the formation of a cluster holding the total number of monomers, n, fixed. There are two stationary points: the unstable critical nucleation size, Nc, and a larger stable equilibrium at Ns. In the lower panel, the free energy is shown for the more usual constraint of constant pressure. This yields only a single stationary point at the nucleation threshold.

solving the relation µcl ) µg for the pressure to obtain the familiar Gibbs-Thomson relation

ΓΩ (rkT )

p(r) ) p∞ exp

where p∞ is the bulk vapor pressure of a flat interface,

(

If the cluster is formed under the conditions of constant (T,p), then the free energy change given by eq 2.1 can be directly integrated to yield

∆G(T,p) ) Nµ∞cl + ΓL - Nµ0g - NkT ln p

(2.5)

≡ Nφ + ΓL where L is the cluster perimeter length, L ) 2πr. Under the conditions appropriate to condensation, that is, φ < 0, the free energy is seen to exhibit a single extremum (a maximum) at the critical size for nucleation

Nc )

πΩ2Γ2 φ2

ΓΩ rc ) φ

)

µ∞cl - µ0g p∞ ) exp kT

(2.9)

Alternatively, the cluster vapor pressure can be defined at a stable equilibrium18 point by constraining the system differently. In particular, if we form the cluster holding constant the temperature and total number of monomers for a fixed area, then gas pressure decreases as the cluster grows forcing a stable equilibrium at a position past the critical size. Hence, we hold the total number of monomers, n ) N + ng, constant when evaluating the free energy change in eq 2.1 yielding

(2.6)

or, equivalently, at the critical radius

(2.8)

∆G ) Nµ∞cl + ΓL - Nµ0g - kT

∫0Nln p dN

(2.10)

where

(2.7)

at which point the two chemical potentials balance. The critical nucleation size is depicted in the schematic for constant pressure presented in Figure 2. The vapor pressure at the unstable equilibrium can be found by (41) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon Press: Oxford, 1982.

∫0Nln p dN ) ∫0Nln(

)

(n - N)kT dN A - NΩ

(2.11)

In addition to the unstable equilibrium at the critical nucleation size, eq 2.5 yields a stable solution at Ns solving

( )

(n - Ns)kT ΓΩ ) p∞ exp A - NsΩ kTrs

(2.12)


Langmuir, Vol. 19, No. 17, 2003 7133

with NsΩ ) πrs2. The two equilibrium sizes are depicted schematically in the upper panel of Figure 2. The equilibrium vapor pressure at this size is of course again given by the Gibbs-Thomson eq 2.6. Numerical simulations carried out near this stable equilibrium offer a practical scheme to empirically determine the vapor pressure for a model system. We are interested in the Gibbs-Thomson formalism for multicomponent systems, specifically for binary clusters containing numbers of species A and B, denoted by NA and NB, where NA + NB ) N. Thus, we wish to derive an expression for the equilibrium partial vapor pressures pA and pB for a cluster of total size N with mole fractions xA ) NA/N and xB ) NB/N of the two species. It was recognized long ago in the context of CNT that the size and composition of a critical cluster could be obtained from balancing the chemical potential for each species present. Thus, instead of eqs 2.2 and 2.4 for the singlecomponent system we have

µicl ) µ∞,i cl +

ΓΩi r

i ) A,B

µig ) µ0,i g + kT ln pi

i ) A,B

(2.13) (2.14)

Equating the chemical potential of the two phases for each species, we obtain a formal expression for the partial pressures,

( )

pi(r) ) p∞i exp

(

p∞i ) exp

kT

Γ ) const

)

µ0,i g

(2.16)

ideal

(2.17) (2.18)

ideal

(2.19)

which yields

( )

0 xi exp pi(r) ) p∞,i

ΓΩi kTr

∂ ln aA ∂ ln aB xA ) xB ∂xA ∂xB

(2.22)

lim ai ) 0 lim ai ) 1

(2.23)

and the limits xif0

xif1

Then a reasonable approach to treating the properties of the bulk real solution is to expand the activities in a power series in x consistent with the constraint eqs 2.22 and 2.23. Thus, the excess free energy, GE ≡ G - Gideal, is

GE ) kTNxAxB[A + B(xA - xB) + C(xA - xB)2 + ...] (2.24) with

N ) N A + NB

γi )

ideal

Ω A ) ΩB ) Ω

(2.21)

subject to the constraint of the Gibbs-Duhem equation

(2.15)

The reason this expression is not immediately applicable is, of course, that each of the system parameters contains an unknown dependence on mole fractions x ) (xA,xB) of the cluster. To make further progress, we require a model that describes the composition dependence of the thermodynamic variables. Since the clusters under consideration are liquidlike in character, we are naturally led to consider the theory of solutions as a framework to understand their composition dependence. The simplest model that one could adopt for the cluster is that of an ideal solution. For that simple case, we have i µ∞,i cl ) µ0 + kT ln xi

i µ∞,i cl ) µ0 + kT ln ai

(2.25)

Using the activity coefficients,

ΓΩi kTr

µ∞,i cl

solutions is provided by the general theory of real solutions.42 The chemical potential for the condensed cluster phase is written in terms of the activities, ai where i ) A,B, as

ideal

(2.20)

0 is the vapor pressure of a bulk liquid of pure where p∞,i species i. It is clear that the cluster composition has a rather trivial effect on critical cluster properties. The composition of the cluster is the same as the composition of the gas phase, and the critical radius is independent of composition. Of course, the ideal solution is an unrealistic approximation for almost any physical system except in extremely dilute cases. A more satisfactory account of

ai xi

(2.26)

we obtain

ln γA ) (A + 3B + 5C + ...)xB2 - 4(B + 4C...)xB3 + ... (2.27) ln γB ) (A - 3B + 5C + ...)xA2 + 4(B - 4C...)xB3 + ... These relations, based on the bulk properties of the 0 in eq 2.15. solution film, can then be used to expand p∞,i The size dependence of the cluster partial pressures is localized in the exponential term of eq 2.15. Thus, we should consider the composition dependence of line tension, Γ, and the partial molar lengths, Ωi. In principle, Γ is expected to exhibit a dependence on x since the surface energy is a byproduct of the monomer-monomer interactions that will, of course, differ for A-A, A-B, and B-B neighboring pairs on the cluster perimeter. Furthermore, one also expects the surface concentrations of the two species to differ from that of the overall composition of the cluster. At the level of the present theory, we assume that surface tension as it appears in eq 2.20 is independent of the cluster size and thus is equivalent to the bulk line tension of a 2D-solution with the same composition. This conjecture, of course, must be tested against numerical results. The Butler equation43 is often used as the basis for describing bulk surface tension of solutions. However, since this expression requires system parameters only obtainable from simulation (or experiment), we adopt a more direct procedure of expanding the bulk surface tension in a series to be fit to bulk simulation, that is,

Γ ) Γ0 + ΓAxA + ΓAAxA2 + ...

(2.28)

where Γ0 represents the line tension for single-component clusters. On the other hand, the monomer areas are

7134


Brunelli and Skodje

expected to show the least composition dependence in that they reflect the underlying lattice spacing. We can often approximate them as constants in order to minimize the number of fitting parameters. In summary, the formalism presented above provides a well-defined fitting procedure to model the equilibrium vapor pressure of binary clusters within the context of the Gibbs-Thomson theory. Numerical simulations can be carried out at a stable equilibrium point to determine the expansion parameters that best fit the data. The 0 , parameters of the bulk solution, required to obtain p∞,i can be obtained from a simulation of a flat interface without the need to model any finite clusters. The parameters of the line tension expansion may be determined from a simulation of the bulk solution, or from a fitting of the size dependence of the vapor pressure itself for very large clusters. III. The Kinetics of Evaporation and Condensation To understand the kinetics of surface clusters, we must go beyond equilibrium thermodynamics and consider the rate laws associated with the underlying atomic processes. In this section, we shall consider two of the most important processes that govern the evolution of clusters, the evaporation and condensation of atoms at the cluster perimeter from the surrounding two-dimensional gas phase. In previous work, these kinetic processes have been studied for single-component clusters. A most important observation made earlier bears reiteration. Namely, it was found14,22,23,26,27 that the evaporation rate from a cluster did not scale with the perimeter length, L, until the cluster size became very large. Thus, the growth rates for small clusters showed very strong deviations from results of the nominal thermodynamic scaling relations. In this section, we generalize this result to multicomponent clusters. In a kinetic description of the growth of a single binary cluster, the number of constituent monomers, (NA,NB), changes in time due to the capture of atoms from the gas phase and the evaporation of surface atoms to the gas. If the densities of gas atoms (number of atoms/available sites) at a radius r from the cluster center are denoted by FA(r) and FB(r), the rate equations for growth in the dilute gas limit are

dNi ) RicFi(r) - Rie dt

(3.1)

The capture and evaporation rate expressions, Ric and Rie, depend on (NA,NB) but not on (FA,FB). The value of r is selected as the instantaneous cluster radius where the capture process occurs. To explicitly solve eq 3.1, a diffusion model is required to obtain the gas densities at the cluster boundary as a function of time. The details of that model will, of course, depend on the nature of the boundary conditions employed to simulate the growth conditions. Since the two-dimensional diffusion problem with sources and sinks has been discussed at length in the literature, we omit a detailed analysis here. We note only that simple expressions for the diffusion current may be constructed under steady-state conditions for the cases of interfacetransfer-controlled and diffusion-controlled transport. (42) Pitzer, K. S. Thermodynamics; McGraw-Hill: New York, 1995; Chapter 12. (43) Butler, J. A. V. Proc. R. Soc. A 1932, 135, 348.

At equilibrium, dNi/dt will vanish, which implies a relation between the rates and equilibrium vapor density, i RicFeq i ) Re

(3.2)

where the density is related to the vapor pressure through the ideal gas law. Thus, if the Gibbs-Thomson expression for the equilibrium partial pressures is assumed at equilibrium, then only one of the two rate expressions is independent. The exact form of the size dependence of the evaporation and capture rates has been the subject of some discussion for unary clusters. Clearly in the limit of extremely large clusters, one expects that the evaporation rate should be proportional to the cluster perimeter, in our case 2πr ∼ N1/2. However, numerical simulations have indicated that this thermodynamic limit is approached quite slowly.14,22,23,27 Indeed, Lo and Skodje have found significant differences from the thermodynamic scaling even for model 2D clusters of size 100 000. Instead, it was found that often a simple power law expression, Re ∝ Nδ with δ ∼ 0.36, was much more accurate for intermediate sized clusters, N ∼ 500. Metiu and Rosenfeld22 have suggested an alternative expression which goes to the thermodynamic limit as N f ∞ yet seems capable of reproducing the power scaling empirically observed for smaller clusters, viz.,

[βr]

Re ) Rr exp

(3.3)

with

β)

ΓΩ kT

(3.4)

Thus the parameter β can be determined from a GibbsThomson analysis of the equilibrium vapor pressure while only the R should be fit to the dynamical evaporation data. Using the Metiu-Rosenfeld expression along with eq 3.2, the capture rate is then exactly a Brownian capture kernel,

Rc )

R r F∞

(3.5)

F∞ ) p∞/kT

(3.6)

which is how the expression was originally derived. Thus, given a prior determination of the equilibrium vapor pressure, the growth kinetics of a cluster may be modeled by using a one-parameter fitting to obtain R. The extension of the kinetic rate expressions to multicomponent systems can proceed, as before, by expanding the required parameters in power series of the mole fractions. The evaporation rate is given by

Rie ) aiRir exp

[ ] ΓΩi rkT

(3.7)

where the appropriate limiting behavior of Rie is enforced by writing the prefactor as the product of the bulk activity and a more slowly varying function, Ri(xA). The functions Ri(xA) reflect the escaping probability of surface atoms and differences between the species concentrations on the cluster surface compared to its interior. They are expanded as

RA ) R0A + R1AxB + R2AxB2 + ... RB ) R0B + R1BxA + R2BxA2 + ...

(3.8)


Langmuir, Vol. 19, No. 17, 2003 7135

where R0A and R0B are the escape parameters for singlecomponent clusters. The Brownian capture rates are then determined by balancing the evaporation rates with the capture rates separately for A and B,

Ric )

aiRi r F∞i

(3.9)

IV. Simulation of Equilibrium Vapor Pressure To investigate the empirical characteristics of the equilibrium partial pressures of binary nanoclusters and to test the validity of the Gibbs-Thomson formalism, we have undertaken a large-scale numerical simulation of a model system. The substrate is assumed to form a square lattice with periodic boundary conditions where the lattice constant l is taken, for simplicity, as unity. The two species, A and B, are structureless atoms that can bind to any open lattice site. Hence, we have ΩA ) ΩB ) l2 ) 1. A surface adlayer is allowed to evolve in time through singleatom hops to vacant nearest-neighbor lattice sites. The time evolution is affected using the dynamical Monte Carlo method12 based on a simple potential energy surface. Nearest neighbors and next-nearest atoms, separated by a distance r, are allowed to interact through a pairwise Lennard-Jones potential

[( ) ( ) ]

r0 V ) 4 r

12

r0 r

6

(4.1)

We select r0 so that the potential minimum is located at the lattice spacing, l. The well depths for the A-A, B-B, and A-B pairs are taken to be

AA ) 301 K

BB ) 319 K

AB ) 325 K

For these parameters, we expect the cluster solution phase to be miscible. The adsorption energy of the atoms to the substrate is assumed to be large so that no atoms can break away from the surface. The barriers for atomic hopping are determined by the simple broken-bond prescription. Thus, for a given hop i f f we have

∆E(iff) ) ∆E0 + E(f) - E(i)

(4.2)

where ∆E0 ) 60 K is the diffusion barrier on the bare substrate, and E(i) and E(f) are the sums of bond energies for the initial and final configurations summed over the local environments. There are thus 2 × 310 combinametrically distinct barriers that contribute. The hopping probability per unit time for a given event, i f f, is then given by

prob )

1 exp(-∆E(iff)/kT) τ

(4.3)

where the time scale τ is set so that the fastest hop occurs with unit probability at the simulation temperature of 200 K. The Monte Carlo simulation employs two random numbers to select an atom and a hop direction and a third to determine the success of the hop. To obtain the equilibrium vapor pressures, a circular cluster of size N′ with a random distribution of A and B in proportions xA′ and xB′ is created at the center of a 51 × 51 periodic lattice. The system is then propagated for millions of MC steps until a stable equilibrium is established. This equilibrium, with the total number of monomers held fixed, corresponds to the minimum free energy configuration depicted in Figure 2. The results are

Figure 3. A snapshot of the configuration of a binary cluster in equilibrium with a two-dimensional gas used in the Monte Carlo simulations. The composition of the cluster, and the gas, for a given set of constraints is determined by averaging several hundred such configurations.

time-averaged after the equilibrium is established, and also averaged over several hundred MC runs to minimize statistical noise. The final cluster size and composition, N and (xA, xB), are monitored as well as the final gas densities FA and FB. Since the ideal gas law is adequate for this case, the vapor pressures are directly proportional to the corresponding densities. This process is repeated for a variety of cluster sizes and compositions to obtain the complete functions pi(N,xA) that can be evaluated at any value of arguments, 0 < xA < 1 and 60 < N, through interpolation and extrapolation of the data. From inspection of the simulation results, the clusters are clearly liquidlike with high atomic mobility and the absence of persistent crystal faces. A snapshot of the lattice configuration of a mixed cluster of size 337 is shown in Figure 3. We consider first the vapor pressures for the two cases of pure clusters, xA ) 1 and xB ) 1. In Figure 4a, we show the size dependence of the vapor pressure versus cluster size plotted logarithmically. We see a clear linear dependence of ln(pi) upon 1/r, consistent with the GibbsThomson formula. Indeed, the Gibbs-Thomson formula appears to be accurate for clusters as small as N ∼ 60. The resulting line tensions, obtained from the slopes in Figure 4a, are quite distinct for A and B, ΓAl ) 443 K and ΓBl ) 503 K. Next, we investigate the vapor pressures for mixed clusters. In Figure 4b, we show the total vapor pressure, pA + pB, versus xA as a function of size for a variety of fixed compositions. Each composition yields a linear plot versus 1/r, again confirming the validity of the generalized GibbsThomson formula, eq 2.15. The slopes of these lines then can be used to determine the value for the quantity ΓΩ/ kT at each composition. In Figure 5, we plot the line tensions obtained from the slopes versus the mole fraction, xA. It is seen that while the absolute variation of Γ is only about 10% over the range of xA, the shape of the function is rather different from that normally observed for the bulk surface tension of normal solutions. We conjecture that this may be a consequence of the fact that the pairbond strength AB was chosen to be larger than either AA or BB. This may tend to keep the surface concentrations more mixed than is usually observed. As seen from the

7136


Brunelli and Skodje

Figure 4. Equilibrium vapor pressure versus cluster size. The vapor pressure is expressed as density using the ideal gas relation. In (a), the logarithm of the density is plotted versus 1/r, the inverse cluster radius, for the cases of clusters of pure A and pure B. The symbols represent the MC simulations, while the solid lines are the best straight-line fit. In (b), a similar plot is made for a number of mixed cluster compositions. All the results are for the total density, F ) FA + FB.

Figure 5. The line tension function, Γ(xA), obtained from the slopes of the curves ln(F) vs 1/r that were plotted in Figure 4. The cubic fitting, Γ ) Γ0 + ΓAxA + ΓAA xA2 + ΓAAAxA3, is shown with a solid line.

figure, the line tension function is very well represented by the cubic expression

Γ ≈ Γ0 + ΓAxA + ΓAAxA2 + ΓAAAxA3

(4.4)

To complete the parametrization of the GibbsThomson formula, we require an expression for the prefactors F∞i . This quantity could be obtained in either of two fashions: first, by the direct simulation of the flat (bulk) interface, or second, by extrapolation of the cluster curves shown in Figure 4 to infinite size, that is, 1/r f 0. We followed the second of these procedures. Despite the relatively small differences introduced into the binding energies, AA, BB, and AB, it was clear that the solution could not be accurately modeled using the ideal-solution relation, eq 2.20. As seen in Figure 6, the bulk vapor pressure versus xA shows significant nonlinearity, a hallmark of nonideal behavior. Instead, we have parametrized the bulk-solution activities up to fourth order, as

Figure 6. The equilibrium vapor density for the bulk interface vs mole fraction, xA. The densities obtained from MC simulations with finite clusters were extrapolated to infinite size assuming a Gibbs-Thomson size dependence. The circles (squares) represent the extrapolated partial density FA (FB), while the diamonds are the total density. The solid lines are the fitted analytical forms using a fourth-order expansion of the activity coefficient.

discussed above, to obtain the following fitted forms:

ln(γA) ) (A + 3B + 5C)xB2 - 4(B + 4C)xB3 + 12CxB4 (4.5) ln(γB) ) (A - 3B + 5C)xA2 + 4(B - 4C)xA3 + 12CxA4 where A ) -0.559, B ) -0.109, and C ) -0.044 at the simulation condition T ) 200 K. As seen in Figure 6 by the agreement between the numerical simulations and the solid curve based on eqs 2.21 and 2.26, this level of expansion seems quite adequate to model the bulk properties. To comprehensively test the fully parametrized expression for the partial pressures, we now show the results of


Langmuir, Vol. 19, No. 17, 2003 7137

Figure 7. The equilibrium vapor density of mixed clusters vs mole fraction, xA, for a variety of cluster sizes. The open circles, open squares, and filled circles represent the MC results for the densities FA, FB, and F ) FA + FB, respectively. The solid lines are the predictions of eq 2.15. The curvature of the partial densities versus xA reflects the nonideal character of the cluster solution.

MC simulations for the composition dependence of finite size mixed clusters along with the model predictions. In Figure 7, the partial and total pressures are shown as a function of xA for a variety of cluster sizes obtained from simulation (symbols) and based on the prediction of eq 2.15. The nonideal behavior previously observed for the bulk is seen to persist to all sizes. The predictions based on the model are quite satisfactory even for clusters smaller than N ) 100. Interestingly, we see a clear size dependence of the partial pressures that is well reproduced by the analytical expression. In addition to an overall increase in vapor pressure as the cluster size decreases, there is also a subtle shift in the shapes of the partial pressures that reflects the different xA dependence of the constituent functions, F∞i and Γ, in eq 2.15. For example, consider the shallow minimum of F(xA), which defines a most stable cluster composition. For very large clusters, N > 100 000, the value of xA(min) approaches the bulk limit of 0.184. However, for N ) 100 xA(min) has shifted considerably to 0.244. Clearly such size-composition correlations will be expected to affect the cluster distribution function observed in thin film coarsening. V. Simulation of Cluster Evaporation Rates The same model described in the previous section was employed to simulate the evaporation kinetics of binary clusters. The appropriate simulation conditions to obtain Re for single-component clusters have been discussed in depth and validated previously.14,27 Briefly, a cluster was created at the center of a large lattice and allowed to come to equilibrium by a long MC “presimulation”. Then, all the 2D gas atoms were removed from the lattice. The simulation was continued, and any atom that became disconnected from the cluster was removed before the next

MC step. The size and composition of the cluster were monitored, and the results are averaged over several hundred MC runs. The evaporation rates were then obtained numerically from ∆Ni/∆t. The number of MC steps taken per atom determined the time step, ∆t. Since the number of atoms was not conserved in the simulation, the number of surviving atoms continuously rescaled the time step. To test whether compositional nonequilibrium effects might develop in the simulation between the surface and interior atoms of the evaporating cluster, the evaporation rates were extracted for clusters of the same composition but with different initial sizes. In Figure 8, we show the size decay curves obtained for several initial cluster sizes with xA ) 0.2 versus time, averaged over many MC runs. After a brief transient period, during which the shape of the cluster equilibrates to the vacuum environment, it is seen that each cluster follows a common decay pathway. Thus, it appears that the evaporation rate can be extracted from the averaged decay of a single large cluster of a given composition, provided that the overall cluster composition is monitored. The evaporation rates versus size obtained by MC simulation of clusters consisting of pure A and pure B are shown in Figure 9. For pure clusters, the analytical expression, eq 3.7, is particularly simple since the activity term is unity. The difference in the rate between A and B is due mainly to the detachment prefactor, Ri, while the line tension has a much smaller influence. It is clear from the figure that the rate expression obtained from a oneparameter fitting to the prefactor is remarkably accurate over a wide range of cluster sizes. The results of the evaporation simulations for mixed clusters are presented in Figures 10 and 11. In Figure 10,

7138


Brunelli and Skodje

Figure 8. The total evaporation rate vs time obtained by MC simulation for four clusters of composition xA ) 0.2 but of different initial sizes. After a brief transient, each cluster yields the same rate law.

Figure 10. The logarithm of the total evaporation rate versus ln(N) for several cluster compositions. While the absolute rates depend on the value of xA, the slopes of the curves are very similar except at very small cluster sizes.

Figure 9. The evaporation rate vs cluster size for pure clusters, xA ) 1 and xB ) 1. The symbols show the results of MC simulations, while the solid curves are the predictions of eq 3.7 with the activities set to unity.

the evaporation rate versus cluster size is plotted on a logarithmic scale for several compositions. It is seen that while the magnitude of the rates shows clear composition dependence, the size scaling is nearly independent of xA thereby yielding roughly parallel curves. This behavior is consistent with the predictions of eq 3.7 with the parameters of our model system. The slope of the rate curves predicted by eq 3.7 is

βi ∂ ln Rie 1 | ) ∂ ln N x 2 2xN/π

(5.1)

where βi ) ΓΩi/kT ranges from 2.2 for pure A to 2.5 for pure B. Thus for even fairly small sized clusters N ) 300, the slopes will range only from 0.40 to 0.42. The dominant influence of cluster composition in this system is then not through the surface tension, but instead is through the single atom escape rates modeled by the prefactor in eq 3.9. The values of the factors Ri are numerically extracted from the simulations by fitting Rie/aireβi/r to the power expansions eq 3.8. It is found that the data are quite well represented using just the first term of the expansion, RA ≈ R0A ) 0.0327/τ and RB ≈ R0B ) 0.0266/τ. No discernible size dependence was observed for Rie/aireβi/r over the range N ) 100-100000, consistent with the model predictions. It is apparent, in this case, that the escape

parameters can actually be extracted from a simulation of the flat bulk interface. The numerical evaporation rates as functions of composition are shown for a number of cluster sizes in Figure 11 along with the predictions of the now fully parametrized rate law, eq 3.7. The agreement between simulation and prediction is seen to be generally quite good. The total evaporation rate, Rtot ) RAe + RBe , is seen to have a e nonlinear composition dependence similar to that of the total vapor pressure, and also a dependence on total size N. Since species A is the more volatile component of the solution, the rate is generally higher for clusters with higher mole fraction xA. However, because of the strong mixing term in the system Hamiltonian, the rate is seen to exhibit a shallow minimum at an intermediate xA value that defines the most stable cluster of size N. At this composition, the cluster evaporates the slowest. Not surprisingly, the most stable value of xA is very close to the point of minimum total vapor pressure noted previously. The xA value of the most stable cluster is seen to shift to larger values with increasing cluster size until the bulk limiting value is achieved for very large sizes. In kinetic simulations of cluster growth and decay, there is a tendency of the cluster composition to shift toward the most stable configuration relative to the composition of the gas phase. Since the rate laws governing cluster evaporation and condensation are now explicitly given through eqs 3.7 and 3.9, it is possible to model the combined growth-decay kinetics of an isolated surface cluster through the use of eq 3.1. To test whether this important kinetic process can be accurately described by this formalism, we have carried out a final set of MC simulations. An initial cluster is created on a fairly small lattice of 51 × 51 and is timeevolved until it equilibrates with the surrounding twodimensional gas. Then the density of the gas phase is uniformly increased, or decreased, so that the cluster is out of equilibrium. Provided that the initial cluster is larger


Langmuir, Vol. 19, No. 17, 2003 7139

Figure 11. The total and partial evaporation rates vs xA for a variety of cluster sizes. The open circles (squares) are the MC rates -dNA/dt (-dNB/dt), while the filled diamonds are the corresponding total rates, -(dNA/dt + dNB/dt). The solid curves represent the predictions of eq 3.7.

expect the simulation conditions to lie in the regime of diffusion-controlled growth/decay, we have adopted a simpler approximate treatment. Using the steady-state solution to the radial diffusion equation, we replace Fi(r) with the approximate expression44

Fi(r) ) FGT i (r) +

Figure 12. Cluster size vs time for several nonequilibrium initial conditions. The system consists of 347 B-atoms on a 51 × 51 lattice. The MC simulations are shown with symbols, and the predictions of eq 3.1 are depicted with solid lines.

in size than the critical nucleation value, Nc, the size of the cluster will grow/shrink to the stable equilibrium value depicted in Figure 2. The total mass is conserved in the simulation. In Figure 12, we show the time evolution of the total size of the cluster for several initial conditions with the same total mass, n ) 347. The cluster sizes are seen to relax to the equilibrium value of N ) 311 on a time scale of roughly 200τ. To apply the rate expression, eq 3.1, we require an explicit model for the monomer density at the cluster perimeter, Fi(r). This could be accomplished by numerically solving the two-dimensional diffusion equation with the appropriate boundary conditions. However, since we

2πD0(F0i - FGT i (r)) Rie ln(R/r)

(5.2)

where FGT i (r) is the Gibbs-Thomson vapor density adapted from eq 2.15. The free monomer diffusion coefficient, D0, is simply l2/τ in the present case. The evaporation rate, Rie, is taken from eq 3.7, and the boundary term ln(R/r) is set to a constant using a phenomenological screening length from the simulation. The mean monomer density, F0i , depends on r through the condition of conservation of mass. Under this constraint, the growth/ decay of the cluster is found to be nearly exponential in time as opposed to the power law expressions familiar from the theory of Ostwald ripening. It is seen in Figure 12 that the resulting predictions are in good agreement with the MC simulations. VI. Conclusions The kinetics of nanoscale surface clusters consisting of more than one chemically distinct species has been studied using a combination of analytic and numerical methods. It has been found that previous results obtained for unary clusters could be successfully generalized and extended to multicomponent clusters. The focus of the present work was on the growth and decay of binary liquidlike clusters (44) Bhakta, A.; Ruckenstein, E. J. Chem. Phys. 1995, 103, 7120.

7140


due to the exchange of monomers with the surrounding mixed two-dimensional gas phase. The fundamental processes are then the capture and emission of atoms (monomers) of types A and B. The four associated rates will depend on the total cluster size, N, as well as the mean cluster composition, measured by the mole fraction xA or xB. When the intracluster relaxation rate is rapid, the surface layer is in quasi-equilibrium with the interior and thus the rates may be parametrized by (N,xA) without the need for additional variables to represent the condition of the cluster surface. At total equilibrium, the capture and evaporation rates balance and produce a background gas with partial pressures pA(xA,N) and pB(xA,N). It was found that these partial pressures could be remarkably well reproduced using an extension of the unary GibbsThomson formula. The size dependence of the partial pressures arises solely due to the Kelvin term and is parametrized only on the line tension, Γ(xA). However, the composition of the cluster also influences the partial pressures through the bulk-solution pressures, p∞i , which are parametrized by the solution activities. Since the xA dependences of Γ and ai are different, a size-composition correlation develops for the pressures that was illustrated numerically in Figure 7. We expect that such correlations may have a significant effect on the kinetics of thin film coarsening. The evaporation and capture rates are related by detailed balance through the equilibrium vapor pressure, and in principle either can be used to fit the cluster kinetics. Motivated by previous work on unary clusters, we used the evaporation rate to model the decay process. Consistent

Brunelli and Skodje

with the unary case, the evaporation rates showed a strong deviation from the nominal scaling behavior, Rthermo ∝ evap xN leading to rates far greater than expected for small cluster sizes, N < 1000. The cluster composition affected the evaporation rate both through the bulk activity of the solution phase and through the surface tension via eq 3.7. For the model system considered, this leads to a most stable composition at each cluster size that possesses a minimum evaporation rate versus xA. It was shown that the growth-decay kinetics of an isolated cluster could be accurately modeled using the fitted expressions. This opens the possibility of using the rate equation approach to model the coarsening kinetics of multicluster ensembles that can consist of more than one distinct species. In a later work, we shall report the results of such a treatment of thin film kinetics. Finally, we note that while in principle the formalism we have developed should apply generally to miscible binary clusters, the confirming numerical simulations have only been carried out for one parameter set. Moreover, the model system chosen was only modestly nonideal. It would be quite interesting to treat the cluster kinetics of more strongly nonideal species. We conjecture that the size-composition correlations that developed in our model may, for such cases, become more profound. Acknowledgment. This work was supported by a grant from the National Science Foundation. LA034435M

Kinetics of Multicomponent Nanosize Clusters on Solid Surfaces

Recommend Documents