Capture zone area distributions for nucleation and growth of islands during submonolayer deposition

Homogeneous nucleation and growth of islands during submonolayer deposition on perfect low-index extended surfaces have been explored for decades as a route for self-assembly of arrays of supported nanoclusters or “islands.”^1–3 This process involves a competition between simultaneous nucleation and growth of stable islands mediated by terrace diffusion of deposited atoms. Nucleation is often characterized by specifying a critical size i (measured in atoms) above which islands are stable and effectively immobile on the timescale of island formation. Thus, i = 1 corresponds to irreversible island formation where dimers and larger clusters are stable. For i > 1, island formation is reversible as sub-stable clusters of i or less atoms can dissociate. For i = 0, one regards isolated diffusing adatoms as being able to spontaneously convert to stable nuclei (e.g., by exchange with the substrate) about which stable islands then grow.⁴ 

Theoretical and experimental analyses have focused on behavior of the mean island density, N_isl, and of the island size distribution (ISD).^1–3 However, there has also been a long-standing appreciation of the significance of capture zones (CZs) associated with individual islands.^5,6 The key idea underlying the CZ concept is that the surface can be tessellated into a distribution of CZs, one per island, where the area of the CZ for each island determines its growth rate.^3,5–11 Recently, attention has turned to analysis of the CZ area distribution (CZD).^12–21 A heuristic formulation by Pimpinelli and Einstein suggested that the CZD has a Gaussian tail for large areas, and that the small area behavior is related to the critical size, i, characterizing nucleation.¹⁴ 

Our strategy for analysis of these surface deposition systems is to develop exact beyond-mean-field evolution equations for the basic quantities of interest: the ISD, the CZD, or related quantities such as the joint probability distribution (JPD) for island sizes and CZ areas.^9,12,13,16 While exact analysis is not possible, these equations generally provide a solid foundation for reasonable approximate treatments. Also from this approach, one immediately recognizes the importance of characterizing spatial aspects of the nucleation process and of the resulting island distribution. This includes description of correlations between island size and island separation (or, equivalently, the dependence of CZ area on island size). Critically, here for analysis of the CZD, it is essential to provide a detailed description of various subtle spatial aspects of nucleation. Indeed, the central importance of the latter has been noted previously both in analytic theories based on the JPD,^12,13 and in the development of an efficient “geometry based simulation” (GBS) procedure to describe island formation.^22–24 Of particular significance is the feature that nucleation is enhanced by a higher adatom density. It therefore occurs preferentially far away from stable islands, the edges of which are sinks for diffusing adatoms, and thus nearby CZ boundaries.^3,22–24

In Sec. II, we describe our model for island formation, develop the fundamental equation for evolution of the CZD, and indicate expected scaling behavior of the solution. In Sec. III, we describe in detail key quantities characterizing the spatial aspects of nucleation. The scaling form of the fundamental evolution equation for CZD is described in Sec. IV, as well as the anticipated asymptotic form of its solutions for small and large CZ area. Scenarios for creating CZs with small areas are described in detail in Sec. V. Simulation results for i = 1 and i = 0 are presented and analyzed in Sec. VI. Further discussion and conclusions are presented in Sec. VII. A brief sketch of the basic formulation presented here and limited results were provided in a previous Comment.¹⁶ 

Our atomistic model for island formation during submonolayer deposition includes the following basic ingredients: (i) Random deposition at rate F per site on a square lattice of adsorption sites. (ii) Hopping of isolated adatoms to empty nearest-neighbor (NN) sites with rate h per direction. (iii) Aggregation of adatoms into two-dimensional (2D) islands, which are stable above a critical size, i (measured in atoms). Nucleation requires aggregation of i + 1 atoms for i ≥ 1. Nucleation is mediated by exchange of single adatoms with the substrate for i = 0. Sub-stable islands can be regarded as in quasi-equilibrium with isolated adatoms. (iv) Periphery diffusion is facile so that islands relax to achieve compact square shapes during growth. (v) Atoms deposited on top of islands are instantaneously transported to the island edge and hop down. See  Appendix A for further details.

We will sometimes comment on behavior for analogous “point-island” models where all islands occupy just a single lattice site but carry a label to indicate their size.²⁵ These models avoid some complications seen in the above more realistic models for which the “footprint” of the 2D compact islands expands during growth. Basic features of the ISD, CZD, and JPD are similar for point and compact islands at lower coverage, although we find that some quantities characterizing the details of nucleation differ significantly.

We define the coverage, θ, from the deposition time, t, via θ = Ft in monolayers (ML). Then, the theory below is applicable to the regime of lower coverage, θ ≤ 0.1 ML, where the typical linear dimension of islands is well below the mean island separation. Let N_s denote the density of islands of s adatoms, where densities are measured per adsorption site. Thus, N₁, N_i, and N_isl = ∑_s>i N_s denote the densities of adatoms, critical clusters, and stable islands with size s > i, respectively. The Walton relation for quasi-equilibrated critical clusters implies that N_i ∝ (N₁)ⁱ.^1,26 For i = 0, one naturally sets N₀ ≈ 1. For i ≥ 1, in terms of the nucleation rate, K_nuc ∝ hN₁N_i, and the rate of aggregation, K_agg ∝ hN₁N_isl, of adatoms with stable islands, one has that^1–3 

Integration of (1) for i ≥ 1 reveals a “transient regime” of rapidly growing N₁ for θ below a crossover value³ θ^∗ ∼ (h/F)^−2/(i+3) and a subsequent steady-state regime where adatom gain due to deposition roughly balances the loss due to aggregation so that F ≈ K_agg. The island density scales like N_isl^∗ ∼ (h/F)^−χ* with χ^∗ = (i + 1)/(i + 3) at crossover and like N_isl ∼ (h/F)^−χ with χ = i/(i + 2) in the steady-state regime.^1–3 Since χ^∗ − χ = 2(i + 2)⁻¹ (i + 3)⁻¹ > 0, nucleation in the steady-state regime dominates that in the transient regime for large enough h/F when i ≥ 1.³ Indeed, this dominance underlies selection of scaling forms for the ISD and CZD.

For i = 0, one must also specify an exchange probability, σ, per hop, so that the nucleation rate becomes K_nuc ∝ σhN₁.⁴ In this case, there are also transient and steady-state regimes, although the characterization of behavior differs.^3,27 Our version of this model is tailored to represent a “pure” i = 0 mechanism with no nucleation by aggregation of multiple diffusing adatoms. See  Appendix A.

Let N_A denote the density of stable islands with CZ area A, where areas will be measured in units of the surface unit cell area (and distances will be measured in units of surface lattice constant). Also let A_av denote the mean CZ area. Then, one has that

Fig. 1 illustrates schematically the terms in this equation. We shall see in Sec. III A that P_A and P⁺_A are related, and one has that ∑_AP_A = ∑_AP⁺_A = M, where M denotes the mean number CZs of pre-existing islands overlapped by the CZ of a just-nucleated island. This equality, together with the normalization condition, ∑_AP^∗_A = 1, on P^∗_A, ensures consistency of (3) after applying ∑_A to both sides, and noting the identity N_isl = ∑_AN_A. We remark that there are some subtleties in the definition of P_A, etc., for smaller A_av where A are strongly discrete, and thus there is a chance that the CZ of a just-nucleated island can overlap two or more existing CZs of the same size. These issues are briefly discussed in  Appendix B.

Note that the mean number of CZs overlapped by that of a just-nucleated island is M ≈ 5.5 for point-island models with i = 1,¹³ and we will find a similar value for compact islands. We emphasize that this result demonstrates the failure of a simple fragmentation picture where each nucleation event is regarded as simply fragmenting an existing CZ into two parts (corresponding to M = 1).²⁸ The reason for the high M-values is that, as noted above, nucleation tends to occur near existing CZ boundaries where overlap with multiple CZs is expected.^22–24 

Our focus will be on the asymptotic scaling regime of large A_av or small N_isl (corresponding to h/F ≫ 1 for i ≥ 1) where N_A is proposed to have the form^12–21 

Unless stated otherwise, ∫dα ranges over 0 ≤ α ≤ ∞. Further comment on two aspects of this scaling are appropriate. First, we should note that deviations from the limiting scaling form can be significant for h/F ∼ 10⁶-10⁷ considered here.^3,29,30 The true limiting or asymptotic shape, g, of the distribution may only be achieved for much higher h/F. Second, the formulation (4) leaves implicit any dependence of g on coverage, θ, given that in this study we just focus on a single θ = 0.1 ML. However, there certainly is a dependence on θ, although we do not expect this dependence to be strong for models incorporating compact islands in the regime for small θ (specifically in the regime θ ≤ 0.1 which we consider). See Ref. 30 for a detailed discussion and analysis of θ-dependence in the context of the ISD. Regarding the CZD, a weak θ-dependence of g will require that key quantities characterizing nucleation also have weak dependence on θ. In Secs. III and V, we will offer some comments related to the origin of the θ-dependence for compact islands. From these observations, it is possible that basic behavior of the CZD for small CZ areas could depend on θ, but this is not expected of behavior for large areas. As an aside, we remark that θ-dependence should strictly only be absent in asymptotic scaling functions for point-island models.^13,25,29,30

Our main goal is to provide insight into the form of the scaling function g(α) based on simulation studies, and also on analysis of (3). The complexity of (3) undoubtedly precludes simple forms for this solution, providing a significant caveat for most existing analyses of simulation data for g(α). However, just as for the ISD, there is a natural appeal in developing simple fitting-forms which can be conveniently used for experimental data analysis. In this context, generalizing a strategy used in the field of stochastic geometry for analysis of tessellations,^31–33 it is natural to consider fits to the generalized gamma (GG) distribution of the form¹⁶ 

where

The value of “b” is determined by the unit mean constraint, and that of “a” from normalization. The special case where ν = 2 corresponds to the so-called Generalized Wigner Surmise (GWS) choice in Ref. 14. Another significant component of Ref. 14 was to suggest that β is determined by the critical size i for nucleation, and specifically that β_GWS = i + 1¹⁴ (later changed to β_GWS = i + 2¹⁷) for island formation in a 2D surface deposition system. It should also be noted that traditional applications to analysis of tessellations in stochastic geometry usually set ν ≈ 1.^31–33 

Application of (3) to determine the CZ area distribution, N_A, or the associated scaling function, g(α), requires characterization of the quantities P⁺_A, P_A, and P^∗_A describing nucleation. For analysis of P⁺_A and P_A and the relationship between them, we introduce a more detailed quantity, P_A(A_sn), giving the probability that the CZ of a just-nucleated island overlaps that of a pre-existing island with CZ area A by an amount A_sn, where A_sn < A.³⁴ See Fig. 1. Then, it follows that

Summing both quantities over A immediately yields the identity $∑ A P + A = ∑ A P A$ (= M) mentioned in Sec. II B. See  Appendix B.

Next, we note that since P_A should be proportional to N_A, it is natural to write

It is also convenient to write P_A(A_sn) in the factorized form

We regard Q_A as an “intrinsic” overlap probability, and E_A(A_sn) as characterizing the distribution of overlap areas for CZs of area A. In addition, it will be appropriate to introduce the quantity A_subnuc(A) giving the mean area of the portion of the CZ of a just-nucleated island which overlaps that of a pre-existing island with CZ area A (or equivalently, the mean area of the portion of the existing CZ of area A which is overlapped by the new CZ). It then follows that

Averaging A_subnuc(A) over all CZ areas, A, yields the quantity

defining μ_av as the average fractional overlap with each existing CZ. In fact, previous analyses for point islands with i = 1 found that μ_av ≈ 0.17 and also indicated that A_subnuc(A) ≈ μ_avA for all A.^12,13 We will find a significantly lower value of μ_av for square islands at ∼0.1 ML.

Finally, it is also appropriate to introduce a quantity A_avnuc = ∑_AAP^∗_A corresponding to the average area for the CZs of just-nucleated islands. Note that α_avnuc = A_avnuc/A_av ≈ 0.97 for point islands with i = 1,¹³ but we will find a significantly smaller value for compact islands at ∼0.1 ML. It is tempting to equate A_avnuc with the product MA_avsubnuc = Mμ_avA_av,^12,13 i.e., to identify α_avnuc with Mμ_av ≈ 0.94 for point islands. However, this is only a (reasonable) approximation as it neglects the correlation between the number of CZs overlapped by a new CZ and the area of the portion of the new CZ overlapping each individual existing CZ.

As indicated in Sec. II B, we will focus on scaling behavior, N_A ≈ (N_isl/A_av) g(A/A_av), of the CZD for large A_av. Correspondingly, for the key quantities described above, we introduce scaling functions

and

and

The simplest anticipated form of the scaling function in (12) for the intrinsic overlap probability, q(α), might be q(α) ≈ c αⁿ. The increase in q(α) with α reflects enhanced nucleation near larger CZs where the adatom density is generally higher. For i = 1, a simplistic analytic estimate might suggest n ≈ 3 (cf. Sec. V), but previous limited simulation data for point islands indicated a value of roughly n ≈ 1.5,¹³ and more complex forms with effective n varying with α were also suggested. In (13), e(μ) represents the scaling function for the distribution of fractional overlaps μ = A_sn/A.

Applying the above results to determine the key quantities in the fundamental evolution Equation (3), one obtains

In  Appendix B, we confirm that these expressions preserve the equality ∑_AP⁺_A = ∑_AP_A. Since A_av is large, one can replace the sum ∑_Asn in (15) with an integral ∫dA_sn. After making the change of variable to the fractional overlap μ_sn = A_sn/(A + A_sn), so that (A + A_sn)/A dμ_sn = dA_sn/(A + A_sn), one obtains

where 〈…〉_sn = ∫μ_sne(μ_sn)... denotes a weighted average over fractional overlaps integrating over 0 ≤ μ_sn ≤ 1.

Given the somewhat obscure nature of the average in (16), it is instructive to consider a simple, although very crude approximation which greatly simplifies its form. If one assumes that the overlap area A_sn with a pre-existing CZ of area A is always close to A_subnuc(A) ≈ μ_avA, then the distribution e(μ) is sharply peaked around μ = μ_av. In the simplest approximation, one could set e(μ) ≈ δ(μ − μ_av), where δ is the Dirac delta function. Then, (16) reduces to

which imparts a simpler but non-local nature to the fundamental evolution equation analogous to that seen in evolution equations for fragmentation and coalescence processes. It is also readily confirmed that (17) preserves the equality ∑_AP⁺_A = ∑_AP_A. However, again we emphasize the crude nature of the approximation yielding (17).

Given the forms in Sec. III of various quantities describing nucleation, we recast (3) as an equation for the scaling function, g(α), of the CZ area distribution. From (4) and using that N_isl = 1/A_av, one obtains

Then substituting (18) and the scaling relations from Sec. III B for P_A, P⁺_A, and P^∗_A into (3), we obtain the scaled evolution equation

Again 〈…〉_sn = ∫dμ_sne(μ_sn)... denotes a weighted average over fractional overlaps μ_sn = α_sn/(α + α_sn) which can be simplified using (17), but this crude approximation is not pursued here. As a consistency check, integration ∫dα of the left-hand side (LHS) of (19) using integration by parts yields ∫dαg(α) = 1. See  Appendix B. Integrating the right-hand side (RHS), the first and second terms cancel exactly, and the last term yields ∫dαp^∗(α) = 1 matching the LHS.

Of course, the solution g(α) to (19) depends on the detailed form of the scaling functions for the intrinsic overlap probability, q(α), the CZ areas of just-nucleated islands, p^∗(α), and the distribution of fractional overlaps, e(μ). The complex form of (19) precludes a simple form such as GWS or GG for g(α). This is analogous to the situation for the scaling function for the ISD where the complex form of its exact evolution equation also strictly precludes commonly adopted simple forms.³ 

Despite the difficulty of exact analysis for (19), it is instructive to explore some approximations which might provide insight into asymptotic behavior of the solution g(α). First, for large α, one might anticipate that the major impact on large CZs is reduction in their size due to overlap by CZs of just-nucleated islands. (Such large CZs are perhaps created early at the end of the transient regime, and surviving into the steady-rate regime.) In this case, it is plausible that the first loss term on the RHS of (19) dominates the gain terms on the RHS, so that

where the integrals are from 0 to α and we utilize the simplified relation q(α) ∼ αⁿ introduced in Sec. III B and neglect any α-dependent prefactor. We emphasize the various assumptions or approximations leading to (21): a dominant loss term on the RHS of (19); the simplified form for q(α) and the approximation of the integral in the exponential. Nonetheless, this analysis suggests that ν ≈ n in (5) relating behavior of the large area tail of g(α) to details of the nucleation process.

Given uncertainties and approximations in the above analysis, it is appropriate to comment on and compare with the quite different heuristic analysis of Pimpinelli and Einstein.¹⁴ This approach anticipates more universal asymptotic behavior for the stochastic process of island and CZ formation, independent of the details of nucleation. Specifically, the CZ areas are assumed to satisfy a strongly damped Langevin equation subject to an effective force field. For large CZ areas, there is an effective “restoring force” which is roughly linear in area and has the effect of reducing area. This force is regarded as reflecting an effective pressure from surrounding CZs. From the corresponding Fokker-Planck equation, one concludes that the CZD has a Gaussian tail for large areas so that ν = 2.³⁵ Since neither (20)-(21) or this heuristic Langevin formulation is rigorous, there remains some uncertainty.

Second, for small α, given that q(α) → 0 as α → 0, the first term on the RHS should be negligible compared to other terms. For simplicity, let us initially also assume that the third p*-term dominates the RHS, i.e., that “direct creation” of new CZs with small area dominates (but see below). Then, terms on the LHS of (19) must balance the p*-term yielding

If p^∗(α) ∼ α^β, it follows that also g(α) ∼ α^β as α → 0, and specifically that¹⁶ 

This analysis is consistent with assumption in previous analysis¹⁴ that g(α) for small α is determined by the form of the distribution of new CZ areas, p^∗(α), for small α.

However, our picture for the formation of small CZs is quite different from Ref. 14 (see Sec. V), and furthermore is actually more complicated than suggested above. New CZs with small area can also be created “indirectly” by the CZs of just-nucleated islands mostly overlapping those of existing islands with typical areas. This corresponds to the second term on the RHS of (18) with small α and α_sn = O(1), and thus does not involve p^∗(α) for small α. More generally, suppose that the overall rate of creation of new CZs with small area α, combining contributions from the second and third terms on the RHS of (19), has the scaling form p_new(α) ∼ α^β. Then, it still follows from (19)) that g(α) ∼ α^β, as α → 0.

The scenario considered by Pimpinelli and Einstein for creating small CZs involves nucleation of new islands at a location on the surface surrounded by multiple pre-existing CZs having small size.^14,17 However, such regions are extremely rare, their population scaling like a higher power of the population of small CZs, a feature not analyzed in these previous studies.³⁶ Instead, we claim the dominant pathway for creating small CZs does not involve pre-existing small CZs, i.e., that new small CZs are primarily created in the absence of pre-existing small CZs.¹⁶ Below, we present and analyze various key scenarios, although a comprehensive analysis of all possibilities is difficult.

One scenario is where a small CZ is created by nucleation of an island in between the sides of two (m = 2) closely spaced near-square islands. See Fig. 2(a). We note that the likelihood of this scenario is dependent on the structure of the islands. Here, the finite extent of islands is key in facilitating creation of a small CZ between a nearby existing pair. Also higher coverage, θ, will facilitate this scenario. Another scenario is nucleation of an island in the middle of a triangle of three (m = 3) closely spaced islands. See Fig. 2(b). In both these cases, the small CZ is associated with the just-nucleated island, and thus with P^∗_A or p^∗(α) described in Sec. III. Of course, one could also consider nucleation in the middle of larger collections of nearby islands (m > 3). However, analysis below indicates that the likelihood of such events is negligible. A semi-quantitative analysis of the p^∗(α) associated with these scenarios has two components. First, we must determine the small population, P_m, of the above “rare” CZ configurations. Second, we must determine the probability, P_nuc, that nucleation occurs in the appropriate restricted sub-region as needed to create a desired small CZ.

For the first part of the analysis, consider the density, N_isl(r), of islands with their centers separated from the center of a specified island by a vector separation r. One has that N_isl(r) ≈ C_isl(r) N_isl, where the island-island correlation function satisfies C_isl(r) ∼ [N₁(r)/N₁]ⁱ⁺¹, and N₁(r) represents the adatom density at position r from the center of an island.^1,3 N₁(r) decreases linearly to zero at the island edge,^1,3 so we expect that N₁(r)/N₁ ∼ r/r_isl for closely spaced islands, where r = |r| and again r_isl is the mean island separation. We conclude that the relative probability to find a pair (m = 2) or trio (m = 3) of closely spaced islands with separations r ≪r_isl should scale like

The factor of (C_isl)^m−1 gives the probability for a specific orientation of the m-tuple of islands, and the additional factor of (r/r_isl) accounts for the volume of phase-space for arbitrary rotations.

The second factor is the relative probability for nucleation in an appropriate restricted subregion, as indicated above, separated by island edges by a distance of order r. First, we note that the total rate of nucleation, K_CZ, for a typical CZ scales like

which follows from the steady-state results from Sec. II B. This type of scaling of nucleation rates with linear size was a key component of previous geometry-based simulations of island nucleation and growth.^22–24 Next, we consider the total nucleation rate for the above-mentioned restricted sub-region of area ∼r² denoted by K_SUB. Let N₁(loc) denote the typical value of the adatom density in this subregion. Then, it follows that K_SUB scales like

Thus, one has that

Combining the results (24) and (27), and using that α ∼ (r/r_isl)² ≪ 1, the likelihood to create a new small CZ by the above mechanisms scales like

Thus, one has β₂ = 3i/2 + 3 and β₃ = 2i + 7/2. For compact islands with finite extent, one expects that contributions from m = 2 (i.e., the smaller β_m) will dominate.

Another scenario can generate island and CZ configurations similar to those described in Sec. V A, but with a different order of island nucleation. This scenario will be most viable for m = 3, i.e., a triangular triple of nearby islands surrounding a fourth, rather than for m = 2. The idea for m = 3 is that islands could first nucleate creating two corners of the triangle, then a third island nucleates which will end up in the middle, and finally a fourth nucleates constituting third corner. See Fig. 3. In this case, the CZ of the last just-nucleated island is not generally small. However, the CZ of the third nucleated island, which was generally large upon formation, is greatly reduced (i.e., mostly covered) upon nucleation of the fourth island. In this case, the formation of a small CZ is associated with the term P⁺_A in dN_A/dt. Specifically, it reflects the small but non-zero weight in the upper end of the fractional overlap distribution, e(μ), with μ close to unity.

For a semi-quantitative analysis of the scenario in Fig. 3, the probability to nucleate the first three islands in the above configurations should scale as P′₃ ∼ (r/r_isl)²ⁱ⁺³ just as does P₃. However, the last “outer” island will nucleate in a region where N₁(r)/N₁ ∼ r/r_isl noting the linear increase of adatom density moving away from the first three islands. This suggests a relative nucleation probability of P′_nuc ∼ (r/r_isl)² [N₁(r)/N₁]ⁱ⁺¹ ∼ (r/r_isl)ⁱ⁺³. Using that α ∼ (r/r_isl)², the overall rate for formation of the small CZ scaling like

(versus β₃(i) = 2i + 7/2 in Sec. V A). Since β′₃ has the same value as β₂ given in Sec. V B, this process should also make a significant contribution to the measured β. In fact for point islands, where the m = 2 pathway of Sec. V A is not effective, one expects this alternative pathway for m = 3 to dominate.

Efficient KMC simulation for large values of h/F utilizes a standard Bortz-algorithm to enable collection of extensive statistics on island configurations. To assess quantities such as g(α), CZ tessellations are constructed which here are taken as edge cells, so that CZ boundaries are equidistant from island edges.^{10,11,22–24} However, to collect adequate statistics at a selected coverage for quantities associated with new islands such as p^∗(α), the following refined approach is implemented:^9,12,13 (i) we perform a standard simulation to generate an island distribution at the desired coverage; (ii) we continue simulation until one new island is nucleated and measure desired quantities associated with this new island; (iii) we remove the new island and continue simulation until another new island is nucleated and perform desired measurements; (iv) we repeat this process for ∼100 nucleations per configuration generated in (i). This approach is far more efficient than collecting information on just one new island from each simulated configuration.

We focus below on analysis of data from simulations with i = 1 and h/F = 10⁶ at θ = 0.1 ML. The total number of islands and thus CZs used to assess g(α) is ∼268 500, with average area A_av = 372.4. With regard to detailed analysis of quantities such as p(α), p^∗(α), e(μ), etc., we consider a total of 50 000 new CZs which have an average area of A_avnuc = 201.6. Key parameters extracted from this analysis include α_avnuc = A_avnuc/A_av = 0.54, μ_av = 0.099, and M = 4.65. Note that the α_avnuc and μ_av values are significantly lower than for point islands. This feature is a natural consequence of the finite extent of the square islands at 0.1 ML reducing the relevant full or partial CZ areas. Also, as for point islands, α_avnuc is somewhat above Mμ_av ≈ 0.46.

We start by performing least-square fits of GG forms to g(α).³⁷ See Fig. 4. First, this is done by selecting the reasonable choices of n-value and fitting β. For h/F = 10⁶, this yields β = 3.05 for n = 2, β = 3.48 for n = 1.75, and β = 4.07 for n = 1.5. One can also simultaneously adjust both β and n. For h/F = 10⁶, this yields β = 3.84 and n = 1.59. As an aside, from more limited simulation data, we obtain β = 3.55 and n = 1.65 for h/F = 10⁷, and β = 3.36 and n = 1.70 for h/F = 10⁸.³⁸ 

The advantage of above fitting approach, which assumes that g(α) is described by a GG form, is that one utilizes simulation data for all α. These data are extensive for CZ areas around the average value (α ∼ 1), thereby compensating for limited data in the upper and lower tails of the distribution. However, there are also significant disadvantages. First, the fit is expected to be sensitive to deviations for finite h/F from asymptotic behavior as h/F → ∞. Variation with h/F in the peak if the CZD has been assessed previously and significantly impacts estimates of β.¹⁵ Second, a more fundamental concern is that g(α) will not be exactly described by a GG form (as indicated above). This feature can introduce some bias and correlation in the resulting values of β and n (e.g., higher n produces smaller β).

The above observations motivate direct unbiased fitting of the lower tail of g(α) to assess β, and thus the proposed relationship to i. Here, one must choose the range of α over which to fit. Lower (and narrower) ranges better describe asymptotic behavior as α → 0, but these suffer from limited data. Thus, we give β-values from fitting over different ranges with results shown in Table I. The decrease in estimated β with broader range of α just reflects non-asymptotic behavior as log g(α) versus log α is convex. Thus, extrapolating values for broader ranges to narrower ranges should give the most reliable estimates. Results from analysis of extensive data for h/F = 10⁶, and more limited data for higher h/F, consistently suggest that β is around 3.5.

As indicated in the theoretical development of Secs. II–V g(α) behavior is controlled by subtle spatial aspects of the nucleation process. Thus, it is appropriate to also examine simulation results for the key scaling functions p^∗(α), p(α), q(α), and e(μ). See Fig. 5. Consistent with the (over)simplified analysis of (22) and (23), the distribution, p^∗(α), of areas of new CZs appears to scale in the same way as g(α) for small α. Also, g(α) is naturally smaller than p^∗(α) in the small-α range. However, direct creation of new small CZs associated with p^∗(α), as described in Sec. V A and as assumed in (23), is not necessarily the dominant pathway. Rather, examination of simulation movies reveals a significant fraction of cases where such small CZs are formed indirectly from a substantial reduction in area of existing CZs, as described in Sec. V B. See Fig. 6. While the associated values of the fractional overlap scaling function, e(μ), for larger μ are small, so are values of p^∗(α) for correspondingly small α ∼ 1 − μ. (We should also note that e(μ) is far from the delta function form corresponding to a narrow distribution of fractional overlaps where typically μ ≈ μ_av. Instead, small fractional overlaps are common.)

With regard to behavior of g(α) for large α, as anticipated, the intrinsic overlap probability, q(α), increases non-linearly with α. The precise form is not clear from the limited data, but a simple fit to q(α) ∼ αⁿ yields n ≈ 1.5 with considerable uncertainty.

For the i = 0 model, the island density scales like N_isl ∼ σ^1/2 independent of h/F.^3,4 Most studies considered the regime of “rapid conversion”^4,19 with σ ≫ (h/F)^−1/2 yielding persistent nucleation in the steady state and a monotonically decreasing ISD.^3,4 A crossover occurs reducing σ to σ = (h/F)^−1/2 and below to growth behavior avoiding persistent nucleation. Thus, in this work, we select σ = (h/F)^−1/2 which leads to relatively low N_isl ∼ (h/F)^−1/4. Our construction and analysis of CZs includes single converted adatoms as well as larger clusters, but not isolated diffusing adatoms. Even at a low θ ∼ 0.1 ML, there is significantly more island coalescence for i = 0 than for i = 1 consistent with the larger island-island correlation function, C_isl(r) ∼ (r/r_isl)ⁱ⁺¹, for small r. Since our theory is geared to the pre-coalescence regime, we also consider behavior for lower θ = 0.04 ML to limit coalescence.

We present a more limited analysis of behavior for i = 0 than for i = 1 focusing on the case h/F = 10⁶ and σ = 10⁻³ for θ = 0.1 ML. The total number of islands and thus CZs used to assess g(α) is ∼346 500 with an average CZ area A_av = 192.71. Least-squares fitting of a GG form adjusting both β and n yields β = 2.80 and n = 1.42 for θ = 0.1 ML. One obtains β = 3.46 and n = 1.42 for θ = 0.04 ML. Given the limitations or bias in least square fitting to GG forms as discussed for i = 1, we do not give a more extensive analysis here. The reader is referred to Fig. 7 where GG forms are compared with simulation data for 0.1 ML for various choices of β and n. Choosing n = 2 requires β values above 2.5 for a reasonable fit. It appears that the data are better fit by lower n values and higher β closer to 3 which is consistent with β ≈ β₂(0) = 3 for i = 0 in (27). As for i = 1, direct fitting of the lower tail of g(α) to assess β avoids bias of an assumed GG form, and is thus most reliable provided there is sufficient simulation data. Results are shown in Table II which suggest again that β is just below 3.0.

Finally, we remark that analysis of simulation movies suggests that as for i = 1, a significant fraction of small CZs is created indirectly rather than directly. See Fig. 8.

Recent extensive interest in the capture zone area distribution (CZD) has been prompted by the proposal of Pimpinelli and Einstein (PE) that the CZD has a GWS form, and furthermore that the exponent, β, describing behavior for small areas is directly related to the critical size.^14,17,21 Considering 2D systems for i ≥ 1, Shi, Shim, and Amar (SSA) utilized point-island models to critique the original proposal of PE that β_WS = i + 1 = 2 for i = 1.¹⁵ Their simulation data for typical h/F = 10⁷ fitted to a GWS form instead suggested that β ≈ 3. Furthermore, SSA demonstrated that key features of the scaled CZD, particularly the peak height, had a significant dependence on h/F. Also, slightly different versions of the point-island model had significantly different peak heights for smaller h/F. However, extrapolation to h/F → ∞ removed discrepancies between different versions of the model, as might be expected. This analysis again indicated that β ≈ 3 assuming that the CZD, and particularly its peak, is described by a GWS form. This analysis appears consistent with the modified proposal of PE that β_WS = i + 2 (= 3 for i = 1).¹⁷ However, as discussed above, there are limitations to an analysis based on fitting to a GWS (or GG) form which does not correspond to the exact form of the CZD.

Additional studies by Oliveira and Aarão Reis^18–20 considered the CZD for point, square, and fractal island models utilizing a modified scaling formulation to reduce dependence of the scaled CZD upon h/F. However, this modified formulation was not designed to highlight or analyze behavior for small CZ areas. Point-island behavior for i = 1 is fit reasonably but not perfectly by β ≈ 3 and n = 2, or by β ≈ 4 and n = 1.5. Square island behavior is also fit reasonably by β ≈ 3 and n = 2, but with significant discrepancies in the tails of the CZD. A key focus of these studies was on the form of the right tail of the CZD for large areas. For i ≥ 1, simulation data for point islands and also for square islands at pre-coalescence coverages appear consistent with a Gaussian tail, i.e., with n = 2, or at least with n close to but probably not above 2. For square islands at higher coverages where coalescence is significant, there was a clear deviation from Gaussian towards exponential tails for large CZ areas.

Our own analysis for 2D systems with i ≥ 1 is motivated by our development of an exact evolution equations for the CZD. This equation suggests (but does not rigorously prove) the possibility of a non-Gaussian tail for large areas. It also provides and incorporates a different picture from PE of behavior for small areas. Specifically for the latter, appropriately identifying and analyzing the dominant pathways for creation of small CZ areas (which in our picture do not require pre-existing small CZs), we predict larger β-values than the refined PE prediction of β_WS = i + 2. Furthermore, we argue that the most reliable assessment of β comes from direct fitting of CZD data for small α. This contrasts the strategy of previous studies of fitting to a GWS or GG form (which can bias the outcome since the CZD is not exactly described by these forms). Such direct fitting does indicate somewhat larger values of β than β_WS = i + 2 = 3 for i = 1, although the values are below those for the pathways which we identified for creation of new small CZs. We note one limitation of our analysis is that data are obtained for relatively small h/F (as in analyses of PE). Data for higher h/F are preferred to avoid possible discrete lattice effects for small α, but its acquisition is computationally more demanding. With regard to behavior for large areas, our simulation data are not sufficiently extensive to make a definitive assessment. Behavior is consistent with Gaussian tails (n = 2) but could also be described well by somewhat smaller n for i = 1. Also taking into consideration the analysis of Refs. 18 and 19, one concludes that n is close to 2 (and that it does not exceed 2).

Considering 2D systems for i = 0, our analysis indicates β values significantly above the refined PE prediction of β_WS = i + 2 = 2 for i = 0. It also seems likely that the large area tail of the CZD is described by n < 2 rather than by a Gaussian tail. Previous studies for i = 0 with larger exchange rates, σ > > (h/F)^−1/2, and thus persistent nucleation, indicated exponential tails for large areas.¹⁹ This is certainly reasonable since persistent random nucleation produces quasi-random island distribution for which one expects a CZD closer to that of a Poisson point process. For the latter, a Gamma function form with exponential tail is expected.^29–31 However, for our choice of lower exchange rate σ = (h/F)^−1/2 which eliminates such persistent nucleation,³ the form of the large area tail is less clear a priori.

In summary, our exact evolution equation for the CZD provides a valuable additional formulation and perspective beyond a heuristic treatment for analysis and understanding of this distribution. In particular, our formulation highlights the importance of subtle spatial details of the nucleation process, and specifically of related quantities such as the intrinsic overlap probabilities and fractional overlap distribution. These quantities have received limited attention and are not considered in heuristic treatments, but their detailed form impacts that of the CZD. Given the complexity of the process of homogeneous island nucleation and growth of islands during deposition, there does not yet exist a definitive characterization and understanding of either the ISD or the CZD. However, this goal is facilitated with the aid of detailed analytic theory such as that developed here, together with more extensive simulation targeting key quantities characterizing the spatial aspects of nucleation. Another challenge is that extensive simulation data are most readily obtained for smaller h/F = 10⁶-10⁷, but one is primarily interested in the scaling limit h/F → ∞, so reliable extrapolation is also required.²⁹ 

We thank Tiago Oliveira for valuable comments on the manuscript. Y.H. and J.W.E. were supported for this work by NSF Grant Nos. CHE-1111500 and CHE-1507223. Computations utilized NSF-supported XSEDE resources. M.L. was supported by NSF of China under Grant No. 51271197. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358.

Our tailored models for island formation with i = 1 involve: (i) Random deposition at rate F on a square lattice of adsorption sites with periodic boundary conditions. (ii) Isolated adatoms hop to a NN site with hop rate h per direction. (iii) Typically, singly coordinated “edge atoms” can hop at rate h to either NN or next-NN edge sites provided they are still coordinated (i.e., they cannot detach from islands). (iv) Multiply coordinated atoms are immobile. (v) Atoms deposited on top of islands immediately undergo a random walk until they reach the edge of the island and hop down to the substrate. Also, atoms landing on isolated adatoms immediately hop down.

Incorporation of facile diffusion of atoms along the island periphery both along straight edges and around corners or kinks suffices to achieve compact islands.³ Implementation of rapid transport of atoms landing on top of islands to the first layer avoids second layer population.

One extra constraint not described above is that singly coordinated adatoms in clusters of size s = 2 or 3 are not allowed to undergo hops to next-NN edge sites. This constraint excludes long-range diffusion of these small clusters, which if present would have the effect of changing the scaling behavior of N_isl.³ For clusters of 4 or more atoms, restructuring will produce a rectangular core of doubly coordinated adatoms which is permanently frozen, thus excluding long-range diffusion.

Our version of a “pure” i = 0 model excludes nucleation by aggregation of multiple diffusing adatoms. This is achieved by preventing isolated diffusing adatoms from reaching NN sites, and also exclude deposition on top of or next to isolated adatoms. Island nucleation occurs exclusively by isolated adatoms converting to permanently immobile adatoms at rate σh. Diffusing adatoms aggregate irreversibly with converted immobile adatoms and with clusters of adatoms which form around these, i.e., detachment from such islands is prohibited. Edge diffusion is included as for i = 1 ensuring the formation of compact islands. Atoms deposited on top of islands are also treated as for i = 1. In our construction and analysis of CZs for i = 0, we will count as islands single converted adatoms in addition to larger clusters, but not isolated diffusing adatoms.

We will not present results for i > 1, however we mention that the model for i = 1 can be naturally extended to i = 2 or 3. For example, for i = 2, singly coordinated edge atoms can only detach at rate h if their unique neighbor is also singly coordinated (i.e., if they are part of an ad-dimer).²² Otherwise detachment is excluded. Doubly coordinated adatoms are again frozen. In these models with i ≥ 1, edge diffusion, as described above, ensures the formation of compact islands. For KMC images of compact islands in such models with i = 1-3, see Refs. 23 and 24.

First, we comment on the regime of smaller A_av so that CZ areas A are discrete. Then, there is a significant probability that the CZ of a just-nucleated island can overlap k = 2 (or more) pre-existing CZs of the same size, A. Then, in determining P_A, one considers a large number, j = 1 − N, of nucleation events and counts the number k_j of overlaps in each case, setting P_A ≈ ∑_j k_j/∑_j 1 = ∑_j k_j/N. A similar treatment is implemented for related quantities P⁺_A and P_A(A_sn).

Second, we show that $∑ A P A + = ∑ A P A$ exploiting a change of variable B = A + A_sn,

The same change of variable B = A + A_sn applies to simplify

Third, we comment on the consistency of the evolution Equation (19) for the scaling function which follows from integrating over α. After integration, the second term on the RHS of (19) yields

which cancels the first term given the normalization of e(μ_sn). Finally, we remark that for the most simple δ-function choice of e(μ_sn) ≈ δ(μ_sn − μ_av) (which again has limitations), (19) reduces to the non-local equation