# Ehrlich–Schwoebel barrier controlled slope selection in epitaxial growth

## 1 Introduction

Molecular beam epitaxy (MBE) has attracted much interest from both, theoretical and experimental physicists. On the one hand it allows the fabrication of high quality crystals with arbitrary composition and modulated structures with atomically controlled thickness . On the other hand it represents a model of nonequilibrium physics which still lacks a general theory . In particular, the appearance and the dynamics of three dimensional (3D) structures (pyramids or mounds) in crystal growth are not well understood in terms of the underlying microscopic processes.

A long time ago Burton, Cabrera and Frank introduced the BCF–theory of crystal growth . Within this theoretical approach the crystal surface is described by steps of single monolayer height. The evolution of the surface is calculated by solving the diffusion equation on each terrace. Within this framework the growth of spirals and the step flow has been investigated. Elkinani and Villain investigated such a model including the nucleation probability of new islands . They found that the resulting structures are unstable. Towers appear which keep their lateral extension and grow in height only. They called this effect the Zeno–effect. The same observation has been made with a “minimal model” of MBE where fast diffusion together with a high Ehrlich–Schwoebel barrier has been implemented .

Even though the Zeno–effect has been observed recently on Pt(111) , quite typically a coarsening process with appearance of slope selection emerges which has been reported for such diverse systems as Fe(001) , Cu(001) , GaAs(001) , and HgTe(001) . In addition, slope selection seems to be the generic case of solid–on–solid computer simulations .

In terms of continuum equations the selection of a stable slope has been related to the compensation of uphill and downhill currents . An uphill current can be generated by an Ehrlich–Schwoebel barrier . The barrier hinders adatoms to jump down a step edge. Hence, more particles attach to the upper step edge which leads to a growth–instability  and 3D–growth.

Another process, which constitutes a downhill current, has been recognized using molecular dynamics simulations . Such diverse mechanisms as downward funneling, transient diffusion or a knockout process at step edges lead to the incorporation of arriving particles at the lower side of the step edge. In addition, it has been suggested that such a process is responsible for reentrant layer–by–layer growth 

Recently we have proposed a simplified model of epitaxial growth quite similar to the “minimal model” of Krug . In particular we found that an incorporation mechanism is crucial to achieve slope selection. However, one simplifying assumption of the model is an infinite Ehrlich–Schwoebel barrier.

In this article we will present in more detail the argument leading to slope selection and we will generalize our results using a continuous step dynamics model analogous to . In sec. 2 we will introduce our extension of the BCF theory and will discuss the relation to existing results (sec. 2 and 3). Typical mound morphologies and the growth dynamics are compared in section 4. Afterwards we will investigate the emergence of slope selection within the framework of this model (sec. 5). We will show that the selected slope has a temperature–dependence which is solely determined by the Ehrlich–Schwoebel barrier. Hence, the determination of selected terrace widths in experiments would give direct insight into microscopic properties such as the Ehrlich–Schwoebel barrier. We confirm the predicted importance of the incorporation mechanism using a kinetic Monte–Carlo simulation of a Solid–On–Solid model in sec. 6. Another effective downward current could be due to detachment from steps and subsequent desorption. We will show in sec. 7 that slope selection cannot be achieved by these two processes alone. In section 8 we will calculate the saturation profile in the limiting case of an infinite Ehrlich–Schwoebel barrier.

## 2 BCF theory

The model is based on the Burton–Cabrera–Frank model in 1+1 dimensions. Within this framework the crystal surface is specified by the position and direction (upward or downward) of steps. Figure 1 shows the crystal surface from the point of view of the BCF–theory. It is a coarse grained view – the detailed positions of atoms are not important. However, the terraces of the height of one atomic monolayer (ML) can still be distinguished. The most fundamental assumption is that at each time $`t`$ the adatom concentration $`\rho `$ is a function of the step positions only. In other words, the diffusion of adatoms is considerably faster than the step velocity. Thus, the diffusion equation becomes

$$\frac{\rho }{t}(x,t)=0=D^2\rho (x,t)+\frac{F}{a}$$
(1)   
where $`D`$ is the diffusion constant and $`F/a`$ is the flux density with $`a`$ denoting the lattice constant. Hence, $`1/F`$ is the time necessary in order to deposit one monolayer. Up to now, this equation was solved with special boundary conditions at $`x=\mathrm{}/2`$ and $`+\mathrm{}/2`$ in the literature. Clearly, the boundary conditions are chosen depending on whether the terrace is a vicinal, a top, or a bottom terrace. In the following we will discuss the typical case of a vicinal terrace. The extension to top and bottom terraces is straightforward.

To include an incorporation mechanism it is necessary to extend the theory. We assume that there exists an incorporation radius such that all particles arriving close to a downward step within this radius immediately jump down the step edge. Hence, one has to split the density of diffusing particles into two regions. The first region close to the upper edge where eq. (1) holds, and the second one given by the incorporation radius close to the downward step where no particles arrive ($`F=0`$). To describe the motion of steps the flux of incorporated particles must be taken into account separately.

In the following we will discuss in detail the situation $`\mathrm{}>R_{\text{inc}}`$ as sketched in fig. 1. For smaller terraces only one region exist and the calculations are much easier. Since our analytical calculations will show that $`\mathrm{}>R_{\text{inc}}`$ is the generic case we concentrate on this situation.

The general one-dimensional solution of eq. (1) is a parabola characterized by three parameters. In addition to the two diffusion equations, four boundary conditions are necessary to determine the two distributions $`\rho _1`$ and $`\rho _2`$ .

$`\rho _1(\mathrm{}/2)`$ $`=`$ $`0`$ (2)   
$`\rho _1(\mathrm{}/2R_{\text{inc}})`$ $`=`$ $`\rho _2(\mathrm{}/2R_{\text{inc}})`$ (3)   
$`\rho _1^{}(\mathrm{}/2R_{\text{inc}})`$ $`=`$ $`\rho _2^{}(\mathrm{}/2R_{\text{inc}})`$ (4)   
$`D\rho _2^{}(\mathrm{}/2)`$ $`=`$ $`{\displaystyle \frac{D}{\mathrm{}_1}}\rho _2(\mathrm{}/2)`$ (5)   
Condition (2) is for the special case of perfect absorbing step edges. (3) and (4) are necessary to obtain a smooth density between region 1 and 2. The left hand side of (5) is the particle current at the step edge. On the right hand side this is reformulated using the number of jump attempts $`D\rho _2(\mathrm{}/2)`$ multiplied by the probability of overcoming the Ehrlich–Schwoebel barrier $`E_S`$. This probability is expressed as the inverse of a typical length $`\mathrm{}_1`$

$$\frac{1}{\mathrm{}_1}=\frac{1}{a}\mathrm{exp}\left(\frac{E_S}{k_\text{B}T}\right)$$
(6)   
where $`a`$ stands for the lattice constant.

The resulting density distribution has the form as indicated in fig. 1: a parabola in the upper region and linear close to the downward step. The detailed expressions of $`\rho _1`$ and $`\rho _2`$ are not of much interest since the evolution of the crystal is determined by the currents at the edges. In the following we will call $`u(\mathrm{})`$ the upward current, i.e.

$`u(\mathrm{})`$ $`=`$ $`D\rho _1^{}(\mathrm{}/2)`$ (7)   
$`=`$ $`{\displaystyle \frac{F}{2a(\mathrm{}+\mathrm{}_1)}}\left(\mathrm{}^2+2\mathrm{}\mathrm{}_12R_{\text{inc}}\mathrm{}_1R_{\text{inc}}^{}{}_{}{}^{2}\right).`$   
The downward current due to diffusion (the contribution of the incorporation mechanism is not included) becomes

$`d(\mathrm{})`$ $`=`$ $`D\rho _2^{}(+\mathrm{}/2)`$ (8)   
$`=`$ $`{\displaystyle \frac{F}{2a(\mathrm{}+\mathrm{}_1)}}\left(\mathrm{}R_{\text{inc}}\right)^2.`$   
Note, that these results are very similar to the corresponding equations (2.2) and (2.3) of reference  where no incorporation was considered. Setting $`R_{\text{inc}}=0`$ we regain their results.

The absence of a dependence on $`D`$ reflects the ansatz of a quasi–stationary distribution. All arriving particles are compensated for by the loss of particles at the borders and hence the currents are proportional to $`F`$. The density itself is proportional to the ratio $`F/D`$ which again is intuitively clear.

## 3 Closure of bottom terraces

In the following, we will reinvestigate the discussion of  concerning the closure of a bottom terrace (c.f. fig. 2). In the limiting case of an infinite Ehrlich–Schwoebel barrier the dynamics of the steps become very simple. We denote by $`x(t)`$ the position of the right step which of course will depend on the time $`t`$. The origin is chosen to be in the middle of the bottom terrace. Due to the infinite Ehrlich–Schwoebel barrier the movement of the right and the left step are symmetric. The evolution is then described by $`\dot{x}(t)=Fx(t)FR_{\text{inc}}`$. The first term corresponds to the particles which do fall on the bottom terrace and diffuse to the right. The second term is the contribution of particles which are incorporated from the step above (which is valid as long as the bottom terrace is more than a distance $`R_{\text{inc}}`$ away from a top terrace). As a result $`x(t)`$ evolves as

$$x(t)=(x_0+R_{\text{inc}})\mathrm{exp}\left(Ft\right)R_{\text{inc}}.$$
(9)   
As long as $`R_{\text{inc}}>0`$ there exists a closure time

$$t_c=\frac{1}{F}\mathrm{ln}\left(\frac{x_0+R_{\text{inc}}}{R_{\text{inc}}}\right).$$
(10)   
Without an incorporation mechanism ($`R_{\text{inc}}`$=0) the bottom terrace will never be closed. This is the reason why Elkinani and Villain called their model the Zeno–model to remind the greek philosopher and his paradox. Even though the situation is changed if the discrete structure of the terraces is considered<sup>1</sup><sup>1</sup>1The currents can be translated into probabilities of placing a particle at the step edge. Hence, a bottom terrace of width one always has a nonvanishing probability to be filled. they showed that this trend still holds which gives rise to the formation of deep cracks. Likewise they found that even finite values of the Ehrlich–Schwoebel barrier do not change this growth scenario which has been investigated in more detail in . Once mounds are built up they remain forever with a fixed lateral size. Our discussion of this limiting case shows that the inclusion of an incorporation mechanism changes the growth in a fundamental manner.

## 4 Growth dynamics

To set up the basic ideas of the behaviour during crystal growth we show two typical surface profiles according to the numerical integration of the step system. In fig. 3 we compare the resulting structure of the Zeno model  without an incorporation mechanism and with the inclusion of such a mechanism.

The simulations were carried out on on a system of $`485a`$ width with parameters corresponding to the model used in sec. 6.

$`D`$ $`=`$ $`10^{12}\mathrm{exp}\left({\displaystyle \frac{0.9\text{eV}}{550\text{K}k_\text{B}}}\right){\displaystyle \frac{a^2}{\text{s}}}5664{\displaystyle \frac{a^2}{\text{s}}}`$   
$`\mathrm{}_1`$ $`=`$ $`\mathrm{exp}\left(+{\displaystyle \frac{0.1\text{eV}}{550\text{K}k_\text{B}}}\right)a8.2a`$   
$`R_{\text{inc}}`$ $`=`$ $`1a`$   
$`F`$ $`=`$ $`\text{1 ML s}^1`$   
As in  the Ehrlich–Schwoebel barrier has been suppressed for bottom terraces of one lattice constant width. Without an additional incorporation mechanism the appearance of trenches is unavoidable in accordance to . The incorporation mechanism gives rise to a well defined slope which does not change with time. Another fundamental difference is the coarsening behaviour. Without an incorporation mechanism the trenches are stable and the number of mounds remains constant. The additional incorporation mechanism leads to a coarsening behaviour.

In lattice models as well as for continuum equations the coarsening is driven by fluctuations  and in 1+1 dimensions the corresponding exponent is 1/3. This is in accordance to Ostwald-ripening which has been predicted from the similarities of the relevant continuum equations . However, since we treat the step evolution in a deterministic manner we do not obtain a scaling behaviour. The only way fluctuations come into play during the simulation is when new islands are nucleated. As a consequence the evolution of e.g. the width of the height distribution $`w`$ is characterized by jumps (data not shown). A jump in $`w`$ appears each time when two mounds merge. These findings are a direct confirmation of the relevance of the fluctuations for the coarsening behaviour.

## 5 Slope selection

The inclusion of an incorporation mechanism leads to slope selection which is apparent from fig. 3b. Siegert and Plischke  required a cancelation of upward and downward currents in the continuum equations. Again, in the case of an infinite Ehrlich–Schwoebel barrier the calculations are straightforward. In this case the downward current on a vicinal terrace of size $`\mathrm{}`$ is solely due to the incorporation mechanism, i.e. proportional to $`FR_{\text{inc}}`$. All the remaining diffusing adatoms will contribute to the upward current and hence the current will be $`F(\mathrm{}R_{\text{inc}})`$. As a consequence the slope selection will be achieved with a mean terrace width of size

$$\mathrm{}^{}=2R_{\text{inc}}$$
(11)   
in accordance to the findings in .

It remains to calculate the terrace widths for arbitrary parameters. Since we know the currents (equations (7), (8) and the incorporation mechanism) we obtain the overall slope (resp. terrace width) dependent current

$`J(\mathrm{})`$ $`=`$ $`u(\mathrm{})+d(\mathrm{})+FR_{\text{inc}}`$ (12)   
$`=`$ $`{\displaystyle \frac{F}{2(\mathrm{}+\mathrm{}_1)}}\left(2R_{\text{inc}}^{}{}_{}{}^{2}+4R_{\text{inc}}\mathrm{}_12\mathrm{}\mathrm{}_1\right)`$   
Note that a positive $`J(\mathrm{})`$ signifies a downward current (to the right in fig. 1). The stable slope, where no net upward or downward current remains, is given by the condition $`J(\mathrm{})|_\mathrm{}=\mathrm{}^{}=0`$ and yields

$`\mathrm{}^{}`$ $`=`$ $`2R_{\text{inc}}+{\displaystyle \frac{R_{\text{inc}}^{}{}_{}{}^{2}}{\mathrm{}_1}}`$ (13)   
$`=`$ $`2R_{\text{inc}}+{\displaystyle \frac{R_{\text{inc}}^{}{}_{}{}^{2}}{a}}\text{e}^{E_S/k_\text{B}T}`$   
As can be seen from expression (12) the current is positive for small values of $`\mathrm{}`$ and becomes negative for $`\mathrm{}>\mathrm{}^{}`$. Hence $`\mathrm{}^{}`$ is stabilized by the current.

The stable slope does not depend on the diffusion constant. However, it should be clear from the derivation, that in order to achieve slope selection the typical diffusion length should be much larger than $`\mathrm{}^{}`$. Otherwise the vicinal terraces would not proceed via step flow. Rather new nucleation events on the terraces would lead to a rugged surface structure.

## 6 Solid–On–Solid model

In order to verify the predicted importance of the incorporation mechanism we use computer simulations of the Solid–On–Solid (SOS) model on a simple cubic lattice. All processes on the surface are (Arrhenius-) activated processes which are described by one common prefactor $`\nu _0`$ and an activation energy which is parameterized as follows: $`E_B`$ is the barrier for surface diffusion, at step edges an Ehrlich–Schwoebel barrier $`E_S`$ is added. However, this barrier is not added for a particle which sits on top of a single particle or a row . Each next neighbour contributes $`E_N`$ to the activation energy. Within this framework the diffusion constant becomes $`D=\nu _0\mathrm{exp}\left(E_B/k_\text{B}T\right)`$.

Here, we concentrate on a particular set of parameters even though other parameter sets were used as well. We choose $`\nu _0=10^{12}s^1`$, $`E_B=0.9\text{eV}`$, $`E_N=0.25\text{eV}`$, and $`E_S=0.1\text{eV}`$. This model was already investigated in  and reproduces some kinetic features of CdTe(001). The deposition of particles occurs with a rate $`F`$. The two simulations shown in fig. 4 are carried out on a $`300\times 300`$ lattice at 560 K and started on a singular (flat) surface.

In fig. 4 the resulting surfaces with and without the inclusion of the incorporation mechanism are shown. Without an incorporation mechanism no slope selection occurs. Clearly, without incorporation the configuration of the towers remains unchanged whereas the inclusion leads to coarsening. The number of mounds diminishes with time. Hence, without an incorporation mechanism no coarsening can be identified. We want to mention that it seems that at higher temperatures the attachment/detachment kinetics of atoms at step edges yields a coarsening effect (data not shown). However, still no slope selection has been observed.

At first glance our findings contradict previous results obtained with a very similar model. Šmilauer and Vvedensky obtained a formation of mounds with slope selection irrespective of the inclusion or exclusion of an incorporation mechanism . However, they implemented the Ehrlich–Schwoebel barrier in a different way. Rather than to hinder the jump over a step edge they impede the jump towards a step edge. Their motivation for this implementation was to allow the adatoms to leave a small line of particles of width one which has been tested as a cause for reentrant layer-by-layer growth . In our simulations the same goal is achieved by suppressing the Ehrlich–Schwoebel barrier in such a situation. However, in their simulations particles arriving directly at a step edge have a probability of 1/4 to jump down the edge, 1/4 to jump away from the edge and 1/2 to jump along the step edge. Effectively this leads to an incorporation radius of length 1/2.

Other simulations of SOS-models used bcc(001)  in order to study the growth of typical metals. In these simulations the SOS–restriction is implemented in such a way that an adatom must be supported by the four underlying atoms. Hence, the downward funneling process is directly implemented. Again, as a result slope selection is achieved, which has already been discussed in great detail in .

## 7 Detachment and desorption

One might assume that other mechanisms could lead to a zero in the slope dependent current. In the following we will carry out an analogous calculation with an adatom–detachment rate from steps and inclusion of desorption . One might assume that both processes generate an effective downward current which can compensate for the Ehrlich–Schwoebel effect. To investigate whether they are sufficient to obtain slope selection (and to simplify notation) we exclude the incorporation mechanism. Thus, the distinction of the two regions on a terrace is not necessary.

The desorption of diffusing adatoms is easily incorporated including a term $`\rho (x)/\tau `$ in the diffusion equation (1) . In order to include detachment from steps we have to replace boundary condition (2) by

$$D\rho ^{}(\mathrm{}/2)=\gamma \frac{D}{a}\rho (\mathrm{}/2)$$
(14)   
where $`\gamma `$ stands for the detachment rate from steps. Accordingly, the boundary condition at downward step has to be corrected and reads now

$$D\rho ^{}(\mathrm{}/2)=\frac{D}{\mathrm{}_1}\rho (\mathrm{}/2)\gamma \frac{a}{\mathrm{}_1}.$$
(15)   
The overall slope dependent current becomes

$$J(\mathrm{})=\frac{\left(\mathrm{\Delta }1\right)\left(\mathrm{}_1a\right)\left(\frac{a\gamma }{\tau }DF\right)}{(\mathrm{}_1+a)\sqrt{\frac{D}{\tau }}(\mathrm{\Delta }+1)+\frac{a\mathrm{}_1}{\tau }(\mathrm{\Delta }1)+D(\mathrm{\Delta }1)}$$
(16)   
where

$$\mathrm{\Delta }=\text{e}^{2\mathrm{}/\sqrt{D\tau }}$$

has been introduced. Note that $`\mathrm{\Delta }`$ is always greater than one.

To discuss the qualitative behaviour it is sufficient to look at the numerator of $`J(\mathrm{})`$ (the denominator is always positive). The first important result is that no slope selection is possible. Only for $`\mathrm{}=0`$ the current is zero (of course, there is no upward and downward current as well).

Even though there is no slope selection, one can discuss whether growth will proceed via layer–by–layer growth ($`J(\mathrm{})>0`$ for all $`\mathrm{}`$, i.e. terraces tend to grow larger) or a growth instability is present ($`J(\mathrm{})<0`$ for all $`\mathrm{}`$, i.e. particles are preferably incorporated at the upper steps).

In the well known limit of negligible detachment or desorption rates  ($`\gamma 0`$ or $`\tau \mathrm{}`$) the Ehrlich–Schwoebel effect alone determines the sign of $`J(\mathrm{})`$. As expected, for positive step edge barriers ($`\mathrm{}_1>a`$) growth becomes instable whereas negative values of $`E_S`$ stabilize layer–by-layer growth.

If we assume a positive $`E_S`$ in the general case we obtain a critical current

$$F_C=\frac{a\gamma }{D\tau }$$
(17)   
where the current changes its sign.

If one expresses the diverse rates as used in the Solid–on–Solid simulations and sets $`a=1`$ one obtains

$$F_C=\nu _d\text{e}^{(E_D+E_{\text{bind}})/k_\text{B}T}.$$
(18)   
where the desorption rate $`\nu _d\text{e}^{E_D/k_\text{B}T}`$ has been introduced. $`E_{\text{bind}}`$ represents the typical binding energy of a detaching adatom. In 2+1 dimensions this should be approximately $`E_{\text{bind}}2E_N`$. Using the model of the previous section, $`\nu _d=\nu _0`$, and $`E_D=1.1\text{eV}`$ (parameters which are a reasonable guess for CdTe(001), ) one obtains a critical flux $`F_C=0.004`$ ML/s. However, this crossover should not be observable in experiments, since at such low external fluxes the step flow of the preexisting steps will dominate the surface evolution.

## 8 Saturation profile with infinite step edge barrier

In this section we will calculate the saturation profile for the model with infinite step edge barrier. The discussion of the closure of the bottom terrace already showed that in this limit the calculations become very simple. As for the bottom terrace the dynamics of higher steps become independent of the above lying terrace. The steps $`x_i(t)`$ evolve according to

$`{\displaystyle \frac{\text{d}x_1}{\text{d}t}}(t)`$ $`=`$ $`F(x_1(t)+R_{\text{inc}})`$   
$`{\displaystyle \frac{\text{d}x_2}{\text{d}t}}(t)`$ $`=`$ $`F(x_2(t)x_1(t))`$   
$`\mathrm{}`$   
$`{\displaystyle \frac{\text{d}x_i}{\text{d}t}}(t)`$ $`=`$ $`F(x_i(t)x_{i1}(t))`$   
$`\mathrm{}`$   
Measuring the time in units of $`1/F`$ (i.e. setting $`F=1`$ in the above equations) the time to grow one monolayer is equal to one. In addition, to simplify notation we will measure all lengths in units of $`R_{\text{inc}}`$.

The solution for the lowest terrace has been given in eq. (9). The solution for all steps is

$$x_n(t)=\underset{i=1}{\overset{n}{}}\frac{t^{ni}}{(ni)!}(1+x_i(0))\text{e}^t1$$
(20)   
as can be easily verified. If we want to calculate the steady state saturation profile we have to require

$$x_{i+1}(1)=x_i(0).$$
(21)   
It should be stressed that this is the only assumption: the surface morphology is a self-reproducing structure. If the upper steps $`x_{i+1}`$ after deposition of one monolayer would be greater than $`x_i(0)`$ this would result in a flattening of the surface. Otherwise the slope would become steeper. Using the solution for the bottom terrace (9) we obtain the initial value

$$x_1(0)=\text{e}1$$
(22)   
when we require that the bottom terrace will be closed at time $`t=1`$.

For the upper terraces eq. (20) yields a recursion relation

$$x_n(1)+1=x_{n1}(0)+1=\underset{i=1}{\overset{n}{}}\frac{1+x_i(0)}{(ni)!\text{e}}$$
(23)   
which can be solved as described in the appendix using the generating function. As a result one obtains the initial positions of the steps on an infinite symmetric step profile. Every time the bottom terrace is closed the steps (with new indices) are located at these positions.

In table 1 we show the analytical expressions for the step positions derived from the generating function. In addition, the numerical values of the terrace widths are shown. With growing index the terrace widths are rapidly approaching two. Even though they oscillate around this value it can be shown that

$$\underset{i\mathrm{}}{lim}\left(x_i(0)x_{i1}(0)\right)=2.$$
(24)   
Note that we measure the lengths in units of $`R_{\text{inc}}`$. Hence we predict a slope selection with slope $`1/(2R_{\text{inc}})`$ as derived from the simple argument of section 5.

The derivation shows that the selected slope is controlled by the closure of the bottom terrace only. Another length scale is the nucleation length. It is defined by the typical length of a top terrace at which nucleation of an island occurs . This length scale is responsible for the rounding of the towers in fig. 3b and causes a perturbation of the steady–state saturation profile.

## 9 Conclusion

We have investigated the effect of an incorporation mechanism on the morphology of growing surfaces. The inclusion of an incorporation mechanism in a 1+1 dimensional BCF–theory as well as in SOS computer simulation in 2+1 dimensions is necessary in order to obtain slope selection and a coarsening process. We were able to derive analytically the temperature dependence of the selected slope. We found that the Ehrlich–Schwoebel barrier alone controls the temperature dependence. In the limit of an infinite step edge barrier we derived the steady state saturation profile. In this case the resulting mound morphology is controlled by the closure of the bottom terrace.

## Appendix A Generating function

To simplify notation we introduce the shifted step positions $`b_j=x_j(0)+1`$. We will try to extract the generating function

$$f(z)=\underset{j=0}{\overset{\mathrm{}}{}}b_jz^j$$
(25)   
for the shifted step positions. Clearly, the $`b_j`$ are only of physical meaning if $`j>0`$ and $`b_0`$ can be chosen arbitrarily. Starting from equation (23)

$`b_{n1}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{b_i}{(ni)!\text{e}}}\text{ for all }n2`$ (26)   
$`\text{e}zb_{n1}z^{n1}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{b_iz^n}{(ni)!}}`$ (27)   
$`\text{e}z{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}b_mz^m`$ $`=`$ $`{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{b_iz^n}{(ni)!}}`$ (28)   
Choosing $`b_0=0`$ and using $`b_1=\text{e}`$ we arrive at

$`\text{e}zf(z)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{i=0}{\overset{n}{}}}{\displaystyle \frac{b_iz^n}{(ni)!}}\text{e}z`$ (29)   
$`\text{e}zf(z)`$ $`=`$ $`f(z)\text{e}^z\text{e}z`$ (30)   
Thus, we finally obtain

$$f(z)=\frac{z}{\text{e}^{z1}z}.$$
(31)   
The lowest coefficients

$$b_j=\frac{1}{j!}\frac{^jf}{z^j}|_{z=0}$$
(32)   
derived from the generating function are shown in table (1). In addition, the generating function can be used to formally prove equation (24).

\***

This work has been supported by the Deutsche Forschungsgemeinschaft DFG through SFB 410.