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 , 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.
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 , with .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
II. THE MODIFIED SIGMUND MODEL
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
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 with the z-axis as their axis of symmetry. Note that the integral of 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 . 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
as shown in Ref. 19.
If the density of deposited energy is integrated over the interior of the solid, the result is less than the ion energy ϵ. This difficulty arises because 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 .
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 to be given by Eq. (2), but we replace the angle θ in Eqs. (3) and (4) by the larger angle ψ to yield
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 .
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 .
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 , and so it yields a value of ER that is greater than the one the SMS gives. However, even if 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 is proportional to the power deposited at the solid surface P, i.e.,
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 and P in the Appendix and demonstrate that the nonlinearity has only a modest effect on the early-time dynamics of the surface.
III. EQUATION OF MOTION
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
where the C's are constants that depend on θ. For later use, we note that if we set , Eq. (8) becomes
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
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 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
where H(x, y) is given by Eq. (10).
We introduce the moments of the crater function
As shown by Harrison and Bradley,
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 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 and for the MSMS by replacing θ by in the corresponding expressions for the Sigmund model that are given in Ref. 26. This yields
Note that by definition ,
A, B1, B2, and C are all positive for . It follows that T11 and T22 are positive for values of θ that have . For , 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 .
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 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, is negative for both values of i. Arguing in the same way, we arrive at the conclusion that is positive for 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 . Conversely, the curvature dependence of the crater function produces a destabilizing effect for θ such that .
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 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 ; this lies between 45° and 90°.
The sputter yield for a flat surface will be denoted by . Equation (8) shows that the erosion rate v0 for a flat surface is . This rate is also equal to , where Ω is the atomic volume. Using Eq. (15), we arrive at the conclusion that
A straightforward analysis shows that the crater function moment is a function of ψ that has a maximum at . From this it follows that has a maximum at , 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 .
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, , attains its maximum value at . 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 that attains a maximum at a value of θ between and , it does not produce a sputter yield that falls to zero for . However, if α∕a and β∕a are reasonably small and ψ is close to for , 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 and are negative for . 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.
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.
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 as θ is increased through a critical value . 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 that does not coincide with θc, whereas C22 is negative and less than C11 for . Replacing the SMS by the MSMS does have the effect of altering the value of , 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 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 is a decreasing function of Kii for i = 1 and 2, both and are positive for , as we found in Sec. III.
and are positive for θ less than a critical value θc for the modified Sigmund model we introduced in this paper. In contrast, and are negative for . This means that in this regime, an ion striking the concave part of the surface 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 , 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 .
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 . If R/a is sufficiently large, we may approximate the surface height of the sphere by the paraboloid
where and .
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 and retain only terms of first order in K11 and K22. This yields
Recalling the definitions of T11 and T22, we see that the average sputter yield is
The results of Nietiadi et al. therefore imply that is positive for and that it is negative for , 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 , 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 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 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 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 . 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 is proportional to the power deposited at the solid surface P. Molecular dynamics simulations have shown that 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.
APPENDIX: NONLINEAR DEPENDENCE OF THE SPUTTER YIELD ON THE DEPOSITED POWER
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.,
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.
As before, we will only retain terms of first order in h. Equation (A3) therefore reduces to
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
and introduce the rescaled time
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
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.
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 .