# A Discrete Model for Nonequilibrium Growth Under Surface Diffusion Bias

\[

## Abstract

A limited mobility nonequilibrium solid-on-solid dynamical model for kinetic surface growth is introduced as a simple description for the morphological evolution of a growing interface under random vapor deposition and surface diffusion bias conditions. Simulations using a local coordination dependent instantaneous relaxation of the deposited atoms produce complex surface mound morphologies whose dynamical evolution is inconsistent with all the proposed continuum surface growth equations. For any finite bias, mound coarsening is found to be only an initial transient which vanishes asymptotically, with the asymptotic growth exponent being $`0.5`$ in both 1+1 and 2+1 dimensions. Possible experimental implications of the proposed limited mobility nonequilibrium model for real interface growth under a surface diffusion bias are critically discussed.

\]

An atom moving on a free surface is known to encounter an additional potential barrier, often called a surface diffusion bias , as it approaches a step from the upper terrace — there is no such extra barrier for an atom approaching the step from the lower terrace (the surface step separates the upper and the lower terrace). Since this diffusion bias makes it preferentially more likely for an atom to attach itself to the upper terrace than the lower one, it leads to mound (or pyramid) - type structures on the surface under growth conditions as deposited atoms are probabilistically less able to come down from upper to lower terraces. This dynamical growth behavior is sometimes called an “instability” because a flat (“singular”) two dimensional surface growing under a surface diffusion bias is unstable toward three dimensional mound/pyramid formation. There has been a great deal of recent interest  in the morphological evolution of growing interfaces under nonequilibrium growth conditions in the presence of such a surface diffusion bias. In this paper we propose a minimal nonequilibrium cellular automata - type atomistic growth model for ideal molecular beam epitaxial - type random vapor deposition growth under a surface diffusion bias. Extensive stochastic simulation results presented in this paper establish the morphological evolution of a surface growing under diffusion bias conditions to be surprisingly complex even for this extremely simple minimal model. Various critical growth exponents, which asymptotically describe the large-scale dynamical evolution of the growing surface in our minimal discrete growth model, are inconsistent with all the proposed continuum theories for nonequilibrium surface growth under diffusion bias conditions. Our results based on our extensive study of this minimal model lead to the conclusion that a continuum description for nonequilibrium growth under a surface diffusion bias does not exist (even for this extremely simple minimal model) and may require a theoretical formulation which is substantially different from the ones currently existing in the literature. Our results in the initial non-asymptotic transient growth regime (lasting upto several hundred or a few thousand layers of growth) do, however, agree with existing theoretical and (many, but not all) experimental findings in the literature.

In Fig. 1(a) we schematically show our solid-on-solid (SOS) nonequilibrium growth model : (1) Atoms are deposited randomly (with an average rate of 1 layer/unit time, which defines the unit of time in the growth problem — the length unit is the lattice spacing taken to be the same along the substrate plane and the growth direction) and sequentially on the surface starting with a flat substrate; (2) a deposited atom is incorporated instantaneously if it has at least one lateral nearest-neighbor atom; (3) singly coordinated deposited atoms (i.e. the ones without any lateral neighbors) could instantaneously relax to a neighboring site within a diffusion length of $`l`$ provided the neighboring site of incorporation has a higher coordination than the original deposition site; (4) the instantaneous relaxation process is constrained by two probabilities $`P_L`$ and $`P_U`$ ($`0P_L,P_U1`$) where $`P_{L(U)}`$ is the probability for the atom to attach itself to the lower(upper) terrace after relaxation (note that a “terrace” here could be just one other atom). The surface diffusion bias is implemented in our model by taking $`P_U>P_L`$, making it more likely for atoms to attach to the upper terrace. Under the surface diffusion bias, therefore, an atom deposited at the top of a step edge feels a barrier (whose strength is controlled by $`P_U/P_L`$) in coming down compared with an atom at the lower terrace attaching itself to the step. Our model is well-defined for any value of the diffusion length $`l`$ including the most commonly studied situation of nearest-neighbor relaxation ($`l=1`$). (We should emphasize, however, that the definition of a surface diffusion bias is not unique even within our extremely simple limited mobility nonequilibrium growth model — what we study in this Letter is the so-called edge diffusion bias .)

We have carried out extensive simulations both in $`1+1`$ and $`2+1`$ dimensions (d) varying $`P_L`$, $`P_U`$ as well as $`l`$, also including in our simulations the inverse situation (the so-called ‘negative’ bias condition) with $`P_L>P_U`$ so that deposited atoms preferentially come down attaching themselves to lower steps producing in the process a smooth growth morphology. Because of lack of space we do not present here our ‘negative’ bias ($`1P_L>P_U`$) results (to be published elsewhere) except to note that the smooth dynamical growth morphology under our negative bias model obeys exactly the expected  linear Edwards - Wilkinson universality . Our growth model is the most obvious finite bias generalization of the well-studied cellular automaton model referred to as the DT model or the 1+ model in the literature .

Before presenting our numerical results we point out two important features of our growth model : (1) For $`P_L=P_U=1`$ our model reduces to the one introduced in ref. 22, often called the DT model, (and studied extensively  in the literature) as a minimal model for molecular beam epitaxy in the absence of any diffusion bias; (2) we find, in complete agreement with earlier findings  in the absence of diffusion bias, that the diffusion length $`l`$ is an irrelevant variable (even in the presence of bias) which does not affect any of our calculated critical exponents (but does affect finite size corrections — increasing $`l`$ requires a concomitant increase in the system size to reduce finite size effects). To demonstrate this, we compare two systems (see Fig. 2(c)), one with diffusion length $`l=1`$ while the other with $`l=5`$. Both systems are on sufficiently large substrates of size $`L=10^4`$ to prevent finite size effect and both have the same bias strength ($`P_U=1.0`$ and $`P_L=0.9`$). Except for the expected layer-by-layer growth seen in the first few layers in the $`l=5`$ system, we see the same critical behavior from the two systems with the coarsening processes slightly faster in the $`l=5`$ system in early time regime. In the rest of this paper (except for Fig. 2(c) where we present some representative $`l=5`$ results), we present our $`l=1`$ simulation results emphasizing that our critical exponents are independent of $`l`$ provided finite size effects are appropriately accounted for. Our calculated exponents are also independent of the precise values of $`P_U`$ and $`P_L`$ ( $`<P_U`$ ) as found  in the unbiased $`P_U=P_L`$ case.

In Fig.1 we show our representative d=1+1 (a and b) and 2+1 (c and d) simulated dynamical growth morphology evolution. The diffusion bias produces mounded structures which are visually statistically scale invariant only on length scales much larger (or smaller) than the typical mound size. Note that the mounding in the growth morphology starts very early during growth and is already prominent in the first 100 monolayers (ML) as is obvious in Fig. 1(b). In producing our final results we utilize a noise reduction  technique which accepts only a fraction of the attempted kinetic events, and in the process produces smoother results (reducing noise effects) without affecting the critical exponents. We have explicitly checked that the noise reduction technique does not change our calculated critical exponents  and has only the cosmetic effect of suppressing noise in our simulated growth morphology.

To proceed quantitatively we now introduce the dynamic scaling ansatz  which seems to describe well all our simulated results.

We have studied the root mean square surface width or surface roughness ($`W`$), the average mound size ($`R`$), the average mound height ($`H`$), and the average mound slope ($`M`$) as functions of growth time. We have also studied the various moments of dynamical height-height correlation function, and these correlation function results (to be reported elsewhere) are consistent with the ones obtained from our study of $`W(t)`$, $`R(t)`$, $`H(t)`$, and $`M(t)`$. The dynamical scaling ansatz in the context of the evolving mound morphologies can be written as power laws in growth time (which is equivalent to power laws in the average film thickness) : $`W(t)t^\beta `$ ; $`R(t)t^n`$ ; $`H(t)t^\kappa `$ ; $`M(t)t^\lambda `$ ; $`\xi (t)t^{1/z}`$, where $`\xi (t)`$ is the lateral correlation length (with $`z`$ as the dynamical exponent) and $`\beta `$, $`n`$, $`\kappa `$, $`\lambda `$, $`z`$ are various growth exponents which are not necessarily independent. We find in all our simulations $`nz^1`$, and thus the coarsening exponent $`n`$, which describes how the individual mound sizes increase in time, is the same as the inverse dynamical exponent in our model. We also find $`\beta =\kappa `$ in all our results, which is understandable in a mound-dominated morphology. In addition, all our results satisfy the expected exponent identity $`\beta =\kappa =n+\lambda `$ because the mound slope $`MH/R`$. The evolving growth morphology is thus completely defined by two independent critical exponents $`\beta `$ (the growth exponent) and $`n`$ (the coarsening exponent), which is similar to the standard (i.e. without any diffusion bias) dynamic scaling situation  where $`\beta `$ and $`z`$ ($`=n^1`$ in the presence of diffusion bias) completely define the scaling properties.

In Figs.2 and 3 we show our representative scaling results in d=1+1 (Fig.2) and 2+1 (Fig.3) for nonequilibrium growth under surface diffusion bias conditions. It is clear that we consistently find the growth exponent $`\beta 0.5`$ in both d=1+1 and 2+1 in the long time asymptotic limit independent of $`P_L`$ and $`P_U`$ as long as $`P_L/P_U<1`$. This $`\beta 0.5`$ is, however, different from the usual Poisson growth under pure random deposition with no relaxation where there are no lateral correlations. Our calculated asymptotic coarsening exponent $`n`$ in both d=1+1 and 2+1 is essentially zero ($`<0.1`$) at long times. In all our results we find the effective coarsening (growth) exponent showing a crossover from $`n`$ ($`\beta `$) $`0.2`$ ($`0.25`$ or $`0.33`$ depending on d=2+1 or 1+1) at early times ($`1<t10^3`$) to a rather small (large) value ( $`n<0.1`$, $`\beta 0.5`$ as $`t\mathrm{}`$) at long times — we believe the asymptotic $`n`$ ($`\beta `$) to be zero (half) in our model. Our calculated steepening exponent $`\lambda `$ satisfies the exponent identity $`\lambda =\beta n`$ rather well, indicating that steepening ($`Mt^\lambda `$) and coarsening ($`Rt^n`$) are competing processes. We find that during the initial transient regime (for $`1<t10^3`$ depending on d and $`P_L/P_U`$), when considerable mound coarsening takes place ($`n0.2`$), $`\lambda `$ is rather small ($`\lambda <0.1`$) and does not change much. We note that the mound formation dominates our growth morphology even during the initial transient — mounding starts early, coarsens rapidly, and then coarsening slows down or almost stops. After the initial transient, however, $`\lambda `$ is finite and large ($`\lambda 0.4`$), indicating significant steepening of the mounds. Based on our simulation results we are compelled to conclude that the evolutionary behavior of our mound dynamics makes a crossover (at $`t=t^{}10^3`$) from a coarsening-dominated preasymptotic regime $`I`$ ($`n_I0.2`$, $`\lambda _I<0.1`$) to a steepening-dominated regime $`II`$ ($`n_{II}<0.1`$, $`\lambda _{II}0.4`$) which we believe to be the asymptotic regime. The initial transient (regime $`I`$: $`t<t^{}`$) can be construed as a “slope selection” regime where the steepening exponent $`\lambda `$ is very small (and the coarsening exponent, $`n0.2`$, is approximately a constant), but the asymptotic long time behavior (regime $`II`$: $`t>t^{}`$) is clearly dominated by slope steepening (large $`\lambda `$) with coarsening essentially dying down ($`n<0.1`$). The crossover time $`t^{}`$ between regime $`I`$ (coarsening) and $`II`$ (steepening) depends on the details of the model (e.g. d, L, $`P_L`$, $`P_U`$, $`l`$), and could be quite large ($`t^{}10^210^4`$).

We emphasize that our initial transient ($`t<t^{}`$) exponents in regime $`I`$ ($`\beta _I0.26`$; $`n_I0.17`$; $`\lambda _I0.09`$ in d=2+1, Fig.3) agree quantitatively with several  non-minimal detailed (temperature-dependent Arrhenius diffusion) growth simulations as well as with direct numerical simulations  of a proposed (empirical) continuum growth equation, which have all been claimed  to be in good agreement with the observed coarsening behavior in various epitaxial growth experiments. Other simulations and theories  find no slope selection, which also agree with some experiments. We believe that the observed slope selection ($`\lambda 0`$) in experiments and simulations is a (long-lasting) transient (our regime $`I`$, $`t<t^{}`$) behavior, which should disappear in the asymptotic regime, at least in models without any artificial transient downward mobility. Theoretically, a distinction has been made  between ‘weak’ and ‘strong’ diffusion bias cases, corresponding respectively to our regime $`I`$ (‘mounding’) and $`II`$ (‘steepening’), respectively. Our results show that this theoretical distinction is not meaningful because the ‘weak’ bias case crosses over to the ‘strong’ bias case for $`t>t^{}`$, and the long time asymptotic regime is invariably the ‘strong’ bias steepening regime for any finite diffusion bias. (There could be accidental slope selections at rather large slopes when crystallographic orientations are taken into account , a process neglected in our minimal growth model.)

In comparing with the existing continuum growth equation results we find that none can quantitatively explain all our findings. Golubovic  predicts $`\beta =0.5`$, which is consistent with our asymptotic result ($`\beta _{II}0.5`$), but his finding of $`n=\lambda =1/4`$ in both d=1+1 and 2+1 is inconsistent with our asymptotic results (asymptotically $`n<0.1`$, and $`\lambda 0.40.5`$) while being approximately consistent with our initial transient results (for $`t<t^{}`$). The analytic results of Rost and Krug  also cannot explain our results, because they predict, in agreement with Golubovic, that if $`\beta =1/2`$, then $`n=\lambda =1/4`$. We also find our asymptotic $`\beta `$ to be essentially $`0.5`$ independent of the actual value of $`P_U/P_L`$, which disagrees with ref. 19. Interestingly our initial transient regime is, in fact, approximately consistent with the theory of ref. 19. We have approximate slope selection ($`\lambda <0.1`$) only in the initial transient regime, which could, however, be of considerable experimental relevance because the pre-asymptotic regime is a long-lasting transient. The only prior work in the literature that has some similarity to our results is that of Villain and collaborators , who in a one dimensional deterministic (i.e. without the deposition beam shot noise) macroscopic continuum description of growth roughness and instabilities under surface diffusion bias, called the Zeno model  by the authors, found, in agreement with our atomistic two dimensional stochastic cellular automaton simulation results, a scenario in which coarsening becomes “extremely slow after the mounds have reached a (characteristic) radius” , which is reminiscent of our crossover from weak bias ($`t<t^{}`$) to strong bias ($`t>t^{}`$). While the precise relationship between our two (and one) dimensional atomistic/stochastic model and their  one dimensional macroscopic/deterministic model is unclear at the present time, it is interesting to note that the authors of refs.  came to a similar negative conclusion as we do about the non-existence of any continuum growth equation describing the strong bias asymptotic regime where the microscopic lattice size may play a crucial role . Finally, we note that a very recent experimental work  reports growth morphology of Ge(001) surface which agrees qualitatively with the scenario predicted in this paper, namely, that even for a very weak diffusion bias, the mound slope continues to increase without any observable coarsening.

Although our diffusion bias model is an extremely simple limited mobility model which may be viewed as unrealistic, our d=1+1 results are remarkably similar to those obtained from a study of a full temperature-dependent Arrhenius hopping model with step edge barrier . This study  offered essentially the same picture as what we present here, i.e. a smaller value of the growth exponent that crosses over to $`\beta 0.5`$ at larger time, and a very large dynamical exponent corresponding to very little coarsening ($`n0`$) after approximately 100 monolayers of deposition. To be specific, the dynamical exponent determined by the growth of the correlation length is $`z16.6`$ in the temperature-dependent edge bias model  while our study yields $`z10`$. Although the dynamical exponents from the two studies are not exactly the same, they are both exceptionally large indicating the mound coarsening process to be negligible in the large time regime. This qualitative agreement between our simple limited mobility results and a d=1+1 full diffusion results argue strongly in favor of our minimal growth model being of reasonably general qualitative validity in experimental situations.

Finally, we note that the introduction of limited mobility models  has opened a whole new way of studying kinetic surface roughening in molecular beam epitaxial growth. These limited mobility nonequilibrium models make it possible to study very large systems in very large time limit when it is impossible to do so in the realistic temperature-dependent full-activated diffusion models. Particularly, we point out the success of the DT model  in providing an excellent zeroth order description of molecular beam epitaxial growth in the absence of any surface diffusion bias. The d=2+1 critical exponents  in the unbiased (DT) model , which belongs to the same universality class as the conserved fourth-order nonlinear continuum MBE growth equation , are $`\beta =0.250.2`$ and $`\alpha 0.60.7`$, which are in quantitative agreement with a number of experimental measurements  where surface diffusion bias is thought to be dynamically unimportant. With this in mind, it is, therefore, conceivable that our study of the limited mobility model which includes surface diffusion bias (i.e. a generalized version of the DT  model) presented in this paper may benefit the subject in a way similar to what the DT model did for the unbiased molecular beam epitaxy growth study. This is particularly significant since there are still many open questions regarding the interface growth under surface diffusion bias condition. Since the model we study here is a generalized DT model  and an approximate continuum description for the original unbiased model  has recently been developed , one could use that as the starting point to construct a continuum growth model for the biased growth situation. Such a continuum description is, however, extremely complex  as it requires the existence of an infinite number of nonlinear terms in the growth equation, and it therefore remains unclear whether a meaningful continuum description for our discrete simulation results is indeed possible .

In conclusion, we want to emphasize the fact that the limited mobility model studied in this paper is an extremely simple model (“the minimal model” in the sense that it is perhaps the simplest nonequilibrium model which captures the minimal features of growth under surface diffusion bias), and realistic growth under experimental conditions should be substantially more complex than this minimal model (we speculate that this minimal model is in the same growth universality class as realistic growth under a surface diffusion bias for reasons discussed above, but we certainly cannot prove that at this incomplete stage of the development of the subject. A word of caution is in order in comparing experimental results with our calculated critical exponents because of the extreme simplicity and the limited mobility nature of our nonequilibrium growth model.) If we do not have a reasonable theoretical understanding (from a continuum equation approach) of even such a simple minimal model, as we have found in this paper, then current efforts  at understanding realistic growth under surface diffusion bias must be quite futile. The main weakness of limited mobility model (of the type presented here) is that they are manifestly nonequilibrium model — this, however, should not be a particularly serious problem in the context of the mound/pyramid formation in the growth morphology, which by definition is a nonequilibrium effect and must disappear in a properly equilibrated surface. Our conclusion based on the results presented in this Letter is that a continuum growth equation for nonequilibrium growth under a surface diffusion bias does not exist at the present time.

This work is supported by the US-ONR and the NSF-DMR-MRSEC.