Ion-Induced Nanoscale Ripple Patterns on Si Surfaces: Theory and Experiment

Nanopatterning of solid surfaces by low-energy ion bombardment has received considerable interest in recent years. This interest was partially motivated by promising applications of nanopatterned substrates in the production of functional surfaces. Especially nanoscale ripple patterns on Si surfaces have attracted attention both from a fundamental and an application related point of view. This paper summarizes the theoretical basics of ion-induced pattern formation and compares the predictions of various continuum models to experimental observations with special emphasis on the morphology development of Si surfaces during sub-keV ion sputtering.


Introduction
Back in the 1960s, Navez et al. studied the morphology of glass surfaces bombarded with a 4 keV ion beam of air [1]. During the sputtering, they found the surface to develop periodic structures with lateral dimensions ranging from 30 to 120 nm depending on the angle of incidence. The orientation of the structures was determined by the direction of the ion beam. For grazing incidence, ripple patterns oriented parallel to the projection of the ion beam were observed whereas the ripples were rotated by 90 • at near-normal incidence. At normal incidence, however, the surface developed dot-like features. In the following years, sputter-induced ripple structures were found on all kinds of amorphous as well as crystalline materials like insulators [2], semiconductors [2,3], and metals [4].
During the 1990s, several in-situ and ex-situ studies investigated the ion-induced formation of nanoripples by means of new techniques for the exact characterization of the eroded surfaces like light scattering [5] and x-ray methods [6], as well as scanning tunneling [7] and atomic force microscopy [8,9]. In 1999, Facsko et al. observed the formation of hexagonally ordered nanodots on GaSb surfaces during normal incidence ion sputtering [10]. Such regular dot patterns have been found on various semiconductor surfaces sputtered at normal incidence [11] as well as off-normal incidence with [12] and without sample rotation [13].
Nowadays, ion-induced nanopatterns become interesting for various technological applications. Recent experiments demonstrate the principal applicability of nanoripples in the fabrication of microelectronic devices [14] and optically active nanostructure arrays [15,16]. Another approach uses nanodot formation under normal incidence sputtering of layer stacks to create isolated magnetic islands for magnetic storage media [17,18]. In addition, rippled substrates are becoming popular as templates for thin film deposition. It was shown that the morphology of the nanorippled substrates modifies the magnetic properties of ultrathin single-crystalline [19] and poly-crystalline [20][21][22] metal films. In a similar manner, arrays of close-packed nanomagnets could recently be obtained by shadow deposition on hexagonally ordered dot patterns [23]. Moreover, the self-organized alignment of physical-vapor deposited metal nanoparticles on nanorippled substrates was recently observed, leading to large arrays of nanoparticle chains exhibiting polarization-dependent plasmon absorption [24,25]. With the same technique, also arrays of metallic nanowires could be produced [26][27][28]. Most of these applications crucially depend on certain properties of the template patterns such as a high degree of order in the case of storage media [17,23] or a well defined ripple wavelength that fits to the growth conditions of the nanoparticles [24]. A precise control of the pattern properties in turn requires detailed knowledge of the pattern formation process and the contributing mechanisms. Up to now, however, this knowledge is still incomplete.
Although several possible origins of the ripple patterns like ion-induced local stresses or initial surface defects have been suggested in the years following their discovery [3], no conclusive explanation could be found until 1988. In this year, Bradley and Harper developed a continuum model [29] to describe the formation of the ripple patterns based on the so-called micro-roughening instability [30]. It was already shown by Sigmund [30] that the local erosion rate of a surface under ion bombardment is higher in depressions than on elevations. This curvature dependence of the sputter yield induces an instability of the surface against periodic disturbances which leads to an amplification of all initial modulations. In the presence of a competing smoothing process like surface self-diffusion, however, a wavelength selection is observed with the most unstable mode growing fastest [29].
The resulting linear continuum equation, the so-called Bradley-Harper (BH) equation, is able to reproduce some of the main experimentally observed features of the formation and early evolution of the patterns like their orientation with respect to the ion beam and the exponential growth of the ripple amplitude. For long sputtering times, however, certain experimental observations such as the saturation of the ripple amplitude cannot be explained within the framework of the linear model. This disagreement was attributed to a growing influence of nonlinear terms that dominate the morphology at later times. Hence, in 1995, Cuerno and Barabási derived a nonlinear continuum equation of the Kuramoto-Sivashinsky (KS) [31,32] type to describe the ion-induced formation of periodic surface structures [33]. In the early time regime, this equation behaves like the linear BH equation. At a certain transition time, however, the nonlinear terms start to control the evolution of the surface [34]. When entering this nonlinear regime, the amplitude of the ripples saturates as found experimentally. However, a transition to kinetic roughening with a loss of lateral order is observed in this regime [34,35]. Whereas such a transition has been observed in a few experiments [36], other studies report a stabilization of the regular patterns at high fluences [37][38][39]. Another feature of the experimental pattern evolution that could not be reproduced by the KS equation is the occasionally observed coarsening of the pattern wavelength [9,11,12,[39][40][41][42][43][44]. In order to overcome these discrepancies, several other nonlinear models based on the KS equation have been proposed [45][46][47][48][49]. These models all show a similar behavior in their linear regime and make different predictions only for the surface evolution in the nonlinear regime corresponding to rather long sputter times [50]. Therefore, a distinct demand for high fluence experiments has evolved which investigate the evolution of the surface morphology in the nonlinear regime in order to identify the continuum model that describes the given experimental system.
In the following section, the theoretical basics of ion-induced pattern formation are summarized and the various continuum equations available at present are discussed. Section 3 shows experimental results on the pattern formation and evolution on Si surfaces and tries to identify a certain continuum equation to describe the surface evolution. In addition, dependencies on experimental parameters are discussed with respect to possible applications. Section 4 provides a summary.

Continuum Theory of Ripple Formation During Low Energy Ion Sputtering
If a solid surface is bombarded with energetic ions, surface material will be removed [51,52]. The theoretical description of this mechanisms called sputtering has already been formulated in the 1960s by Sigmund [53]. The ions penetrating into the target surface are slowed down and lose their kinetic energy and momentum in elastic and inelastic collisions with target nuclei and electrons, respectively. For kinetic energies of the order of some keV and below, however, the momentum and kinetic energy of the ions are transferred to the target atoms in nuclear collisions mainly and inelastic collisions play only a minor role [54]. A target atom taking part in one of these collisions receives some of the ion's kinetic energy and momentum and can, therefore, be set in motion. If such an atom obtains sufficient energy, it can induce further collisions with other target atoms, thus increasing the number of moving atoms. This situation is then called collision cascade [54]. For typical ion fluxes, the collision cascades do not overlap in space and time and can therefore be treated independently. Within one collision cascade, it may happen that a target atom receives momentum directed towards the surface. If the kinetic energy of such an atom is high enough to overcome the surface binding energy, it will leave the surface and be sputtered away. Under continuous irradiation, the surface will be eroded as a whole. Additional effects that also might cause the removal of target material such as the deposition of potential energy during the impact of slow multiply-charged ions [55] will not be treated in this review.
When bombarding a crystalline non-metallic surface, e.g., a semiconductor, one can observe an additional effect. The number of generated defects in the crystal increases with the number of ion impacts. Therefore, for a large number of ion impacts, the crystal structure of the surface becomes unstable and the whole surface gets amorphized [54]. For single crystalline Si surfaces bombarded at energies of a few hundred eV at room temperature, this amorphisation is observed already after the impact of about 10 15 ions per cm 2 [56]. For higher fluences, the surface can be treated as fully amorphous.

Sigmund's theory of sputtering
A keV ion penetrating a solid surface loses its kinetic energy mainly in nuclear collisions with target atoms. The energy loss per unit path length, or stopping power, is then given by with the atomic density N of the solid and the nuclear stopping cross section S n (E). E is the initial kinetic energy of the penetrating ion. The nuclear stopping cross section S n (E) depends on the interaction potential used to model the collision between ion and target atom. With the power approximation of the Thomas-Fermi potential as a common choice, S n (E) reads [53] Here, m accounts for the Coulomb screening of the nuclei due to the electrons in the solid and ranges from 0 to 1. In the lower-keV and upper eV region, m = 1/3 is commonly assumed, whereas m should be close to zero in the eV region [53]. C m and ω are constants that incorporate the atomic parameters of the projectile and target species: M p,t is the atomic mass and Z p,t the atomic number of the projectile and the target atom, respectively. λ m is a dimensionless function of m with values ranging from λ 1 = 0.5 to λ 0 ∼ 24 and a T F is the Thomas-Fermi screening length.
The average number of sputtered atoms per incident ion is given by the sputtering yield Y . For linear collision cascades, i.e., for a sufficiently small number and isotropic distribution of binary collisions within one cascade [54], the sputtering yield Y is proportional to the energy F D (z) deposited per unit depth in the surface at z = h by a certain ion at the lateral position (x, y), with the ion energy E and the angle of incidence θ. Λ is given by Here, E sb is the surface binding energy and Γ m a function of m given by Because the majority of the sputtered particles originates from secondary collisions with low energy (< 50 eV) recoils, Sigmund suggested m = 0 for Equation (5) [53], resulting in Γ 0 = 6/π 2 . Therefore, Equation (5) becomes with C 0 = 0.0181 nm 2 [53]. For a plane and homogeneous surface, the deposited energy does not depend on the lateral position of the ion impact and is given by with α being a dimensionless function of the angle of incidence θ and the mass ratio M t /M p [53]. Then, the sputtering yield becomes According to Equation (9), the sputter yield depends on the surface binding energy and due to Equations (2) and (3) also on the atomic species. Therefore, for a multicomponent material, different sputter yields for individual atomic species i might be observed. In a first approximation, the total sputter yield can be treated as the sum of the different components according to their surface concentration. For this purpose, so-called "component" sputtering yields Y c i are defined such that the partial sputtering yields Y i follow the relation with the surface atomic fractions q s i . Then, the total sputtering yield is given by Different component sputtering yields then lead to i.e., one or more components are sputtered preferentially. Due to this preferential sputtering, the surface concentrations are altered at increasing fluence even in a homogeneous material. For a two-component material with the components A and B, preferential sputtering of A leads to a decrease of the surface concentration and thus also to a decrease of the partial sputtering yield of A. Prolonged sputtering will then lead to a stationary state described by which is characterized by the stationary partial sputtering yields Y ∞ i and the bulk atomic fractions q i . In the stationary state, the altered composition profiles remain constant but are moved into the bulk due to sputter erosion, so that atoms sputtered at the surface must be balanced by atoms fed from the bulk into the altered surface layer. From Equation (10) and (13), for the stationary surface composition is obtained with the stationary surface atomic fractions (q s i ) ∞ and the initial partial sputtering yields Y 0 i .

The Bradley-Harper model
If a surface is bombarded with a homogeneous flux of ions j, then the over-all energy deposited in a given point A of the surface is the sum of the energy deposited in this point due to all surrounding ion impacts. Therefore, with Equation (4), the local erosion rate in point A is given by the integral over all contributing events [30] v where ϕ(r) is the flux of incoming ions j corrected for the local angle of incidence and E D (r) is the energy deposited per unit volume at r = (x, y, z).
The spatial distribution of the deposited energy E D (r) can be approximated by a Gaussian, Here, µ and σ represent the lateral and longitudinal width of the distribution, respectively, and a is the mean penetration depth of the ion. A contour plot of the energy distribution is shown in Figure 1. For a rough surface sputtered with an uniform flux of ions, the energy deposited in the surface is not constant but rather depends on the lateral position r. To some extent, this is caused by the angular dependence of the ion flux at the surface. In addition, however, the energy deposition into the surface depends on the local shape of the surface. This lateral variation of the energy deposition causes a lateral variation of the local erosion rate and, therefore, a change of the surface morphology with sputtering time [30]. A closer inspection of the underlying mechanisms reveals that the local erosion rate is higher in troughs than on crests. This is demonstrated in Figure 2 where ions penetrate into a surface region with positive ( Figure 2 left) and negative ( Figure 2 right) curvature, respectively. The Gaussian distribution of the deposited energy is centered at the mean penetration depth a of the ions and indicated by the (broken) lines of constant energy. From Figure 2 it is obvious that the distance from the surface point A where the sputtering occurs to the contributing impact at B is shorter than the distance A * − B * . Therefore, the over-all deposited energy and also the erosion rate is larger in points with positive curvature (A) than in those with negative curvature (A * ). Obviously, the surface becomes unstable and the initial surface roughness gets amplified. This mechanism is called surface micro-roughening [30].  (17) with a = 3 nm, σ = 0.9 nm, µ = 0.5 nm, and E = 500 eV. The surface at z = 0 is indicated by the broken line. Figure 2. Schematic drawing of the energy deposition in rough surfaces, see text.
In order to explain the formation of periodic ripple patterns during sputtering, Bradley and Harper have calculated the integral (15) under the assumption of large radii of curvature R x and R y [29]. Then, the time evolution of the continuous surface height function h(x, y, t) is given by with φ being the angle between the direction of the ion beam and the local surface normal [33]. The projected direction of the ion beam is parallel to the x axis. Equation (18) can then be expanded in terms of derivatives of the surface height [33]. To first order in the surface curvature, Bradley and Harper obtained ∂h Here, v 0 is the erosion velocity of the planar surface, γ causes a lateral movement of the structures, and the micro-roughening instability is incorporated by the coefficients ν x,y . These coefficients are given by the following relations [45]: When sputtering a surface at finite temperature, atoms will diffuse on the surface leading to a relaxation of the surface. This effect, the so-called Herring-Mullins surface diffusion [57,58], can be introduced by adding a term proportional to the fourth derivative of the surface height to Equation (19), resulting in [29] ∂h In the Bradley-Harper (BH) equation (25), K is the relaxation rate due to thermally activated surface self-diffusion [29], with the surface self-diffusivity D s , the surface free energy per unit area ϱ, the areal density of diffusing atoms n d , the Boltzmann constant k B and the temperature T . The behavior of Equation (25) shall be analyzed by calculating its Fourier transform. Beh(k, t) the Fourier transform of the surface height function h(r, t) with the wave vector k = k x e x + k y e y and r = (x, y). Then, Equation (25) can be written as Integration of Equation (27) yieldsh with the growth rate . Therefore, spatial frequencies k with positive R k grow exponentially in amplitude, whereas those with negative R k decay exponentially with time. Because of the positive value of K, surface roughening occurs only for negative ν x,y . The maximum value of R k is reached for Therefore, the Fourier component of the initial roughness spectrum with the wave number k c will grow fastest, resulting in a wavelike surface pattern with a periodicity For ν x < ν y and ν x > ν y , the wave vector of the observed pattern is k c = k c e x and k c = k c e y , respectively. The angular dependence of ν x,y for a certain set of microscopic parameters is shown in Figure 3. At an angle of θ ∼ 73 • , one observes a change from ν x < ν y to ν x > ν y what corresponds to a rotation of the observed ripple pattern from normal to parallel with respect to the projected direction of the ion beam. This is demonstrated in Figure 4 which depicts numerical integrations [59] of Equation (25) at θ = 65 • (upper row) and θ = 75 • (lower row) at different times t. This type of pattern rotation with increasing incident angle has been observed in several experiments [1,36,40,44,[60][61][62]. Some other predictions of the BH equation, however, are at variance with certain experimental observations: • The amplitude of the ripples should grow exponentially without saturation. In experiments, however, saturation of the ripple amplitude at a constant value is observed after an initial exponential increase [63,64].
• Furthermore, from the same equations λ follows to be a function of the ion energy E and the penetration depth a, which again is a function of E. Therefore, one expects the ripple wavelength to decrease with the ion energy as λ ∝ E p with the negative exponent p [66]. However, this behavior is in general only observed at relatively high temperatures [67]. At low and moderate temperatures, several studies report the ripple wavelength to increase with energy [38,39,44,68,69].
• Equations (30) and (26) indicate a dependence of λ on the sample temperature. However, in the case of GaAs and InP, such a dependence of the wavelength was only observed at elevated temperatures whereas λ was found to be constant at room temperature and below [70]. Another study on SiO 2 surfaces found λ to be relatively constant with temperature even up to about 200 • C [71].
Several attempts have been made in order to overcome these deficiencies of the BH equation and shall be discussed in the following.

Kuramoto-Sivashinsky equation
In the series expansion of Equation (18), Bradley and Harper considered only linear terms. Cuerno and Barabási, however, took the expansion to lowest nonlinear order resulting in [33] ∂h ∂t The additional nonlinear terms in this equation are non-conserved Kardar-Parisi-Zhang (KPZ) nonlinearities [73,74] that incorporate the dependence of the local erosion velocity on the absolute value of the surface slopes. Their coefficients are given by [45] In order to account for the stochastic arrival of the ions, the Gaussian white noise term η, defined as was added. Here, D η is the strength of the noise and d the dimension of the surface. Equation (31) is an anisotropic stochastic generalization of the so-called Kuramoto-Sivashinsky (KS) equation which was originally proposed to describe chemical waves [31] and the propagation of flame fronts [32]. For short sputtering times, this equation behaves like the linear BH equation with an exponential increase of the ripple amplitude and constant ripple wavelength. Then, at a certain transition time the surface enters a nonlinear regime and a saturation of the ripple amplitude as in the experiments is observed [34]. However, numerical analyses of the noisy KS equation in 1 + 1 and 2 + 1 dimensions show that the saturation of the ripple amplitude is accompanied by a transition to kinetic roughening [34,35]. In this regime, the surface does not exhibit any lateral order. Although such a transition has been observed in few experiments [36], it is at variance with several other experimental reports of a pattern conservation at high fluences [37][38][39].

Damped Kuramoto-Sivashinsky equation
Inspired by the observation of stationary patterns in numerical simulations of the isotropic damped KS (dKS) equation by Paniconi and Elder [75], Facsko et al. adopted this equation for normal incidence ion sputtering [46]. The isotropic dKS equation is frequently used to describe different processes like compact electrodeposition growth [76] or directional solidification [75]. For oblique ion sputtering, however, the anisotropic dKS equation must be applied: This equation differs from the undamped KS equation (31) just by the additional damping term −αh with α being a damping coefficient that enters the effective growth rate of the ripple amplitude R * kc = R kc − α. This damping term induces smoothing of all spatial frequencies and, therefore, prevents kinetic roughening.
In the case of sputter erosion, the damping term in Equation (36) violates the translational invariance of the surface in the erosion direction. However, translational invariance can be restored by replacing the term −αh by −α(h −h) withh being the mean height of the surface and thus transforming Equation (36) into a nonlocal dKS equation [46] which again, as has been demonstrated [77], can be exactly mapped to a local dKS equation. The physical meaning of α, however, is still not clear in the case of sputter erosion.

General continuum equation
Although Equation (31) includes KPZ-like nonlinearities, other higher order terms are neglected [33]. The most general nonlinear equation that results from the expansion of Equation (18) is given by [45] ∂h ∂t The coefficients of the additional linear and nonlinear terms then read Actually, the ξ and Ω terms in Equation (37) have already been derived in reference [33] but were neglected since their influence on the asymptotic scaling of the surface was assumed to be of minor importance. The terms with the coefficients D ij enter Equation (37) in the form of diffusion-like terms proportional to the fourth derivative of the height function and thus lead to an additional anisotropic smoothing of the surface. Therefore, this relaxation mechanism is usually called effective or ion-induced surface diffusion (ISD) [81]. However, it is important to note that ISD results from preferential erosion during the sputtering which appears as a reorganization of the surface and does not involve any mass transport along the surface. Thus, ISD is strictly speaking no diffusion mechanism. This is also displayed by the fact that the coefficient D xx might even become negative at large incident angles, leading to an additional instability of the surface [81].
Since ISD does not depend on the temperature (cf. equations (42) -(44)), this smoothing mechanism is able to explain the temperature independence of the wavelength at low temperatures where thermal diffusion can be neglected. In this case, the ripple wavelength is given by From Equations (22), (23), (42), (43) and (45) it follows that the wavelength at low temperatures does no longer depend on the ion flux. Moreover, with a, µ, and σ being proportional to E 2m [54], we find λ ISD ∝ E 2m and, therefore, an increase of λ ISD with the ion energy. At high temperatures, however, thermal diffusion becomes the dominating smoothing mechanism and the wavelength follows from Equation (30). Hence, with the incorporation of ISD into Equation (37), one is able to explain the experimentally observed flux and temperature independence of the wavelength, as well as its increase with ion energy. However, the fluence dependence of the ripple wavelength as observed in some experiments [9,[39][40][41][42][43][44]72] still cannot be explained by the general continuum equation.
In the special case of normal ion incidence, the general continuum equation (37) is reduced to the isotropic stochastic KS equation with γ = ξ x = ξ y = Ω x = Ω y = 0, ν x = ν y , ζ x = ζ y , and D xx = D yy = D xy /2. For off-normal incidence, however, Equation (37) has a highly nonlinear character with a rich parameter space which might lead to rather complex morphologies and dynamic behaviors. Although some general features of Equation (37) have been studied [45], its detailed behavior, and especially the role of the additional nonlinearities with the coefficients ξ x,y , is still to be investigated.

Coupled two-field model
In order to overcome the inability of the KS-type Equations (31), (36), and (37) to predict ripple coarsening, Muñoz-García and co-workers recently developed a new nonlinear model following a hydrodynamic approach [49]. In this approach, Muñoz-García et al. considered two coupled fields where h and R represent the surface height function and the thickness of the mobile surface adatom layer, respectively. Here,φ = (1 − ϕ) is the fraction of eroded adatoms that become mobile, Γ ex is the curvature dependent erosion rate and Γ ad is the rate of addition to the immobile bulk. Γ ad is given by with the mean nucleation rate for a flat surface γ 0 , the variation in the nucleation rate with the surface curvatures γ 2x,y , and the thickness of the layer of mobile atoms generated thermally without bombardment R eq . Γ ex follows from microscopic derivations [50,82], The coefficients α i of Equation (49) . In the framework of Sigmund's theory of sputtering, these coefficients can be related to those of the general Equation (37) Equations (46)-(49) can be approximated by performing a multiple scale expansion with a subsequent adiabatic elimination of R. This results in an equation similar to the general continuum equation (37) but with additional conserved KPZ nonlinearities [49]: The coefficients of the coupled two-field (C2F) model differ from those of the general equation and are given by [82] The main novelty of the C2F model is the incorporation of redeposition of eroded material to the surface with the parameter ϕ controlling the amount of redeposited atoms. A key feature of this model is the presence of ripple coarsening which is probably induced by the conserved KPZ nonlinearity [49,50,83]. Depending on the ratio between the coefficients of the conserved and the nonconserved KPZ terms, i.e., ζ (1) i and ζ (2) ij , very different time dependencies of the ripple wavelength have been observed, ranging from marginal logarithmic to strong power-law coarsening. Moreover, in agreement with some experiments [9,11,12,39,41,44,84], the observed coarsening is interrupted at a certain time and the wavelength saturates at a constant value [49,50].

Morphology of Ion-sputtered Si Surfaces
Because of its great technological relevance, e.g., in micro-and nanoelectronics, silicon has attracted considerable attention during the last decades as an interesting material for nanopatterning by ion erosion [3,5,9,11,14,38,42,43,72,[85][86][87][88][89][90][91]. Thus, pattern formation on Si surfaces under various experimental conditions is well studied. However, the morphology of ion-sputtered Si surfaces exhibits some rather peculiar features and thus represents an interesting challenge for comparison with continuum theories. In this section, the morphology development of the Si surface during sub-keV ion sputtering will be summarized and discussed in the context of the different continuum models and in view of potential applications in thin film growth.   (Figure 5(a,b)), the Si surface remains flat. At a slightly larger incident angle of θ = 55 • (Figure 5(c)), however, the formation of shallow and rather disordered ripples that are oriented normal to the direction of the ion beam is observed. The wavelength of these ripples is about 50 nm. A further increase of the incident angle to θ = 67 • leads to a well ordered pattern of long homogeneous ripples with a periodicity of about 35 nm.
The observation that the Si surface remains flat at small incident angles is at variance with the BH model and most of its nonlinear extensions which predict an instability of the surface during ion sputtering independent of the experimental parameters. Carter and Vishnyakov explained a similar observation on Si surfaces bombarded with Xe ions of 10 to 40 keV energy as caused by an additional ion-induced mass transport along the surface that acts mainly at normal and near-normal incidence but is of minor importance at larger incident angles [9]. This so-called ballistic diffusion can also be introduced into the BH equation where it results in an additional smoothing term proportional to ∇ 2 h [9,92]. A similar mechanism has also been proposed for lower ion energies [92]. On the other hand, other experimental studies report dot and ripple pattern formation on Si surfaces also under normal and near-normal ion incidence [11,38,93,94]. However, recent experiments indicate that pattern formation under these low incidence conditions requires the presence of metal contaminations on the surface that may originate from the ion source [95] or the sample holder [96,97]. It has also been demonstrated that the resulting morphology of the Si surface can be tuned by varying the amount of metal contaminations during the sputtering [95,97]. A possible explanation for this so-called seeding effect invokes local variations of the sputter yield along the surface due to the segregation of deposited metal atoms that have a different component yield than Si [96]. A similar mechanism could also be responsible for the formation of dot patterns on compound semiconductors since there, preferential sputtering induces a form of "internal seeding" due to the enrichment and segregation of one atomic species on the surface. It has been shown theoretically that preferential sputtering can lead to a compositional modulation of the rippled surfaces of compound materials with the ripple crests having a different chemical composition than the valleys [98]. Since ion bombardment leads to an increase of the number of free bonds on the Si surface, also silicide formation could occur which would again alter the surface chemistry and thus also lead to a variation of the local sputter yield [96,99]. However, the presence of silicides on the sputtered Si surface could not be verified yet [95,96]. Also an increase of surface stress due to the seed atoms has been suggested as a possible origin of the dot patterns, a hypothesis that is supported by the experimental observation of tensile stress development in the presence of seeding [96,100].
With increasing angle of incidence, the BH model and the resulting linear and nonlinear continuum equations predict a rotation of the ripple pattern from normal to parallel with respect to the ion beam. Although this ripple rotation has been confirmed on various materials like metals [40,61], SiO 2 [1,44,60,62], and graphite [36], the formation of ripple patterns oriented parallel to the direction of the ion beam at grazing incidence seems to be suppressed on Si surfaces at room temperature, so that only shallow anisotropic structures have been observed [101,102] that do not resemble the well ordered patterns obtained at elevated sample temperature [5,72]. However, recent experiments by Mollick and Ghose [103] showed that the formation of a clearly developed rotated ripple pattern under 80 • incidence can be induced also at room temperature by a chemical pre-roughening of the Si surface which is known to influence the dynamics of the pattern development [45,91].

Evolution of the surface morphology
The various continuum models discussed in section 2.3 make different predictions for the temporal evolution of the surface morphology especially in the limit of long times where nonlinearities dominate. Therefore, the fluence dependence of certain parameters that characterize the surface morphology, e.g., the ripple amplitude and wavelength, is of particular importance for identifying a potential continuum description of the given experimental system. In addition, as will be shown below, the ion fluence is also a crucial parameter for the optimization of the pattern quality which therefore directly affects possible applications of the nanopatterned surfaces. Thus, in this section, the morphology evolution of Si surfaces will be discussed in detail for the example of sub-keV sputtering under 67 • incidence. At this incident angle, the formed ripple patterns exhibit the highest quality, a fact that might be correlated with the maximum of the sputter yield in this angular region.  13 nm (e), and 28 nm (f). The size of the images is 1 × 1 µm 2 (a-c) and 5 × 5 µm 2 (d-f), respectively; the ion beam was entering from the left. Insets: corresponding FFT ranging from -75 to +75 µm −1 (a-c) and from -4 to +4 µm −1 (d-f) [39]. Larger area AFM scans ( Figure 6(d-f)) reveal that the corrugations overlaying the normal pattern become anisotropic with increasing fluence and finally form a quasi-periodic pattern at high fluences, which is oriented parallel to the beam direction ( Figure 6(f)). This pattern is referred to as parallel pattern. Although the parallel pattern exhibits a much lower degree of order, side peaks can be identified (indicated by the white arrows) in the FFT, as shown in the inset of Figure 6(f). The side peaks indicate Figure 7. Evolution of (a) normal wavelength λ n , (b) parallel periodicity λ p , and (c) ratio of parallel to normal periodicity λ p /λ n over fluence for 300 eV and 500 eV. The solid lines in (a) represent power law fits, yielding coarsening exponents of n = 0.085 ± 0.006 and n = 0.084 ± 0.007 for 500 eV and 300 eV, respectively. The dotted lines represent logarithmic fits [39]. the quasi-periodicity of the parallel pattern and their position yields a much larger spatial periodicity of λ p ∼ 900 nm.
In Figure 7(a) the fluence dependence of the normal wavelength λ n , determined from the FFT of each AFM image, is depicted. Interrupted wavelength coarsening following a power law or logarithmic dependence is observed as soon as the ripple pattern is formed. Since wavelength coarsening is a nonlinear phenomenon, this indicates that nonlinearities start to dominate the surface evolution so early that no purely linear regime can be observed in the current experiments. In addition, λ n is found to increase with ion energy, indicating that ion-induced diffusion is the dominating smoothing process (cf. Section 2.3). This is also in agreement with the observed independence of λ n on the ion flux. The evolution of λ p is shown in Figure 7(b). Again, coarsening is observed. Figure 7(c) depicts the ratio of the wavelengths λ p /λ n . This ratio is quite constant in the investigated fluence range, indicating that both ripple modes exhibit similar coarsening behavior.  The evolution of the root-mean-square (rms) surface roughness w which describes the fluctuations of surface heights around the mean height and was calculated from the AFM images is shown in Figure 8. For both ion energies, w increases following a power law until it saturates at high fluences. One should note, however, that the rms roughness is not determined by the amplitude of the normal ripple pattern but rather by the larger corrugations and the parallel pattern, respectively. This is shown in Figure 8(b) that depicts the evolution of the ripple amplitude A, defined as the half of the average peak-to-peak height of the ripples, for the case of 500 eV sputtering. In the low fluence regime, the amplitude A is increasing from initially 0.4 nm to a maximum value of about 0.8 nm at Φ ≃ 5 × 10 17 cm −2 . For higher fluences, the amplitude decreases again and finally saturates at a value of A sat ≃ 0.6 nm. A similar overshooting before saturation has already been observed in previous experiments under normal ion incidence [37] and simulations of the anisotropic KS equation [34]. However, in contrast to the surface evolution in the KS equation, the experimentally observed saturation of the ripple amplitude is not accompanied by a loss of lateral order as is evident from Figure 6 which clearly shows a conservation of the pattern even at highest fluences. In combination with the observed interrupted wavelength coarsening, this suggests the C2F model as a potential description of the ripple formation and evolution on Si surfaces under these experimental conditions.

Dynamic scaling behavior
In the C2F model, with the interruption of the coarsening the surface enters a long-time regime that exhibits kinetic roughening at large lateral scales and a preservation of the ripple pattern at small scales [83]. Such a behavior is also seen in the experimental results presented in Figure 6. A kinetically rough surface is invariant under appropriate rescaling of its lateral and vertical dimensions and the time t [74]. This results in a certain behavior of its surface roughness w( , t) is the surface height function, l is the size of the observation window over which w has been calculated, and the angular brackets denote spatial averaging. In the case of Family-Vicsek (FV) dynamic scaling [105], the roughness should scale as w(l, t) ∼ t β until the correlation length ξ(t) ∼ t 1/z has reached the window size l. Then, the roughness will saturate with the saturation value depending on the window size, w(l) ∼ l α . The roughness exponent α, the growth exponent β, and the dynamic exponent z = α/β characterize the surface in space and time and can be used to attribute the system to a certain universality class and, therefore, to a certain continuum equation [74]. With this intention, the dynamic scaling behavior of the ion-sputtered Si surface has been analyzed by evaluating its one-dimensional structure factor. According to the dynamic scaling hypothesis [74], the one-dimensional structure factor should obey the relation with the scaling function s(u) ∼ u 2α+1 and s(u) ∼ const. for u ≪ 1 and u ≫ 1, respectively. In the case of anisotropic surfaces, this behavior is modified and the surface is characterized in the normal and parallel direction in real and momentum space by four different roughness exponents [106]. However, for kt 1/z ≫ 1, the dynamic scaling behavior of the one-dimensional structure factor can still be described by Equation (59) [104,106]. The structure factor S p (k p ) calculated in the direction parallel to the ion beam is given in Figure 9(a). For Φ ≥ 5 × 10 16 cm −2 , a peak appears at the spatial frequency k * p corresponding to the wavelength λ of the ripple pattern. For k p ≫ k * p , the S p curves all collapse. The slope m (in the log-log plot) of the curves in this regime is about −4, corresponding to a roughness exponent of 1.5. With increasing fluence, the ripples coarsen and the position of the peak is shifting to smaller k p values. Also the structure factor increases with fluence for k p ≪ k * p and a second scaling regime develops at high fluences. Here, the roughness exponent was determined to be α p = 0.41 ± 0.04. In Figure 9(b), the structure factor curves in the direction normal to the ion beam, S n (k n ), are depicted for different fluences. At large values of k n , the data is consistent with a slope m = −4. At small k n values, S n (k n ) again increases with fluence and until a power-law behavior with a roughness exponent α n = 0.76 ± 0.04 appears at high fluences.
The observed peak in the structure factor S p in the direction parallel to the ion beam with the -4 slope at large k p values (cf. Figure 9(a)) indicates the presence of a KS like instability in this direction [35]. The orientation of the ripples with respect to the incident ion beam is determined by the signs of the linear coefficients: the wave vector of the observed ripple structure is parallel to the direction with the smallest negative ν (cf. section 2.2). Therefore, for the here presented experiment ν p < ν n . In the direction normal to the ion beam, the experimental S n curves shown in Figure 9(b) do not exhibit a local maximum. The determined low-fluence behavior for the n direction S n (k n ) ∼ k −4 n corresponds to the scaling behavior of the one-dimensional linear molecular beam epitaxy (lMBE) equation with α lM BE = 3/2 [74]. This indicates that the very short-distance behavior of the sputtered Si surface is dominated by the diffusion term. This behavior holds even at the highest applied fluence of Φ = 1 × 10 20 cm −2 without any noticeable crossover. This indicates that |ν n | ≈ 0 [104].  In the limit of high fluences Φ ≥ 10 19 cm −2 , the morphology of the Si surface exhibits anisotropic algebraic scaling at large lateral scales with α n = 0.76 and α p = 0.41. The KS equation (31) is not able to reproduce such an anisotropic scaling behavior since the only term breaking the x → −x symmetry is the one with the coefficient γ which acts only at rather short length scales. On the other hand, the dispersive nonlinearities with the coefficients ξ x,y that appear both in the general continuum equation (37) and in the C2F model (50) have been found to induce anisotropic scaling under certain conditions [106]. Therefore, the appearance of anisotropic scaling supports above assumption of the C2F model being a suitable description of the Si surface during sub-keV ion sputtering [104].

Dynamics of topological pattern defects
In view of possible applications of the nanorippled Si surfaces, the appearance of kinetic roughening, i.e., of a disordered state, at high fluences is not favorable since most of these applications rely on a well-ordered homogeneous pattern. Therefore, the applied ion fluence is an important parameter in the fabrication of nanopatterned surfaces and vital for optimizing their quality. The quality of the ripple patterns can be quantified by calculating a normalized density of topological pattern defects from the AFM images [86,103,107]. In this context, topological pattern defects means either a bifurcation (B) of a ripple, i.e., a Y junction of ripples, or an interstitial (I), i.e., a discontinuous or broken ripple. Figure 10(a) shows an AFM image of the rippled Si surface in which these defect types are indicated.
The procedure of determining the normalized density of topological pattern defects is demonstrated in Figure 11. In order to determine the total number of defects of a given AFM image (Figure 11(a)), the image is Fourier-filtered to remove the long-wavelength surface morphology (Figure 11(b)). The filtered image is then converted into a binary image by applying Otsu's threshold [108] (Figure 11(c)). Finally, the ripples in the binary image are thinned to lines of one pixel width (Figure 11(d)). Then, every black pixel with more or less than two black neighboring pixels is counted as a defect. The normalized density of defects is then calculated as N D = N λ 2 /A S with the total number of defects N of the image, the ripple periodicity λ, and the scan area A S . A value of N D = 0 then corresponds to a perfect pattern without any defects and N D = 1 to a pattern in which each ripple contains one defect per length λ.
Following this approach, the normalized density of pattern defects N D has been calculated for different fluences in order to monitor the evolution of the pattern quality. The result is shown in Figure  10(b). The N D values are comparable for both energies, although in average N D appears to be slightly lower for 300 eV than for 500 eV. At the lowest fluence Φ = 5×10 16 cm −2 , the normalized defect density is around 0.3. With increasing fluence, N D decreases until it reaches a minimum value of N D ∼ 0.1 (500 eV) and 0.07 (300 eV) at a fluence of Φ ∼ 2 × 10 18 cm −2 . This decrease of N D is caused by the growth of the ripple length and the annihilation of pattern defects due to a complex interplay of different "annealing" processes [107]. At higher fluences, N D increases again until it saturates at Φ ∼ 10 19 cm −2 at a value of N D ∼ 0.28. This increase results from the appearance of kinetic roughening which induces a certain disorder in the pattern that leads to the formation of "defect clusters" [107]. Interestingly, the coarsening of the ripple wavelength does not seem to be related to the evolution of the pattern defects (cf. Figures 7(a) and 10(b)). This is in contrast to other experimental systems like Pt(111) surfaces under grazing incidence sputtering where rapid coarsening proceeds due to the annihilation of defects [109]. These results demonstrate the influence of the applied fluence not only on the ripple amplitude and wavelength but also on the pattern quality. Therefore, in order to fabricate patterns of a certain periodicity in the highest quality possible, the interplay between fluence, energy, pattern wavelength, and pattern quality needs to be known.

Summary
We have presented an overview of the continuum approach to ion-induced pattern formation on amorphous surfaces. The predictions of the various linear and nonlinear continuum models have been discussed and compared to experimental observations with a special focus on sub-keV ion sputtering of Si surfaces. Because of its potential applications, pattern formation on Si surfaces induced by low-energy sputtering has been investigated excessively during the last two decades. These studies revealed several peculiarities of the morphology of sputtered Si surfaces such as the stability of the flat surface at near-normal incidence, (interrupted) wavelength coarsening or the absence of a pattern rotation with increasing angle of incidence. In addition, contradictory observations have been reported, e.g., the occurrence of smoothing and roughening at small incident angles, respectively.
Recent experimental findings such as the importance of metal contaminations during the sputtering, delivered further insight into the basic mechanisms of ion-induced pattern formation on Si surfaces. In addition, novel and elaborated theoretical models provided new explanations for certain experimental observations, e.g., wavelength coarsening or the occurrence of anisotropic scaling. Therefore, a rather coherent picture of the morphology of ion-sputtered Si surfaces has developed during the last few years. However, at the same time new challenges, both experimental and theoretical ones, have appeared, among them the control of surface contaminations and the investigation of its detailed effects on the morphology development which might enable the fabrication of novel nanopatterned surfaces.
On the other hand, the application of nanorippled Si substrates in various fields of modern materials science, especially in nanoscale magnetism and plasmonics, is developing tremendously and demands for a precise control over the fabricated patterns. Besides the tuning of the wavelength and the amplitude, the quality, order, and regularity of the patterns is becoming more and more important since the order has a strong effect on the degree of the induced functional anisotropy. Providing nanopatterned substrates of high quality and with tailored properties has thus become a major experimental issue.