Impurity-induced step interactions: a kinetic Monte-Carlo study

A one-dimensional continuum description of growth on vicinal surfaces in the presence of immobile impurities predicts that the impurities can induce step bunching when they suppress the diffusion of adatoms on the surface. In the present communication we verify this prediction by kinetic Monte-Carlo simulations of a two-dimensional solid-on-solid model. We identify the conditions where quasi one-dimensional step flow is stable against island formation or step meandering, and analyse in detail the statistics of the impurity concentration profile. The sign and strength of the impurity-induced step interactions is determined by monitoring the motion of pairs of steps. Assemblies containing up to 20 steps turn out to be unstable towards the emission of single steps. This behavior is traced back to the small value of the effective, impurity-induced attachment asymmetry for adatoms. An analytic estimate for the critical number of steps needed to stabilize a bunch is derived and confirmed by simulations of a one-dimensional model.


I. INTRODUCTION
Impurities and adsorbates affect the growth of crystals and thin films in a variety of ways. Small amounts of CO strongly enhance the nucleation density in homoepitaxial growth of Pt on Pt(111) and reverse the orientation of the resulting triangular islands [1,2]; a floating monolayer of a suitably chosen surfactant species induces layer-by-layer growth in many growth systems [2,3]; and the adsorption of antifreeze proteins on the growth surface prevents the formation of macroscopic ice crystals in the blood of fish living in polar waters [4].
On a vicinal surface growing by step propagation, impurities are generally expected to slow down the steps by pinning. This can lead to the formation of step bunches [5,6,7]. A different mechanism for impurity-induced step bunching not related to step pinning was recently proposed in the context of SiC growth on Si(100) where C plays the role of a codeposited impurity [8,9] (see also [10]). Due to the motion of the steps, at any given time different parts of a terrace have been exposed to the impurity flux for different durations, which leads to a gradient in the impurity concentration directed towards the ascending step. The coupling of the impurity concentration profile to the diffusion of the growth units on the terrace may destabilise the equidistant step train. Originally [8] a lower binding energy to the impurities was suggested and confirmed by kinetic Monte Carlo (KMC) simulations as a possible physical source for the experi- * Electronic address: juergen.vollmer@ds.mpg.de mentally observed behavior. Another scenario was discussed in [9]: Impurities slow down the adatom diffusion without affecting adatom binding energies (random barriers [11]) which also leads to instability.
The linear stability analysis of [9] was based on a onedimensional model of straight steps. The impurity and adatom concentration fields were treated in a continuum approximation and assumed to take on their stationary profiles instantaneously on the time scale of step motion. In the present communication we revisit the problem within a fully microscopic KMC simulation, taking explicit account of non-stationarity and fluctuations.
We identify a range of parameters where the assumption of a one-dimensional array of straight steps is applicable. We then consider systems of two, three and more steps with periodic boundary conditions in order to follow the loss of stability of the equidistant arrangement, and characterise the long-term evolution of the system. The basic stability properties predicted by the continuum theory are confirmed. Simulations with two steps show bound pairs and equidistant steps with only small fluctuations of the terrace width in the appropriate parameter regimes. However, we also find that the impurity-induced step interactions are unable to stabilise larger assemblies of steps. Simulations with 3, 4, 8, and 20 steps approach highly dynamic states with many closely adjacent pairs of steps that frequently exchange partners. This behavior can be explained within the framework of a deterministic step-dynamical model [12]: in spite of attractive interaction between the steps, bunches that contain less than a critical number of steps decay by step emission.
The paper is organized as follows. In the next section the KMC model is introduced and the fundamental growth modes are identified as a function of the system parameter description suppression of diffusion by bonds change of diffusion by impurities ν = g/Φ hopping rate / deposition rate ρ fraction of impurities in deposited atoms parameters. Section III contains a detailed analysis of the spatial distribution of impurities on the terraces, focusing in particular on the fluctuations around the mean impurity concentration gradient. The dynamics of step pairs, triplets and bunches is described and analyzed in Sec.IV, and conclusions are given in Sec.V.

II. MICROSCOPIC GROWTH MODEL
We consider an SOS system with a simple cubic lattice, where the surface has an extension L i × L j . There are periodic boundary conditions along i, and Lees-Edwards boundary conditions along j. The latter are also periodic on the surface, but they induce a change of height of N steps when transversing the system in j-direction. By this topological constraint one enforces the existence of N steps on the surface which are aligned parallel to i. Depending on the growth parameters these steps may be fairly straight, they may meander, or there may be additional islands on the surface.
The kinetic Monte Carlo simulation is constructed to be close to the continuum description [9]. In the following the technical details are described with special emphasis on the differences to previous simulations [8]. The transition rate of a thermally activated hopping process from site A to the neighboring site B is given in transition state theory by [2] The preexponential factor Γ 0 is taken usually to be 10 13 s −1 , but other values are possible also [13]. The activation energy is given by the difference of transition energy E T and binding energy E B . E T may depend on initial and final state, whereas E B only depends on the initial state. From now on these dependencies are suppressed. In the case of a simple cubic lattice usually the binding energy is described by a next neighbor counting model [14] with n being the number of in-plane next neighbors. In the present work the impurities are taken to influence the transition energy E T only, with E t being the transition energy of free adatoms. Previous simulations [8] had only considered the influence on the binding energy E B , which leads to a strong step bunching effect. Inserting Eq. (3) and (4) into Eq. (1) results in the hopping rate of an adatom residing on an impurity The factors J and J im are defined in Tab. I. Therefore, free adatoms on the surface perform a random walk with hopping rate g. Their hopping rate is modified by the existence of neighbors and impurities.
In the simulation adatoms are added randomly to the surface with a flux of Φ atoms per lattice site, and the ratio of diffusion and incoming flux will henceforth be characterized by ν ≡ g/Φ. An adatom on the surface will perform on average ν/(L i L j ) free steps before another atom is added anywhere on the surface.
Impurities disappear from the surface by burying them under normal atoms, and they are created by a flux onto the surface. A fraction of ρ of impinging atoms are impurities. When they hit the surface they immediately exchange position with a "normal" atom in the surface. Subsequently, there is a new impurity in the surface, and an additional adatom diffusing on the surface.
Altogether the dynamics of surface growth is hence characterized by four dimensionless parameters summarized in Table I. Kinetic Monte Carlo simulation of this model (cf [15] for details on the algorithm) show that the model is capable of reproducing the crossover from step flow at large ν, where adatoms mostly attach to step edges, to island nucleation for smaller ν, where adatoms merge into islands before reaching the step edges. Moreover, the simulations also show the formation of step pairs for appropriately chosen binding strength J im and density of impurities ρ (see Sec.IV). The different growth regimes for a system with N = 2 steps are summarized in Fig. 1. Fig. 2 shows snapshots of the time evolution of the surface height and the distribution of impurities for a system of two equidistant steps which evolves into a steady state where the two steps form a pair.

III. DISTRIBUTION OF IMPURITIES
Before further discussing the numerical results it is illuminating to calculate the distribution of impurities on the terraces. To this end we consider two steps which are roughly aligned in parallel such that they enclose a terrace of width w. In a system of two parallelly aligned steps the width w of the terrace does not change in a steady state where the steps form a bound pair, such as in the lower-most panels of Fig. 2. A small surface element is created in this situation when the edge of the step is formed at that position and height, and it is buried after deposition of two ML, when the two steps have reached the same position again. In order to study the distribution of impurities we use a comoving coordinate frame where i denotes the position in lateral direction, and d the distance from the position of the steps measured in the direction of growth. The latter distance will be measured in units of L j , such that there are d L j lattice sites between the considered site and the step edge. The age τ (d) of a position on the surface is proportional to dL j . The age will be measured in units of atoms added to the surface since the site has last been visited by a step. For N = 2 bunched steps moving together it amounts to The probability that there is no impurity at that site can be calculated as follows: The probability to turn a site into an impurity when a single atom is added to the surface is p = ρ/L i L j . Hence, the probability to be changed to an impurity after τ (d) atoms have been added to the surface is The latter approximation applies provided that L i L j /ρ ≫ 1. For the simulations shown in Fig. 2 one has L i L j /ρ = 25 000 such that this approximation is well justified. Moreover, in this situation the argument of the exponential function is small, such that the Note that impurities at the terrace are not uniformly distributed any more. There are more impurities ahead of the steps than behind the steps.
distribution of impurities is approximately linear, Fig. 3 demonstrates that this prediction holds to a very good approximation for the time average over 8 ML. Plotting the ratio of the numerical values and the theoretical prediction (lower panel) shows that the agreement is better than 5% along the full width of the terrace. In contrast to time averages the distribution is very noisy, however, for snapshots of the distribution of impurities as shown by three different examples in the lower panel of Fig. 3. Indeed, for the considered deposition pro- cess the number of impurities in a row d is distributed according to a binomial distribution such that one expects to find N d = L i P i (d) impurities in a row at a normalized distance d from the steps. The standard deviation should be L i P i (d) (1 − P i (d)) such that the relative error, which is indicated by a thick dashed green line in the lower panel of Fig. 3, takes the value For small d hardly any atoms impinged on that part of the surface such that the probability P i is small. In that case the relative error exceeds 100%. From the perspective of understanding the transport problem this is not so problematic, however, because the impurities hardly influence the diffusion on the surface in that case. It is much more important to observe that the relative error remains fairly large even far away from the steps, where the coverage of impurities is large. In the case of Fig. 3 where we consider expectation values for L i = 25 and the final coverage with impurities is close to 25% the relative error still amounts to (0.8/(25/5)) 1/2 = √ 0.16 = 0.4, and even for a coverage of 50% it only drops down to 0.2. Fluctuations of this magnitude are troublesome, when trying to describe the transport by a continuum model which does not take into account fluctuations of the distribution of impurities. The problem can not be resolved by considering averages of larger L i because the continuum description has to be based on local averages, and a one-dimensional description will only apply when the two-dimensional equations obtained in this manner are invariant under translation parallel to the steps. It is this latter property, however, which is lost when the fluctuations in the density of impurities are noticeable.
The observation that the variance in the number of impurities in small neighborhoods of the lattice is large probably applies in general. This poses a major challenge to continuum models of step bunching where the presence of impurities is only taken into account as a modification of the diffusion coefficient, which itself depends on the expectation value for the number of impurities. In view of the large fluctuations of the distribution this might very well be a non-admissible oversimplification, which deserves a close inspection by comparison to numerical results in the following section.

A. Diffusion bias
The key ingredient of the continuum theory developed in [9] is the dependence of the effective adatom diffusion coefficient D(θ) on the local impurity coverage θ. As we are concerned here with impurity concentrations of θ ≃ 0.1 or less, the leading term in an expansion in θ is expected to suffice. We therefore write For completely blocking barriers (J im → ∞) the coefficient is given by α = π − 1 ≃ 2.14 [16]. In effective medium approximation [11] one obtains the simple expression For J im → ∞ it yields α = 2, which is close to the exact result.
Using the general formulae derived in [9], the adatom currents j ± to the ascending (j + ) and descending (j − ) steps bordering a terrace of width w can be computed from Eq. (9). For perfectly absorbing steps one finds that j ± = Φwp ± , where the attachment probabilities p ± are independent of the terrace width w, and given by where the relation p + ≡ 1 − p − ensures flux conservation. For small αρ the attachment probabilities p ± amount to such that the effect of the impurities can be quantified by the diffusion bias parameter The sign is chosen such that step bunching results for b > 0 [9,17].
The following considerations will be based on a simple one-dimensional, deterministic dynamical model for the step positions x k , k = 1, .., N measured along the j-direction. The speed of the k th step is the sum of the fraction p + of the flux incident on the (leading) terrace in front of the step, of width x k+1 − x k , and the fraction p − of the flux incident on the (trailing) terrace behind the step, of width x k − x k−1 . Thus we have with the additional constraint that x k > x k−1 at all times. For convenience we use here dimensionless units where time is measured in units of the time scale Φ −1 needed to deposit a monolayer of new material, and length still in units of L j .

B. Step-step interactions
Before turning to the discussion of our simulation results, we need to address the role of repulsive step-step interactions that are usually added to the right hand side of (14) [12,18]. As no direct step-step interactions are included in our KMC model, we only consider the well-known entropic interactions induced by collisions between neighboring steps, which in turn are a consequence of thermal step meandering. Following [17], we estimate the distance between two such collisions along the transverse (i-) direction to be of the order of whereδ is the step stiffness and w denotes the distance between the two steps. Clearly step collisions are irrelevant as long as L c is larger than the lattice size L i parallel to the steps. We therefore conclude that the range of the repulsive step-step interactions in our simulations is limited to step distances smaller than w c ∼ (L i k B T /δ) 1/2 . Using the expression [2,17] for the step stiffness in the SOS model at low temperatures, we find that w c ≃ 2.8 for L i = 25 and J = 40. Thus the step-step interactions in our simulations are a purely local effect which merely ensures that steps cannot overtake each other. In this sense the situation is similar to that considered in [19] within a one-dimensional model with hard core step-step interactions.

C. Stability of step pairs
Consider first a system of two steps, with normalized step distances d(t) and 1 − d(t) (Fig. 4). It follows from (14) that d(t) evolves according tȯ The equidistant fixed point d = 1/2 is unstable (stable) for b > 0 (b < 0). For b > 0 the steps collide (d →  Fig. 4 for two steps) was needed to clearly show this effect. In the marginally stable case (b) the distances fluctuate showing regularly looking oscillations. They arise due to coupling of step velocities because of shared neighboring terraces. In this case the steps approach each other much closer than in the repulsive case. However, still they always remain well-separated due to entropic repulsion, which is always present. Finally, in the case of impurity-induced attraction between the steps (a), one clearly sees bunches of two steps, which separate however when the third step approaches. 0) in finite time, while for b < 0 a pair of close steps separates and approaches the fixed point exponentially at rate 2b. Using the expression (13) with α given by Eq. (10) the half-time of the decay into the stable states, t 1/2 ≡ ln 2/2b is about 16 ML. This value is consistent with but somewhat larger than the time scale observed in the simulations Fig. 4(a) and (c). A possible source of this deviation is the fact that the impurity concentration profile may not have reached stationarity on the time scale of step motion.

D. Stability of small bunches
Consider next a system of three steps, where d 2 and d 3 denote the normalized distance from the first to the second and the first to the third step, respectively. Fig. 5 shows the evolution of d 2 and d 3 as a function of time t measured in units of deposited ML.
In the absence of impurity-induced step interactions (J im = 1, Fig. 5(b)) and for repulsive impurity-induced step interactions (J im = 0.2, Fig. 5(c)) the system behaves very similar to the one with only two steps. Indeed, histograms for the probability to find a certain distance w between adjacent steps in Fig. 6 consistently show sharply peaked distributions around the the average terrace size w = 50, while the distribution is broad with a maximum at this value in the neutral case. As expected for a system with attractive interactions between the steps, the distribution P (w) has considerably more weight for small w. However, surprisingly, it does not decay but saturates at fairly large constant background extending till w = 100.
The origin of this background becomes clear from inspection of the time traces of d 2 and d 3 shown in Fig. 5(a). The simulation shows the transient formation of step pairs which exchange partners at regular intervals. However, no stable step triplets are formed. To see how this follows from the dynamical equations (14), suppose a triplet of three nearby steps has formed, such that the step positions satisfy x 2 − x 1 , x 3 − x 2 ≪ 1, i.e., they are both much smaller than the system size L j . Then step 3 has a large terrace of size ≃ 1 in front of it, and it moves at speed (1 − b)/2. Step 2 is surrounded by small terraces and moves very slowly, and step 1 is constrained by the no-passing condition to move at the same speed. Thus for every b < 1 step 3 will detach from the triplet.
We conclude that, for any b < 1, step triplets and larger bunches are unstable against the emission of single steps. Bound states of steps moving together at constant speed [20] can form when b > 1 [12], but this condition obviously cannot be reached in the growth model considered here. A detailed analysis of the equations (14) shows that step emission will continue until the number of free steps between two bunches (or, equivalently, between one bunch and its periodic image) has reached the steady state value [12] where N denotes the total number of steps in the bunch and on the terrace. As N f cannot exceed N , we conclude that a stable bunch can form only if the number of steps satisfies the condition 3b N/ ln N > 1.
With the value of b ≃ 0.02 obtained for J im = 4 and ρ = 0.1 this implies N > 65, which (given the constraints on the systems parameters described above) exceeds our computational capacities. Indeed, simulations conducted with systems containing up to 20 steps confirm that step bunching remains a transient phenomenon, even when a bunched initial configuration is chosen (Fig. 7).

E. Stability of large bunches
The agreement between the predictions and the results of the kinetic Monte-Carlo simulations indicates that the stability of bunches can be captured based on the onedimensional growth model. To illustrate the formation of stable bunches with increasing N we therefore relied on simulations of a one-dimensional stochastic growth model described in detail in [19]. In this model particles are deposited onto a one-dimensional vicinal surface and transferred instantaneously (without explicit diffusion) to the ascending or descending step with probabilities p + and p − , respectively. Since steps are allowed to coalesce (but not to pass each other), the formation of bunches can easily be followed by monitoring the maximal step height h max (the largest nearest neighbor height difference) in the system. In Fig. 8 we show results obtained for b = 0.02 on a lattice of L j = 1000 sites. For N = 20, the situation corresponding to Fig. 7(b), the maximal step height remains below 2, showing that only step pairs exist, whereas for N = 100 a stable bunch forms that contains almost half of the steps in the system. and N = 100 steps, respectively. The respective data points in these plots correspond to initial conditions with a single large step (red) and an equidistant step train (green), which were averaged over 100 independent runs. Beware of the different time scales of these two plots. The inset in the upper graph shows the ratio of the data of the 1D model (red curve) to the maximal height of a bunch in the corresponding KMC simulation Fig. 7. The bunch size was calculated by determining the maximal number of steps in a window of fixed width d = 0.03, d = 0.06, and d = 0.09 (from bottom to top) roughly corresponding to the range of the effective hard-core repulsion N wc/Lj . The lowermost panel shows the deviation ∆hj ≡ hj − 10 5 of the height hj from the average height after deposition of 10 5 ML. It is the final configuration of one of the simulations run to generate the graph in the middle panel.

V. CONCLUSIONS
We summarize the main achievements of this work: (i) We have verified by KMC simulations that an impurity-induced increase of the transition energy between neighboring sites, a purely kinetic effect, is a possible source of step bunching. By focusing on the behavior of pairs of steps we have explicitly determined the sign and strength of the impurity-induced step interactions, as quantified by the asymmetry parameter b. The adatoms can make use of massive fluctuations of the spatial and temporal distribution of the impurities (Fig. 3) in order to find optimal paths to the binding sites at the steps. As a consequence an estimate of the order of magnitude of b properly has to account for the diffusion of the adatoms in a 2d disordered arrangement of impurities, which is very different from the 1d setting used in [9].
(ii) The step bunching observed in our work is substantially weaker than that found in simulations of the SiC system [8], where the impurities were assumed to affect the adatom binding energies. This probably indicates that the effective asymmetry b is larger for energetic impurities. Unfortunately, although the theory of [9] correctly predicts step bunching for impurities that lower the binding energy, the magnitude of the effect depends on the precise boundary conditions at the steps, which do not easily translate into the two-dimensional KMC setting.
(iii) Despite the presence of impurity-induced attractive step interactions, triplets and larger assemblies of steps are not necessarily stable when b is small: Instead of agglomerating into macroscopic step bunches, the steps display a peculiar dynamical pattern characterized by transient step pairs that exchange partners much like in a folk dance. Up to now this effect has gone unnoticed, because previous KMC simulations of step bunching during growth have generally considered situations where the effective attachment asymmetry b is of order unity [21,22,23]. The lack of stability of the step assemblies was explained based on a recently developed determinstic theory [12]. It allows us to predict that bunches can form only when the number of steps exceeds the bound (19).
(iv) For real surfaces there is no restriction on the total number of steps. Nevertheless, it is highly improbable to observe bunching in systems with small b. Our simulations for the one-dimensional model (Fig. 8) show that for small b bunches evolve only after exceedingly long times even when the bound (19) is satisfied: In order to see step bunching one has to wait for a fluctuation nucleating a bunch with a minimal size given by (19).
For applications the most noticeable consequence of our study is that the rapid formation of large step bunches seen experimentally in the growth of SiC [8] cannot be explained only in terms of kinetic impurities. Some coupling to the adatom binding energy must also be involved.