Low coverage surface diffusion in complex energy landscapes: Analytical solution and application to intercalation in topological insulators

A general expression is introduced for the tracer diffusivity in complex periodic energy landscapes with more than one distinct hop rate in two- and three-dimensional diluted systems (low coverage, single-tracer limit). For diffusion in two dimensions, a number of formulas are presented for complex combinations of hop rates in systems with triangular, rectangular and square symmetry. The formulas provide values in excellent agreement with Kinetic Monte Carlo simulations, concluding that the diffusion coefficient can be directly determined from the proposed expressions without performing such simulations. Based on the diffusion barriers obtained from first principles calculations and a physically-meaningful estimate of the attempt frequencies, the proposed formulas are used to analyze the diffusion of Cu, Ag and Rb adatoms on the surface and within the van der Waals (vdW) gap of a model topological insulator, Bi$_{2}$Se$_{3}$. Considering the possibility for adsorbate intercalation from the terraces to the vdW gaps at morphological steps, we infer that, at low coverage and room temperature: (i) a majority of the Rb atoms bounce back at the steps and remain on the terraces, (ii) Cu atoms mostly intercalate into the vdW gap, the remaining fraction staying at the steps, and (iii) Ag atoms essentially accumulate at the steps and gradually intercalate into the vdW gap. These conclusions are in good qualitative agreement with previous experiments.

In this study we are interested in the surface diffusion of an adparticle on a complex periodic energy landscape, such as the one shown in figure 1(a). One can discern the presence of 4 different adsorption sites (labeled as f , h, t and b), corresponding to the locations where the energy has a local/global minimum. A diffusing particle proceeds by hopping between the sites, as shown schematically in figure 1(b). If one focuses on a particular site, as illustrated in figure 1(c)-(f) for h, b, f and t, respectively, one may notice that each hop requires surpassing a different energy barrier. As further emphasized in figure 1(g), the different energy barriers and dissimilar shapes of the energy wells (= dissimilar attempt frequencies) result in different hop rates for the different jumps (ν ht , ν th , ν hb , ν bh , etc...), including the forward and backward directions. As a result, the random walk between points A and B in figure 1(b) involves as many as 10 different hop rates for a total of 14 performed hops. In this study we focus on describing analytically the average distance travelled by the adparticle as a function of the hop rates ν ij when the number of hops grows very large, i.e. the diffusion time becomes arbitrarily large.
The average squared distance covered by a single particle per unit time is a well-defined quantity, known as the tracer diffusion coefficient (or tracer diffusivity) [1]: where α = 1, 2, 3 is the number of dimensions, n is the number of adparticles simultaneously present on the surface, r i (t) designates the position of adparticle i at time t, and · is the ensemble average. Not surprisingly, D T is a function of the number of adparticles n or, equivalently, of the coverage θ = n/m, where m is the number of adsorption sites that may be occupied by the adparticles. The larger the number of adparticles the smaller the number of available empty sites where any chosen adparticle can jump to, thus leading to correlation effects between consecutive hops, also known as memory effects [1,2]. This is specially relevant for systems with strong adsorbateadsorbate interactions [36,37].
(a) Example of a complex potential energy landscape for a diffusing particle. Four adsorption site types are indicated: f , h, t and b. (b) A possible diffusion track (random walk) involving 14 performed jumps with 10 different hop rates ν ij . (c)-(f) Possible hops from each site type (h, b, f and t, respectively). Grayshaded landscapes in (b)-(f) are used to highlight the colored arrows/hops. (g) Energy paths for hops starting/ending in h sites. A longer arrow assigned to a hop rate (ν th , ν ht , etc...) indicates a larger rate.
We are interested in the low coverage regime, where the adparticle density is so low that the chance of affecting each other's motion is negligible: The low coverage limit is an important measure as it provides a simple procedure to compare the typical distances covered by different adsorbates across different substrates [1,2,3,4,5]. Previous analytical work on diluted systems has focused on the determination of the center-of-mass diffusivity for small 2D islands and clusters on metal surfaces based on the master-equation [38,39,40] or the continuous time random walk formalism [41,42]. In the framework of bulk-mediated surface diffusion, Revelli et al. described the average motion of the adsorbed molecules, including both Markovian and non-Markovian desorption, by using the generalized master-equation approach [43]. Birnie et al. [44] and Condit et al. [45] derived the overall jump rate for complex, sequential diffusion paths in three-dimensional (3D) crystals, where typical vacancyinterstitial complexes evolve by repeating a particular sequence of hops. The present study generalizes such sequential analysis by providing a universal expression for the diffusivity in complex hopping networks where both parallel and sequential diffusion routes are available between the different adsorption sites. Computationally, adsorbate diffusion is traditionally studied by [1,5]: (i) first principles calculations, typically involving the use of density functional theory (DFT) for the determination of the activation barriers and attempt frequencies, which are then passed to other methods; (ii) molecular dynamics simulations, which numerically solve Newton's equations for the substrate and adsorbate atoms based on effective interaction potentials, enabling the analysis at the picosecond time scale for systems with ∼10 5 atoms; (iii) Langevin models, which describe the adparticle in an effective periodic force field, restricting the analysis to general trends for time scales of picoseconds; and (iv) Kinetic Monte Carlo (KMC) simulations, which focus on describing the hops between adjacent basins (rare events) while disregarding all other vibrations, this way enabling long simulated times (seconds and minutes) with affordable computational resources. In this study we validate the proposed formulas for the diffusivity by direct comparison to KMC simulations in relevant energy landscapes, finally using the formulas to discuss the relative mobility and intercalation of various adsorbates in the context of topological insulators.

Theory
Let us consider a system with S different site types, such as the one shown in figures 2(a) and 2(b) for S = 4, or figure 2(c) for S = 3. Although we assume diffusion in two dimensions, the underlying mathematical treatment and main result are also valid for three dimensions. For a generic hop from site type i to site type j (i, j = 1, ..., S) we consider that the hop distance l ij , hop rate ν ij and hop multiplicity n ij are known, and we define a new variable, the rateplicity µ ij , as the product of the rate and the multiplicity: We then propose that the low coverage diffusivity D θ≈0 T , as defined in equation 2, can be written as a weighted sum of partial diffusivities: where (Σ j µ ij l 2 ij ) is the partial diffusivity from site i, which contains the rateplicities and hop lengths for all jumps from site type i to any accessible site type j, and the dimensionless coefficient w i (or normalized weight for site i) stands for the probability to find the adatom at site type i: where the B i coefficients consist of sums and products of rateplicities, their particular form depending on the number of adsorption sites S. As an example, for S = 2 and S = 3 we have: while a slightly more complex definition can be used for S = 4: In practice, to determine each B i it is convenient to regard i as the end site of a jump from a neighboring site k, so that B i = B ki , which is then determined by applying equations 9-11 (S = 4) or equation 7 (S = 3) or equation 6 (S = 2). Any value of k different from i can be used since the B ij coefficients depend only on the second index . This is easily demonstrated by writing out B ki and B k i according to equation 9 and confirming their correspondence. For S ≥ 5, a general procedure to determine the B i coefficients is described in the Supplementary Data, where also a rigorous derivation of equation 4 is provided. Valid for any value of S in two and more dimensions, equation 4 is our central result.
Although equation 4 considers all S ×S hops between all possible pairs from a set of S different site types, non-occurring jumps can be eliminated from equation 4 by setting their rateplicities to zero (µ ij = 0). This may lead, however, to undetermined values, e.g. b i = 0 0 , and it is better to substitute the zeroed rateplicities by a small value ( ) and take the limit → 0. As an example for the system shown in figure 1, if we consider only the f h, f t, hf , ht and tf hops (with multiplicity equal to 3) while disregarding all other hops, we may set µ f f = µ f b = µ hh = µ hb = µ th = µ tt = µ tb = µ bf = µ bh = µ bt = µ bb = in order to calculate the B i coefficients and then take the limit → 0 to determine the b i factors. We have: Thus, using Equations 4 and 5 the diffusivity is: Examples of low coverage tracer diffusivities (D θ≈0 T ) for different combinations of hop rates between standard adsorption sites on square, rectangular and triangular lattices. Site labels: f = 4-fold hollow (square) / 2-fold hollow (rectangular) / fcc hollow (triangular), h = hcp hollow (triangular), t = on-top, b = bridge (square) / short-bridge (rectangular) and B = long-bridge (rectangular). See the Supplementary Data for additional formulas.  4 2α Sequential hops presentation of diffusivity expressions for 2D landscapes, equation 4 is completely general and can be applied to 3D problems as well. In fact, the last two equations of Table 1 for three and four sequential jumps, respectively, are identical to those derived by Condit et al. [45] and Birnie et al. [44], respectively, for vacancy diffusion in three dimensions. Our procedure, however, is more general, taking into account any number of competing diffusion paths from any given site (parallel processes) in addition to any number of successive hops along any given diffusion path (sequential processes).

Numerical validation
Based on the popularity of KMC simulations to determine tracer diffusivities [1,2,3,4,5,36,38,46], we now compare the values obtained from the previous formulas and those determined by KMC simulations. As described in figure 3(a)-(b), we perform two types of simulations. In the first type (KMC-1) the tracer is followed until it hits the perimeter of a circle of radius R o l ij , repeating the process for N RW different random walks (RWs) in order to obtain an ensemble average of the time < t > required to cover that distance, thus determining the diffusivity as D θ≈0 In the second type of simulations (KMC-2) the tracer is followed until it performs a desired number of hops N H , repeating the process for N RW different RWs to determine the average squared distance < X 2 + Y 2 > covered by the tracer and the corresponding average time < t >, obtaining the diffusivity by using D θ≈0 Since the goal is to check the validity of the analytical expressions for the diffusivity, different hop rate values are used to probe situations where the rates have similar values / differ by several orders of magnitude. We also use realistic hop rates for several adsorbates, including Cu, Ag, Rb and Se on the Bi 2 Se 3 (0001) surface and in the vdW gap of this material. In this case the hop rates are expressed as ν = ν 0 e −Ea/k B T , where the Boltzmann factor e −Ea/k B T reflects the probability to perform the jump at temperature T if the energy barrier is E a , and ν 0 is the attempt frequency, which reflects the dynamical coupling between the substrate phonons and the adparticle vibrations [1]. The actual values for these hop rates are obtained by determining the energy barriers through labor-intensive DFT calculations (see the Supplementary Data) and estimating the attempt frequencies by the method described in figure 3(c)-(e), leading to: Here, E a is the energy barrier for the hop, d is the separation between the initial site A and the saddle point T (transition state) and m is the adatom mass. This estimate can be considered as an alternative to (i) the typical assumption of equal prefactors for all hop rates [1,47,48,49,46], and (ii) the large computational cost to determine all the vibrational mode frequencies ν A i and ν T i for the adatom and substrate at the initial and saddle configurations by using DFT methods [1,47,48,49] (see expression 'U' for the attempt frequency in figure 3(d)). Due to typical cancellations of the vibrational modes of the substrate [49], the prefactor ν o is usually approximated by the Vineyard equation (expression 'V' in figure 3(d)). Owing to compensation effects between the surface-parallel (ν A and ν T ) and surface-normal (ν A ⊥ and ν T ⊥ ) vibrational frequencies of the diffusing atom [49], the Vineyard formula is approximated by just keeping the frequency of the longitudinal vibrations (along the diffusion path) of the atom at the initial site (ν A L , see expression 'W' in figure 3(d)). Our estimate consists on approximating the longitudinal path by a sinusoidal path, resulting in a simple and physically meaningful expression for ν A L in terms of the energy barrier, hop distance and adatom mass, as described in figure 3(f) and contained in equation 22. Although the actual energy path can be asymmetric with respect to the saddle point T, the part after this point is irrelevant for the rate calculation and in our approximation it is considered to be a reflection of the part before it.
For all considered systems the simulated and calculated diffusivities agree extremely well with each other (see Table S5 of the Supplementary Data), thus concluding that the proposed formulas are suitable to discuss the low coverage tracer diffusivity of typical diffusion species in complex energy landscapes without any need to perform the corresponding KMC simulations.

Application to topological insulators
Encouraged by the validation of the diffusivity formulas we now consider the temperature dependence of the diffusion of Cu, Ag and Rb adatoms on the Bi 2 Se 3 (0001) surface and the corresponding intercalation of these adsorbates in the Bi 2 Se 3 vdW gap.  Table S5 of the Supplementary Data). For diffusion with a single barrier E a , the underlying assumption that the hops can be treated as rare events (as compared to the fast vibrations around the adsorption sites) is valid if E a > 4k B T (see Ref. [1]). Correspondingly, each displayed curve in figure 4 terminates at T max = E a,max /4k B , where E a,max is the maximum barrier experienced by the corresponding adatom. As an example, the diffusion of Rb on the surface experiences two barriers: E f h = 0.127 eV and E hf = 0.104 eV. Thus, E a,max = 0.127 eV and T max = 368 K. Similarly, for Cu in the vdW gap we have T max = 342 K. Nevertheless, the total temperature range is restricted to 600 K since the desorption of Cu is reported to start at ∼550 K [50].  It can be seen from figure 4 that the hierarchy of the diffusion length is Λ Rb > Λ Cu > Λ Ag on the surface and Λ Cu > Λ Ag > Λ Rb in the vdW gap. The Rb (Cu) atoms are the most mobile species on the surface (in the vdW gap), capable of covering more than 1 µm within one minute even at 100 K. However, the vdW (surface) diffusion length of the Rb (Cu) atoms is much lower than that on the surface (in the vdW), which is due to significantly higher diffusion barriers. Interestingly, the Ag atoms travel with almost equal rates on the surface and in the vdW gap.
Let us now bring the discussion closer to the available experiments [21,23,30].  These indicate indirectly the occurrence of partial [21] or almost complete [23,30] intercalation of the metal adatoms inside the Bi 2 Se 3 vdW gap at room and higher temperatures. Recently, it has been argued [28] that the intercalation of the metal atoms in the Bi 2 Se 3 vdW gap is step-mediated, in the sense that they penetrate into the vdW gap after reaching the steps, which are typically present at the Bi 2 Se 3 surface [23]. This is favored by the geometrical alignment between the terrace and the vdW gap in a stepped surface, as schematically shown in 5(a). At the same time, penetration of Cu [23] and Ag [28] in the vdW gap via interstitials and/or vacancies of the topmost Bi 2 Se 3 QL is significantly less probable due to high energy barriers. Therefore, one may also rule out the possibility of vertical penetration of the Rb atom, whose covalent radius is 1.5 (1.66) times larger than that of Ag (Cu), whereupon the diffusion barriers are expected to be even higher.
By using large and small double-head arrows, figure 5(b)-(d) presents a graphical description of the relative diffusivity of the three types of adatoms on the terraces and within the vdW gaps. In addition, the figure also provides the relative rates to enter the vdW gap and to bounce back to the terrace by assigning them large/small unidirectional arrows. Moreover, the figure collects all the diffusion barriers determined by our DFT calculations for the three adatoms on the terrace and in the vdW gap, as well as for their penetration into the vdW gap and re-entry into the terrace for two different step orientations, namely, [1120] and [0110], the latter having two possible atomic terminations [28]. To ease the discussion, the three barriers for vdW penetration (terrace re-entry) for each adatom are algebraically averaged and used to estimate the vdW penetration (terrace re-entry) rate of that atom, accordingly assigning them large/small unidirectional arrows. We center the discussion at room temperature (295 K) and long diffusion times.
The Rb atoms have the largest terrace and smallest vdW diffusivities, with a very low vdW-penetration rate and a rather large terrace-reentry rate. Thus, the Rb atoms are expected to quickly diffuse across the terraces and hit the steps, where a small fraction will remain trapped while the majority will bounce back to the terraces. This is described schematically in figure 5(b) by drawing a large number of Rb atoms in the terrace region, with additional atoms at the step, while hardly any atoms in the vdW gap. In comparison, the Cu atoms are characterized by moderate terrace and largest vdW diffusivities, with the largest vdW-penetration rate and a low terrace-reentry rate. Accordingly, the Cu atoms need a longer time to arrive at the steps but they eventually penetrate with relative ease into the vdW gap, where they diffuse rather fast. Thus, the Cu atoms are expected to mostly intercalate in the vdW gap, although a notable fraction will remain at the steps, as sketched in figure 5(c). Finally, the Ag atoms display slightly lower terrace diffusivity as compared to Cu and a medium diffusivity in the vdW gap among the three species under consideration. Having reached the steps, the Ag atoms are forced to linger along them due to the large terrace-reentry barrier and significant vdW-penetration barrier ( figure 5(d)). Since the latter is smaller, the Ag atoms are expected to gradually intercalate into the vdW gap. The behavior is similar to that of Cu, although intercalation is slower for Ag. Valid only for the low coverage limit, these trends are in good qualitative agreement with those from available experiments [21,23,30].

Conclusion
Focusing on the analysis of adsorbate diffusion on surfaces and within the twodimensional gap of layered materials, we present a combination of formulas for future reference and their application to intercalation in a topological insulator (Bi 2 Se 3 ). We start by presenting a general expression to determine the average motion of the diffusing particles at low densities in complex, periodic energy landscapes consisting of various energy barriers located between distinct adsorption sites in any number of dimensions. For adsorbate diffusion in two dimensions, formulas are provided for the low coverage tracer diffusivity for complex combinations of hop rates in systems with triangular, rectangular and square symmetry. The analytical expressions are validated against Kinetic Monte Carlo simulations, obtaining an excellent agreement between the calculated and simulated diffusivities. Thus, one can determine the overall diffusion coefficient without performing the KMC simulations. Based on diffusion rates from energy barriers obtained by labor-intensive density functional theory calculations, we analyze the temperature dependence of the diffusion of Cu, Ag and Rb on the Bi 2 Se 3 (0001) surface and within the van der Waals (vdW) gap inside the (layered) crystal. We also analyze the occurrence of adsorbate intercalation due to the alignment between the vdW gaps and terraces on [1120]-and two types of [0110]-stepped surfaces of this topological insulator. At room temperature and low coverage, we conclude that the Rb atoms quickly diffuse across the terraces and hit the steps, where a small number remain trapped while the rest bounce back into the terraces. Thus, Rb is expected to partially decorate the steps while remaining present on the terraces. In comparison, the Cu atoms take a longer time to arrive to the steps, eventually penetrating with relative ease into the vdW gap. Thus, Cu atoms are expected to mostly intercalate into the vdW gap while partially remaining at the steps. Ag atoms are expected to take an even longer time to diffuse across the terraces, eventually penetrating into the vdW gap with time but meanwhile remaining around the steps.