Point island models for nucleation and growth of supported nanoclusters during surface deposition

Analysis of nucleation and growth in gas and condensed phases is a classic problem in chemical physics,^1–3 the kinetics of which is traditionally analyzed in terms of Becker-Döring-type rate equations. Studies continue both exploring the fundamental theory⁴ and also including refinements to incorporate a more accurate molecular-level description of the process.⁵ Already by the 1950’s, the above treatments were refined to describe nucleation and growth during deposition on surfaces and to account for the feature that often one finds small critical sizes for stable nuclei in these systems.^6,7 Such mean-field (MF) level treatments for surface deposition were developed and applied extensively over the next few decades.⁸ However, only much later were certain fundamental shortcomings recognized.⁹ These shortcomings derive from a failure of MF treatments to capture subtle spatial aspects of the nucleation and growth process on surfaces. Avrami-type models also fail to capture these features.¹⁰ 

We consider here the nucleation and growth of nanoclusters (NCs) or islands during deposition on flat crystalline surfaces involving continuous random adsorption of atoms at a periodic array of adsorption sites together with rapid terrace diffusion of deposited atoms by hopping between neighboring adsorption sites.^8,9 The latter facilitates the nucleation of new stable NCs in competition with growth of existing NCs. Characterization of these far-from-equilibrium systems involves two basic classes of questions. The first relates to the spatial and size distributions of NCs, i.e., NC ensemble-level properties. The second relates to the structure of individual NCs which depends sensitively on the details of the relaxation or rearrangement within NCs of aggregated atoms. Here, we focus on the first class of questions where the detailed structure of the individual NCs is not a significant factor.

To this end, we will consider tailored deposition models which simplify the treatment of NC structure by regarding these as “point islands” occupying a single adsorption site, but carrying a label to track their size (in atoms).^9,11 See Fig. 1. Naturally, these models are applicable only for coverages in the pre-coalescence regime, i.e., below about 0.2 ML for two-dimensional (2D) islands, and for a significantly broader range of coverage for three-dimensional (3D) islands. These so-called “point island models (PIMs)” are particularly valuable for two reasons. First, they are effective in elucidating fundamental aspects of the nucleation process, e.g., removing the complicating effect of the expanding island footprint. We emphasize that the simplicity of PIMs does not detract from their capability to capture basic features of the NC distribution. Second, this type of modeling is particularly versatile as it can be readily modified to incorporate a diverse range features and alternative pathways impacting nucleation. Kinetic Monte Carlo (KMC) simulation of these stochastic non-equilibrium lattice-gas models is also readily implemented with PIMs.

In Sec. II, we describe the PIM for basic diffusion-limited homogeneous nucleation and growth processes and review the analytic theory for such processes. We describe key features of the transition from transient to steady-state behavior which are important for our subsequent analysis. We also emphasize open challenges for a precise beyond-MF theory of the NC size distribution. In Sec. III, we present a detailed analysis and new insights into scaling behavior for high surface diffusivity, as well as finite-size effects. Attachment-limited nucleation and growth is analyzed in Sec. IV, both for a standard PIM and also for refinements mimicking formation of 2D and 3D islands. Various other refinements of the PIM to treat modified nucleation and growth processes are briefly described in Sec. V. Discussion of applications of various versions of the PIM to a variety of specific systems, as well as conclusions, is presented in Sec. VI.

First, we specify the details of the PIM for diffusion-limited homogeneous nucleation and growth of NCs (also described as islands) during deposition.^9,11 Below, T will denote the surface temperature, and we set β = 1/(k_BT) where k_B is Boltzmann’s constant. The surface is described by a periodic square lattice of adsorption sites (although basic behavior is preserved for other lattices). See again Fig. 1. Deposition occurs randomly at rate F per site. Terrace diffusion of deposited adatoms involves hopping to nearest-neighbor (NN) sites at rate h per direction. We will write h = ν exp(−βE_d), where E_d denotes the terrace diffusion barrier. As indicated above, islands occupy a single adsorption site, but carry a label to track their size, s, in units of atoms. With regard to nucleation and aggregation, there are two reasonable choices. We can specify that these processes are implemented when a hopping adatom reaches either (A) the actual site on which another adatom or an island resides or (B) a NN site to another atom or island. We should also note that nucleation and aggregation can also occur via direct deposition: (A) on the same site as the adatom or island or (B) on a NN site. However, contributions from direct-deposition pathways are small, and thus are neglected in the analytic theory below which is focused on low island densities. Basic behavior will not depend strongly on the version (A) versus (B), and the difference between the models will vanish in the scaling limit h/F → ∞ for fixed coverage θ = Ft, where t denotes the deposition time. Results presented in Sec. III correspond to version (B) and in Sec. IV to version (A). In the simplest implementation of this model, there is no diffusion of stable islands, no extra barrier for adatom attachment to islands, etc. These and other modifications are discussed in Secs. IV and V.

Another key feature of the model is the specification of a critical size, i, such that only islands of s > i atoms are stable, meaning that adatoms cannot detach from such islands.^7–9 Let E_s < 0 denote the binding energy for adatoms within a sub-stable island of size s ≤ i, where E₁ = 0. One naturally chooses a value for E_s so that |E_s+1| > |E_s| since larger clusters should have stronger overall binding than smaller clusters. Then, in the absence of attachment barriers (so h gives the attachment rate per direction), the rate of detachment per direction from an island of size s is given by h exp[ − β(E_s−1 − E_s)] according to detailed-balance. Actually, critical size i is realized only if |βE_i| is not too large, so that detachment from islands of size i is facile on the time scale of aggregation. Otherwise, one has a smaller effective critical size. Furthermore, it is then expected that the population of sub-stable islands is quasi-equilibrated relative to the current density of diffusing adatoms (which reflects both the deposition flux and rates of nucleation and aggregation processes).^7–9 In this case, behavior is strongly dependent on the binding energy, E_i, for the critical cluster, and effectively independent of E_s for smaller s < i. We will focus on the behavior for i = 1 which corresponds to irreversible island formation where behavior at a specified coverage, θ, depends only on the ratio h/F. As an aside, we note that benchmark studies for i > 1 have often set E_s = 0 for s ≤ i, so that detachment from sub-stable islands occurs at rate h.⁹ 

The quantities of primary relevance in this study are the densities (per adsorption site) of adatoms, N₁, and of stable islands given by N_isl = Σ_s>iN_s, where N_s is the density of islands of s atoms. The coverage satisfies θ = Σ_s≥1sN_s in units of monolayers (MLs).

Figure 2 illustrates the two distinct regimes occurring during nucleation and growth of NCs during deposition.^9,12,13 Both regimes will feature in our analytic treatment and can be described as follows:

• In the initial transient regime, N₁ ≈ Ft builds up uniformly over most of the surface leading to nucleation of far-separated islands at near-random locations. Depletion zones (DZs) form around each island. N₁ is reduced within the DZ due to aggregation with that island which can be regarded as a sink for diffusion adatoms. DZ radii increase like R_DZ ∼ (hδt)^1/2, where δt is the time since nucleation.¹² Later, growing DZs collide to form boundaries of capture zones (CZs) which surround each island, as described below. These boundaries are constructed to be roughly equidistant from the relevant point islands (a non-trivial construction process, the details of which are not critical here).

• In the subsequent steady-state regime, DZs have expanded sufficiently so that the surface is completely covered or tessellated by the corresponding CZs surrounding each island. The basic feature of CZs is that most deposited atoms aggregate with the island associated with the CZ in which they are deposited, and ideally the growth rate of each island should be exactly proportional to the corresponding CZ area. Consequently, CZs must strictly be constructed by solving a boundary value problem for deposition-diffusion equation, but in practice CZs for point islands are reasonably described by Voronoi cells based on the island distribution.⁹ We emphasize that island nucleation persists and is dominant in the steady-state regime especially for small i.

In Sec. II A, we briefly review MF rate equation analysis for the behavior of the average island density.^8,9 We follow in Sec. II B with an analysis of the island size distribution (ISD), where MF treatments fail qualitatively to predict the correct behavior, and where open challenges for reliable beyond-MF treatments remain.⁹ 

Traditional MF rate equations are effective for describing the behavior of average island and adatom densities in 2D surface systems^8,9 despite ignoring logarithmic corrections discussed in  Appendix A. Refinement is needed for 1D systems as discussed in Sec. V and also in  Appendix A. We first describe the two key rates. The nucleation rate is given by K_nuc = σ_ihN₁N_i, where the capture number, σ_i, describes the propensity of critical clusters for the capture of diffusing adatoms. The aggregation rate is given by K_agg = σ_avhN₁N_isl, where σ_av is the average capture number for stable islands. The mean adatom and island densities satisfy

and

These equations are closed by invoking a quasi-equilibrium Walton relation⁷ for the density of critical clusters N_i ≈ exp(−βE_i)(N₁)ⁱ. Integration of (1) reveals the transient regime where dN₁/dt ≈ F so N₁ ≈ Ft = θ, followed by crossover to a steady-state regime. The latter reflects a balance between a gain of adatoms due to deposition and loss primarily due to aggregation, so that F ≈ K_agg, and N₁ ≈ F/(σ_avhN_isl). Substituting this form for N₁ into (1) for N_isl and integrating yields (3a) and (3c) below. Matching the transient and steady-state solutions for N_isl reveals crossover at (low) coverage of

Specifically, one finds that

in the transient regime,

at crossover where θ ≈ θ^∗, and

for θ = O(1), the latter being the steady-state regime. The scaling exponents satisfy

As an aside, one can immediately calculate the average island separation at crossover from L_isl^∗ = (N_isl^∗)^−1/2 and in the steady state from L_isl = (N_isl)^−1/2. Also see Fig. 2 for L_isl and L_isl^∗.

These results immediately yield several fundamental insights into the behavior of N_isl and of related quantities. First, consider the key result (3c) for steady-state N_isl. The F-dependence of N_isl determines χ and thus i. The Arrhenius energy of N_isl is given by E = (iE_d − E_i)/(i + 2), so that E = E_d/3 for i = 1 determining E_d. Second, a particularly significant observation is that the steady-state N_isl far exceeds its crossover value N_isl^∗ for large h/F since χ^∗ − χ > 0, dramatically so for smaller i. Correspondingly, L_isl^∗ far exceeds L_isl. The predominance of nucleation in the steady-state regime for large h/F is quantified by the ratio

Third, we consider the behavior of N₁. One can use the steady-state relation N₁ ≈ F/(σ_avhN_isl) together with (3c) to obtain

in the steady state. By matching (6) to N₁ ≈ θ in the transient regime, one immediately recovers the result (2) for the crossover θ^∗. From (6), it is also clear that N₁ < < N_isl for i = 1, N₁ ∼ N_isl for i = 2, and N₁ > > N_isl for i > 2 when θ = O(1) for large h/F.¹⁴ Fourth, we consider the mean island size, s_av. If most adatoms (or even a finite fraction of them) are incorporated into islands, then one has s_av ∼ θ/N_isl, so that

in the steady state. It is common to write s_av ∼ θ^z and N_isl ∼ θ^1−z with dynamic scaling exponent z = (i + 1)/(i + 2). One can show that a finite fraction of adatoms are incorporated into islands at crossover, so that s_av^∗ ∼ θ/N_isl^∗ applies there also. It then follows that the corresponding s_av^∗ is far below s_av for smaller i and large h/F, with s_av^∗ = O(1) for i = 1, so that island growth up to crossover is “negligible.” Fifth, one can show that the nucleation rate has the form

where k(u) ∼ uⁱ⁺¹ (correcting a typographic error in Ref. 9) in the transient regime, and k(u) ∼ u^{−(i+1)/(i+2)} in the steady-state regime. Thus, k(u) and the nucleation rate are strongly peaked around crossover u ≈ 1. The universal nature of the prefactor in front of k(θ/θ^∗) for both transient and steady-state regimes follows since the distinct forms of N₁ for these regimes match at θ ≈ θ^∗.

Of primary importance in analysis of the ISD is the rate at which diffusing adatoms aggregate with an island of size s, denoted by K_agg(s) = σ_shN_sN₁ where σ_s is the “capture number” for islands of size s. One then obtains^8,9

for s > i. To recover the simpler reduced Equation (1) for N_isl, one sums (9) over s > i to obtain dN_isl/dt ≈ K_agg(i) = K_nuc. We do not give details here, but one can also show that

where σ_av = Σ_s>iσ_sN_s/Σ_s>iN_s. We remark that these rate equations are effective in describing the mean density of islands (and adatoms), as confirmed by comparison with KMC simulation. These equations require as input a reasonable choice of capture numbers. This is relatively simple for point islands where there is a unique size-independent capture number for stable islands. For compact 2D islands, a self-consistent mean-field analysis of the σ_s is required, and then σ_av follows with knowledge of the ISD.¹⁵ We emphasize, however, that rate equations utilizing even sophisticated self-consistent mean-field analysis of capture numbers fail to describe the ISD. See below.

We will focus on the analysis of scaling solutions to (9) for large h/F or large s_av in the steady-state regime. Introducing a scaled island size variable, x = s/s_av, we anticipate that these solutions have the form^9,11,16

for large s_av. Here, one has that $∫0∞f(x)dx=∫0∞xf(x)dx=1$⁠. The PIM is special in that f(x) has no explicit θ-dependence, a feature which is lost for islands of finite extent. Our analysis also introduces a scaling function for the capture numbers, σ_s/σ_av ≈ c(x).¹⁷ Then, we perform a quasi-hydrodynamic coarse-graining analysis of (9) treating rescaled island size as a continuous variable. Also exploiting the temporal scaling, s_av ∼ (Ft)^z, and the steady-state relation for N₁, one obtains¹⁶ 

with the auxiliary relation c(x = 0) f(x = 0) = 1 − z > 0. See also Ref. 17. The auxiliary relation shows that f(x = 0) > 0 is in qualitative contrast to the most commonly assumed heuristic Amar-Family form for the scaling function

where¹⁸ 

One advantage of the PIM is that in a traditional MF theory, islands of any size, s, have the same capture number, as they have the same spatial extent. This in turn implies that c_MF(x) ≡ 1. In this case, (12) exhibits a singularity in the integral yielding the explicit non-analytic form^9,11

where

and H(u) = 0 (1) for u < 0 (u > 0) is the Heaviside step function. In fact, (13) is also qualitatively incorrect since the exact c(x) increases in a way which avoids this type of MF singularity.

It is clear that the above formulation (12) for the ISD provides an exact but incomplete theory since the form of c(x) is not yet specified. Certainly, adopting a simple MF form c_MF(x) ≡ 1 is inadequate. Perhaps the greatest open challenge is to provide a quantitative beyond-MF theory for c(x). We emphasize that any such formulation must account in detail for the subtle spatial aspects of nucleation during the steady-state regime. These spatial aspects, together with the persistent nature of nucleation in the steady-state regime, control the form of c(x) and thus the form of f(x).

Here, we just briefly indicate the current status of this effort. To this end, it is useful to adopt a geometric picture of adatom capture by islands, as already implied in the CZ picture of Fig. 2. The idea is that the CZ area measures the rate of aggregation with an island, so that K_agg(s)/N_s = σ_shN₁ = FA_s, where A_s denotes the mean area of CZs for islands of size s.^8,9,16 It is also expected that CZs for PIM constructed as Voronoi cells will reasonably describe the true CZs for which CZ area exactly describes island growth rate.^9,19 Then, introducing a scaling function for A_s via A_s = A_ava(x), where A_av = 1/N_isl is the mean CZ area, it follows that a(x) ≡ c(x). Thus, the most effective theories for beyond-MF c(x) are actually formulated as rate-equation-type theories for A_s and thus for a(x). Indeed, one can write down a heuristic rate equation for the fractional CZ area, A_sN_s, associated with islands of size s accounting readily for gain and loss due to island growth, and less easily for loss associated with nucleation of new islands whose CZs “cut into” this area.²⁰ The most sophisticated approach starts with the joint probability distribution, N_s,A, for islands of size s and CZ area A.^21–25 A moment analysis leads to a refinement of the above-mentioned heuristic rate equation for A_s which in turn leads to an equation for the scaling function a(x). However, this equation requires as input information on subtle spatial aspects of nucleation such as the probability that the CZ of a new island overlaps that of an existing island of size s and also the typical amount of overlap.²⁴ 

The predictions of the MF rate equation theory of Sec. II A for the scaling behavior of N_isl with h/F and with θ have been confirmed by precise KMC simulation analysis at least for i = 1 and i = 2.^9,11 We note, however, that there are corrections to ideal asymptotic scaling due to finite h/F and due to logarithmic corrections.^9,11 (Logarithmic factors were not included in the capture numbers in our analysis, although such refinement is readily implemented. See  Appendix A.) Thus, here we focus on more delicate scaling issues for the ISD and related quantities.

Deviations from asymptotic large-h/F scaling for the ISD are more severe than for N_isl, not surprisingly since this quantity is more sensitive to the details of the nucleation and growth process and is not described by a simple MF theory. Previous studies have tended to explore the approach to the scaling limit for fixed rather low θ ∼ 0.1 ML and increasing h/F. These studies suggested that the ISD achieves the limiting scaling form somewhat slowly only when h/F is above 10⁹ (corresponding to rather demanding simulations for large system sizes). This feature is illustrated in Fig. 3(a). For an assessment of convergence to the scaling limit of various quantities describing the shape f(x) of the scaled ISD, in Fig. 4 we show the coefficient of variation, C = (〈s²〉_c)^1/2/〈s〉, the skewness, S = 〈s³〉_c/(〈s²〉_c)^3/2, and the (excess) kurtosis, K = 〈s⁴〉_c/(〈s²〉_c)². Here, 〈sⁿ〉_c denotes the nth cumulant of the ISD, where 〈s²〉_c = 〈(s − 〈s〉)²〉, 〈s³〉_c = 〈(s − 〈s〉)³〉, 〈s⁴〉_c = 〈(s − 〈s〉)⁴〉 − 3(〈s²〉_c)². C gives the standard deviation of f(x) in the scaling limit, S describes lack of reflection symmetry about s = 〈s〉, and K measures the weight of the distribution in the tails relative to a Gaussian distribution. Simple extrapolation to h/F → ∞ indicates that C_∞ ≈ 0.46 to 0.47, S_∞ ≈ − 0.32 to −0.41, and K_∞ ≈ − 0.83 to −0.89. For comparison, for the MF scaling function f_MF(x), one has C_∞ = 1/√5 ≈ 0.447, S_∞ = − 2√5/7 ≈ − 0.639, and K_∞ = − 6/7 ≈ − 0.857.

To provide a framework to understand the convergence to the scaling limit (and also to potentially facilitate alternative analysis of this limit), we suggest that deviations from the asymptotic form of the ISD are controlled by the extent of persistent nucleation in the steady-state regime (beyond the transient regime). This extent is naturally quantified by the ratio N_isl/N_isl^∗ ∼ (θ/θ^∗)^1/(i+2), and thus equivalently by the ratio θ/θ^*. Thus, since θ^∗ ∼ (h/F)^−2/(i+3), proximity to the scaling limit is reflected by the magnitude of the combination ℜ = θ^(i+3)/2(h/F). Asymptotic behavior is usually probed for smaller θ = 0.1 ML, say, by increasing h/F to at least 10⁹. The above observations suggest that the asymptotic behavior can instead be probed retaining moderate h/F ∼ 10⁶, say, by sufficiently increasing θ. Support for this idea is provided by results in Fig. 3(b) which show that the shape of the ISD for i = 1 as θ increases from 0.05 to 0.5 ML for fixed h/F = 10⁶ appears to approach the asymptotic scaling form for h/F → ∞. Specifically the shape for θ = 0.5 ML and h/F = 10⁶ corresponds closely to that for θ = 0.1 ML and h/F = 2.5 × 10⁷ where for i = 1 one has ℜ = θ²(h/F) = 2.5 × 10⁵ in both cases.

Finally, we mention that one can of course perform more detailed analysis of the approach to the scaling limit for the ISD and various cumulant ratios indicated in Fig. 4 describing its shape. In addition to scaling with increasing h/F, one can consider scaling with increasing θ. However, a definitive characterization of behavior is difficult even with our current extensive data. Analysis for K either increasing h/F or θ yields reasonably consistent results, and a well-defined limiting value K_∞. Analysis of S is more complex due to a slower approach to limiting behavior as h/F → ∞. We leave further analysis and discussion of these subtleties for a separate publication.

In practice in performing simulations, one analyzes behavior for finite L × L site systems with periodic boundary conditions. One hopes that the side length, L (in units of surface lattice constant), is large enough to avoid finite-size effects. In general, one expects such effects to become significant when the system size, L, becomes comparable to an intrinsic characteristic length. One potential complication for the considered systems here is that in addition to the natural characteristic length, L_isl, for θ = O(1), there is a significantly larger length, L_isl^∗. Certainly, choosing system size L < < L_isl^∗ will greatly impact system evolution in the transient and crossover regimes. The initial transient regime will no longer involve island nucleation at essentially random locations. However, we will show that for such sizes, behavior of at least basic quantities in the steady-state regime is not affected, provided that L > > L_isl. Results shown in Fig. 5 for the island density in a finite system, N_isl(L), relative to N_isl(L = ∞) do reveal clear scaling with L/L_isl. This contrasts previous speculation.^9,13 Furthermore, limiting behavior is achieved quite quickly already by L ≈ 4L_isl. As an aside, defining L_isl based on either finite (L_isl = [N_isl(L)]^−1/2) or infinite (L_isl = [N_isl(L = ∞)]^−1/2) systems, or using the asymptotic form (L_isl = θ^−1/6R^1/6), yields similar data collapses. See Fig. 5(b). The last choice of L_isl has practical utility, since from the results in Fig. 5(b), it allows ready assessment of the minimal L which must be used in simulations, for given θ and h/F, to avoid finite-size effects. One might anticipate more complex behavior for more “delicate” quantities characterizing the ISD shape, but in  Appendix B we show that these also exhibit simple scaling with L/L_isl.

To analyze the behavior for L/L_isl = O(1) in more detail, we introduce the probability P_n(θ) that there are at least n islands at coverage θ for a finite L × L site system. Equivalently, one can consider the probability, p_n = P_n − P_n+1, for exactly n islands (with p₀ = 1 − P₁). Then, one has that

where

For small L < L_isl, one expects a “single island regime” with P₁ ≈ 1 and P_n≥2 ≈ 0, so that 〈n〉 ≈ 1 and N_isl(L)/N_isl(∞) ≈ (L/L_isl)⁻², as confirmed in Fig. 5. Moreover, one anticipates that as L increases above L_isl, the populations of P_n≥2 increase quickly. This behavior is supported by the simulation results for P_n for various L/L_isl shown in Fig. 6(a).

Further elucidation of the behavior shown in Fig. 5 comes from making analytic estimates for P_n, at least for small n, and thus for N_isl(L) for L ∼ L_isl and θ = O(1). To this end, we start by considering nucleation of the first island. The probability, Q₁, that no island exists at coverage θ (where Q₁ = p₀) is given by

where

Here, K_Tnuc0(θ) denotes the total nucleation rate in the L × L system (with periodic boundary conditions) at coverage θ in the absence of any island. This result shows that choosing L to correspond to L_isl for θ = O(1), the first island is nucleated at a coverage

For large h/F, θ_nuc(1) is just above θ^*, and far below θ = O(1). Thus, for L ∼ L_isl, the first island is nucleated very early, and P₁ = 1 − Q₁ ≈ 1 for most θ above small θ_nuc(1).

After nucleation of the first stable island, a steady-state adatom density, N₁(ss), is quickly established in the finite system. The basic features of N₁(ss) follow from solving F + h∇²N₁(ss) = 0 in an annular domain with inner radius r_in = 1 and outer radius r_out ≈ L/2.²⁶ From this solution, one finds a typical magnitude for N₁(ss) ∼ (h/F)⁻¹L², revealing that the total nucleation rate, K_Tnuc1, for this single-island system scales like

where this quantity is independent of time or θ. Then, the probability that the second island has not been nucleated at coverage θ is

and

with b = O(1). The existence of the third island will depend on the time-independent total nucleation rate, K_Tnuc2, for the third island in the presence of two existing islands. Actually, K_Tnuc2 depends on the position of the second island relative to the first. However, it is clear that K_Tnuc2 scales in the same way as K_Tnuc1 with respect to L/L_isl. The same applies for nucleation of subsequent islands. See  Appendix B.

For L close to L_isl with significant probability of a small number of islands, then one has that 〈n〉 ≈ P₁ + P₂ + …, and consequently that

where

The corresponding behavior with b = 1.1 (consistent with behavior of P₂ in Fig. 6(a) for L/L_isl just above unity) is shown in Fig. 6(b) also including the contribution from P₃.

Analysis of the crossover regime for higher L/L_isl ≈ 1.5–4 is challenging. Nucleation rates become more complex as they depend on the locations of the multiple previously nucleated islands (which are not random as new islands prefer to be nucleated far from existing islands). However, from a more general perspective, just as convergence to the scaling limit in Sec. III A is controlled by N_isl/N_isl^∗, one might expect that convergence to L = ∞ behavior reflects the ratio of the number of islands for L > L_isl to that when L = L_isl which is unity. This ratio is given by L²N_isl = (L/L_isl)².

Our treatment of the effect of additional barriers on nucleation and growth^27,28 refines the standard PIM. Here, for simplicity, we just describe these changes for the case of irreversible island formation i = 1, and below we only consider this case. We now include separate hopping rates h₀ = h for terrace diffusion, h₁ = hexp(−βδ^∗) for nucleation, and h₂ = hexp(−βδ) for aggregation. Specifically, h₁ and h₂ describe the rate of the last hop leading to nucleation and aggregation, respectively. Since δ ≥ 0 and δ^∗ ≥ 0, h₁ and h₂ are generally reduced below h₀ = h, corresponding to inhibited nucleation and aggregation, and we set r_j ≡ h_j/h ≤ 1, for j = 1 or 2. It is sometimes useful to introduce L_δ^* = exp(βδ^∗) = 1/r₁ and L_δ = exp(βδ) = 1/r₂, which denote “attachment lengths” (in units of the surface lattice constant) for nucleation and aggregation, respectively, larger attachment lengths implying more difficult attachment processes.²⁹ Note that r(=r₁ = r₂) = 1 corresponds to the classic diffusion-mediated PIM considered in Secs. II and III. Here and below, we always use r for r₁ = r₂. A key feature of these models, discussed further below, will be that increasing attachment barriers (or reducing r_j) makes the adatom density more spatially uniform relative to the diffusion-limited regime r = 1.

One benchmark case of the model sets δ = δ^∗ so that h₁ = h₂, and r = exp(−βδ) ≤ 1. KMC simulation results for i = 1 shown in Fig. 7 indicate that the scaling behavior N_isl ∼ (h/F)^−1/3 familiar in the classic PIM for r = 1 is actually preserved even for r < < 1. This behavior is explained below in Sec. IV B. The other key feature of the simulation results is that N_isl increases with increasing barrier δ (or with decreasing r) at fixed h/F. Note that the increase in N_isl with increasing δ or decreasing r helps avoid corruption due to finite-size effects in simulation for large h/F. Just increasing the barrier for nucleation would reduce N_isl, but increasing the barrier for aggregation will have the opposite effect boosting the adatom density and thus the nucleation rate. For δ = δ^∗, the latter effect dominates. This feature has been observed in previous experimental and simulation studies mentioned in Sec. VI. A quantitative analysis of behavior for fixed h/F and increasing δ based on the results described below in Sec. IV B suggests that N_isl ∼ exp(βδ/3) = r^−1/3 for δ = δ^∗. Our simulation data are consistent with this behavior.

In addition, we explore evolution of the ISD for increasing the attachment barrier, i.e., for decreasing r, where we simultaneously increase h/F to maintain a roughly constant N_isl. Results shown in Fig. 8 suggest that the shape evolves towards the singular MF form (13) for i = 1 described in Sec. II A. This behavior is naturally understood from the feature that the more uniform adatom density leads to adatom capture which is more independent of island size s. This in turn implies that appropriate rate equation analysis will naturally produce MF-type behavior. See Sec. IV B for further discussion.

Finally, we discuss a second benchmark case of the model setting δ^∗ = 0 so that h₁ = h, but δ ≥ 0 so that h₂ ≤ h and r₂ = exp(−βδ) ≤ 1, i.e., we inhibit aggregation but not nucleation. In this case, it is obvious that N_isl will increase with increasing δ, and the prediction of the formulation in Sec. IV B is that N_isl ∼ exp(2βδ/3) = (r₂)^−2/3 for δ^∗ = 0. We have checked that our simulation results are consistent with this behavior. For this type of model, there is a general and reasonable perception that inhibiting aggregation by inclusion of the additional barrier will lead to more persistent nucleation relative to the standard model with r₂ = 1. This, in turn, suggests that the variation of N_isl with coverage, θ, should become closer to linear in contrast to the strongly sublinear behavior quantified by (3c) for r₂ = 1. Clearly, an infinite barrier, δ = ∞, allows formation only of dimers, so one has that N_isl ≈ θ/2. Indeed, KMC simulation results for i = 1 shown in Fig. 9 reveal that increasing δ (or decreasing r₂) does ultimately produce linear behavior, but this transition only occurs “very slowly.” For expected typical physical magnitudes of barriers where βδ ≤ 4 (e.g., δ ≤ 0.1 eV at 300 K), it is clear that strongly non-linear variation of N_isl with θ should still be observed.

In this subsection, we focus on the regime of large additional barriers for nucleation and growth, so that these processes are attachment-limited rather than diffusion-limited. This, in turn, means that the adatom density, N₁, will be relatively uniform across the surface. We will exploit this feature to develop an analytic theory for island nucleation and growth which is appropriate for real 2D (d = 2) and 3D (d = 3) island systems, as well as for standard point islands. To this end, we note that if P_s denotes the perimeter length for an island of size s, then the aggregation rate for that island for quasi-uniform adatom density, N₁, is given by

It is clear that $P̃s∼1$ for standard point islands and that $P̃s∼s1/d$ for 2D and 3D islands. Note that one can also use the latter form for standard point islands by assigning d = ∞.

The form (20) suggests a simple refinement of the prescription of aggregation with large extra barriers within the PIM framework in order to mimic behavior in real 2D and 3D island systems. This refined PIM retains a hop rate h₀ = h for terrace diffusion, and h₁ = hexp(−βδ^∗) for nucleation, but now assign h₂ = h₂(s) = hexp(−βδ) s^1/d for aggregation. This assignment reflects the enhancement of aggregation rate with increasing island size and thus perimeter length. Certainly for exp(−βδ) s^1/d < 1, this formulation with d = 2 or d = 3 will incorporate the enhanced propensity for adatom capture at larger islands.³⁰ As an aside, we note that this refined PIM is even better suited for analysis of 3D islands than that of 2D islands since the “footprint” of the former is smaller, or more point-like, at a fixed coverage. Thus, the coalescence regime occurs at higher coverage for 3D islands, and the refined PIM accurately describes behavior for a broader range of coverage.

Before presenting simulation results for these refined PIM, we develop an appropriate analytic treatment for their behavior. To this end, we consider the average aggregation rate per island, $kagg=exp(−βδ)hP̃avN1$⁠, where the average perimeter length satisfies

In the steady-state regime, N₁ is determined by a balance between the rate of gain per CZ due to deposition, k_dep ∼ F(L_isl)², and the rate of loss, k_agg, due to aggregation with the island in that CZ. This implies that

With this result, we can determine the nucleation rate,

Then, integration of dN_isl/dt = K_nuc yields the expression for the island density

For the integral determining the θ-dependence, one should just integrate over the steady-state regime for θ ≥ θ^∗ for some suitable initial conditions.

Focusing on the scaling with h/F, we obtain

The results for the standard PIM (d = ∞) and for the refined PIM mimicking 3D islands (d = 3) are new, and those for the refined PIM mimicking 2D islands (d = 2) are found in Ref. 27. The standard PIM reveals the same scaling with h/F and with θ as for diffusion-mediated growth with δ = δ^∗ = 0. The scaling N_isl ∼ (h/F)^−1/3 for PIM with i = 1 is clearly illustrated in Sec. IV A.

The study in Ref. 27 in fact performed a more comprehensive analysis for 2D islands assessing the crossover from diffusion-limited to attachment-limited growth with increasing magnitude of the additional barriers. The basic idea, which extends to other island structures, is that crossover should occur from the diffusion-limited to the attachment-limited regime as the relevant attachment length, L_attach, increases beyond the mean island separation, L_isl. Another perspective is that for fixed high barriers, and thus fixed r_i < < 1, one expects a transition from attachment-limited to diffusion-limited behavior with increasing h/F as L_isl grows to exceed L_attach. For the regime of modified scaling where L_isl < < L_attach, we do not expect significant finite-size effects in an L × L system for L < L_attach, provided that L is significantly larger than L_isl.

Next, we comment briefly on the expected form of the ISD for large attachment barriers. For standard point islands (d = ∞) where $P̃s∼1$⁠, as noted above, uniformity of the adatom density, N₁, implies that islands will capture adatoms at the same rate independent of size. Thus, the rate equation (9) for the densities, N_s, of islands of size s will apply with K_agg(s) = σ_shN_sN₁ and σ_s = σexp[ − βδ] independent of s. This in turn indicates that the MF results (13) described in Sec. II B will be applicable. For 2D and 3D islands, one instead has that $Kagg(s)∝hexp(−βδ)P̃sNsN1$⁠, where $P̃s∼s1/d$⁠. This feature modifies the form of the ISD as explored in previous analyses.^31,32

Finally, we describe selected KMC simulation results for our modified PIM mimicking 2D and 3D islands for i = 1 and with large δ = δ^∗, so that r (=r₁ = r₂) < < 1. See Fig. 10. The goal is to probe modified scaling with h/F relative to behavior for point islands. For r = 0.01, the attachment length, L_δ ≈ 100, is well above the mean island separation, L_isl, for h/F up to about 10¹² from Fig. 7. Thus, modified scaling should be found for such h/F. From Fig. 10(a), we estimate that χ ≈ 0.42 for h/F from 10⁷ to 10¹² for 2D islands [significantly below χ = 1/2 from (25) for d = 2]. Similar analysis for 3D islands yields χ ≈ 0.39 [versus χ = 3/7 ≈ 0.428 from (25) for d = 3]. Plausibly, scaling exponents will be closer to analytic estimates for smaller r = 0.001 corresponding to L_δ = 1000. This L_δ is far above the average island separation for all h/F accessed by simulation. From the data in Fig. 10(b), we find that χ ≈ 0.44 for 2D islands, and χ ≈ 0.41 for 3D islands (versus χ = 1/3 ≈ 0.333 from (25) for d = ∞).

For the above small r, it is challenging to precisely quantify crossover from modified scaling for 2D or 3D islands for smaller h/F to scaling with χ ≈ 1/3 for large h/F (just because it is difficult to obtain precise data for such large h/F). One possibility is to select a larger r = 0.04, say, so that crossover should occur for smaller h/F likely around 10⁸ to 10⁹. In this case for 2D islands with r = 0.04, we do find a transition from higher values of χ ≈ 0.39 for h/F from 10⁶ to 10⁹, say, to lower values of χ ≈ 0.34 for h/F from 10¹¹ to 10¹⁴.

In this section, we highlight the versatility of PIM for treating a diverse variety of mechanisms as well as other possible features of nucleation and growth systems.

It is straightforward to modify the standard model to include anisotropic or even 1D diffusion.¹¹ In the latter extreme 1D case, one must specify whether nucleation and aggregation only occur with adatoms and clusters in the same row (version A) or also in the adjacent rows (version B) along which diffusion occurs. For strong anisotropy, the basic scaling behavior of the island density,^11,33

with χ = 2/(2i + 1), for critical size i, is distinct from the isotropic case for both versions A and B (which exhibit similar behavior for large h/F). This behavior is obtained from modified rate equations accounting for the distinct nature of diffusion in 1D as indicated in  Appendix A. Crossover from isotropic to strongly anisotropic behavior in PIM can be analyzed and quantified.³⁴ 

Diffusion of small stable clusters also impacts scaling of N_isl in nucleation and growth if diffusivity is comparable to that on single adatoms. Here one introduces a new parameter, i^∗, wherein clusters of size i^∗ and smaller are mobile. Thus, the cases i = 1 and i^∗ = 2 correspond to irreversible island formation with mobile dimers. For i = 1 and general i^∗ where all mobile clusters have hop rates comparable to h, a rate equation analysis reveals that N_isl ∼ (h/F)^−χ with χ = i^∗/(2i^∗ + 1).^35,36 The basic trend in behavior is that increasing i^∗ allows greater incorporation of small stable mobile islands with larger stable immobile islands, and thus N_isl decreases. These models can be naturally implemented within a PIM framework specifying desired hop rates for clusters with sizes s ≤ i^∗. We have implemented these models for i = 1 and i^∗ = 2 and 3 (in addition to the standard model for irreversible island formation with i = i^∗ = 1). Results shown in Fig. 11, for the case where all mobile clusters have the same hop rate, h, as that of adatoms, are consistent with the above scaling law, but also illustrate the extent of decrease in N_isl produced by small cluster mobility.

In some systems, interactions between adatoms are not just dominated by very short-ranged attractions responsible for NC formation^28,37 (as assumed above). Also, any additional barriers inhibiting attachment may not just impact the last hop before nucleation or aggregation, but also hopping rates for further separated adatoms. The simplest scenario for implementation of long-range interactions is to modify PIM model type A (where implicitly there are zero range attractions since nucleation and aggregation occur when an adatom hops onto the same site as another adatom or island). In this case, we specify an additional longer-range interaction, ϕ₁(j, k) between adatoms at adsorption sites with separation r_s = (j, k) in units of surface lattice constant for r_s = |r_s| > 0, e.g., repulsive ϕ₁(r_s) = ε/[1 + (r_s/r_c)ⁿ] for ε > 0, where r_c measures the interaction range.²⁸ Distinct interactions ϕ_s(j, k) are prescribed between an adatom and island of size s.

The key ingredient for models of nucleation and growth is to specify barriers for all hops. For a diffusing adatom, one can determine the total interaction Φ_tot with other islands and adatoms by summing the above pair interactions, ϕ, for both the initial (init) and final (fnl) states before and after hopping. In practice for large h/F and small i, one expects the dominant contribution to come from one ϕ associated with a single nearby island or adatom. One can write the diffusion barrier as E_act = E_d + Φ_tot(TS) − Φ_tot(init). Here, Φ_tot(TS) is the total interaction at the transition state for hopping which is in principle not determined by the above ϕ. For slowly varying ϕ, one can reasonably utilize a symmetric Brønsted–Evans–Polanyi-type approximation, Φ_tot(TS) = ½[Φ_tot(init) + Φ_tot(fnl)].^38,39

For heteroepitaxial growth, a lattice mismatch between the substrate and the crystalline form of the deposited material leads to strain buildup, and thus a significant energetic penalty for larger islands inhibiting their growth.⁴⁰ A simple strategy to model this behavior is to reduce the rate of the last adatom hop leading to aggregation in the PIM for larger islands by a factor of q(s) ≤ 1. Specifically, one chooses a threshold size s^∗ for growth inhibition and sets q(s) ≈ 1 for s < < s^∗ and q(s) ≈ 0 for s > > s^∗. A simple default choice is q(s) = tanh[(s − s^∗)/w]. For this model, one can also include an additional barrier, δ, for nucleation and aggregation so those rates are reduced by an additional factor r = exp(−βδ). In Fig. 12, we show simulation results for a modified PIM with i = 1, h/F = 3 × 10¹³, s^∗ = 1000 for various w for both r = 1 and r = 0.01. Compared to models with no growth inhibition factor q(s), there is the “ready development” of a linear increase of N_isl with θ, noting that s_av = θ/N_isl quickly exceeds s^∗ for the above choice of model parameters.

A distinct nucleation mechanism has been suggested for some heteroepitaxial metal-on-metal systems.^41,42 Here, deposited adatoms can exchange with atoms in the top surface layer thereafter remaining essentially immobile and embedded in that layer. Furthermore, these embedded adatoms provide nucleation centers for subsequently deposited adatoms. This nucleation scenario has been labeled as i = 0 since the single embedded adatom forms a stable nucleus for subsequent island growth. A PIM for this process has been implemented previously. In addition, rate equations were also developed to describe a potentially unusual variation of island density with temperature.⁴² 

For sequential codeposition of species with strongly differing surface mobilities, the nucleation and growth processes depend strongly on the order of deposition.^43,44 If the less mobile species is deposited first, a higher density of NCs of that species is formed which provide nucleation centers during the second stage of deposition. This results almost exclusively in two-component (possibly core-shell) NCs after the second deposition. If the more mobile species is deposited first, then a lower density of pure NCs of that species is formed which allows nucleation of new pure NCs of the less mobile species in the second stage of deposition. In addition, NCs formed in the first stage of deposition are mainly converted to two-component NCs. PIM has been implemented for these processes, also incorporating additional features described in Sec. V I.⁴⁴ PIM can also readily treat the case of simultaneous codeposition.

Naturally, PIM can be adapted to describe heterogeneous (rather than homogeneous) nucleation on surfaces where preexisting immobile defects provide nucleation centers. However, behavior for these processes is relatively trivial compared to homogeneous nucleation. It is natural to tessellate the surfaces into CZs surrounding each defect, where in the simplest case tessellation cells can be approximated by Voronoi cells. Then, it is clear that island growth rates are controlled by the corresponding CZ areas, so that the ISD should trivially match the CZ area distribution⁴⁵ (in contrast to the more complicated and subtle behavior for homogeneous nucleation).⁹ 

For deposition using an e-beam evaporator, some typically small fraction of the depositing species can be ions rather than neutral atoms. These ions can damage more fragile substrates, thereby creating nucleation centers during the deposition process.⁴⁶ Such nucleation centers will be created at constant rate and formed at random locations on the surface. Thus, the nucleation process is effectively described by a Johnson-Mehl-Avrami-Kolmogorov model¹⁰ potentially modified to account for a variable rate for radial island growth in this model.⁴⁶ PIM has been recently modified to treat this type of nucleation process.⁴⁶ 

For either homogeneous nucleation or heterogeneous nucleation at preexisting defects, there is effectively no control of nucleation locations. This shortcoming is addressed in directed assembly, where nucleation locations are guided by a templated substrate.⁴⁷ Consider a substrate which is periodically templated or modulated in the sense that some aspect of the adsorption energetics varies periodically on some coarse length scale, L_c (i.e., L_c is significantly larger than the surface lattice constant). Here, there is the potential to induce the formation of periodic arrays of NCs. PIM has been implemented for such assembly which can be thermodynamically directed by modulating the adsorption energy or kinetically directed by modulating the diffusion barrier (or both).⁴⁸ 

The above presentation is intended to highlight into the value and utility of PIMs for elucidation of fundamental issues in homogeneous diffusion-limited nucleation and growth (Secs. II and III); distinct behavior in the regime of attachment-limited nucleation and growth (Sec. IV); the effect of a diverse variety of modified mechanisms and features of nucleation and growth processes during deposition on surfaces (Sec. V).

We have emphasized model development and analysis rather than applications. However, PIM can be applied directly to elucidate and interpret behavior in a host of specific systems, some of which are indicated above. An expanded list of examples includes (a) homogenous nucleation and growth of metal NCs during deposition on metal, semiconductor, oxide substrates,^8,9 with particular utility for Volmer-Weber growth of 3D islands on weakly binding substrates; (b) nucleation inhibited by attachment barriers as observed in metal (111) homoepitaxial systems which exhibit enhanced long-range adatom interactions with a “repulsive ring” due to surface states;^38,49 inhibited attachment was also recently suggested for Fe deposition on graphene;⁵⁰ (c) nucleation and growth with strongly anisotropic diffusion as observed for homoepitaxy on dimer-row reconstructed Si(100) surfaces;⁵¹ (d) significant effects on nucleation of small mobile clusters as anticipated for homoepitaxy on metal (100) and (111) surfaces;⁵² (e) growth inhibition in strained-layer heteroepitaxy with large mismatch;⁵³ (f) codeposition to form bimetallic NCs, e.g., of Pt and Au on TiO₂(110),⁴³ and Pt and Ru on graphene;⁴⁴ (g) exchange-mediated nucleation in Fe on Cu(100),⁴¹ and Ni on Ag(111)⁴² systems; (h) nucleation of metal NCs on graphite is often facilitated by sputtering to create surface damage and heterogeneous nucleation centers, but even a small fraction of Cu ions from an e-beam evaporator can create sufficient damage that heterogeneous nucleation dominates for Cu deposition on HOPG;⁴⁶ (i) deposition of metals on metal-supported graphene, with periodically rumpled morié structure due to lattice mismatch, often results in directed assembly of 3D NCs,⁵⁴ and the PIM framework is ideally suited to modeling of this complex process.⁴⁸ 

Finally, all of the above discussion of PIM analysis pertains to nucleation and growth of NCs during deposition. We should mention that the PIM approach can be extended to treat post-deposition coarsening of NC arrays.⁵⁵ Coarsening can occur either by cluster diffusion and coalescence (where one would incorporate appropriate size-dependent diffusion coefficients for NCs) or via Ostwald ripening (where one would incorporate an appropriate size-dependent detachment rate of adatoms from NCs).

Y.H. and J.W.E. were supported for this work by NSF Grant Nos. CHE-1111500 and CHE-1507223. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358. Computations utilized USDOE NERSC, OLCF, and NSF-supported XSEDE resources. E.G. was granted access for this work to the HPC resources of GENCI (Grand Equipement National de Calcul Intensif) under the allocation 96339. T.J.O. acknowledges the support from CNPq and FAPEMIG (Brazilian agencies).

Rather than the standard MF expressions for diffusion-mediated nucleation and aggregation rates in Sec. II, an alternative more fundamental and flexible formulation is based on the lifetime, τ, of deposited adatoms.^9,11,33,34 Let N_t be the density of “traps” for a diffusing adatom, so N_t = N_isl + N_i ≈ N_isl in the steady state, where again N_isl is the density of stable islands and N_i is the density of critical clusters. Then, the mean number of hops, m_t, of a deposited atom before trapping through nucleation or aggregation satisfies m_t = hτ ∼ (N_t)⁻¹ln[(N_t)⁻¹] for isotropic diffusion in 2D, and m_t = hτ ∼ (N_t)⁻² for 1D diffusion. Then, the rate for aggregation with NCs of size s are given by K_agg(s) ≈ N₁q_s/τ where q_s = N_s/N_t is the probability to aggregate with an island of size s. Similarly, K_nuc = K_agg(i) and K_agg ≈ N₁/τ. Matching these results with MF expressions, we conclude that σ_s ∼ σ_i ∼ σ_av ∼ 1/ln[(N_t)⁻¹] for isotropic 2D diffusion, but σ_s ∼ σ_i ∼ σ_av ∼ N_t for 1D diffusion. The former result yields just log-corrections to the previously derived results for the 2D case. The latter yields fundamentally different scaling behaviors for 1D diffusion.

For the case of diffusion-limited homogeneous nucleation and growth (without additional barriers for nucleation and aggregation), we further quantify finite-size effects on the ISDs by analyzing the behavior of the first four cumulants of ISDs versus L/L_isl . Results in Fig. 13 show clear scaling with L/L_isl, but slower convergence to limiting behavior at L/L_isl ≈ 20 − 30 than observed for N_isl. This is reasonable since the ISD is expected to be more sensitive to finite-size effects than the mean density, N_isl.

Next, we expand the analysis in Sec. III B of the few-island regime where L/L_isl is of order unity. We have already assessed the probability for nucleation of the first and second islands as a function of coverage. To assess the existence of the third island, one should convolute the probability density dP₂/dθ′ = F⁻¹K_Tnuc1exp(−θ′K_Tnuc1/F) for the second island to be formed at coverage θ′ < θ with the probability 1 − exp(−δθK_Tnuc2/F) that the third island has formed after an additional coverage δθ = θ − θ′. Here, K_Tnuc2 is the time-independent total nucleation rate for the third island in the presence of two existing islands. Actually, K_Tnuc2 depends on the position of the second island relative to the first. Thus, one should really account for these island-position-dependent rates, but we ignore this complication. Evaluating the convolution integral, we obtain

Note that P₃ initially increases quadratically with θ (in contrast to P₂ which initially increases linearly). Below, we will exploit the feature that K_Tnuc2 has the same scaling as K_tnuc1 but is just reduced by a factor of p₂ < 1 since the steady-state N₁ is reduced on the more crowded surface. Setting K_Tnuc2 = p₂K_Tnuc1 with p₂ < 1. Then

where

The data in Fig. 5(a) with θ = 0.1 ML around L/L_isl ≈ 1 are best fit by choosing b = 1.1 and p₂ = 0.75.

As a final aside, we note that analysis of finite-size effects is particularly relevant for deposition on vicinal surfaces, where a transition from island formation on terraces to step flow occurs for terrace widths L_terr ∼ L_isl.⁵⁶ We expect that a similar crossover behavior occurs to that described above.