Continuum modeling of particle redeposition during ion-beam erosion

We present a detailed analysis of a laterally one-dimensional continuum model for target particle redeposition during ion-beam erosion under normal incidence. We show that the redeposition mechanism can be approximated by a Taylor expansion around the mean surface height and clarify the impact of redeposition on nano-scale patterning. Most importantly, we show that stable well-ordered periodic patterns can result in a common Kuramoto-Sivashinsky-like continuum model for ion-beam erosion in cooperation with the model for particle redeposition and identify parameter ranges where they appear.


Introduction
Although studied for many decades [1][2][3] the theoretical understanding of ion-beam erosion processes on solid surfaces still constitutes a highly active research topic. This is basically caused by the intriguing structure forming processes on the surface of the target that originate from the microscopic complexity of the erosion process on the atomic scale. Generally speaking, in ion-beam erosion an incoming flux of eroding particles, typically rare gas ions, penetrates the target surface, transfers its energy to target particles and eventually causes surface particles of the target to be detached. This gives rise to a persistent outflow of eroded particles and in cooperation with various induced surface restructuring effects, e.g. surface diffusion, to a spatio-temporal evolution of the surface morphology on the nano-or meso-scale.
The primary focus of our study aims for the understanding of the potential impact of outgoing eroded particles on the pattern forming processes at the surface. At first sight, one might argue that all eroded particles will also leave the system. This, however, is only true if the target is perfectly flat and the incoming beam hits the surface in normal direction. As soon as the target surface develops the typical mountainous structure, the trajectories of eroded particles that do not stem from the uppermost crests will be potentially obstructed by the already build-up surface structure. The outcome of this can be a redeposition of at least a part of the eroded par-Supplementary material in the form of one mp4 file available from the Journal web page at http://dx.doi.org/10.1140/epjb/e2013-40555-7 a e-mail: c didd01@uni-muenster.de b e-mail: slinz@uni-muenster.de ticles at generally distant, i.e. not neighboring, parts of the surface or additionally induced secondary sputtering processes. This elucidates the nonlocal character of the redeposition process. Redeposition processes rarely receive any consideration in most of the common theoretical approaches such as continuum models for the evolution of the surface morphology [4] and more microscopic inspired solid-on-solid (SOS) models [5]. Only in the context of the highly regular surface patterns observed on semiconductors [6][7][8][9][10][11], Facsko et al. [12] have suspected that redeposition might be of physical relevance.
A common approach to describe the appearing nanopatterns is a continuum model where the surface height is given by a continuous function H(x, t). Assuming normal ion-incidence, a minimal expansion in terms of spatial derivatives compatible with the symmetry [4] results in a Kuramoto-Sivashinsky equation (KSE) Each term on the r.h.s is related to a physical mechanism: while a 0 < 0 is a constant erosion velocity, the third term a 2 ∇ 4 H (with a 2 < 0) models the surface smoothing due to Mullins diffusion [13] and the nonlinear term a 3 (∇H) 2 stands for the tilt-dependence of the sputter yield to the lowest order [14]. The second term a 1 ∇ 2 H (with a 1 < 0) represents a roughening mechanism which is mandatory for surface patterns to emerge. Regarding the roughening mechanism at normal incidence, there are different experimental conditions to distinguish. On binary compound targets a morphological instability, the Bradley-Harper mechanism [15], is a possible candidate, even if it is questioned by recent MD simulations [16]. Roughening due to chemical interactions has also been proposed [17]. Metallic impurities are a prerequisite for nanodot patterns on Si [18][19][20] and, therefore, there must be an impurity induced roughening mechanism which could be attributed to the deflection of approaching metal atoms. This roughening effect has been successfully identified in the field of metallic thin film vapor deposition on Si [21,22]. For nanopatterns formed by cluster-like ion bombardment on Ge [23], the underlying roughening mechanism is still speculative. In our model, we represent the roughening mechanism by the simple anti-diffusion term a 1 ∇ 2 H (with a 1 < 0).
However, equation (1) is not able to reproduce wellordered dot structures, as it generically exhibits spatiotemporal chaotic solutions [24] when starting from a plain, slightly perturbed initial surface. Therefore, there must be a further mechanism that triggers a long-range order on the surface.
Besides other proposals [12,[25][26][27][28][29][30][31][32][33], the inclusion of a damping term proportional to H−H (withH denoting the mean height), as proposed by Facsko et al. [12], can lead to stable and long-range arranged hexagonal dot structures during the spatio-temporal evolution of this damped variant of equation (1). The theoretical basis for such an extended Kuramoto-Sivashinsky equation, the exact mapping to a local, damped Kuramoto-Sivashinky equation and anisotropic extensions of that equation have been discussed in [29][30][31][32][33]. Physically, the damping term has been interpreted by Facsko et al. as a lowest order approximation of the height evolution stemming from redeposition of eroded material [12]. Of course, the damping term can only be understood as an approximation since there are some differences between the modeling term and the physical process of redeposition: on the one hand, redeposition can only take place between two points of the surface which are connected by a line of sight, whereas the nonlocalityH in the damping term does not take this aspect into account. On the other hand, the damping term is independent of the wavelength of the structures on the surface, while redeposition of particles is expected to vary in some way with the lateral extent of the structures.
Taking for granted at the moment that redeposition processes generically act as a surface smoothing process, an in-depth investigation of redeposition requires the study of two facets: (i) an 'in vitro' examination of the mesoscopic features and effects of redeposition on the nano-scale without invoking other physically present roughening and smoothing mechanisms on the surface. In that case, the transient evolution of the target surface is of fundamental interest and, therefore, indispensably requires to investigate appropriately large initial surface profiles. Restricting the investigation to the case of small surface gradients is insufficient since, on a microscopic level, one has to expect a sort of stochastic roughening [34]. (ii) The other aspect is how redeposition cooperates with the other physically present roughening and smoothing processes and whether redeposition can act as a stabilizing mechanism for periodic structures. Although in that case small rough initial surface profiles can be considered, a restriction to a small gradient approximation of the redeposition mechanism is still inadequate due to the appearance of roughening mechanisms. In our present investigation we will address both aspects (i) and (ii).
To our knowledge, the detailed theoretical modeling and analysis of the redeposition mechanism seems to be a relatively recent subject. On the basis of a one-dimensional solid-on-solid (SOS) block model adapted specifically for an 'in vitro' study of the redeposition process, Anspach and Linz [34] could numerically show that redeposition processes generically share significant features of the damping term in equation (1), the linear decay of the mean height profile and the exponential decay of the rms roughness with time. This holds true for a variety of emission direction probabilities and redeposition probabilities. In a subsequent investigation by these authors [35], focusing on a semi-analytically manageable special case, they showed that the evolution of the surface profile subject to pure redeposition can be represented by In contrast to that, Bradley [36] argued on the basis of a continuum formulation of the pure redeposition process for small height gradients of the surface that redeposition cannot give rise to a linear damping term and cannot lead to periodic patterns with short range hexagonal order if combined with the afore-mentioned KSE. The apparent dissimilarity of the SOS model results and the continuum formulation and the general interest in a detailed understanding of redeposition necessitates a thorough investigation of the continuum formulation. In this study, we develop a generalized version of the continuum formulation of Bradley [36] that is not restricted to small gradients of the surface morphology and particularly focus on the laterally one-dimensional case. This is because we can make contact with the SOS model approach [34,35], are able to derive analytical solutions for special cases and can lay out a general methodology to extract the impact of redeposition on nano-scale patterning. Moreover, an in-depth investigation of the impact of various angular erosion distributions can be achieved. The laterally two-dimensional case will be presented in depth in a subsequent study [37], a brief account is given in [38]. In short, we are able to show that (i) a parametrization of the dynamics of 'in vitro' redeposition leads to the same functional form as in the SOS model 2 and (ii) combining the continuum formulation for redeposition with the standard Kuramoto-Sivashinsky equation can stabilize periodic surface patterns. For the latter case, we are also able to determine parameter constellations where stable periodic patterns can emerge.

Modeling
As mentioned above, we limit ourselves to one lateral dimension x. At a time t the surface height at x is given by H(x, t). In this nomenclature, a general continuum model for the spatio-temporal evolution of H(x, t) for ion-beam erosion under normal incidence that includes the effect of redeposition can be written as In equation (3), the functional F E describes the surface evolution stemming from pure erosion and surface relaxation processes, whereas the contribution from redepositing particles solely enters into the functional F R . ∂ x H is used here as shorthand notation for all spatial derivatives of H which are compatible with the symmetry. Since equation (1) is an established ansatz for a continuum model for the pure erosion process F E , i.e. F R = 0, and we are primarily interested in the redeposition process, we temporarily ignore F E in equation (3) in order to extract the effect of redeposition only. The functional form of F R will be derived in a similar way as it was done by Bradley [36], however, generalized for arbitrary height gradients, and by Smith et al. [39,40] in the context of ion-beam etching. The geometry of the redeposition process is sketched in Figure 1 and illustrates the subsequent derivation of the model for F R . Assuming a homogeneous and temporal constant irradiation of the target, the ion flux J, i.e. the number of ions impinging on the surface per lateral element and per time is constant. Within a time step dt on an infinitesimal interval [x E , x E + dx E ] there are J dx E dt ions hitting the surface. Upon impact, each ion can potentially erode target material. The average number of eroded target atoms per incident ion is given by the sputter yield Y . This quantity at least depends on the local slope of the surface [41]. By assuming normal ion-incidence, the local angle of incidence α ∈ (− π 2 , π 2 ) measured with respect to the surface normal n E can be calculated by using the shorthand notation H E = H(x E , t) and e x as unit vector in x-direction. The exact formula for Y (α) is not relevant at this stage. However, by taking the Bradley-Harper effect into account, Y also depends on the local curvature ∂ 2 x H E of the surface. For simplicity, we neglect this influence for now and investigate its impact subsequently in Section 5.3.
In an infinitesimal time dt at a lateral element dx E there are J Y (α) dx E dt surface atoms that will be eroded. Each of them can leave the surface under a different angle θ being measured with respect to the normal n E . Let f (α, θ)dθ be the probability for an eroded atom to leave the surface between an angle θ and θ + dθ. The angular ejection distribution f is for now an arbitrary function of the angle of incidence α and the leaving angle θ of the sputtered atom. The only requirement on f is that −π/2 f (α, θ)dθ = 1 holds for all incident angles α since the total number of eroded atoms per impinging ion is already entirely given by the sputter yield Y (α).
In order to obtain the number n(x R , x E ) dx R dx E dt of atoms which are eroded at a position in the interval between x E and x E + dx E and redeposited between the redeposition points x R and x R + dx R within a time interval dt, we just have to determine the angular interval [θ, θ + dθ] which encloses the surface element dA R at x R with the vertex x E . Using the abbreviation and dθ can be directly calculated yielding In the case, that there is a line of sight between (x E , H E ) and (x R , H R ) we can calculate the redeposition density n(x R , x E ) from atoms eroded at x E and redepositing at x R at the time t by using the relation n(x R , x E ) dx R = JY (α)f (α, θ)dθ resulting from equations (4) and (5). If, however, there is no line of sight between these two points, i.e. there is another part of the surface between them, the target material eroded at (x E , H E ) cannot reach the point (x R , H R ) and, therefore, n(x R , x E ) = 0 holds. For this purpose, we introduce a visibility function v(x R , x E ) → {0, 1} which returns one if both points see each other and zero when there is another blocking part of the surface between them. At least at this point we must assume that during the time of flight of a target atom between erosion and redeposition the surface height H does not significantly evolve, so that the redeposition happens basically instantaneously after the erosion. Otherwise the temporal evolution of the surface height would influence every single redeposition process and retardation effects would factor into the redeposition functional F R . Additionally, we have neglected effects of erosion caused by the redeposited particles.
The last step of the derivation of the redeposition functional F R is to integrate n(x R , x E ) over all lateral elements dx E to sum up the total amount of redeposited atoms at x R . With the volume of a target atom being V a , the height evolution arising from pure redeposition reads with α and θ given by equations (4) and (5). The resulting functional (6) is the one-dimensional generalization of the model proposed by Bradley [36] for the case of arbitrary, not necessary small surface gradients ∂ x H. It is obvious that equation (6) is highly nonlocal and nonlinear and, therefore, it is far too complicated for a minimal continuum representation of the redeposition mechanism that can easily plugged in a general continuum equation. In the following sections we will analyze F R and show that it can be approximated by a Taylor expansion around the mean heightH: Another interesting quantity is the redeposition ratio η which indicates the number of redeposited atoms in relation to the number of eroded atoms. The total number of eroded atoms per time is determined by whereas the number of redeposited atoms per time is given by Consequently, the redeposition ratio η is simply defined by and can only take on values between zero and one.

Model characteristics 3.1 The visibility function v
In the last chapter we have introduced the visibility func- i.e. the target, and the connecting line the visibility function can be formally expressed as Consequently, the visibility function v fulfils a symmetry holds. This fact can be exploited to increase the performance in the subsequently following numerical analysis. Furthermore, for all s > 0 the visibility v(x R , x E ) is invariant under the substitution of the surface H(x, t) by the scaled surface sH(x, t), implying the invariance of v on the surface amplitude as long as it is not inverted or set to zero. An inversion of the amplitude, i.e. H → −H, causes that all parts of the connecting line L previously located inside the solid afterwards lie above the surface and vice versa. Furthermore, a dilatation or contraction of the whole surface along the lateral dimension has no impact on the values of v as long as x R and x E are transformed as well.

Scale invariance of F R
There is one key property of the redeposition functional F R in equation (6): it is invariant under an aspect ratio conserving transformation where the height and the lateral dimension are scaled by the same factor s, i.e. H → sH and x → sx. While Y (α) and f (α, θ) only depend on angles being conserved by this transformation, dθ scales with s −1 . On the other hand, the integration v(x R , x E ) × . . . × dx E scales with s, so the resulting F R is invariant. This means that the redeposition process acts on a large-scale surface morphology with the same magnitude as on a small-scale one as long as both have same shape and aspect ratio. Numerical evidence for this symmetry has also been noticed in the SOS-model [35]. From the physical point of view, the invariance can be seen as follows: the larger the extent of the surface the more target atoms are overall eroded, but on the other hand their spots of redeposition are wider distributed over the surface. Both factors cancel each other resulting in the invariance of F R . This invariance also holds for the redeposition ratio η introduced in equation (10). When considering the invariance of F R in the evolution equation for pure redeposition, the time has to be transformed as well by t → st to absorb the scaling factor on the l.h.s. In this way, a small-scale morphology relaxes faster under redeposition than the corresponding large-scale one. However, this only holds true for the same absolute redeposition rate acting on both morphologies equally, leading to a faster relative relaxation of the small-scale surface.
The scale invariance has a significant impact on gradient expansions of F R . Especially, when assuming a small relative surface height h = H −H ≈ 0, one is interested in a linearization of F R for small amplitudes: With the coefficients c 2n = const., such an approximation can only fulfill the scale invariance if c 2n = 0 holds for all n. Thus, the redeposition functional F R cannot have a universally valid linear portion with constant coefficients. In contrast, there are several (nonlinear) terms like (∂ x h) 2n , . . satisfying the invariance. However, as it is subsequently shown, for physically reasonable surfaces it is far more accurate to approximate F R by an expansion around the mean height like it is done in equation (7), than an approximation consisting of a combination of terms preserving the scale invariance.
In the derivation of F R the sputter yield Y is restricted to be a function of the incident angle α only. However, by taking the Bradley-Harper effect [15] into account, Y also depends on the local curvature κ = ∂ 2 x H E of the surface. At this point, Y (α, κ) -and with it the whole redeposition mechanism F R -loses its scale invariance, since κ scales with s −1 under the aforementioned aspect ratio conserving transformation. The impact of the curvature-dependence of sputter yield is discussed subsequently in Section 5.3.

Approximation for small gradients
To investigate the limit of an almost flat surface, i.e. h = H −H ≈ 0, we rewrite the relative surface height h as h = A ·ĥ. Hereby,ĥ represents the same surface with unity amplitude, whereas A > 0 denotes the surface amplitude. Together with the invariance of v on A, this allows an expansion of the integrand of equation (6) for small amplitudes: In this limit, α ≈ 0 and θ ≈ ± π 2 holds. Therefore, by using symmetry considerations for Y (α) and f (α, θ) and defining the constants As a consequence, a linear increase of the redeposition yield with the amplitude A can only be expected when eroded atoms can graze the surface tangentially with a non-vanishing probability (i.e. especially f 0 = 0). Otherwise, for small amplitudes the influence of redeposition goes with A 2 and higher order terms. Due to the scale invariance, the approximation can be expressed in terms of the aspect ratio instead of the amplitude. Thus, the approximation is valid for small gradients of the surface. The approximation (15) with f 0 = 0 and f α = 0 (i.e. no erosion along the grazing direction θ = ±π/2) can be cast to a one-dimensional variant of the redeposition model proposed by Bradley [36]: Unlike the two-dimensional model, the denominator in the integrand of (18) is raised to the third instead of the fourth power. However, when plugging a two-dimensional surface of the form H(x, y) = H(x) into the two-dimensional small-gradient approximation by Bradley, the resultant redeposition functional is the r.h.s of (18) times the factor π/2. The latter factor is the consequence of the redeposition contribution stemming from eroded particles at all visible surface positions (

Analytical investigation
Due to the complexity of F R , it is impossible to find a general analytical solution for F R for any given arbitrary surface H(x, t). In order to keep it as simple as possible, we are not interested in the complete temporal evolution of a surface, but only in a formula to calculate F R for a given surface H(x, t) at a fixed time t. In the following, two surfaces H(x) and their corresponding, analytically solvable redeposition distributions F R are discussed. Here, A denotes the surface amplitude given as difference of maximum and minimum height, whereas L is the lateral period length, i.e. (19) Due to the periodicity, the investigation can be restricted to x ∈ − L 2 , L 2 . Furthermore, we introduce the aspect ratio m, the normalized surface heightĥ(x) and the normalized spatial redeposition coordinatex R by The model used for the angular erosion distribution f (α, θ) reads

Triangle wave
One of the most elementary nontrivial surfaces is the following symmetric triangle wave The incident angle is given by the constant For this surface, redeposition can take place between x E and x R (i.e. v(x R , x E ) = 1) if and only if the signs of x R and x E are different. Using the non-approximated redeposition model (6) one obtains As already stated in Section 3.2, the solution is only dependent on the aspect ratio m. Using the inverse function of (22), equation (24) can be regarded as a function of the normalized relative heightĥ. This allows for a Taylor expansion aroundĥ = 0. The three terms entering on the r.h.s of equation (25) are determined by Furthermore, the redeposition ratio, equation (10), yields For comparison, we have also calculated the result of the small gradient approximation (18) which reads The graphs of the coefficients (26) are plotted as function of the aspect ratio m in Figure 2. As can be also seen in the following sections where we analyze arbitrary height profiles, the saturating monotonic increase of the F 0 -curve and the shape of the b-curve are characteristic for the redeposition mechanism. It is clearly evident from the indicated errors that the Taylor expansion around the mean height (25) maintains its accuracy for high aspect  (24), shows the quality of the Taylor approximation even for high aspect ratios. In comparison to this, the small gradient approximation (28) is only valid for shallow surfaces (coefficients indicated by corresponding thin lines). The factor VaJY was set to unity.
ratios, whereas the approximation for small gradients (18) quickly loses its validity. The triangle wave surface (22) is a prime example for the impracticality to approximate F R by terms preserving the scale invariance. All scale invariant terms like (∂ x h) 2n , (22) yield a constant or zero. Approximations in these terms can thereby only produce homogeneous redeposition distributions on this surface. This fact is obviously contradictory to the results in Figure 2 and, thus, a scale invariant approximation by terms consisting of local derivatives cannot successfully reproduce the redeposition mechanism.

Parabolic wave
Another simple but more physically realistic surface is a trigonometric wave. However, the evaluation of the visibility function v for this surface yields a transcendental equation. To avoid this problem, an approximation by a piecewise-defined parabolic surface is considered: .
(29) Since this surface is already too complicated to solve it in the general redeposition model (6), the analysis is confined to the approximated model (18). In the limit of small aspect ratios m one obtains the redeposition distribution The surface (29) and the corresponding redeposition function (30) are depicted in Figure 3. The analytical solution qualitatively shows the typical redeposition distribution, i.e. no redeposition at hilltops, maximum redeposition at the slopes and a local minimum in the sink. The numerical solutions for the general model (6) are of similar shape. Sharp peaks in the redeposition curve are only present for the parabolic surface (29) and are a result of the discontinuous change of the surface curvature. The trigonometric variant reveals a smooth curve with redeposition maxima shifted towards the valley.

Numerical investigation
Although it was possible to solve F R for special cases in the previous section, for most surfaces H(x, t) one is confined to numerical investigations. By discretizing the surface H into N grid points (assuming periodic boundary conditions) with a lateral discretization Δx, the slope ∂ x H and angles α and θ can be calculated by finite differences. The spatial integration in F R is substituted by the corresponding Riemann sum, whereas the value of the visibility function v can determined by checking each height point against the connecting line between (x R , H R ) and (x E , H E ). The accuracy of the numerical routine has been verified by comparison with the analytical solutions (24) and (30). The resulting deviation of below 0.1% (for N = 1000, Δx = 2L N and different aspect ratios m = A L between 0.01 and 100) is definitely within the bounds of numerical accuracy. The factor JV a has been set to unity. This corresponds to a rescaling of time in the equation for pure redeposition ∂ t H = F R .

F R for t = const
First of all, we focus our investigation on a single time step. We are interested in the redeposition yield F R (x R ), for the case that the surface H(x), the angular erosion distribution f (α, θ) and the sputter yield Y (α) are given. We start by setting the surface to a sinusoidal profile In this way, we are able to investigate the impact of amplitude and mode variations on the redeposition mechanism. The scale invariance predicts that F R does not change with A and L as long as the aspect ratio m = A L is held constant. Furthermore, we have to specify Y (α) and f (α, θ). For the latter, we introduce four reasonable versions of the angular erosion distribution f : f skew (α, θ) = 1 2 cos(θ) (1 + sin(θ) sin(α)). (32d) While f const represents an isotropic erosion direction distribution, f cos θ centers the preferred erosion direction around the local surface normal and ensures that eroded atoms cannot graze the surface. f spec models an angular distribution focused around the specular reflection given by α = θ. The model f skew is a cosine distribution skewed to the direction of specular reflection with a vanishing erosion probability for grazing angles θ. Thus, the latter unifies the experimentally observed enhanced erosion peaked into the specular direction [42] with the property  f (α, θ). The inset shows the validity of the Taylor expansion (15) for small gradients, which predicts a linear increase of η for fconst and a quadratic growth for the other angular distributions. For physically reasonable aspect ratios m ≈ 1 more than 50% of all eroded particles are redeposited on the surface. f skew (α, θ = ±π/2) = 0. In each case To get first insights into the numerical results for F R and the influence of the choice of f on it, we set the sputter yield Y (α) = Y const = 1 for the moment.
In Figure 4 the redeposition ratio η is plotted as function of the aspect ratio m. As predicted by the Taylor expansion for small gradients (cf. Sect. 3.3), for the angular distributions with f (0, ±π/2) = 0 the redeposition sets in with m 2 , whereas it rises with m for f const only. In addition, the angular distributions f spec and f skew , with maximum erosion peaked into the specular direction α = θ, exhibit steeper increases of η than the others, since the ejection of eroded atoms is directed more downwards to the surface. We would like to point out that for physically reasonable surfaces the redeposition ratio η can reach values above 50%. As a consequence, the influence of this mechanism should be considered in theoretical approaches to the pattern formation under ion-beam erosion.
However, the redeposition ratio η does not provide any information about the spatial distribution of the redeposited atoms. Some representative numerical results for F R (x R ) with different values of m are shown in the upper graph of Figure 5. Since F R determines directly the rate of change of the height H, one can read off from this figure that maximum redeposition happens in the vicinity of the valleys whereas no redepositing particles arrive at the crests.
It is possible to rearrange the redeposition functional F R as a function of the surface height H. This allows for a fit by a Taylor expansion around the normalized mean heightĥ = 1 A (H −H). As long as m is approximately below 1.5, i.e. especially for   physically reasonable surfaces, this fit (bottom graph in Fig. 5) reproduces the numerical data almost exactly.
We have plotted the fitted coefficients as function of m for the different angular distributions in Figure 6. The qualitative behavior of the coefficient curves is only marginally influenced by the explicit choice of the angular erosion distribution f . While the coefficient F 0 of zeroth order exhibits a monotonic but saturating increase, the linear coefficient b drops down to a minimum at m between 1 and 2, which is followed by a slow asymptotic convergence to zero. The shown fit error is the rms of the deviation, which is quite small for realistic values of m.
Another possibility to get insight into the redeposition mechanism is a cosine decomposition of F R : The coefficientsF n (m) can be calculated from the numerical data by Fourier transformation. Beside the mean rede-positionF 0 , the results in Figure 7 show a strong feedback to the basic modeF 1 , which is proportional to the surface height. Higher harmonics aboveF 2 can be neglected for realistic aspect ratios m. Due to the orthogonality of the modes, the Fourier analysis is another evidence for term proportional toĥ. In addition, the cosine decomposition (34) up to the second harmonicF 2 can be transformed   (α, θ). The plots show that FR mainly consists of the modesF0,F1 andF2. Thus, at least for a cosine surface, the Taylor expansion (33) is a good approximation for the redeposition mechanism.
to the functional form of the Taylor expansion (33) by

Tilt-dependent sputter yield
So far, we have restricted the discussion to a constant sputter yield given by its value for α = 0. To generalize it to a tilt-dependent sputter yield Y (α) = Y tilt (α), we use a sputter yield model derived by Wei et al. [41]: The maximum sputter yield occurs at α = α max . According to experiments, α max typically lies in the interval between 50 • and 70 • . For the comparison with the case Y const. = 1, the constant Y 0 is chosen in a way that Y const. and Y tilt have the same mean value with respect to α. Representative numerical results for the sinusoidal surface (31), f = f cos θ and α max = 65 • are shown in Figure 8. For small aspect ratios m, there is obviously no qualitative influence on Y tilt observable. Furthermore, the graph for the coefficient b is almost independent of the choice of the sputter yield, whereas the graph of F 0 and the modulus of c first increase monotonically until they reach an extremum when m ≈ 1.1 holds. Only for larger m, both curves decay monotonically in contrast to the case of constant sputter yield. We conclude that in the limit of reasonable surface morphologies, as e.g. observed for semiconductor surfaces under low energy ion-beam erosion [2], the tilt-dependent sputter yield has only a quantitative impact on the redeposition mechanism.

Inclusion of the Bradley-Harper effect
Due to the Bradley-Harper effect, the sputter yield Y is not only a function of the local incident angle α, but also of the local surface curvature κ = ∂ 2 x H E . The Bradley-Harper effect can be modeled by separating the sputter yield Y (α, κ) into where the function γ(α) denotes the strength of the Bradley-Harper mechanism [15], i.e. the enhanced sputter yield in sinks with respect to mounds. The advantage of the ansatz (37) is that the redeposition functional now decomposes in the same manner: with describing the curvature-dependent portion of the redeposition mechanism. Since we have already investigated equation (6), we only have to examine F R,κ . We set Y (α) γ(α) to unity to get a qualitative result which does not depend on additional parameters. This simplification, i.e. the evaluation of Y (α) γ(α) at the global and not at the local incidence angle, is valid for sufficiently smooth surfaces [4,15]. In contrast to equation (6), F R,κ may be negative and it violates the aforementioned scale invariance. When m = A/L is held constant, F R,κ is inversely proportional to L and to A, respectively. Therefore, we have plotted representative coefficients of the fitted Taylor expansion F R,κ ≈ F 0,κ + b κĥ + c κĥ 2 multiplied by L against the aspect-ratio m in Figure 9. With growing m, the coefficient b κ decreases, whereas the coefficient c κ increases almost linearly. Thus, the consideration of a curvaturedependent sputter yield enhances the negative b-coefficient and leads to a greater effective c-coefficient. interested in the spatio-temporal evolution of the equation ∂ t H(x, t) = F R (H, x, t), i.e. the relaxation induced by pure redeposition. To that end, we solve this partial integro-differential equation using an explicit Euler method with respect to time t starting from different initial surfaces H(x, 0). Representative results of our investigations are shown in Figure 10, where two generic surface evolutions (upper panels) and the corresponding surface roughness (lower panels) given by are plotted as function of t. From these figures it can be inferred that redeposition is manifestly a relaxation mechanism as it was also observed in the solid-on-solid block model by Anspach and Linz [34]. The decay of the surface roughness w(t) with time is not exactly exponentially as it is in the simple damping term model This deviation, however, is caused by the additional quadratic term c(H −H) 2 and the variation of b and c with the surface amplitude.

Competition of erosion and redeposition
A central question, already posed in the introduction, is whether erosion and redeposition processes can conspire in a way that well-ordered stable periodic structures can emerge. For that purpose we investigate here the onedimensional Kuramoto-Sivashinsky equation combined with the redeposition functional F R , i.e.
Setting F R [H, x, t] to zero, equation (41) exhibits surface profiles that are spatio-temporally chaotic [24] and, therefore, cannot account for steady periodic solutions.
Since we are temporarily interested in a qualitative result, we set a 1 = a 2 = −1 (for example by an appropriate rescaling) and a 0 = 0 (since it is irrelevant for the dynamics). With this parameters, the linear portion − ∂ 2 x + ∂ 4 x H of equation (41) selects the lateral extent of the surface with a dominating wave number k = 1 √ 2 . The typically arising surface height and thereby the aspect ratio can be controlled by the choice of a 3 . The redeposition yield F R is calculated in the same manner as in the previous subsection, i.e. by a full numerical simulation of F R using equation (6). By the choice of the constant of proportionality C = JV a the intensity of the redeposition mechanism can be adjusted.
The upper right panel of Figure 11 shows from top to bottom for six representative times the spatio-temporal evolution of the height (full lines) generated by equation (41) without redeposition, F R = 0, and started at time t = 0 from a weakly rough initial surface. In contrast, the upper left panel of Figure 11 represents the surface evolution for exactly the same parameters and initial conditions except that now redeposition processes have been added, F R = 0. The dashed blue lines represent F R and determine the spatial variations of the magnitude of redeposition. The difference between spatio-temporally chaotic dynamics for F R = 0 and the evolution into a periodic and steady surface profile is striking. Moreover one can also read off the long-time evolution for F R = 0 that (i) there is no redeposition taking place on the top of the surface profile, (ii) in a wide range around the bottom of the surface profile an almost constant redeposition represented by the plateaus of F R is achieved, (iii) right at the minima of the surface profile there are even small local minima of F R indicating a slightly reduced redeposition there. To substantiate our findings, a spatio-temporal visualization of the full numerical simulations leading to the results in Figure 11 is attached as supplementary material.
In the lower graph of Figure 11, we have plotted the temporal evolution of the Taylor coefficients F 0 , b and c during the surface evolution with redeposition. Since the amplitude A of the evolving non-converged surfaces at early times cannot be defined in a proper way, we have fitted the coefficients by i.e. by the non-normalized relative surface height H−H instead of the normalized relative heightĥ = (H −H)/A. As predicted by the small-gradient approximation (18), there is only a quadratic damping term (c -term) present for the flat surfaces at early times. With growing roughness w, the coefficient of the linear damping term b and the quantity F 0 set in and converge with saturating roughness. Moreover, c changes the sign at early times but tends to converge slower than F 0 and b . As a result, redeposition can be neglected at early times but it has to be considered with increasing and saturating surface roughness. In this regime, it is valid to approximate the coefficients F 0 , b and c by constants. The latter can be extracted by fitting (42) with the converged surface (A = 1.179, L = 7.937) at t = 500: These results are in accordance with the coefficients F 0 = 0.0923, b = −0.0698 and c = −0.2389 with respect to the normalized heightĥ for the perfect cosine surface in Figure 6. The primed coefficients can be obtained from Figure 6 by using C = V a J (here: C = 4.0).

Stability maps
The previous section raises the question in which parameter ranges the combined model (41) exhibits ordered patterns. By transforming into a coordinate system comoving with a 0 and taking advantage of the scale invariance of F R , the evolution equation (41) can be rescaled to with the effective parameters The definition of β already contains the factor J V a that has been set to unity in the redeposition functional F R of (45). Starting with a slightly perturbed prestructured wave pattern (16 periods with the critical wave number k c = 1/ √ 2) we have performed a parameter scan in the (β, κ)-plane. After a simulation time of T = 10 000 for each parameter tuple, we have analyzed the resulting surface whether it has retained the ordered pattern or changed over to spatio-temporal chaos. By that, we are able to determine the weakly attracting regions where the effect of redeposition can stabilize prestructured patterns. This gives insight into the dynamics of equation (45) and might be of importance for experiments on prepatterned surfaces. In a second stage, an almost flat random surface was used as initial height field. By determining whether the surface converges to a stable pattern, the globally attracting patterning region can be obtained. As criterion to distinguish a stable ordered pattern from a spatiotemporally chaotic surface evolution, we have applied the 0-1 test for chaos [43] on the surface roughness w(t). Since the 0-1 test cannot directly distinguish between stationary and periodically oscillating solutions, we have additionally checked if the rms of w(t) and all subharmonic modes converge to zero. Periodically oscillating solutions were, however, not found in our investigation.
The resulting stability maps for two representative sets of models for Y and f are depicted in Figure 12. There is typically a weakly attracting region for both positive and negative κ and sufficiently large β 1. Additionally, for κ < 0, a globally attracting region can be found. This region is located at higher redeposition strengths β 1 and small negative κ. In the limit κ → 0, (45) tends to diverge. If f const is used as erosion distribution, only flat solutions can be found at larger β.
The most important result is that, independently of the specific sputter yield and angular erosion model, an extended region in the parameter space can be found where flat surfaces evolve to stable periodic patterns due to the impact of redeposition in (45). The same applies if the small-gradient approximation (18) is used or the Bradley-Harper effect is considered within the redeposition model. In white areas, no pattern formation can be observed. Prestructured wavepatterns are stable against perturbations in dark gray shaded regions, but no ordered patterns arise from a initially flat surface. In the light gray region stable periodic wave patterns evolve from flat initial surfaces. For κ → 0, the surface diverges. Different models for Y and f yield similar stability maps.

Conclusion
The results of our detailed analysis of the continuum model for the redeposition mechanism under ion-beam erosion can be summarized as follows.
The model exhibits a scale invariance which predicts the same absolute redeposition yield on large and small scaled surfaces as long as the aspect ratio is held constant. This fact simplifies the analysis, but also restricts the compatible terms of an approximation. However, as it can be seen from the triangle wave surface (cf. Sect. 4.1), a local approximation in terms preserving the scale invariance cannot be successful. Approximating the redeposition in the limit of small gradients (cf. Eq. (18)), as proposed by Bradley [36], is only valid for very flat surfaces and the resulting model is still very complicated. Of course, redeposition is a complex, nonlocal and nonlinear process, but precisely because its complexity one needs a simple approximation for the inclusion into a theoretical approach. The most basic nonlocality in the model is the mean heightH, suggesting an approximation by a Taylor expansion around this quantity. Even if it is a crude approximation, the expansion up to (H −H) 2 is able to reproduce the redeposition mechanism for generic surfaces. The coefficients are not constant, but depend on the aspect ratio. The explicit models for the sputter yield Y and the angular ejection distribution f , as well as the consideration of the Bradley-Harper effect do not alter the curves for the characteristic coefficients in a qualitative manner. Moreover, this fact, which was already observed in the one-dimensional solid-on-solid redeposition model of Anspach and Linz [34], also holds for the extension to two lateral dimensions (as briefly discussed in [38]): hexagonal dots or holes exhibit similar coefficient curves as in the one-dimensional case, by what our simplification to the one-dimensional problem is justified.
The central point we have shown is that the redeposition model can trigger the emergence of stable periodic structures for a wide range of the entering parameters (cf. Fig. 12) if it is combined with a Kuramoto-Sivashinsky type erosion model. Thus, the redeposition model shows the same stabilizing feature as the inclusion of a damping term proportional to H −H into the Kuramoto-Sivashinsky equation. There is, however, a key difference between the damping term and the presented redeposition model: the damping term already factors into the linearized Kuramoto-Sivashinsky equation, whereas the redeposition mechanism has no linear contribution in the limit of small gradients (as already stated by Bradley [36], except for the case f const ). However, when the surface roughens during erosion, the redeposition mechanism becomes more and more important for pattern formation. The coefficients F 0 , b and c saturate with the roughness after a short transient time. At this stage, the redeposition mechanism can be approximated by damping terms with constant coefficients (which can be read off from the converged final surface) and stabilizes well-ordered patterns. The damping term approximation could be improved by blending the coefficients b and c in with growing roughness w. However, this has no effect on the resulting final surfaces. To conclude, and since there is no sufficiently simple alternative to it, the linear (and the quadratic) damping term can be used to approximatively model the emergence of stable periodic patterns on the surface caused by redeposition during ion-beam erosion.