Communication: Diverse nanoscale cluster dynamics: Diffusion of 2D epitaxial clusters

Nanoscale dynamics of crystalline clusters,^1–4 epitaxial structures,^5–7 and even liquid jets⁸ can differ qualitatively from mesoscale or macroscale behavior described by continuum formalisms. This applies for reshaping of 3D crystalline nanoclusters^1–3 and 2D epitaxial nanoclusters.⁶ Differences occur when linear nanocluster sizes, L, decrease below various characteristic lengths.^1–3,9–12 For L < L_k, the characteristic separation between kinks on close-packed step edges, 2D clusters become effectively faceted and reshaping is controlled by inhibited nucleation of new edges.⁹ Similar anomalies have been discussed for 3D nanoclusters, but here facets can be macroscopic and nucleation of new layers becomes prohibitive above nanoscale.¹³ Also relevant are characteristic Ehrlich-Schwoebel type lengths associated with additional barriers to round kinks and corners in 2D and to cross step edges and move between facets in 3D.^10,14,15 In general, key exponents describing both temporal scaling and size-scaling of relaxation times differ from macroscale theory.^1–3,5,6 Note also that reshaping of 2D epitaxial clusters and pits, or 3D clusters and voids in bulk crystals, is equivalent on the macroscale,¹⁶ but not on the nanoscale.^11,12 The size-dependence of diffusivity of supported metal nanoclusters was also found to deviate from macroscale predictions.^17–22 Refined continuum treatments with multiple order parameters might describe such anomalous scaling of diffusion.²³ However, the current study reveals a lack of recognition of the full diversity of possible nanoscale dynamics, which undoubtedly cannot be captured even with refined continuum treatments.

To explore this issue, we provide a comprehensive analysis of a canonical model for diffusion of 2D epitaxial clusters on metal(100) surfaces mediated by periphery diffusion.¹⁰ Cluster diffusion impacts kinetics of coarsening of arrays of supported nanoclusters via Smoluchowski ripening, i.e., cluster diffusion and coalescence.^5,21,24–28 Thus, elucidation of coarsening kinetics has importance for stabilizing functional nanomaterials, e.g., in limiting catalyst degradation. In our stochastic lattice-gas model, clusters of N adatoms on a square lattice of adsorption sites (lattice constant a) interact with just the nearest-neighbor (NN) attractions of strength ϕ > 0. Edge atoms can hop to NN and also to the second NN empty sites, provided that hopping retains at least one NN adatom in the cluster. This preserves cluster connectivity and size. All hop rates have the Arrhenius form h = ν exp[−E_act/(k_BT)] for surface temperature T and Boltzmann constant k_B. In our simplest model, the activation barrier, E_act, satisfies E_act = E_e + (n_NN − 1)ϕ, where n_NN denotes the number of in-plane NN adatoms prior to hopping. Thus, E_act = E_e applies for diffusion of isolated edge atoms along close-packed steps and around kinks (or corners), which occurs at rate h_e = ν exp[−E_e/(k_BT)]. For reference, choosing E_e = 0.29 eV and ν = 10^12.5 s⁻¹ mimicking Ag/Ag(100) yields h_e = 10^7.6 s⁻¹ at 300 K. The above formulation of E_act could be refined to include an additional kink rounding barrier, δ, for second NN hops around kinks, but we concentrate on δ = 0.

Model behavior is precisely determined by Kinetic Monte Carlo (KMC) simulation. Our focus is on the cluster diffusion coefficient, D_N = lim _δt→∞ D_N(δt), where D_N(δt) = ⟨[$δr$(δt)]²⟩/(4δt). Here, $δr$(δt) is the displacement in the cluster center-of-mass (CM) in a time interval δt, and ⟨⟩ denotes an average over a CM trajectory. See Fig. 1. Our model has D_N(δt) ∝ a²h_e so that D_N/(a²h_e) is independent of E_e and ν. Due to backward correlations in the walk of the CM, D_N(δt) is not constant but decays to a plateau value for sufficiently long δt > δt_c, where ⟨[$δr$(δt_c)]²⟩ ∼ a². Thus, precise determination of D_N from ⟨[$δr$(δt)]²⟩ for δt ≥ δt_c requires a CM trajectory with length, t_tot, of at least O(10³ δt_c). We choose t_tot ∼ 35 000δt_c. Macroscale theory predicts simple monotonic size scaling of the diffusion coefficient, D_N ∼ σ_PD N^−β, with β = 3/2 where the step mobility σ_PD has Arrhenius energy²⁹ E_PD = E_e + ϕ, but we find distinct behavior. Insight beyond KMC simulation comes from analytic treatments. If Ω_N denotes the number of distinct configurations of a cluster of size N, then D_N follows exactly from analysis of the “acoustic” eigenmode of an Ω_N × Ω_N matrix encoding allowed transitions between cluster configurations. This matrix is extracted by Fourier transformation of the linear master equations.^30,31 However, this analysis is only viable for small N, as Ω_N grows rapidly with N. Thus, we instead apply alternative combinatorial analyses.

Prior to presenting KMC results, we characterize nucleation-mediated versus “facile” cluster diffusion for moderate sizes with N ≥ 9. We identify “perfect” sizes N_p = L² and L(L+1) with L = 3,4,… for which clusters have unique square or near-square ground state shapes. This uniqueness does not apply for sizes N = L(L + n) with n ≥ 2. For perfect sizes, diffusion is nucleation-mediated in the sense that after an atom is extracted from a corner in the ground state to an edge (raising the energy by ΔE = +ϕ), another atom must quickly detach from a corner or kink to join the first atom before the first atom returns to the corner. Such a pair of atoms is regarded as nucleating a new edge. The most direct pathway to create another perfect configuration with displaced CM shifts atoms from kinks and corners of the opposite edge to complete this new edge. Shifting the second and subsequent atoms from one kink to another does not change the energy after each reattachment, so the system evolves through a series of first excited state configurations with energy ΔE = +ϕ above the perfect configuration. Only when the last (isolated) edge atom is shifted to recover the perfect configuration is the energy lowered by ΔE = −ϕ. See Fig. 2(a). Diffusion of clusters with sizes N_p+3, N_p+4, etc., is also nucleation-mediated requiring formation of a dimer on an outer edge to facilitate long-range diffusion of the cluster CM. Nucleation-mediated diffusion involves an atom breaking out of a kink (or corner) site with rate h_k = ν exp[−(E_e + ϕ)/(k_BT)] and aggregating with an isolated edge atom which has low quasi-equilibrium density n_eq = exp[−βϕ/(k_BT)]. Thus, D_N ∼ n_eq h_k ∼ exp[−E_eff/(k_BT)] with effective barrier^9,10,18 E_eff = E_e + 2ϕ.

For “facile” sizes N_p+1 and N_p+2, edge nucleation is not necessary. For N_p+1, an isolated edge adatom on a perfect core (a “special” ground state configuration) can readily diffuse around the cluster. For N_p+2, a NN pair of edge atoms or edge dimers on a perfect core (a ground state) can dissociate and reform on another edge. Neither process results in a net change of energy. After the isolated edge atom or dimer is transferred to a new edge, atoms can be transferred from the opposite edge of the core to complete the edge to which the isolated adatom or dimer was moved. This again leaves an isolated adatom or dimer on the edge of a displaced perfect core. Shifting of atoms from one kink to another does not change the energy after reattachment, and in this case the system evolves through a set of configurations iso-energetic with the initial ground state configuration. See Fig. 2(b). Facile diffusion just involves breaking atoms out of kink sites and subsequent edge diffusion so that E_eff = E_e + ϕ.

In addition to the above direct pathways for long-range diffusion, there are numerous less efficient indirect pathways with atoms removed from multiple corners. Here, the cluster wanders through a large phase space of configurations iso-energetic with the first excited state for perfect sizes or with the ground state for facile sizes. See Fig. 2(c). However, for long-range diffusion, most of these eroded corners must be rebuilt as the cluster must repeatedly pass through the unique ground state for perfect sizes and the “special” ground state for facile sizes (unless accessing higher excited states).

Next, we present an overview of KMC results illustrating various regimes and branches of D_N behavior for ϕ = 0.24 eV and δ = 0 at 300 K. See Fig. 3. Small sizes N = 4-8 all have the form N_p+1 or N_p+2 and exhibit facile diffusion with high D_N. (Even higher D_N values for N = 2 or 3 where E_eff = E_e are not shown.) For moderate sizes, N = 9 to O(10²), for clarity we just show four distinct branches: facile N_p+1, facile N_p+2, perfect N_p, and slow N_p+3. Key features are as follows: (i) initially high values and rapid decay of D_N ∼ N^−βf for facile sizes N_p+1 with large β_f ≈ 2.6 up to N ∼ 100 (see the supplementary material) and similarly high D_N but less regular decay for facile sizes N_p+2; (ii) lowest values of D_N for sizes N_p+3 (not N_p) with slow decay D_N ∼ N^−βs, where β_s ≈ 0.53 for N ∼ 67-200; (iii) very weak size-dependence of D_N for perfect sizes up to N_p ≈ 81; also perfect and N_p+3 branches merge for small N = 12; (iv) intermingling of D_N for perfect N_p with facile branches for N_mingle ≈ 81 (dashed yellow arrow in Fig. 3), and subsequent transition to a rapid decrease of D_N for perfect sizes; (v) near-merging of all branches for N ≈ N_merge ≈ 250 (solid brown arrow in Fig. 3). For larger sizes N > N_merge, if we write D_N ∼ N^−βeff, the effective exponent slowly varies from β_eff ≈ 0.75 for N just above N_merge, to β_eff ≈ 1.12 for N from 500 to 1000, and to β = 1.5 (the macroscopic value) for N → ∞. Here, the kink separation is L_k = ½ exp[½ϕ/(k_BT)] ≈ 52, and the asymptotic regime N >> (L_k)² ≈ 2700 is not achieved in our simulations.

We have also assessed behavior for ϕ = 0.20 eV at 300 K. All the features for ϕ = 0.24 eV are preserved qualitatively, but now β_f ≈ 2.3 up to N ∼ 100, β_s ≈ 0.84 for N ∼ 67-200, N_mingle ≈ 49, N_merge ≈ 150, with β_eff ≈ 1.06 just above N_merge, β_eff ≈ 1.33 for N from 500 to 1000, and β_eff ≈ 1.48 for N from 2000 to 3600. Here L_k = 24 consistent with achieving asymptotic scaling N >> (L_k)² ≈ 570.

Next, to provide fundamental insight into D_N behavior, we exploit instructive combinational analysis of cluster configurations and also apply first-passage concepts:

Fast decay of D_N for facile cases. For facile N_p+1 clusters, we find scaling of D_N ∼ N^−βf exhibits exponents β_f ≈ 2.3-2.6 far exceeding any identified previously for 2D cluster diffusion. We attribute this behavior to the feature that the cluster can wander between a large number of iso-energetic ground state configurations, Ω_N(0) ∼ N^α, with α ≈ 2.8 for N = O(10²). Most configurations are far removed from the special configuration with one edge atom on a perfect core through which the cluster must pass for long-range diffusion. If ⟨t_N⟩ denotes the first-passage time for the walk through configuration space to return to the special configuration, one expects that D_N ∝ 1/⟨t_N⟩. Analysis for random walks in any dimension suggest that³² ⟨t_N⟩ ∝ Ω_N(0) so that D_N ∼ 1/Ω_N(0) implying that β_f ≈ α.

Intermingling of D_N for perfect and facile branches. The distinction between perfect and other sizes of clusters is predicated by the former primarily existing in square or near-square ground state shapes. However, the number of configurations, Ω_Np(1), corresponding to first excited states with energy ΔE = +ϕ above the ground state increases strongly with N_p, e.g., Ω_Np(1) = 1140, 2472, 5152, 10 352, 20 208,… for N_p = 36, 49, 64, 81, 100,…, respectively. Perfect sized clusters have a significant probability of being in the 1st excited state when¹⁷ Ω_Np(1) ≈ exp[ϕ/(k_BT)] = 2208 (10 764) for ϕ = 0.20 (0.24) eV at 300 K corresponding to N_p ≈ 49 (81). This roughly matches the size where perfect and facile branches intermingle (dashed yellow arrow in Fig. 3 for ϕ = 0.24 eV).

Behavior for perfect sizes. When N is not too large, diffusion is largely controlled by the nucleation step, which depends weakly on N and not so much on the subsequent transfer of atoms to complete the new edge. This explains the weak dependence of D_Np on N_p. The property that D_Np for larger N_p = L² actually exceeds that for smaller N_p = L(L − 1) reflects the feature that nucleation for perfect rectangular shapes is slightly more likely to occur on the longer side which increases the chance that subsequent nucleation will return the cluster to its original configuration.

Variation of D_N within each cycle. D_N varies quasi-periodically within each cycle between consecutive perfect sizes. The local minimum always occurs for size N = N_p+3. Nucleation-mediated diffusion corresponds to sizes N = N_p+n with n = 3, 4, …, n*, where n* corresponds to the next largest perfect size (e.g., n* = 4 for N_p = 7 × 7 and N_p+n* = 7 × 8). D_N increases smoothly with N = N_p+n in this range from a minimum for n = 3 to a maximum for n = n*. A relevant observation is that the degeneracy of the ground state decreases strongly with n from a maximum for n = 3 to a minimum of 1 or 2 for n = n*. These degenerate states include multiple atoms shifted from various corners of the cluster and thus multiple kinks which trap atoms. This makes it more difficult to efficiently nucleate a new outer edge (as the lifetime of isolated atoms is reduced) and also to transfer atoms to complete that new outer edge (which is required for long-range diffusion). Consequently, D_Np+n unexpectedly increases with n. A local maximum within each cycle in D_N occurs for the facile case N = N_p+1 (or n = n* + 1).

“Oscillations” were observed in previous simulation studies,^17,33 but limited analysis can provide misimpressions, e.g., that perfect sizes N_p = L² diffuse slowest. Note that D_N values for n = n* (perfect sizes) can be quite close to those for the local maximum for facile n = n* + 1 and especially n = n* + 2. Thus, one might question the assignment of nucleation-mediated versus facile. However, an Arrhenius-type plot of D_N does show clearly the distinction between E_eff for these classes. See Fig. 4.

Weak variation of D_N versus N for N = N_p+3. While the high ground state degeneracy for facile sizes N_p+1 produces rapid decay of D_N versus N, this is not the case for N = N_p+3. Why? The fundamental difference is that long-range diffusion of clusters for sizes N_p+3 does not require repeatedly passing through a single special configuration, unlike for N_p+1.

Merging of branches and post-merging behavior. Near-merging of all branches of diffusivity, and in particular of facile N_p+1 and slow N_p+3 branches, occurs for N ≈ N_merge. Here, we argue that the distinctive nature of N_p+3 clusters (relative to N_p+1) is lost when the ratio of the number of 1st excited states Ω_Np+3(1) to the number of ground states Ω_Np+3(0) satisfies Ω_Np+3(1)/Ω_Np+3(0) ≈ exp[ϕ/(k_BT)]. We find that Ω_Np+3(1)/Ω_Np+3(0) ≈ 2565, 7002, and 11 116 for N_p = 196, 324, and 400, respectively, which implies that merging should occur for N_merge ≈ 400 for ϕ = 0.24 eV. This estimate is above N_merge ≈ 250 (solid brown arrow in Fig. 3) cited earlier. However, closer examination of data (see the supplementary material) reveals that D_N values for N_p+3 and N_p+1 branches actually slightly cross at N ≈ 250 and more properly merge at N ≈ 400. After merging, the effective scaling exponent, β_eff, describing the size-dependence of diffusivity is well below the asymptotic value of β = 1.5 from continuum theory. However, our simulations for ϕ = 0.20 eV indicate that β_eff increases with N to approach β = 1.5 for N above ∼6(L_k)².

In summary, KMC simulation of a model for diffusion of 2D epitaxial clusters reveals an extraordinary diversity of behavior with multiple distinct branches in the nanoscale size regime. We elucidate this behavior by exploiting combinatorial analysis. Surprisingly, perfect “closed shell” sizes are not the slowest, a feature plausibly extending to 3D supported clusters.⁷ Experimentally observed anomalous scaling^5,22 corresponds to behavior in the merged regime where the scaling exponent varies continuously towards the asymptotic value. Additional analysis (see the supplementary material) for non-zero kink rounding barrier, δ > 0, reveals that behavior for δ = 0 is preserved qualitatively, although β_eff is lower.

See supplementary material for simulation results for size-scaling of D_N for facile clusters, merging of facile and also branches, and behavior for non-zero kink rounding barrier.

K.C.L. and J.W.E. were supported for this work by NSF Grant No. CHE-1507223. They developed theoretical interpretations of the data and performed detailed KMC analysis. D.-J.L. was supported by the USDOE, Office of Science, Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biological Sciences, and his work was performed at Ames Laboratory which is operated by Iowa State University for the USDOE under Contract No. DE-AC02-07CH11358. He performed the exploratory KMC simulations revealing diverse behavior together with preliminary analysis.