Ripples and Ripples: from Sandy Deserts to Ion-Sputtered Surfaces

We study the morphological evolution of surfaces during ion sputtering and we compare their dynamical roughening with aeolian ripple formation in sandy deserts. We show that the two phenomena can be described within the same theoretical framework. This approach explains the different dynamical behaviors experimentally observed in metals or in semiconductors and amorphous systems. In the case of ion erosion, we find exponential growth at constant wavelength up to a critical roughness $W_c$. Whereas, in metals, by introducing the contribution of the Erlich-Schwoebel barrier, we find a transition from an exponential growth to a power law evolution.


I. INTRODUCTION
When an ion hits a surface liberates locally a large amount of energy that melts a region of the solid immediately below. For geometrical reasons, the sputtering effect depends on the surface-curvature: the energy concentrates on regions of positive curvature and this favorites the excavation of valleys and the growth of hills. On the other hand, thermal diffusion and surface tension tend to smoothen the irregularities by flattening the surface. It has been observed that under the combined action of these mechanisms the surface tends to create spontaneously ripples [1,2,3,4,5,6,7,8]. In nature, ripples are commonly observed in sandy-deserts as the result of dynamical instability of the sand surface under the action of a sufficiently strong wind [9]. In this case, the formation of ripples is commonly associated with the effect produced by some grains that are lifted from the sand-bed and accelerated by the wind. These grains, when re-impact with the bed, splash up a number of other grains. Most of these grains return to the bed leading to a local rearrangement, whereas some other are accelerated by the wind and impact again on the bed after a certain 'saltation' length. In the literature, many studies have been devoted to understanding the mechanism of ripple formation [10,11,12,13,14,15]. In particular, an hydrodynamical model for aeolian ripple formation, based on a continuum dynamical description with two species of grains (immobile and rolling grains), was proposed with success by Bouchaud et al. [16,17,18,19]. The main ingredient of such a model is a bilinear differential equation, for the population of the two species of grains, which shows the instability of a flat bed against ripple formation.
In this paper we show that the same reasoning which has been used to describe the sand ripples formation in deserts, applied to the studies of dynamical surface roughening, leads to an accurate description of the mor- * Electronic address: tomaso.aste@anu.edu.au phogenesis and evolution of ripples on crystal and amorphous surfaces during ion sputtering. The present approach contains not only the Bradley-Harper approach [20,21] (based on a Kardar-Parisi-Zhang type equation [22]), but it is also able to describe some of the crucial experimental features recently observed in these systems [6]. In particular, by means of this approach we can explain the two distinct dynamical behaviors experimentally observed in amorphous/semiconductors systems and in metals [4]. In the first case (amorphous/semiconductors) we find that the ripples growth exponentially fast at constant wavelengthλ up to a critical roughening W c at which the growing process interrupts. On the other hand, in metals (when the Erlich-Schwoebel barrier is active), we find a transition between an initial exponential to a slower power-law growth of the roughness (the root mean square of the height profile). In this regime the ripple-wavelength tends also to growth with time.

II. PARTICLE MOBILITY AND RIPPLE DYNAMICS
When the surface of a solid is taken under ion sputtering some atoms in the proximity of the surface receive energy from the sputtered ions and pass from a bounded -'immobile' -solid state to a 'mobile' melted unbounded state. The opposite mechanism is also allowed: some mobile atoms can gain in energy by becoming immobile and bounding in a given position in the solid. A certain fraction of atoms might also be dispersed into the atmosphere. Let us call h(r, t) the height of surface profile made of immobile -bounded-atoms and call R(r, t) the height of mobile -melted-atoms. In analogy with the theory developed to explain the dynamical evolutions of dunes in deserts [16,17,18,19], we describe the mechanisms of excavation, exchange between mobile and immobile atoms and surface displacement of mobile atoms in term of the following differential equation: Where Γ(R, h) ex is the rate of atoms that are excavated under the action of the sputtering, and (1 − φ) is the part of them that pass from immobile to mobile, whereas φ is the fraction that is dispersed into the atmosphere.
Let us now write in details the various terms contained in Eq.1.

A. Excavation
The excavation effect must clearly depend on the number and velocity of the sputtered ions (i.e. its flux), but also the local shape and orientation of the surface might play an important role. Indeed, the energy transmitted by the impacting ions concentrate more in regions of the surface with positive curvature. Moreover, part of the surface facing the flux are likely to experience a different erosion respect to others which are less exposed to the flux. Crystalline orientation and anisotropies might be also taken into account. We can write: here η is the sputtering flux, whereas a and b are respectively associated with the flux-direction-dependent and with the curvature-dependent sputtering erosion.

B. Adsorption
The rate of adsorption of mobile atoms into immobile solid positions must be dependent on the quantity of mobile atoms in a given spatial position. Similarly to the excavation process, the adsorption is also dependent on the local curvature and orientation. We can write: where the parameter γ is the recombination rate and c and d are associated to the different probabilities of recombination in relation with the local orientation and shape of the surface.
Note that Eqs. 2, 3 contains the same terms as the ones proposed in the literature for the formation of aeolian dunes in the so-called hydrodynamical model [16,17,18,19,23]. Indeed, in deserts, sand grains are lifted from the sand-bed and readsorbed into it with a probability which is dependent on the local shape and orientation of the dunes. Eqs. 2, 3 represent the simplest analytical expressions which formally take into account these shape and orientation dependences. In the search for simple explanations, such equations are therefore rather universal.

C. Mobility
Mobile atoms will move on the surface, and the quantity J(R, h) in Eq.1 is the 'current' of these atoms. In surface growth, there are two main mechanisms that are commonly indicated as the responsible for the surface mobility of atoms [24]. The first is a current, driven by the variations of the local chemical potential, which tends to smoothen the surface asperity moving atoms from hills to valleys. The second is the current induced by the Erlich-Schwoebel barrier which -on the contrarymoves atoms uphill. In addiction to these main mechanisms we might also have to take into account a drift velocity and a random thermal diffusion, obtaining: (4) In this equation, the first term describes a deterministic diffusion driven by the variations of the chemical potential which depends on the local shape of the surface; the second term is associated with the uphill current due to the Erlich-Schwoebel barrier and α d is a constant associated with the characteristic length. The quantity v is a drift velocity of the mobile atoms on the surface, whereas D is the dispersion constant associated with the random thermal motion.
Note that Eq.4 is substantially different from the one proposed in the literature to describe ripples in granular media [16,17,18,25,26]. Here the current is supposed to be dependent on the local shape and orientation of the surface (the h(r, t) profile). The equations describing sand deserts can be retrieved from Eq.4 by imposing K = 0 and s = 0, but -on the contrary-in surface growth these two parameters are the leading terms of the equation and play the role of control parameters in the dynamics of ripple formation. Nonetheless, these terms describe a rather simple dependence of the dynamic of particles on a surface on the geometrical shape of the surface itself. Again, in our seek for universality, we expect that similar terms can be profitably introduced in the context of aeolian sand ripples in order to describe specific phenomena (associated, for instance, with packing properties [27] or granular flow [28]) which relate the current of grains with the dune-shapes. It should be noted that the factors a, c and v in Eqs.2,3 and 4 are vectors (i.e. they have -in general-different components in the two horizontal directions). Indeed, crystal surfaces are in general anisotropic and therefore one must take into account the dependence of the parameters on the relative orientation of the crystal-surface and the sputtering direction.

III. DISPERSION RELATION
A trivial solution of Eq.1 can be written for a completely flat surface: h(r, t) = h 0 (t) and R(r, t) = R 0 . In this case, we obtain R 0 = (1 − φ)η/γ and h 0 (t) = −φηt+const.. This describes a surface that rests flat and it is eroded with a speed equal to φη. But this behavior is only hypothetical since -in general-the dynamics of the surface-profile presents instabilities against spontaneous roughening and therefore its evolution is more complex. For instance, a numerical solution of Eq.1, is shown in Fig.2 (for the 1-dimensional case). We observe that, in a certain range of the parameters, the surface is unstable and periodic ripples are formed spontaneously.

A. Stability analysis
In order to infer indications about the amplification or the smoothing of small perturbations and to deduce an analytical expression for the ripples wave-length at their beginning, we performe a stability analysis on Eq.1. For this purpose we assume that the surface-profile is made by the combination of a flat term plus a rough part: with R 1 (r, t) =R 1 exp(iωt + ikr) and h 1 (r, t) = h 1 exp(iωt + ikr). We substitute these quantities into Eq.1 and linearize the equation by taking only the first order in R 1 and h 1 . A Fourier analysis (see Appendix A) shows that such a linearized equation admits solutions when the frequencies ω and the wave vectors k satisfy: where, to simplify the equations, we have introduced the following notation: establishes a dispersion relation ω(k) that is a complex function with two branches corresponding to the solutions of the quadratic Equation 6.

IV. SURFACE INSTABILITIES
The kinetic growth of the surface instability is related to the immaginary part of ω(k). Indeed, Im(ω(k)) corresponds to modes with amplitudes that change exponentially fast in time, and negative values correspond to unstable modes that increase with the time. We can therefore study Im(ω) from the solution of Eq.6 and search for the range of k in which Im(ω) is negative. The most unstable mode is the one that grows faster and it corresponds to the value of k at which Im(ω) reaches its most negative value (see Fig.3).
The solution of Eq.6 for Im(ω), is where we have and can assume negative values which are associated with the surface instability (arbitrary units) [29]. The amplitude of modes with wavelengths λ > 2π/k * will grow exponentially fast. The tick line is the imaginary part of the analytical solution of Eq.6, whereas the tinny-gray line is the approximated expression (at the fourth order in k) obtained for small ion flux (η small).
Let us first observe that in absence of sputtering (i.e. when η = 0 and therefore, v 1 = 0, v 2 = 0, D 1 = 0, D 2 = 0, s 1 = 0, K 1 = 0) the solutions of Eq.6 are ω(k) = 0 and ω(k) = −kv + i(γ + k 2 D). In this case, the imaginary part of ω(k) is non-negative, therefore we -correctly-expect no spontaneous corrugation of the surface. On the contrary, when the sputtering is active (η = 0), the immaginary part of ω(k) can assume negative values as shown in Fig.3 where a plot of Im(ω) − is reported (along one direction of the vector k). As one can see, typically the branch Im(ω) − takes negative values for k between 0 and a critical value k * at which it passes the zero [29]. The critical point k * , fixes the minimal unstable wavelength. We therefore expect to find unstable solutions associated with the formation and evolution of ripples with wavelengths λ ≥ λ * = 2π/k * .

V. RIPPLE WAVELENGTH
Several analytical solutions of Eq.6 can be found in special cases which are discussed in Appendix B. But the study of the surface instabilities can be highly simplified if we consider the first order effects when the sputtering flux η is small.

A. Approximate equation
In the case of small spattering fluxes, the branch of Im(ω(k)), with negative values can be approximated to: Im(ω) − ≃ P 1 k 6 + P 2 k 4 + P 3 (k)k 2 + P 4 k 2 + P 5 (k) D 2 k 4 + 2γDk 2 + (vk) 2 + γ 2 (10) with When k is sufficiently small ( k ≪ γ/η ), we can develop Eq.10 at the 4 th order obtaining: with Here v, v 1 and v 2 are respectivelly the components of v, v 1 and v 2 in the direction parallel to k). (In Fig.3 a comparison between this approximate solution and the exact one is given.)

B. Solutions
The expected wavelength of the ripples is associated with the fastest growing mode, which corresponds to the value of k at which Im(ω) − reaches its most negative point. Here the minimum of Im(ω) is at Therefore, at the beginning, the roughness will grow exponentially fast as W ∼ exp(B 2 t/(4A)) with associated ripple-wavelength at:λ C. Some special cases Let us first observe that, when K 1 , s 1 and φ = 0, the ripple wavelength, given by Eq.15, coincides with the one found for sand dunes in deserts (see for instance [18]). In our notation the 'reptation length' is l 0 = v/γ, the 'cut-off length' is l c = (D 2 − D 1 )/v, whereas v 1 − v 2 is the collective drift velocity of the dunes. The approximations usually applied in this context [17,18], imply: l c ≫ D/γ, and γl c ≫ v 1 − v 2 . Giving, from Eq.15 Let us now consider the dynamical evolution of a surface under ion sputtering and in particular the case when the effect of the Erlich-Schwoebel barrier is not present (as for semiconductors and glasses). In this case, s = 0, s 1 = 0 and we also expect that the drift velocity v and the dispersion constant D are equal to zero or infinitesimally small. Indeed, here the current of mobile atoms on the surface is mainly induced by the differences in the chemical potential. Under these assumptions, from Eq.15, the wavelength of the most unstable ripple is: where we called ν = γbφ/(1 − φ), a quantity which plays the role of an effective surface tension. Note that Equation 17 is the same result as from the Bradley and Harper theory [20,21,24,30,31]. When the Erlich-Schwoebel barriers are active (s, s 1 = 0), effects can be observed on the ripple-wavelength at their beginning, which becomes:

VI. EXPONENTIAL/POWER-LAW GROWTH AND CRITICAL ROUGHNESS
In metals, when the Erlich-Schwoebel barrier is active, there is an important non-linear contribution in the current of mobile atoms which becomes sizable when the roughness becomes sufficiently large and therefore (α d ∇h) 2 ∼ 1 (see Eq.4) where the average is over all the surface positions. We observe numerically that this changes the law of growth of ripples: from an exponential to a power-law behavior. This effect is shown in Fig.4. Numerically, all the computed exponents follow in the range between 0.65 and 0.85. The theoretical evaluation of this exponent is under current investigation. In this regime the ripple-wavelength tends also to grow with time. Experimentally, power law growth of the roughness and growth of characteristic wavelengths were observed in erosive sputtering on Ag(001) [4].
In semiconductors or glasses, when no Erlich-Schwoebel barrier is present, it is physically intuitive that the exponential growth of the surface roughness (which is a characteristic of the beginning of the surface instability) cannot continue indefinitely. Indeed, from the expression R(r, t) = R 0 +R 1 exp(iωt + ikr), which we used to derive Eq.15, we can immediately observe that whenR 1 > R 0 = (1 − φ)η/γ, the amount of mobile atoms might become negative. Since a negative amount of atoms is physically impossible, the process of exponential roughness-growth described above must necessarily finish around a critical roughness given by: This behavior is confirmed by numerical solutions of Eq.1 and it is expected to be observable in semiconductors and glasses after sufficiently long times.

VII. CONCLUSIONS
We have shown that the same theoretical approach introduced to describe the formation of aeolian sand ripples can be conveniently applied to the study of the formation of periodic structures on surfaces under ion sputterning. Although the two phenomena are rather different, they can be described within the same conceptual framework by using rather general ideas that relate mobility, excavation and adsorption rates with the surface shape and orientation.
We have obtained general expressions for the ripples wavelength in term of the system parameters. It has been shown that in some particular cases such a solution coincide with the ones already known in the literature for sand dunes and surface instability [17,18,20,21,22]. We have discussed the effect of the Erlich-Schwoebel barriers and compared the result with numerical solutions. We pointed out the Erlich-Schwoebel barrier can be responsible for a dramatic change in the system dynamics: from the exponential growth to a power law. Finally, we have demonstrated that the occurrence of a critical roughness is predictable within the present theoretical framework.
It should be noted that the main purpose of this paper is to point out a relevant example of universality: two processes which have completely different scales present a dynamical evolution which obeys to the same geometrical constraints and thus can be described by using the same phenomenological model. On the other hand, we must observe that the class of solutions of Eq.1 is rich and complex -even in the linear approximation. Exhaustive, systematic studies of the classes of solutions of this equation and their dependence on the set of parameters will be the subject of future studies and publications. By substituting Eqs.5, 4, 3 and 2 into Eq.1 and by neglecting the second order terms (in R 1 and h 1 ), we obtain the following linearized equation: A Fourier analysis of Eq.A1 leads to withR 1 andĥ 1 the Fourier components of R 1 and h 1 respectively. This equation is a simple linear equation in two variables. It admits a non-trivial solution when the determinant of the coefficients is equal to zero. This leads to Eq.6.

APPENDIX B: EXACT SOLUTIONS
Analytical expressions for the value of k at which Im(ω) = 0 (k * ) can be calculated from Eq.7 in some special cases.
In particular, when φ = 0, s 1 = 0, K 1 = 0 and D = 0, we obtain where v, v 1 and v 2 are the components of v, v 1 and v 2 in the direction of k * . On the other hand when, K 1 , s 1 , D 1 and D 2 = 0, we find The effect of the deterministic diffusion induced by the chemical potential can be studied from the solution which holds when v = 0, v 1 = 0, v 2 = 0, s 1 = 0 and D − D 1 + (1 − φ)D 2 > 0. Whereas, when v 1 , v 2 , D 1 and D 2 = 0, we find which implies that the uphill current due to the Erlich-Schwoebel barrier can generate instability even when the shape-dependent erosion and recombination terms are inactive.

APPENDIX C: NUMERICAL SOLUTIONS
The numerical solutions of Eq.1 presented in this paper and in particular the ones shown in Figs.2 and 4 have been performed as follows. We considered a one-dimensional flat substrate (h(x, 0) = h0) of length L, with periodic boundary conditions. An infinitesimal quantity of mobile atoms were added randomly to the substrate (with 0 < r(x, 0) < L10 −10 ). We then computed the profile-evolution using Eq.1 with the derivative substituted with finite differences. To this purpose, substrate has been divided into N discrete points. The -adimensional-time indicated in Fig.4 is the number of numerical steps. The height is in unit of L/N and the roughness is defined as w(t, L) = [h(x ′ , t) − h(x, t) x ] 2 1/2 x ′ (see, for instance, [25]). Several computations with a number of points equal to N = 100, 200 and 300 (the one published have N = 300) have been performed to verify the effect of boundary and discretization. Moreover, simulations with no periodic boundary conditions and with the sputtering term (Eq.2) applied only to a central mask, have also been performed obtaining very similar results. The robustness of the present approach has been verified varying the parameters, the time steps, the initial roughness of the substrate, etc. Comparable results have been always found but, we must stress that, under some conditions, numerical instabilities (in particular small surface-deformations with λ ∼ L/N ) can be trigged on depending on the protocol utilized.