Hydrodynamic approach to surface pattern formation by ion beams
Abstract
On the proper timescale, amorphous solids can flow. Solid flow can be observed macroscopically in glaciers or lead pipes, but it can also be artificially enhanced by creating defects. Ion Beam Sputtering (IBS) is a technique in which ions with energies in the 0.1 to 10 keV range impact against a solid target inducing defect creation and dynamics, and eroding its surface leading to formation of ordered nanostructures. Despite its technological interest, a basic understanding of nanopattern formation processes occurring under IBS of amorphizable targets has not been clearly established, recent experiments on Si having largely questioned knowledge accumulated during the last two decades. A number of interfacial equations have been proposed in the past to describe these phenomena, typically by adding together different contributions coming from surface diffusion, ion sputtering or mass redistribution, etc. in a nonsystematic way. Here, we exploit the general idea of solids flowing due to ion impacts in order to establish a general framework into which different mechanisms (such as viscous flow, stress, diffusion, or sputtering) can be incorporated, under generic physical conservation laws. As opposed to formulating phenomenological interfacial equations, this approach allows to assess systematically the relevance and interplay of different physical mechanisms influencing surface pattern formation by IBS.
keywords:
Ion Beam Sputtering, Hydrodynamics, Solid flow, Pattern formation, Stability, Viscous flow, Stress, Erosion1 Introduction
Observations of nanoscale patterns on the surfaces of solid targets that undergo ionbeam sputtering (IBS) by 100 eV to 10 keV ions, date back at least to the early 1960’s cunningham:1960 ; navez:1962 . They correspond either to amorphous materials like glass, or to targets that are amorphized by this type of irradiation, like semiconductors gnaser:1998 . The fascinating resemblance with macroscopic structures, like ripples on water or sandy dunes,^{1}^{1}1In Fig. 1 we show the nice nakedeye similarity between water and Silicon ripples at very different orders of magnitude. emphasized in navez:1962 , goes beyond a mere visual analogy munozgarcia:2010 ; chan:2007_et_al , underscoring the basic interest of this type of structures. Moreover, the high efficiency of IBS to produce surface nanopatterns (ripples and dots) over large areas on top of a wide variety of targets including metals and insulators chan:2007_et_al ; munozgarcia:2009 , furnishes the technique with a large potential for applications.
Specifically, in IBS both the erosive process and the appearance of the patterns are due to the complex iontarget interactions: thus, ions create permanent defects in a surface layer whose thickness is of the order of the average penetration depth, trigger material ejection (sputtering) from its surface and induce material redistribution and flow volkert:1993 ; mayr:2003 . A conceptual framework explaining pattern formation in amorphous or amorphizable systems is available since the late 80’s, through the work by Bradley and Harper (BH) bradley:1988 that was based on Sigmund’s theory of linear collision cascades. Physically, a morphological instability cross:2009 originates in the higher erosion rate (yield) at surface minima as compared with surface maxima, in such a way that a pattern forms for any value of the angle of incidence and ion energy. This has become the theoretical paradigm supporting most experimental observations. Minor discrepancies with this view have been attributed to secondary effects, leading to refinements of BH’s theory (see reviews in chan:2007_et_al ; munozgarcia:2009 ).
However, many of the experiments done on elemental targets have been shown to be affected by contamination ozaydin:2005 , suggesting a correlation between pattern formation and a sizeable concentration of undesired species. In view of this, very recently clean experiments have been specifically tailored in order to avoid this feature madi:2008 ; madi:2009 ; macko:2010 ; macko:2011 , taking Si targets as representative cases for the large class of substrates that become amorphous under irradiation. Remarkably, the outcome of these works invalidates the BH paradigm, since various morphological transitions are observed between unstructured and patterned surfaces as a function of both incidence angle and ion energy , confirming earlier observations at higher energies carter:1996 and contradicting one of the main predictions of BH theory. Thus, ironically, while for the case of multiplecomponent systems differential sputtering and species segregation may have been recently identified shenoy_et_al ; leroy:2010 ; bradley:2010 ; zhou:2010 as the main physical mechanisms behind the morphological instability, the a priori simpler case of singlecomponent systems still remains to be understood.
In addition, none of the current phenomenological models that elaborate on the BH paradigm makeev:2002 ; facsko:2004 ; davidovitch:2007 are able to capture the complex morphological features elucidated by the recent experiments on Si. In these theories, the dynamics of the target surface is described through an effective evolution equation in which the various physical mechanisms are combined adhoc, and defect dynamics are neglected. The same limitations when confronted with experiments occur even for improved physical models in which the evolution of the surface height is explicitly coupled to that of the density of defects whose transport is confined to a thin surface layer, akin to pattern formation on e.g. aeolian sand dunes, see e.g. castro:2005 ; munoz:2006 ; munoz:2008 ; munoz:2009b , and cuerno:2011 for a recent overview.
In order to overcome the limitations of previous continuum approaches to surface dynamics by IBS, we propose a generic framework based on general conservation laws (such as those governing mass or momentum), that allows to study systematically the role of the different mechanisms involved in the process. This approach, reminiscent of classical hydrodynamics, will allow us to identify the physics underlying the patterning of the eroded target surface. In summary, the main contribution of this work is to present a general consistent physical framework to study this system, that inherits the benefits of traditional fluid mechanics problems for which it has been proved highly successful in the past.
2 Governing equations
In the energy range considered, an impinging ion interacts mostly with the nuclear structure of the target, generating vacancies and interstitials. Ions and defects can diffuse at different rates, annihilate by recombination or form clusters. As borne out from Molecular Dynamics (MD) simulations, the overall effect of these processes is threefold: amorphizing pelaz:2004 , stressing the target kalyanasaundaram:2008 , and displacing the mean atom position inside the material moseler:2005 . Moreover, as also seen in experiments gnaser:1998 , an amorphous layer forms rapidly moore:2004 ; kalyanasaundaram:2008 , its thickness and other mechanical properties like density, , becoming stationary after a fluence of order ions cm (corresponding to a few seconds of irradiation for a typical ion flux of to ions cm s). Such a small ion flux drives dynamics slowly, inducing (as in macroscopic solid flow for e.g. a glacier) a time scale separation between atomistic relaxation rates ( ps for collision cascades or adatom hopping), and the s times in which significant morphological changes occur. This fact legitimates a continuum description.
2.1 Conservation of mass
We are interested in the stages of the bombarding process in which the amorphous layer has already been formed. Thus, although its density value changes with respect to that for the pristine target, it does not evolve further after the initial transient just described. Mathematically, this fact can be expressed in the form
(1) 
where is the velocity field of the fluidized layer, whose components are given by
(2) 
In the language of fluid mechanics, the amorphous layer is said to be incompressible in the steady state.
2.2 Conservation of momentum
Conservation of momentum can be universally expressed as oron:1997
(3) 
where is the stress tensor. Indeed, the physics of the IBS problem enters the constitutive equation for this tensor.
Although a complete characterization of the stress tensor in the bombarded material is still lacking, there are strong evidences that the amorphous layer can be considered as a highly viscous fluid, hence can be written as
(4) 
where contains all the terms that do not come from hydrostatic pressure, , or viscous flow, and are the components of the velocity field (2).
Moreover, because the radiation induced viscosity, , is around orders of magnitude larger than the viscosity of water mayr:2003 , we can safely drop the left hand side of Eq. (3). Thus, combining the conservation of mass with this assumption for the amorphous layer, we can write Eq. (3) simply as
(5) 
where we define as a body force acting in the bulk of the fluid layer. Actually, this force contains the relevant physics besides viscous flow. It is effective in the sense that detailed knowledge of the stress induced beneath the surface is not yet known experimentally. However, it is reasonable to expect that this force can be split into an amplitude, , and an angular contribution, , which is a function of the local angle of incidence. Mathematically,
(6) 
where is the local slope of the surface and is the incidence angle of the ion beam. Here, contains the coarsegrained information about the effect of the residual stress created in the target, due to ioninduced mass redistribution, and has dimensions of a gradient of stress.
In three dimensions, and using standard spherical coordinates, we can write
(7)  
(8)  
(9) 
where is the azimuthal angle (the angle in the plane of the target with respect to the plane of incidence of the ions). Possible anisotropies can enter the dynamics from the difference in the values of the amplitudes and . We are assuming that the collision cascade that produces the damage and, ultimately, the strain/stress that drives the amorphous layer, is oriented in the direction of the incoming ion (see Fig. 2). The fact that is the same in the equations for and comes from the fact that there are, actually, just two important directions: parallel and perpendicular to the ion flux. When we project the parallel contribution in those directions ( and ), the prefactor remains the same.
This body force can be understood as an effective gravity term. Thus, in Fig. 2 we provide a cartoon to illustrate the analogy between ion induced solid flow and a genuine fluid flowing down an inclined plane. The two systems are mathematically equivalent once a rotation of angle is performed in the laboratory reference frame. The main difference between both systems comes from the corresponding gravity fields. While for the standard fluid problem gravity is a constant, for the erosion system the effective gravity changes from point to point on the surface, its value depending strongly on the local surface orientation with respect to the ion beam. Thus, unlike the fluid problem, in the erosive case one can obtain a morphological instability for appropriate angles of incidence.
2.3 Boundary conditions
Equations (1)(9) need to be supplemented with proper boundary conditions involving the erosive terms and the stress at the interfaces. In the most general case, we have to specify the evolution for two important boundaries of the system: the amorphousvacuum interface (the free surface), , and the amorphouscrystalline interface, .
At the free surface, the balance of stress takes the form oron:1997
(10)  
(11)  
(12) 
where is the unit outward normal vector, and are two mutually perpendicular vectors in the tangent plane to the surface, is the mean surface curvature, and is surface tension. Finally, are the projections along the corresponding directions of an external stress that is applied over the surface.
Here we will consider that the surface tension is constant. However, for multicomponent materials where one species segregates at the surface, one can assume that there is a surface tension gradient that may affect the dynamics of the surface (the socalled Marangoni effect oron:1997 ; Kree ).
Finally, a kinematic condition oron:1997 at each interface leads to evolution equations for both and as
(13) 
and
(14) 
where, as defined above, is the vertical component of the velocity field. The terms and stand for the rates of amorphization and erosion, respectively. Thus, Equations (13) and (14) express mathematically that the difference between the vertical velocity field and the actual motion of the interface is equal to the rate at which material is removed from both phases.
For instance, the erosive (BHtype) mechanisms enter directly through , so that their contribution to the full dispersion relation will be additive, see section 3. For simplicity, we assume that the rate of amorphization is equal to the rate of erosion, since in experiments both interfaces evolve at the same rates in the steady state. Moreover, for mathematical tractability, we assume that the crystallineamorphous interface is flat. Under these assumptions, a final boundary condition is required, by which we specify that, at , the fluid moves with the same velocity as the crystalline phase; this is often referred to as the noslip boundary condition oron:1997 .
3 Linear analysis
The standard procedure to study the linear stability of a surface is developed in three steps:

Obtaining the flat (zeroth order) interface solution.

Making a periodic infinitesimal perturbation of the solution found in the previous step, with wavelength ( being the wavenumber).

Calculating the amplification rate, , for every wavenumber (also known as linear dispersion relation).
For the sake of simplicity, we will detail these steps for the two dimensional case (namely, the interface is a line and not a proper surface). At the end of our derivation, we will quote the dispersion relation for the full three dimensional case.
3.1 Flat interface solution
We fix our reference frame on the free surface, so that the crystallineamorphous interface is located at , where we denote by the average thickness of the amorphous layer and by the fluctuations of the free surface. The flat interface condition is given simply by . Introducing it in the equations for the velocity field, we find that the flat solution is further characterized by
(15)  
(16)  
(17) 
3.2 Linear perturbation
To study the morphological stability of the flat interface, we look for solutions of the governing equations of the form
(18) 
for an arbitrary wavenumber . Similarly, we assume that the pressure and the velocity fields can be written as
(19)  
(20)  
(21) 
and solve the governing equations perturbatively to linear order in the small parameter . The precise form of the ensuing corrections to the flat interface solution, , and , is given by Eqs. (54)(56) in A, respectively.
3.3 Dispersion relation
Finally, the dispersion relation can be extracted from the kinematic condition. Specifically, from Eqs. (14) and (18) we get
(22)  
where comes directly from the erosive contribution ; here is where, for instance, the classical theory of Bradley and Harper can enter bradley:1988 .
Using Eqs. (15)(21) and denoting real part with a single prime and imaginary part with a double prime, we find:
(23) 
and
(24) 
3.3.1 Real part: Stability of the flat interface
The real part of controls the stability of the flat interface: if is positive/negative, a periodic perturbation to the flat solution with wavenumber will grow/decay exponentially as given by Eq. (18). Equation (23) generalizes important results from fluid dynamics, and has already been reported in the context of the IBS problem for more specific choices of the function arxiv ; cuerno:2011 . For instance, for , it reduces to Orchard’s classic result orchard:1962 on the flow of a viscous layer of arbitrary thickness on top of an immobile substrate. In particular, then leads to Mullins’ mullins:1959 celebrated rate of flattening when relaxation occurs through bulk viscous flow, , already invoked for some IBS experiments chason:1994 ; frost:2009 .
Generally, equation (23) is nonpolynomial, implying nonlocal effects in the dynamics of the amorphous layer. This nonlocality may be important when other nonlocal effects are present in the dynamics (as redeposition, shadowing, or secondary ion scattering). Moreover, in our case the sign of may depend on the incidence angle, , thus explaining morphological transitions from flat to patterned structures at different values of madi:2008 ; madi:2009 ; macko:2010 ; macko:2011 .
3.3.2 Imaginary part: inplane propagation of the pattern
In those cases in which an unstable wavelength is selected (namely, for some ), the imaginary part of provides information about the velocity of lateral inplane propagation of the pattern. Specifically, if the most unstable wavenumber is (namely, the one with the largest positive value of ), then the pattern propagates with velocity
(25) 
Again, the sign of may depend on the incidence angle; thus, depending on the choice of , the ripples can move in the same direction as or opposite to the ion beam.
3.3.3 “Shallow water” approximation ()
For large enough energies, and because the typical wavelengths of the patterns reported in the literature for energies around 1 keV are in the nm range, thus an order of magnitude larger than the amorphous layer thickness ( nm), we can assume that and Taylor expand . In fluid mechanics this is known as the shallowwater approximation oron:1997 . Thus, we get
(26) 
where it is understood that a similar expansion has been done in the erosive contribution . Note that, if we consider solid flow as the only pattern forming mechanism, we can simply drop this term. Since we are interested in probing the morphological consequences of viscous flow, we shall do this henceforth unless otherwise stated.
Equation (26) predicts that, if the function is positive, the interface is stable (flat or rough, depending on the intensity of the thermal and beam flux fluctuations); on the contrary, if it is negative, a pattern appears with an angle dependent wavelength.
At this point, note that our hydrodynamical formulation is a coarse grained description of the microscopic phenomena. Thus, it has been recently shown that the ion impacts produce a response in the target referred to as a prompt regime norris:2011 . In this reference, starting from atomistic molecular dynamics (MD) simulations, the overall effect of the ion impacts on the evolution of the surface height is summarized into the moments, , of a socalled crater function, that eventually contribute to a continuum evolution equation for the surface height. It is worth noting that we can map the approach norris:2011 with our present formalism through
(27) 
In this respect, we could incorporate MD information into our hydrodynamical description through our function . Let us note incidentally that, as suggested by the simulations in norris:2011 ; kim:2011 , purely erosive contributions to the evolution of the surface height are numerically negligible, justifying our neglect of above.
Patterns generated in the unstable parameter region in which Eq. (26) has a positive maximum for , have a characteristic wavelength given by
(28) 
For instance, if one assumes that the angular correction is (namely, that the body force is proportional to the local flux at the surface), then Eq. (26) predicts a transition from flat to patterned surfaces at for all ion energies, as has been observed in experiments of Ar on Si at energies between eV and keV madi:2009 . Similarly, for this choice of , the velocity of ripple is given by
(29) 
To be more quantitative, we need to determine the values of the main parameters of the theory. As mentioned above, the new parameter is related to the stress gradient created across the amorphous layer, induced by the collisional processes. There is a large experimental uncertainty in the measurement of stress, a property that depends moreover on the energy of the incident ions. Different authors report different values for ioninduced stress on Si, in the range MPa to GPa madi:2009 ; kalyanasaundaram:2008 . We are interested in the order of magnitude, hence we take the geometrical mean (which accounts for the average order of magnitude) of both extreme values, thus, approximately 569 MPa. This corresponds to kg nms for a layer of nm depth. In Fig. 3, this value corresponds to the blue solid line, while the limiting values of MPa and GPa are used to plot the red dashed lines, which we take as our confidence interval. More accurate measurements of stress would provide improved estimations. Additional parameters have not been measured experimentally yet. However, we can provide indirect estimates for their values when these are not explicitly available, as shown in Table 1. With these parameters at hand, Eq. (28) and experimental data are seen to agree quite closely in Fig. 3.
Parameter  Value  Reference 
Surface tension ()  J/m  vauth:2007 
Amorphous layer thickness ( eV)  nm  madi:2008 ; madi:2009 
Viscosity ()  Pa s  norris:2011 
Stress range  MPa– GPa  madi:2009 ; kalyanasaundaram:2008 
The dependence of the pattern properties on other parameters such as energy, flux, etc., that enter the amplitude , will be the subject of future work.
3.4 Dispersion relation for the three dimensional case
For completeness, we include here the dispersion relation for the three dimensional case. Now, the flat interface is a plane, and deviations to linear order take the form
(30) 
with and . In the “shallowwater” limit, the dispersion relation is now given by
(31) 
for a general angular dependence of the body force and general (anisotropic) forms for the strain induced by the ion through the amplitudes and .
4 Connection with other continuum frameworks
4.1 Twofield “hydrodynamic” theory
Despite the fact that we have not considered the dynamics of the amorphouscrystalline interface, here we want to note that, in the shallowwater limit, one can reduce the present full hydrodynamical formulation, making contact with previous theories of erosion, specifically with the twofield model approach pioneered in aste and further developed in castro:2005 ; munoz:2006 ; munoz:2008 ; munoz:2009b , see an overview in cuerno:2011 .
Within this approach (see Ref. munoz:2008 for details), the dynamics of the density of moving species and that of the free surface height, , are explicitly described, within the implicit assumption that the thickness of the surface layer where material transport takes place is negligible. The evolution of these two fields can be expressed through proper rate equations, borrowed from the field of pattern formation on aeolian sand dunes. Specifically,
(32)  
(33) 
where measures the fraction of atoms that, having been dislodged from their equilibrium positions, are actually sputtered. Here, and are, respectively, rates of atom excavation from and addition to the (crystalline) immobile bulk.
4.2 Lubrication theory
The hydrodynamical nature of the model considered in Sections 2, 3 suggests a direct comparison with the socalled lubrication theory of fluids. This theory is valid in the shallowwater approximation which is, in fact, the one relevant to IBS experiments.
Following Oron et al. oron:1997 , we define a small parameter . Additionally, we assume that not only the slopes are small but also the timescales related to the motion of the fluid are slow (with respect to the motion of the very same surface). Mathematically, the following change of variables is defined:
(38) 
In order to capture the essential ingredients that affect the motion of the fluid, we rescale the stress tensor and the corresponding values at the boundary. Thus,
(39) 
and
(40) 
Finally, the surface tension in the proper units is redefined as oron:1997 . Introducing this change of variables into the governing equations and boundary conditions, and keeping the dominant contribution to order , we find
(41)  
(42)  
(43)  
(44)  
(45) 
Integrating this system, one arrives to
(46) 
where is the contribution from erosion in the scaled variables.
In order to solve (42)(45), one needs some closure relation between the stress and the velocity. In our problem, we have split the stress into a body force and a (Newtonian) viscous tensor. For such a case, one arrives to a closed nonlinear equation of the form
with
(47)  
(48)  
(49) 
Hence, with and ,
(51)  
Within our framework, an alternative equation can also be derived that features the same real part of the linear dispersion relation, setting , , and , with as derived from Ref. spencer for a twodimensional coordinate system aligned with the ion beam,
(52) 
where for a twodimensional incompressible fluid. Thus,
(53)  
with for compressive stress (the experimental case).
Linearizing Eq. (4.2) and neglecting the effects of external stresses and the Marangoni term (), we reproduce the shallow water dispersion relation given by Eq. (26). The interest of this lubrication approach is that, provided an appropriate closure relation between the stress tensor and the velocity field is put into Eq. (42), we obtain a nonlinear equation for the surface height through Eq. (46), given Eq. (42). Thus, nonlinear effects can be assessed that go beyond the analysis performed in Section 3.
5 Conclusions
In this work we have introduced a physical framework for surface pattern formation by IBS that is inspired in classic Fluid Mechanics. The strength of the theory is that it is based on general conservation laws, and that all the physics of the erosive process can be included into appropriate constitutive equations for the stress (body force) and boundary conditions. This approach allows to properly account for different mechanisms, such as viscous flow, erosion, or ioninduced strain (damage).
We have performed a linear stability analysis of the governing equations and shown that the forces related to the stress produced by the ion impacts indeed influence the morphological stability of the interface. For a specific choice of the angular dependence of the body force, , we have seen that some of the complexities of the morphological diagram for Si that have been unveiled recently, can be accounted for by our description. In this connection, purely erosive contributions indeed seem to play a minor role in the surface dynamics as compared with viscous flow.
Finally, the proposed hydrodynamic framework can be seen to generalize previous twofield formulations of IBS, to the case of a flowing layer of arbitrary thickness. Such phenomenological models have proved successful to obtain an effective equation describing the nonlinear dynamics of the surface, reaching in some cases quantitative munozgarcia:2010 and semiquantitative kim:2011b agreement with experiments. On the other hand, our present framework does allow for a systematic study of nonlinear effects through a lubrication theory approach.
The specific details of the constitutive equations that are appropriate for different experimental conditions will be the aim of future work. Also the general nonlinear theory will be addressed through standard procedures oron:1997 that should be able to account for nonlinear phenomena such as pattern coarsening, saturation of the roughness, or defect creation and/or annihilation.
Acknowledgments
This work has been partially supported by MICINN (Spain) Grants No. FIS200912964C0501 and No. FIS200912964C0503.
Appendix A Solutions of the linearized governing equations
(55) 
and
(56) 
References
References
 (1) R. L. Cunningham et al., J. Appl. Phys. 31, 839 (1960).
 (2) M. Navez, C. Sella, and D. Chaperot, C. R. Acad. Sci. Paris 254, 240 (1962).
 (3) H. Gnaser, Low Energy Ion Irradiation of Solid Surfaces (Springer, New York, 1998).
 (4) J. MuñozGarcía, R. Gago, L Vázquez, J. A. SánchezGarcıa, and Rodolfo Cuerno, Phys. Rev. Lett. 104, 026101 (2010).
 (5) W. L. Chan and E. Chason, J. Appl. Phys. 101, 121301 (2007).
 (6) J. MuñozGarcía, L. Vázquez, R. Cuerno, José A. SánchezGarcía, M. Castro, and R. Gago, in Towards Functional Nanomaterials, edited by Z. M. Wang (Springer, New York, 2009).
 (7) R. Cuerno, M. Castro, J. MuñozGarcía, R. Gago, and L. Vázquez, Nucl. Instrum. Meth. Phys. Res. B 269, 894 (2011).
 (8) C. A. Volkert, J. Appl. Phys. 74, 7107 (1993).
 (9) S. G. Mayr, Y. Ashkenazy, K. Albe, and R. S. Averback, Phys. Rev. Lett. 90, 055505 (2003).
 (10) R. M. Bradley and J. M. E. Harper, J. Vac. Sci. Technol. A 6, 2390 (1988).
 (11) M. Cross and H. Greenside, Pattern Formation and Dynamics in Nonequilibrium Systems (Cambridge University Press, Cambridge, 2009).
 (12) G. Ozaydin, A. S. Ozcan, Y. Wang, K. F. Ludwig, H. Zhou, R. L. Headrick, and D. P. Siddons, Appl. Phys. Lett. 87 (2005) 163104.
 (13) C. S. Madi, B. Davidovitch, H. B. George, S. A. Norris, M. P. Brenner, and M. J. Aziz, Phys. Rev. Lett. 101, 246102 (2008).
 (14) C. S. Madi, H. B. George, and M. J. Aziz, J. Phys.: Condens. Matter 21, 224010 (2009).
 (15) S. Macko, F. Frost, B. Ziberi, D. F. Förster, Th. Michely, Nanotechnology 21, 085301 (2010).
 (16) S. Macko, F. Frost, M. Engler, D. Hirsch, T. Höche, J. Grenzer, and T. Michely, New J. Phys. 13, 073017 (2011).
 (17) G. Carter and V. Vishnyakov, Phys. Rev. B 54, 17647 (1996).
 (18) V. Shenoy, W. L. Chan, and E. Chason, Phys. Rev. Lett. 98, 256101 (2007).
 (19) S. Le Roy, E. Sondergard, I. S. Nerbo, M. Kildemo, and M. Plapp, Phys. Rev. B 81, 161401 (2010).
 (20) R. M. Bradley and P. D. Shipman, Phys. Rev. Lett. 105, 145501 (2010).
 (21) J. Zhou and M. Lu, Phys. Rev. B 82, 125404 (2010).
 (22) M. Makeev, R. Cuerno, and A.L. Barabási, Nucl. Inst. Meth. Phys. Res. B 197, 185 (2002).
 (23) S. Facsko, T. Bobek, A. Stahl, H. Kurz and T. Dekorsy, Phys. Rev. B 69 153412 (2004).
 (24) B. Davidovich, M. Aziz, and M. Brenner, Phys. Rev. B 76, 205420 (2007).
 (25) M. Castro, R. Cuerno, L. Vázquez, and R. Gago, Phys. Rev. Lett. 94, 016102 (2005).
 (26) J. MuñozGarcía, M. Castro, R. Cuerno, Phys. Rev. Lett. 96, 086101 (2006).
 (27) J. MuñozGarcía, R. Cuerno, and M. Castro, Phys. Rev. B 78, 205408 (2008).
 (28) J. MuñozGarcía, R. Cuerno, and M. Castro, J. Phys. Condens. Matter 21, 224020 (2009).
 (29) L. Pelaz, L. A. Marqués, and J. Barbolla, J. Appl. Phys. 96, 5947 (2004).
 (30) N. Kalyanasundaram, M. Wood, J. B. Freund, and H. T. Johnson, Mech. Res. Comm. 35, 50 (2008).
 (31) M. Moseler, P. Gumbsch, C. Casiraghi, A. C. Ferrari, and J. Robertson, Science 309, 1545 (2005).
 (32) M. C. Moore, N. Kalyanasundaram, J. B. Freund, and H. T. Johnson, Nucl. Inst. Meth. Phys. Res. B 225, 241 (2004).
 (33) A. Oron, S. H. Davis, and S. G. Bankoff, Rev. Mod. Phys. 69, 931 (1997).
 (34) R. Kree, oral communication at International Conference on IonBeam Induced Nanopatterning of Materials (IINM2011), 0610 February, 2011.
 (35) M. Castro and R. Cuerno, arXiv:1007.2144v1 (2011).
 (36) S. E. Orchard, Appl. Sci. Res. 11A, 451 (1962).
 (37) W. W. Mullins, J. Appl. Phys. 30, 77 (1959).
 (38) E. Chason, T. M. Mayer, B. K. Kellerman, D. T. Mcilroy, and A. J. Howard, Phys. Rev. Lett. 72, 3040 (1994).
 (39) F. Frost, R. Fechner, B. Ziberi, J. Völlner, D. Flamm, and A. Schindler, J. Phys. Condens. Matter 21, 224026 (2009).
 (40) S. A. Norris, J. Samela, L. Bukonte, M. Backman, F. Djurabekova, K. Nordlund, C. S. Madi, M. P. Brenner, and M. J. Aziz, Nat. Comm. 2, 276 (2011).
 (41) S.P Kim, B.H. Kim, H. Kim, K.R. Lee, Y.C. Chung, J. Seo, and J.S. Kim, Nucl. Instrum. Meth. Phys. Res. B 269, 2605 (2011).
 (42) S. Vauth and S. G. Mayr, Phys. Rev. B 75, 224107 (2007).
 (43) T. Aste and U. Valbusa, Physica A 332, 548 (2004); New J. Phys. 7, 122 (2005).
 (44) B.J. Spencer, S.H. Davis and P.W. Voorhees, Phys. Rev. B 47 (1993) 9760.
 (45) JH. Kim, N.B. Ha, J.S. Kim, M. Joe, K.R. Lee, and R. Cuerno, Nanotechnology 22, 285301 (2011).