We modify the Sigmund model of ion sputtering so that the distribution of deposited energy is in better accord with the results of molecular dynamics simulations. For a flat target surface, the model gives a sputter yield that displays a maximum as the angle of ion incidence is increased, as observed experimentally. The model also predicts that for sufficiently high angles of incidence, the curvature dependence of the crater function tends to destabilize the solid surface.

Bombarding a solid surface with a broad ion beam can produce a remarkable variety of self-assembled nanoscale patterns with feature sizes that can be as small as 10 nm.1,2 These patterns include periodic height modulations or “ripples”3 and mounds arranged in hexagonal arrays of astonishing regularity.4,5 The spontaneous emergence of these patterns is not just fascinating in its own right: Ion bombardment has the potential to become a cost-effective method to rapidly fabricate large-area nanostructures at length scales beyond the limits of conventional optical lithography. However, before the full potential of this nanofabrication method can be realized, a more complete understanding of it will be needed.

The underpinning of many theories of nanoscale pattern formation induced by ion bombardment is the Sigmund model of sputtering (SMS).6 The Bradley-Harper (BH) theory, for example, provides an explanation for why oblique-incidence ion bombardment can lead to the formation of ripples on solid surfaces.7 In that theory, a partial differential equation is derived that gives the time evolution of a nearly flat solid surface due to the combined effects of sputtering and surface diffusion; the theory applies to long wavelength disturbances at sufficiently early times. The SMS is used to describe the effects of sputtering in the BH theory. Since the BH theory was introduced in 1988, many extensions of it have been developed in which its limitations are to some extent surmounted.8–19 All of these theories are based on the Sigmund model.

In the SMS, a collimated ion beam produces ellipsoidal contours of equal power deposition within the solid, as shown in Fig. 1(a). Let the origin O be the point where the beam first makes contact with the solid and let ê denote the ion beam direction. The common centroid of the ellipsoids is located at aê, i.e., it is a distance a along the beam direction from the point O. The ellipsoids are rotationally symmetric about the beam direction. Finally, Sigmund made the natural assumption that the erosion rate at an arbitrary point on the solid surface is proportional to the power deposited there.

FIG. 1.

(a) The Sigmund model of ion sputtering. Ions that strike the solid surface at the point O produce ellipsoidal contours of equal power deposition within the solid with their centroid at aê. The direction of ion incidence is ê, where êx̂sinθ+ẑcosθ and θ is the angle of incidence. (b) The modified Sigmund model of ion sputtering. The contours of equal power deposition are the same as those the Sigmund theory would give for an angle of incidence ψ that is greater than θ.

FIG. 1.

(a) The Sigmund model of ion sputtering. Ions that strike the solid surface at the point O produce ellipsoidal contours of equal power deposition within the solid with their centroid at aê. The direction of ion incidence is ê, where êx̂sinθ+ẑcosθ and θ is the angle of incidence. (b) The modified Sigmund model of ion sputtering. The contours of equal power deposition are the same as those the Sigmund theory would give for an angle of incidence ψ that is greater than θ.

Close modal

Although Sigmund's model has been remarkably successful, it does not agree with all observed aspects of sputtering. According to the theory, the sputter yield from a flat surface increases with the angle of incidence θ for all θ > 0. Many experiments have shown that this sputter yield actually increases up until a critical angle and then begins to decrease. The yield drops to zero for θ = 90°, i.e., for grazing-incidence bombardment. This discrepancy occurs in part because the Sigmund model underestimates the effect of backscattering of ions, an effect that becomes increasingly important as θ approaches 90°.

Molecular dynamics simulations have confirmed that the contours of equal power deposition are to a good approximation ellipsoidal for several choices of target material and ion beams,20 as in the SMS. However, as the angle of incidence is increased from zero, the ellipsoids' common centroid deviates significantly from the ray given by r=lê, with l0.20 The contours of equal power deposition in fact strongly resemble those the Sigmund theory would give for a higher angle of incidence than θ. These contours are therefore nearly the same as those we would have if the ions were suddenly deflected as they entered the solid.

In this paper, we introduce a modified Sigmund model of sputtering (MSMS) that includes the apparent deflection that occurs as an ion passes from the vacuum to the target material. We derive a continuum equation of motion for the solid surface that applies for early times if at t = 0 the surface is nominally flat and its height varies slowly with position. For a flat target surface, the model gives a sputter yield that displays a maximum as θ is increased, as observed experimentally. The MSMS also leads to the prediction that for sufficiently high angles of incidence, an ion that strikes a concave portion of the surface sputters a greater number of atoms on average than an ion incident on a convex surface region. Recent Monte Carlo simulations of the sputtering of amorphous silicon with 20 keV argon ions provide support for the validity of this remarkable and counterintuitive prediction.21,22

This paper is organized as follows. In Sec. II, we introduce our modified Sigmund model. The resulting equation of motion for a solid surface with slowly varying height is determined in Sec. III. We discuss our findings in Sec. IV and conclude the paper in Sec. V.

Consider an ion of energy ϵ that impinges on the surface of an elemental solid at the point P. By shifting the location of the origin O if necessary, we can arrange for O and P to coincide. The height of the surface above the x–y plane, h = h(x,y,t), then vanishes for x = y = 0.

In the Sigmund theory of sputtering, the rate with which material is sputtered from a point on the solid surface is proportional to the power deposited there by the random slowing down of ions. The average energy density deposited at a point (x, y, z) within the solid by an ion which travels along the z axis until striking the surface is taken to be

ED(x,y,z)=ϵ(2π)3/2αβ2exp((z+a)22α2x2+y22β2).
(1)

Here a is the average energy deposition depth and α and β are the longitudinal and transverse straggling lengths, respectively. The contours of equal energy deposition are ellipsoids of revolution centered at the point aẑ with the z-axis as their axis of symmetry. Note that the integral of ED(x,y,z) over all space is the ion energy ϵ.

If the angle of incidence θ is positive rather than zero, the direction of the incident beam is ê, where êx̂sinθ+ẑcosθ. In the Sigmund model, the average distribution of deposited energy is obtained by rotating the distribution (1) through the angle θ about the y-axis, as shown in Fig. 1(a). Explicitly, the density of deposited energy at an arbitrary point r within the solid is given by

ED(r)=ϵ(2π)3/2αβ2exp(ρ22α2ρ22β2),
(2)

where

ρ=axsinθ+zcosθ
(3)

and

ρ=[(xcosθ+zsinθ)2+y2]1/2,
(4)

as shown in Ref. 19.

If the density of deposited energy ED(r) is integrated over the interior of the solid, the result ED,tot is less than the ion energy ϵ. This difficulty arises because ED(r) is nonzero in the vacuum. The Sigmund model makes sense physically, however, if we take the average energy that is carried off by a backscattered ion to be ER=ϵED,tot.

In our modified Sigmund model of ion sputtering, we take into account the apparent deflection of the ion as it enters the solid, as shown in Fig. 1(b). Thus, we take ED(r) to be given by Eq. (2), but we replace the angle θ in Eqs. (3) and (4) by the larger angle ψ to yield

ρ=axsinψ+zcosψ
(5)

and

ρ=[(xcosψ+zsinψ)2+y2]1/2.
(6)

The distribution of deposited energy in the MSMS is therefore the same as it would be in the Sigmund theory if the angle of incidence were ψ rather than θ. The apparent deflection angle is ψθ.

The angle ψ will be taken to be a given function of θ. Clearly, ψ = 0 for θ = 0. In addition, ψ(θ) must exceed θ for θ > 0. We will not specify an explicit form for the function ψ(θ) for the time being. We assume without loss of generality that θ0.

For a given choice of target material, ion species, energy ϵ and nonzero angle of incidence θ, the angle ψ could be determined using the results of molecular dynamics or Monte Carlo simulations of the impact of an ion on a flat surface. The distribution of deposited energy for normal-incidence bombardment would first be averaged over a number of ion impacts and fitted to the form (1) to yield estimates for the parameters a, α and β. Simulations of multiple ion impacts for the angle of incidence θ would then be carried out. Once the average distribution of deposited energy had been computed, it would be fitted to the form given by Eqs. (2), (5), and (6) to give a value for the parameter ψ. After repeating this procedure for a range of θ values, an empirical form for the function ψ(θ) could be obtained.

The point of maximum energy deposition was on the surface of the solid for sufficiently large values of θ in the molecular dynamics simulations of Hossain, Freund, and Johnson.20 In addition, the maximum density of energy deposition decreased as θ increased in this regime. Therefore, if the fitting procedure just described were applied to their data, the angle ψ would exceed 90° for sufficiently large θ. The common centroid of the ellipsoids is in the vacuum, above the surface of the solid, for ψ>90°.

For ion bombardment of a flat surface, the SMS gives a value of ER approaches ϵ∕2 as the angle of incidence θ tends to 90°. This is problematic because in this limit, all of the ions should be backscattered and ER should tend to ϵ. Our MSMS partially remedies this difficulty because ψ(θ)>θ for θ>0°, and so it yields a value of ER that is greater than the one the SMS gives. However, even if ψ(90°) were equal to 180°, the MSMS value of ER for glancing-incidence bombardment would be still be less than ϵ.

In the SMS, it is assumed that the erosion rate h/t is proportional to the power deposited at the solid surface P, i.e.,

ht=ΛP,
(7)

where the value of the constant of proportionality Λ depends on the choice of target material and the ion species and energy. Hossain, Freund, and Johnson tested the validity of this assumption in their molecular dynamics simulations and found that in fact the erosion rate is not simply proportional to the power deposited.20,23 In our MSMS, we will assume that Eq. (7) holds for the sake of simplicity. We will, however, consider a nonlinear relationship between h/t and P in the Appendix and demonstrate that the nonlinearity has only a modest effect on the early-time dynamics of the surface.

We will now develop an equation of motion for the solid surface that applies for our MSMS. Suppose that the initial surface height h(x,y,0) is small and varies slowly with x and y. We will only take into account the effect of sputtering, and so will neglect mass redistribution,9,14,24 surface diffusion and ion-induced viscous flow.25 The modified Sigmund model will be taken to describe the effects of sputtering and we will confine our attention to times early enough that we may work to first order in h. We will also work to second order in the space derivatives. The goal of our analysis is to derive an equation of motion of the form

1Jht=C0+C1hx+C112hx2+C222hy2,
(8)

where the C's are constants that depend on θ. For later use, we note that if we set u(x,y,t)=h(x,y,t)JC0t, Eq. (8) becomes

1Jut=C1ux+C112ux2+C222uy2.
(9)

To find the constant coefficients C0, C1, C11, and C22 in Eq. (8), we will apply the crater function formalism of Harrison and Bradley.26 Consider a single ion incident at the origin on a solid with the surface height

H(x,y)=12K11x2+K12xy+12K22y2
(10)

above the point (x, y) in the x–y plane. The ion species, energy ϵ and angle of incidence θ as well as the composition of the solid are the same as in our original problem. By definition, the crater function F(x,y,θ,K11,K12,K22) is minus the average change in the height H above the point (x, y) in the x–y plane that results from the impact of the ion. F is not just a function of the angle of incidence—it also depends on the curvature of the surface at the origin, i.e., on the constants K11, K12, and K22. The crater function for the MSMS is

F(x,y,θ,K11,K12,K22)=Λϵ(2π)3/2αβ2exp(12α2[axsinψ+H(x,y)cosψ]212β2[xcosψ+H(x,y)sinψ]212β2y2),
(11)

where H(x, y) is given by Eq. (10).

We introduce the moments of the crater function

MK11F(x,y,θ,K11,0,0)dxdy,
(12)
MK22F(x,y,θ,0,0,K22)dxdy,
(13)

and

Mx(n)F(x,y,θ,0,0,0)xndxdy.
(14)

As shown by Harrison and Bradley,

C0(θ)=Mx(0)cosθ,
(15)
C1(θ)=θ(Mx(0)cosθ),
(16)
C11(θ)=S11(θ)+T11(θ),
(17)

and

C22(θ)=S22(θ)+T22(θ),
(18)

where

S11(θ)θ(Mx(1)cosθ),
(19)
S22(θ)Mx(1)cosθcotθ,
(20)

and

Tii(θ)cosθKiiMKii|Kii=0,
(21)

for i = 1 and 2.

Tii(θ) is the contribution to Cii(θ) that comes from the curvature dependence of the crater function. If this dependence is neglected, then we obtain Cii(θ)=Sii(θ) for i = 1 and 2. This truncated crater function approximation has been shown to be quite inaccurate in the case of the Sigmund model.26 For the SMS with θ = 0, for example, the values of C11 and C22 obtained by neglecting the curvature dependence of the crater function are a factor of two larger in magnitude than their exact values.

We can obtain Mx(0),Mx(1) and MKii/Kii|Kii=0 for the MSMS by replacing θ by ψ=ψ(θ) in the corresponding expressions for the Sigmund model that are given in Ref. 26. This yields

Mx(0)=2πDaβB1eaα2/2exp(A22B1),
(22)

where

DΛϵa2(2π)3/2αβ2.
(23)

In addition,

Mx(1)=aAB1Mx(0),
(24)
K11MK11|K11=0=aMx(0)(A2B22B12+A3CB13+B22B1+3ACB12),
(25)

and

K22MK22|K22=0=Mx(0)aaβ2(B22+ACB1).
(26)

Note that by definition aα=a/α,aβ=a/β,

A=aα2sinψ,
(27)
B1=aα2sin2ψ+aβ2cos2ψ,
(28)
B2=aα2cosψ,
(29)

and

C=12(aβ2aα2)cosψsinψ.
(30)

Inserting these results in Eqs. (15)–(21) yields the constant coefficients in the equation of motion (8).

A, B1, B2, and C are all positive for 0<ψ<90°. It follows that T11 and T22 are positive for values of θ that have ψ(θ)<90°. For 90°<ψ(θ)<180°, on the other hand, A and B1 are positive but B2 and C are negative. As a consequence, T11 and T22 are negative for values of θ with 90°<ψ(θ)<180°.

It is actually easy to see why the signs of T11 and T22 show this behavior. The ellipsoids' common centroid lies within the solid for ψ(θ)<90° and so the average sputter yield for the impact of a single ion is a decreasing function of Kii for i = 1 and 2. Thus, MKii/Kii|Kii=0 is negative for both values of i. Arguing in the same way, we arrive at the conclusion that MKii/Kii|Kii=0 is positive for 90°<ψ(θ)<180° for both i = 1 and 2.

This analysis leads us to the conclusion that the curvature dependence of the crater function has a smoothing influence for values of θ that have ψ(θ)<90°. Conversely, the curvature dependence of the crater function produces a destabilizing effect for θ such that 90°<ψ(θ)<180°.

The molecular dynamics simulations of Hossain, Freund and Johnson strongly suggest that ψ is an increasing function of θ that does not exceed 180°.20 Let us assume that this is the case and let θc denote the value of θ that separates the regimes in which ψ(θ) is greater or less than 90°. The analysis given in the preceding paragraphs shows that the Tii's are positive for 0θ<θc and are negative for θ > θc.

To further investigate the behavior predicted by our modified Sigmund model, we need a specific form for the function ψ=ψ(θ). As a simple model that is consistent with all of the requirements we have placed on this function, we will assume that ψ(θ)=ρθ, where ρ is a constant that satisfies the inequality 1 < ρ < 2. This choice of ψ vanishes for θ = 0, is an increasing function of θ, is greater than θ for θ > 0, and is less than 180°. The critical value θc is π/(2ρ); this lies between 45° and 90°.

The sputter yield for a flat surface will be denoted by Yflat(θ). Equation (8) shows that the erosion rate v0 for a flat surface is JC0. This rate is also equal to ΩJYflatcosθ, where Ω is the atomic volume. Using Eq. (15), we arrive at the conclusion that

Yflat=Mx(0)/Ω.
(31)

A straightforward analysis shows that the crater function moment Mx(0) is a function of ψ that has a maximum at ψ=π/2. From this it follows that Yflat(θ) has a maximum at θc=π/(2ρ), as illustrated in Fig. 2 for the parameter values appropriate for bombardment of silicon with a 1.0 keV beam of argon ions.27 Thus, in our MSMS, the sputter yield for a flat surface has a maximum for a value of θ that lies between 45° and 90°, as observed experimentally. An estimate of ρ could be deduced from an experiment or a simulation simply by determining the angle θc that gives the maximum sputter yield Yflat(θ).

FIG. 2.

The sputter yield for a flat surface versus the angle of incidence θ in deg. The sputter yield has been normalized by its value for normal incidence. The parameters were taken to have the values α/a=0.75,β/a=0.69, and ρ=1.4.

FIG. 2.

The sputter yield for a flat surface versus the angle of incidence θ in deg. The sputter yield has been normalized by its value for normal incidence. The parameters were taken to have the values α/a=0.75,β/a=0.69, and ρ=1.4.

Close modal

As we have seen, in the MSMS, T11 and T22 change sign from positive to negative at the critical angle θc as θ is increased. In addition, the sputter yield for a flat surface, Yflat(θ), attains its maximum value at θ=θc. It is important to note that θc is not the critical angle for the onset of ripple formation or for ripple rotation.

Although the MSMS gives a sputter yield Yflat(θ) that attains a maximum at a value of θ between 45° and 90°, it does not produce a sputter yield that falls to zero for θ=90°. However, if α∕a and β∕a are reasonably small and ψ is close to 180° for θ=90°, the MSMS sputter yield will be quite small for grazing incidence.

T11 and T22 are plotted versus θ in Fig. 3 for the same set of parameter values as we employed in Fig. 2. As we have already discussed, T11 and T22 are positive for 0<θ<θc=π/(2ρ) and are negative for θc<θ<π/2. We also see that T11 has a greater magnitude than T22 for the parameters used to make the plot except when both Tii's are zero. This is in fact true so long as β<α, a condition that is always satisfied in practice.

FIG. 3.

T11 (black solid curve) and T22 (red dashed curve) in units of 2πβaα2Deaα2/2 plotted versus the angle of incidence θ in deg. The parameters were taken to have the same values as in Fig. 2.

FIG. 3.

T11 (black solid curve) and T22 (red dashed curve) in units of 2πβaα2Deaα2/2 plotted versus the angle of incidence θ in deg. The parameters were taken to have the same values as in Fig. 2.

Close modal

To get an idea of how accurate the truncated crater function approximation is for the MSMS, we plotted C11, S11, C22, and S22 in Fig. 4. We again used the same parameter values as we employed in producing Fig. 2. As in the case of the SMS, this approximation incurs the largest errors for relatively small values of θ, but there are significant discrepancies between its predictions and the exact values of the coefficients C11 and C22 for nearly all angles of incidence.

FIG. 4.

C11 (black solid curve), S11 (red long dashes), C22 (blue dashes), and S22 (purple dashed-dotted) in units of 2πβaα2Deaα2/2 plotted versus the angle of incidence θ in deg. S11 and S22 are the values that the truncated crater function approximation yields for C11 and C22. The parameters were taken to have the same values as in Figs. 2 and 3.

FIG. 4.

C11 (black solid curve), S11 (red long dashes), C22 (blue dashes), and S22 (purple dashed-dotted) in units of 2πβaα2Deaα2/2 plotted versus the angle of incidence θ in deg. S11 and S22 are the values that the truncated crater function approximation yields for C11 and C22. The parameters were taken to have the same values as in Figs. 2 and 3.

Close modal

The BH theory is based upon the Sigmund model of ion sputtering, as we have already noted. The theory predicts that ripples form for all nonzero angles of incidence θ and that the ripple orientation changes by 90° as θ is increased through a critical value θ̃c. If the modified Sigmund model of ion sputtering is instead adopted as the basis of the theory, these predictions are unchanged: as illustrated in Fig. 4, C11 is negative and less than C22 for θ less than a critical value θ̃c that does not coincide with θc, whereas C22 is negative and less than C11 for θ̃c<θ<90°. Replacing the SMS by the MSMS does have the effect of altering the value of θ̃c, however.

Ion-induced mass redistribution plays an important role in determining whether ripples form when silicon is bombarded with a broad ion beam for a range of ion energies.27 It also influences the wavelength and orientation of the ripples if they do form. These conclusions carry over to Ge targets28 and most likely to other materials as well. For materials of this kind, the erosive contribution to C11 and C22 could be well described by the MSMS, but the redistributive contribution to these coefficients cannot be neglected. Thus, whether or not ripples form and how they are oriented will be strongly influenced by mass redistribution and the predictions of the modified BH theory described in the preceding paragraph will not apply.

Suppose that an ion strikes a convex region on a solid surface S and a second ion is incident on a concave portion of the same surface. We will consider the case in which the local angles of incidence of the two ions are equal. In the Sigmund model of ion sputtering, the first ion sputters a greater number of atoms on average than the second ion.21,26,29 It is intuitively clear why this is so: the convex portion of the solid surface is closer to the centroid of the distribution of deposited energy than the concave portion is. Because the crater function moment MKii is a decreasing function of Kii for i = 1 and 2, both T11(θ) and T22(θ) are positive for 0°θ<90°, as we found in Sec. III.

T11(θ) and T22(θ) are positive for θ less than a critical value θc for the modified Sigmund model we introduced in this paper. In contrast, T11(θ) and T22(θ) are negative for θ>θc. This means that in this regime, an ion striking the concave part of the surface S sputters more atoms on average than an ion incident on the convex area. As we will now discuss, Monte Carlo simulations of the sputtering of a sphere support the validity of this surprising prediction.

Nietiadi et al. have recently carried out Monte Carlo simulations of the sputtering of a sphere composed of amorphous silicon by a single 20 keV argon ion.21 Their results show that for angles of incidence below a critical value θc, the average sputter yield is a decreasing function of the sphere's radius R for large R/a. Their data also show that for θ>θc, the average sputter yield is an increasing function of R if R/a is sufficiently large, although Nietiadi et al. did not comment on this aspect of their data in their published paper.21,22 The critical angle θc is approximately 83°.

To relate these observations to the results of the present paper, consider a single ion incident upon a solid sphere of radius R. Let θ denote the local angle of incidence. We place the origin at the point of impact and put the z-axis along the outward-pointing local surface normal. The x-axis will be oriented so that the direction of the incident ion ê lies in the x–z plane and ê=x̂sinθ+ẑcosθ. If R/a is sufficiently large, we may approximate the surface height of the sphere by the paraboloid

h(x,y)=12K11x2+K12xy+12K22y2,
(32)

where K11=K22=1/R and K12=0.

Let

M(θ,K11,K12,K22)F(x,y,θ,K11,K12,K22)dxdy.
(33)

For a sphere of large radius, K11 and K22 are small and K12 is zero. We may therefore carry out a Taylor series expansion of M(θ,K11,K12,K22) and retain only terms of first order in K11 and K22. This yields

M=M(θ,1/R,0,1/R)M(θ,0,0,0)(K11MK11|K11=0+K22MK22|K22=0)1R.
(34)

Recalling the definitions of T11 and T22, we see that the average sputter yield is

Y=Yflat+(T11+T22)secθΩR.
(35)

The results of Nietiadi et al. therefore imply that T11+T22 is positive for θ<θc and that it is negative for θc<θ<90°, in accord with the prediction of the MSMS.

The curvature dependence of the crater function has a smoothing effect for the SMS for arbitrary angles of incidence. Based on this finding, Harrison and Bradley implicitly assumed that this is so for an arbitrary choice of ion beam and target material.26 For the MSMS, however, T11 and T22 are negative for θc<θ<90°, and so the curvature dependence of the crater function tends to destabilize the solid surface for sufficiently large angles of incidence. The validity of this prediction is supported by the simulations of Nietiadi et al. Thus, although Harrison and Bradley's crater function formalism is perfectly correct, the curvature dependence of the crater function does not always have a stabilizing influence on the surface dynamics.

As we noted in Sec. II, in formulating his model of ion sputtering, Sigmund made the natural assumption that the erosion rate h/t is proportional to the power deposited at the solid surface P, i.e., Eq. (7) holds.6 The molecular dynamics simulations of Hossain, Freund, and Johnson, however, have revealed that the erosion rate is in fact a nonlinear function of P.20,23 In the MSMS studied in this paper, we assumed that Eq. (7) is valid for the sake of simplicity. In the Appendix, we analyze a variant of our MSMS in which the power deposited is still given by Eqs. (2), (5), and (6) but h/t is an arbitrary nonlinear function of P. The overall rate of erosion JC0 is modified by this change, as we would expect a priori. At early times, the linearized equation of motion is a good approximation if the surface is initially nearly flat. In this regime, the only change in the dynamics of the surface morphology that is brought about by the nonlinear relationship between h/t and P is a change in the rate of the time evolution. This change can also be thought of as the replacement of the constant Λ by an effective value Λeff. These observations carry over to the SMS since it is a special case of our MSMS.

Modifications of the Sigmund model of sputtering have been considered previously. Feix et al. carried out simulations of the impact of Cu ions on a copper surface using the binary-collision approximation and found that the distribution of deposited energy was not fit well by Sigmund's Gaussian form.30 The distribution that best fit their simulations instead has a minimum near the position where the ion penetrates the surface, and the decay of energy deposition with the distance from the ion trajectory is exponential rather than Gaussian.31 In addition, in order to gauge how robust the predictions of the BH theory are, Davidovitch, Aziz, and Brenner found C11 and C22 for a variety of non-Gaussian energy distributions.14 However, all of the distributions of deposited energy considered in these two studies were rotationally symmetric about the direction of ion incidence ê, in contrast to the model studied in this paper.

In this paper, we modified the Sigmund model of ion sputtering so that the distribution of deposited energy is in better accord with the results of the molecular dynamics simulations of Hossain, Freund, and Johnson.20 Using the crater function formalism of Harrison and Bradley,26 we derived a continuum equation of motion for the solid surface that applies for early times if at t = 0 the surface is nominally flat and its height varies slowly with position. For a flat target surface, the model gives a sputter yield that displays a maximum as the angle of incidence is increased, in accord with experiment. The modified Sigmund model also leads to the prediction that the curvature dependence of the crater function tends to destabilize the solid surface for sufficiently high angles of incidence. Recent Monte Carlo simulations of the sputtering of amorphous silicon with 20 keV argon ions provide support for the validity of this prediction.21,22

In the Sigmund model of ion sputtering, it is assumed that the erosion rate h/t is proportional to the power deposited at the solid surface P. Molecular dynamics simulations have shown that h/t is actually a nonlinear function of P, as one would expect for sufficiently large values of P.20 We demonstrated that this nonlinearity has only a modest effect on the dynamics of the surface morphology at early times: it simply alters the rate of the time evolution.

We have benefited from helpful discussions and correspondence with Matt Harrison, Wolfhard Möller, Patrick Shipman, Peter Sigmund, and Herbert Urbassek. R.M.B. would like to thank the National Science Foundation for its support through Grant DMR-1305449.

Our analysis of the MSMS in which Eq. (7) holds led to the equation of motion (8) with coefficients given by Eqs. (15)–(20). As a consequence

P=JΛ(C0+C1hx+C112hx2+C222hy2).
(A1)

This expression does not depend on the form of the relationship between the erosion rate and the power deposited.

Let us now suppose that Eq. (7) does not apply and that instead the rate of erosion is a nonlinear function of the deposited power, i.e.,

ht=V(P),
(A2)

where V is a given function of P. In this Appendix, we will determine how the equation of motion obtained in Sec. III is altered by this change to the MSMS.

Combining Eqs. (A1) and (A2), we obtain the equation of motion

ht=V(JΛ(C0+C1hx+C112hx2+C222hy2)).
(A3)

As before, we will only retain terms of first order in h. Equation (A3) therefore reduces to

1Jht=1JV(JΛC0)+1ΛV(JΛC0)(C1hx+C112hx2+C222hy2).
(A4)

We see that the erosion rate for a flat surface is altered in the nonlinear version of the MSMS. To see the effect on the time evolution of the surface morphology, we let

u(x,y,t)h(x,y,t)V(JΛC0)t
(A5)

and introduce the rescaled time

t̃V(JΛC0)tΛ.
(A6)

We obtain

1Jut̃=C1ux+C112ux2+C222uy2.
(A7)

A comparison of Eqs. (9) and (A7) shows that the only change in the dynamics of the surface morphology is a change in the rate of the time evolution. Stated in another way, the time evolution of the surface morphology in the nonlinear MSMS is obtained by replacing Λ by the new, effective value given by

Λeff=V(JC0/Λ).
(A8)
1.
For example, see
B.
Ziberi
,
M.
Cornejo
,
F.
Frost
, and
B.
Rauschenbach
,
J. Phys.: Condens. Matter
21
,
224003
(
2009
).
2.
Q.
Wei
,
J.
Lian
,
S.
Zhu
,
W.
Li
,
K.
Sun
, and
L.
Wang
,
Chem. Phys. Lett.
452
,
124
(
2008
).
3.
For a review, see
J.
Muñoz-García
,
L.
Vázquez
,
M.
Castro
,
R.
Gago
,
A.
Redondo-Cubero
,
A.
Moreno-Barradoc
, and
R.
Cuerno
,
Mater. Sci. Eng., R
86
,
1
(
2014
).
4.
S.
Facsko
,
T.
Dekorsy
,
C.
Koerdt
,
C.
Trappe
,
H.
Kurz
,
A.
Vogt
, and
H. L.
Hartnagel
,
Science
285
,
1551
(
1999
).
5.
F.
Frost
,
A.
Schindler
, and
F.
Bigl
,
Phys. Rev. Lett.
85
,
4116
(
2000
).
6.
P.
Sigmund
,
J. Mater. Sci.
8
,
1545
(
1973
).
7.
R. M.
Bradley
and
J. M. E.
Harper
,
J. Vac. Sci. Technol., A
6
,
2390
(
1988
).
8.
R.
Cuerno
and
A.-L.
Barabási
,
Phys. Rev. Lett.
74
,
4746
(
1995
).
9.
G.
Carter
and
V.
Vishnyakov
,
Phys. Rev. B
54
,
17647
(
1996
).
10.
M. A.
Makeev
and
A.-L.
Barabási
,
Appl. Phys. Lett.
71
,
2800
(
1997
).
11.
M. A.
Makeev
,
R.
Cuerno
, and
A.-L.
Barabási
,
Nucl. Instrum. Methods Phys. Res., Sect. B
197
,
185
(
2002
).
12.
M.
Castro
,
R.
Cuerno
,
L.
Vázquez
, and
R.
Gago
,
Phys. Rev. Lett.
94
,
016102
(
2005
).
13.
J.
Muñoz-García
,
M.
Castro
, and
R.
Cuerno
,
Phys. Rev. Lett.
96
,
086101
(
2006
).
14.
B.
Davidovitch
,
M. J.
Aziz
, and
M. P.
Brenner
,
Phys. Rev. B
76
,
205420
(
2007
).
15.
V. B.
Shenoy
,
W. L.
Chan
, and
E.
Chason
,
Phys. Rev. Lett.
98
,
256101
(
2007
).
16.
J.
Muñoz-García
,
R.
Cuerno
, and
M.
Castro
,
Phys. Rev. B
78
,
205408
(
2008
).
17.
R. M.
Bradley
and
P. D.
Shipman
,
Phys. Rev. Lett.
105
,
145501
(
2010
).
18.
P. D.
Shipman
and
R. M.
Bradley
,
Phys. Rev. B
84
,
085420
(
2011
).
19.
R. M.
Bradley
,
Phys. Rev. B
84
,
075413
(
2011
).
20.
M. Z.
Hossain
,
J. B.
Freund
, and
H. T.
Johnson
,
J. Appl. Phys.
111
,
103513
(
2012
).
21.
M. L.
Nietiadi
,
L.
Sandoval
,
H. M.
Urbassek
, and
W.
Möller
,
Phys. Rev. B
90
,
045417
(
2014
).
22.
W.
Möller
, private communication (
2014
).
23.

The simulations carried out in Ref. 20 show that the relation between the sputter yield and the deposited energy changes during the course of a single ion impact. Such an impact occurs over a time on the order of a few hundred femtoseconds. For the purposes of determining the behavior of the macroscopic surface over time intervals on the order of seconds or longer, however, only the net sputter yield per ion impact is needed. A nontrivial dependence of the sputter yield on the surface mass density ρs was also found in Ref. 20, even though the spatial variation of ρs was small. This is an interesting topic that is beyond the scope of the present paper.

24.
M.
Moseler
,
P.
Gumbsch
,
C.
Casiraghi
,
A. C.
Ferrari
, and
J.
Robertson
,
Science
309
,
1545
(
2005
).
25.
C. C.
Umbach
,
R. L.
Headrick
, and
K.-C.
Chang
,
Phys. Rev. Lett.
87
,
246104
(
2001
).
26.
M. P.
Harrison
and
R. M.
Bradley
,
Phys. Rev. B
89
,
245401
(
2014
).
27.
C. S.
Madi
,
E.
Anzenberg
,
K. F.
Ludwig
, Jr.
, and
M. J.
Aziz
,
Phys. Rev. Lett.
106
,
066101
(
2011
).
28.
E.
Anzenberg
,
J. C.
Perkinson
,
C. S.
Madi
,
M. J.
Aziz
, and
K. F.
Ludwig
, Jr.
,
Phys. Rev. B
86
,
245412
(
2012
).
29.
M. L.
Nietiadi
and
H. M.
Urbassek
,
Appl. Phys. Lett.
103
,
113108
(
2013
).
30.
M.
Feix
,
A. K.
Hartmann
,
R.
Kree
,
J.
Muñoz-García
, and
R.
Cuerno
,
Phys. Rev. B
71
,
125407
(
2005
).
31.

Interestingly, the sputter yield for a flat surface given by the model of Feix et al. has a maximum for an angle of incidence smaller than 90°.