Tests of the modified Sigmund model of ion sputtering using Monte Carlo simulations

Irradiating a solid surface with energetic ions at a given incidence angle can lead to the formation of self-organized nanoscale surface patterns, including ripple patterns of high regularity and dot patterns.^1–3 Most of the theoretical work performed in analyzing these patterns has been based on the continuum Bradley-Harper (BH) theory,⁴ which itself is based on the Sigmund model of ion sputtering.⁵ BH showed that for the Sigmund model the sputter yield at a point on the surface not only depends on the local angle of incidence but also on the local curvature of the surface.

In a recent paper, we introduced the Modified Sigmund Model of Ion Sputtering (MSMS).⁶ In this model, the energy deposited by an ion incident on a solid surface is described by Sigmund's ellipsoidal energy deposition distribution⁵ but the apparent deflection of the ion towards the surface when entering the solid is taken into account. The center of the ellipsoid is shifted towards the surface and may even be outside the surface at larger ion incidence angles. In this way, the model accounts for the increasing fraction of the incident energy that is carried away by backscattered ions. The resulting distribution of deposited energy within the solid is in better accord with the results of the molecular dynamics simulations of Hossain et al.⁷ than the unmodified Sigmund model.

To analyze the MSMS, the crater function formalism of Harrison and Bradley⁸ was used, since it takes into account the dependence of the crater function on the curvature of the surface around the point of ion impact. In contrast, in previous studies of pattern formation that used simulations as input to the crater function formalism, this curvature dependence was neglected and only crater functions of flat surfaces were considered.^3,7,9–13 Using Harrison and Bradley's formalism, we derived a continuum equation of motion for the solid surface that applies for early times if the surface is nominally flat and its height varies slowly with position at time t = 0. 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 experimental data.

The MSMS leads to the surprising prediction that the curvature dependence of the crater function tends to destabilize the solid surface for sufficiently high angles of incidence. For incidence angles smaller than the angle of the sputter yield maximum, the MSMS predicts a lower sputter yield for a surface with positive curvature (i.e., for a concave surface) and a higher sputter yield for negative curvature. At larger angles, i.e., when the center of Sigmund's ellipsoid is positioned outside the solid, the behavior reverses and the higher sputter yield is obtained for the surface with positive curvature. Recent Monte Carlo simulations of the sputtering of amorphous silicon spheres with 20 keV argon ions provide support for the validity of this prediction.¹⁴ 

In this work, we carry out Monte Carlo simulations using the program SDTrimSP¹⁵ to evaluate the angular and curvature dependence of the sputter yield. In addition, we compute the moments of the erosion crater functions for curved surfaces for different choices of target materials and ion beams. We compare the Monte Carlo simulation results with the MSMS for two model cases: 1.5 keV Ar ions incident on Si and 10 keV Xe ions incident on Si. For surfaces curved in the longitudinal direction, the Monte Carlo simulations reproduce the curvature coefficients given by the MSMS quite well, including the predicted destabilizing contribution at larger ion incidence angles θ due to the curvature dependence of the crater function itself. The apparent deflection of ions introduced in the MSMS is described by the positive deflection angle Ψ − θ. We find that the simple linear relation between Ψ and θ proposed in Ref. 6 is a rather good approximation. In contrast, except for low ion energies, the MSMS predictions for surfaces that are curved in the transverse direction are not in satisfactory agreement with the simulation results.

Our motivation for our choice of two model ion/target combinations is as follows. Silicon is a commonly used target material because wafers with extremely flat surfaces are inexpensive and readily available. In several experimental studies, Ar ion irradiation produced patterns on Si for ion energies between 250 eV and 1 keV (Refs. 16 and 17) and for energies above about 20 keV.^18,19 Published data for the energy regime between 1 and 20 keV are sparse, however. Carter et al. reported no patterns on Si for irradiation with 5 keV and 10 keV (Ref. 18) Ar ions. Similarly, Bhattacharjee et al. found no patterns for an 8 keV argon ion beam.²⁰ Recent experimental results obtained by one of the authors (Hofsäss) indicate that patterns are indeed absent between 1.5 and 10 keV.²¹ Because these observations are not currently understood, the case of a 1.5 keV argon ion beam incident on Si is of considerable interest. In contrast to the findings for Ar ions, for Xe ions, irradiation patterns are observed experimentally for many materials and for a broad range of ion energies. In the case of a Si target, experiments using Xe ion beams have produced surface ripples for a variety of energies ranging from 600 eV to 40 keV.^22–25 

SDTrimSP Version 5.06 is a Monte Carlo program that uses the binary collision approximation to simulate the impact of an ion on a planar solid surface and the ensuing collision cascade.¹⁵ We extended the program so that it applies to an ion impact on a curved surface given by the height function

where the ion strikes the surface at the origin and the positive x-direction is the projected direction of the incident ions.

Simulations for a curved surface given by Eq. (1) are implemented in SDTrimSP in the following way: SDTrimSP offers the possibility of setting the ion's starting point at (x, y, z) = (0, 0, −z₀) below the planar surface at z = 0. (Note: SDTrimSP uses a different coordinate system (y,-z,-x)_SDTrimSP= (x,y,z). x is the inward surface normal direction in the program.) Adding the input command line “x0 = 100.0” to the SDTrimSP input script file tri.inp sets the ion impact point to (0, 0, −z₀) inside the target at the depth z₀ = 100 Å. We now define an artificial internal surface

The global coordinate system surface normal is $ẑ$ = (0,0,1). The local surface normal on the surface $h(x,y)$ is

Let the direction of flight of a recoil atom be denoted by the velocity unit vector $r̂=(rx,ry,rz)$⁠. The fraction E_n of the recoil kinetic energy E_recoil that corresponds to the velocity component along the local surface normal $n̂$ is given by $En=Erecoil cos2β$ with $cos β=n̂⋅r̂$⁠. Atoms passing through the internal surface in the positive z-direction are considered to be sputtered if E_n exceeds the surface binding energy E_SBE. The component E_n is calculated for each atom that passes through the internal surface. The criteria for sputtering from the curved surface $h(x,y)$ are therefore $En=Erecoil cos2β≥ESBE and cos β>0$⁠. A noble gas ion reaching the surface is emitted if β > 0.

SDTrimSP considers only one projectile and one target atom at a time. Following a collision, the position of the next target atom to undergo a collision is determined by a random number generator using the mean atomic density of the material. Thus, target atoms are only generated as they are needed. If a recoil atom passes through the internal surface, a new target atom is not generated; instead the procedure to determine whether sputtering or emission occurs is executed. Atoms above the internal surface are therefore not generated and so cannot influence the collision cascade. SDTrimSP does not generate target atoms above the plane z = 0. In the case of a positive surface curvature, employing an internal surface ensures that the program can create target atoms at larger lateral positions (x,y).

Atoms reaching the surface with a low normal velocity component are not sputtered but are instead reflected back into the target. We also implemented a correct treatment of back-reflection from a curved surface in SDTrimSP. In the SDTrimSP program for flat surfaces, the r_z component of the vector $r̂=(rx,ry,rz)$ changes its sign and the back-reflected atom has the velocity unit vector $(rx,ry,−rz)$⁠. For a curved surface, the back-reflection of $r̂=(rx,ry,rz)$ will be denoted by $b̂$ with $|b̂|=1$⁠. At a surface with surface normal $n̂={−K11x ,−K22y , 1}1+(∇h)2$⁠, $b̂$ is determined by the conditions $(r̂+b̂)⋅n̂=0$ and $r̂×n̂=b̂×n̂$⁠. The unit velocity vector $b̂=(bx,by,bz)$ of a back-reflected atom is then given by $b̂=r̂−2(r̂⋅n̂)n̂$⁠. It should be mentioned that back-reflected atoms usually have very low kinetic energies, and so they typically stop inside the target near its surface. Therefore, these atoms do not play a significant role in determining the sputter yield.

The output of the modified SDTrimSP program yields, among many other data, the sputter yield as function of the ion incidence angle and the curvatures K₁₁, K₁₂, and K₂₂. For the present work, we set K₁₂ = 0 and vary K_ii from −0.02 nm⁻¹ to 0.02 nm⁻¹ for either i = 1 or 2. If K₁₁ is varied, K₂₂ is set to zero and vice versa.

The moments $MK11$ and $MK22$ introduced in the Harrison-Bradley crater function formalism⁸ are given by

where Ω is the atomic volume. The value of $dYdKii|Kii=0$ is taken to be average of four to six slopes $Y(Kii≠0)−Y(0)ΔKii$ for the values K_ii = ±0.02 nm⁻¹ and ±0.005 nm⁻¹ for angles up to about 60°. For larger angles, only the two slopes for K₁₁ = ±0.005 nm⁻¹ can be used, because the sputter yield depends non-linearly on K_ii for larger values of K_ii. In order to achieve sufficient accuracy of the slopes $Y(K22≠0)−Y(0)ΔK22$⁠, the sputter yields on curved surfaces with curvature K₂₂ ≠ 0 require simulations with at least 2 × 10⁵ incident ions.

The statistical error for the sputter yield is <0.5% for N_ion = 5 × 10⁴ incident ions. Since $MKii$ depends on the difference of sputter yields for different curvatures, the corresponding statistical error may be large, especially for transverse curvatures K₂₂. The statistical error $σ(MKii)≈ΩY2ΔKiiNion$ can be made as small as desired, however, by choosing N_ion to be sufficiently large. The error in $MKii$ shown in the following figures is calculated for up to n = 6 different values of $ΔKii$⁠. The error of the average of $MKii¯$ due to the individual statistical errors is then $σ(MKii)≈1nΩY2Nion∑i=1n1(ΔKii)2$⁠. Finally, we add the statistical error $s(MKii)≈1n(n−1)∑i=1m(MKii,n−MKii¯)2$ derived from the standard deviation of the slopes for different curvature differences ΔK_ii. The latter takes into account the possibility that the sputter yield could be a nonlinear function of K_ii.

The calculation of the first order moments of the erosion and redistribution crater functions using SDTrimSP has already been described in detail in Ref. 13. The program has now been extended so that higher order moments of crater functions may be computed. This calculation is not done using the standard SDTrimSP output data files which contain data for only a limited number (∼10⁴–10⁵) of sputtered or recoiled atoms (such as “partic_back_r.dat” or “partic_stop_r.dat”), but is instead carried out within the program. All stopped (∼10⁸) and sputtered (∼10⁶) recoils generated in a typical simulation with ∼10⁵ incident ions are taken into account.

The m-th order erosive moments of the crater function for a given ion incidence angle θ are given by

Here, Z_sp = Y(θ)·N_ion is the total number of atoms sputtered during a simulation with N_ion incident ions for a given incidence angle θ. Y is the sputter yield, and $xi,start$ is the x-position where ith sputtered recoil atom is generated. For the comparison with the MSMS, we will only need the first order erosive moment

The redistributive moments of the crater function are given by

Here, N_ion is the number of incident ions, and Z_R = N_R(θ)·N_ion is the total number of recoils which stop in the target rather than being sputtered or transmitted. N_R is the number of recoil atoms created per ion, and $xi,start(θ)$ and $xi,end(θ)$ are the positions where the ith recoil atom is generated and stopped, respectively. In the same way, the program also calculates the higher order moments

and

To compare the MSMS with SDTrimSP Monte Carlo simulations, we have to calculate the parameters a, α, and β of Sigmund's ellipsoidal distribution of deposited energy for a given ion-target combination. Here, a is the average depth of energy deposition measured from the surface along the beam direction, and α and β are the widths of the distribution parallel and perpendicular to the beam direction. To obtain a, α, and β, we apply the procedure used in Ref. 13. For ion incidence angle θ = 0°, we calculate the spatial distribution of recoil end positions within the target using SDTrimSP. We assume that this distribution reflects the spatial distribution of deposited energy, following common practice. From a Gaussian fit to the distribution along the z-direction and the x- or y-direction, respectively, we directly obtain a, α, and β. The distributions of recoil end positions and the corresponding Gaussian fits are shown in Fig. 1 for our two model cases. The quality of the fits is similar for other ion–target combinations and ion energies. For the calculation of curvature coefficients using the MSMS in this work, we will employ the values listed in Table I.

The apparent deflection of ions introduced in the MSMS is described by the deflection angle Ψ. In the simplest version of the MSMS, Ψ is given by Ψ = ρθ, where the constant ρ is not less than 1. In the unmodified Sigmund theory and also in the MSMS, the erosion rate at a surface point is assumed to be proportional to the power deposited there by the incident ions. This proportionality constant is denoted by Λ. We chose the factor ρ so that the sputter yield maximum occurs at the same ion incidence angle in the MSMS as it does in the simulations. The factor Λ is chosen so that the maximum sputter yield obtained from SDTrimSP and from the MSMS are identical for a given target material and ion energy ε. The values of ρ and Λ are listed in Table I.

The fraction of the incident energy f_E that is carried away by backscattered ions is plotted versus the incidence angle in Figures 2 and 3. The targets are amorphous carbon (Figure 2) and silicon (Figure 3) and the projectiles are the noble gas ions He, Ne, Ar, and Xe. The dependence of f_E on the angle of incidence can be described by the relation

The positive constant c₁ depends weakly on the ion energy but strongly on the chosen ion species and target material, as shown by Figures 2 and 3.

The fraction of backscattered energy also depends on the curvature of the surface. This is shown in Figs. 4 and 5 for the case of 1500 eV Ar on Si and 10 keV Xe on Si, respectively. In both figures, K₂₂ is zero and results for several different curvatures K₁₁ are shown. For negative curvature, f_E is significantly larger than for positive curvatures. In the semi-logarithmic plot of f_E, the differences appear as a vertical shift upwards for negative curvature and downwards for positive curvature and can be taken into account by setting

where the constant c₂ is positive. Equation (11) can be used as an approximate expression for f_E as long as K₁₁ is sufficiently small.

For a non-zero curvature K₂₂, we find that f_E is independent of K₂₂ to an excellent approximation. This is shown in Fig. 6 for 1500 eV Ar on Si. The differing dependence of f_E on the curvatures K₁₁ and K₂₂ can be explained by considering the backscattering process itself. Backscattering of ions is a consequence of a single collision or a small number of collisions between ion and target atoms. These collisions usually occur very close to the target surface at locations with small y-values, and so the value of K₂₂ does not significantly affect the backscattering of the ions.

Fig. 7 shows the average sputter yield for an ion that impinges at the origin on a surface with its form given by Eq. (2). The results given are for our two model cases, 1.5 keV Ar ions incident on Si and 10 keV Xe ions incident on Si. The sputter yield depends significantly on the curvature K₁₁ (Fig. 7(a)). Compared to the yield for a flat surface, a negative (convex) curvature K₁₁ generates a higher yield and a positive (concave) curvature a lower yield, but only for angles below a critical angle of about 70°. At higher angles, the situation is reversed. These observations are in agreement with the prediction of the MSMS. The sputter yield's dependence on K₂₂ is shown in Fig. 7(b). To obtain this plot, we simulated the impact of 2 × 10⁵ incident ions per angular step. These simulations reveal that the sputter yield is nearly independent of K₂₂. For all angles of incidence, we find that the sputter yield for a convex surface with K₂₂ < 0 is slightly higher than its value for a concave surface with K₂₂ > 0. This behavior does not agree with the prediction of the MSMS for high angles of incidence.

In order to obtain the weak dependence of the sputter yield on K₂₂ with sufficient statistical significance, one has to perform SDTrimSP simulations with at least 10⁵ ions per angular step. Therefore, it would be quite difficult and very time consuming to extract this curvature dependence from molecular dynamics simulations. In contrast, the SDTrimSP simulations that yielded Fig. 7 only required about an hour of computation time on a standard personal computer.

The goal of the extended crater function formalism introduced by Harrison and Bradley⁸ is to derive an equation of motion of the surface h(x,y,t) of the form

where the coefficients C₀, C₁, C₁₁, and C₂₂ are expressed in terms of moments of the crater function. The curvature coefficients C₁₁ and C₂₂ are obtained by summing the erosive and redistributive contributions. The erosive contributions to these coefficients will be denoted by $C11erosion$ and $C22erosion$⁠. They themselves are the sum of two contributions

and

The erosion coefficients $Siierosion$ take into account the fact that a non-zero surface curvature gives rise to a local angle of incidence that depends on the point of ion impact. The coefficients $Tii$ arise as a result of the curvature dependence of the crater function itself. The curvature coefficients are given by⁶ 

and

where $MKii/Ω$ is the sputter yield for ion impact at the origin on a surface with non-zero curvature K_ii, and the first order erosive moment $Mx,erosion(1)$ can be obtained from SDTrimSP simulations using Eq. (6).

We obtain the moments $Mx(0)$⁠, $Mx,erosion(1)$⁠, and $∂MK11/∂K11|K11=0$ for the MSMS by replacing θ with $Ψ=ρ θ$ in the corresponding expressions for the Sigmund model that are given in Ref. 8. In the following, we compare the sputter yield Y = ΩM_x⁽⁰⁾, the curvature coefficients T₁₁ and T₂₂ (see Eq. (17)), as well as $S11erosion$ (Eq. (15)) and $S22erosion$ (Eq. (16)) given by the MSMS with the sputter yield and the corresponding curvature coefficients obtained from SDTrimSP simulations.

The derivative $(dY/dK11)|K11=0=1Ω(dMK11/dK11)|K11=0$ is plotted in Fig. 8 for our model cases. The derivative is negative up to about 70° and consequently gives rise to a positive (stabilizing) curvature coefficient T₁₁ for θ < 70°. Above 70°, the derivative becomes large and positive and gives rise to a strongly destabilizing negative curvature coefficient T₁₁.

Fig. 9 shows the comparison of the sputter yield and the curvature coefficients T₁₁ and $S11erosion$ for 1.5 keV Ar on Si. The corresponding comparison for 10 keV Xe on Si is shown in Fig. 10. For both of these model cases, we find acceptable quantitative agreement for the sputter yield and the coefficients T₁₁ and $S11erosion$ up to the ion incidence angle where the sputter yield reaches a maximum. For larger angles, the sputter yield decreases and T₁₁ becomes negative in both cases. However, due to significant back-scattering of ions, f_E increases and the sputter yield rapidly drops to zero in the SDTrimSP simulations. This strong decrease of the sputter yield is underestimated by the MSMS. In addition, T₁₁ obtained from SDTrimSP becomes much more strongly negative and leads to a stronger destabilization of the surface than is predicted by the MSMS.

The MSMS treats the backscattering of ions in a simple, approximate fashion. Our results show that this approximation does quite well in the case of 10 keV Xe incident on Si for all but the highest angles of incidence. In contrast, the MSMS seriously underestimates the effect of backscattering for 1.5 keV Ar on Si for angles of incidence beyond the sputter yield maximum, leading to discrepancies in the sputter yield and $S11erosion$⁠. This discrepancy is more pronounced for the lower ion energy because the maximum of the sputter yield appears at a smaller angle. The MSMS also appears to underestimate the longitudinal width α of the distribution of deposited energy for high angles of incidence in the case of 1.5 keV Ar on Si. This leads to the disagreement between the MSMS and SDTrimSP values for T₁₁.

Also shown in Figs. 9 and 10 are the predictions of the SMS (ρ = 1). In this case, the parameter Λ is chosen so that at the angle where the SDTrimSP simulation gives the maximum sputter yield, the SMS and SDTrimSP sputter yields are identical. The values Λ used for the SMS are indicated in Figs. 9 and 10. The SMS gives the highest sputter yield for glancing-incidence bombardment and so does not yield a maximum in Y(θ) for θ < 90° as both SDTrimSP and the MSMS do. Unlike the MSMS, the SMS also fails to reproduce the negative values of T₁₁ that are seen in the simulation results for sufficiently high angles of incidence.

The coefficients $S22erosion$ and $T22$ are plotted in Fig. 11. For 1500 eV Ar on Si, we find satisfactory quantitative agreement between the MSMS and the SDTrimSP simulations. For 10 keV Xe on Si, on the other hand, the coefficients from SDTrimSP are a factor of 2–3 larger than the MSMS predictions, except at large angles. The coefficient T₂₂ from SDTrimSP is always positive for all angles and contributes to a weak stabilization of the surface in the y-direction. In contrast, the MSMS predicts slightly negative T₂₂ values at large angles which would add to the negative $S22erosion$ coefficient and thus would slightly destabilize the surface in the y-direction.

In the case of 10 keV Xe on Si, the MSMS seems to give a distribution of deposited energy that has a lateral width β that is too small. This results in an underestimate of T₂₂ for all but the highest angles of incidence. The large discrepancy between the MSMS and SDTrimSP values for $S22erosion$ stem from the appearance of the factor of cot θ in Eq. (16). Because this factor diverges as θ tends to zero, small deviations in the MSMS values for the crater function moment $Mx,erosion(1)$ from the SDTrimSP values are amplified for intermediate to small angles of incidence.

All in all, for our two model examples, the MSMS describes the angular dependence of the sputter yield and the curvature coefficients T₁₁ and $S11erosion$ reasonably well. However, the MSMS predicts that T₂₂ becomes negative at high angles of incidence, contrary to the results of the simulations. In addition, the MSMS values of T₂₂ and $S22erosion$ for the irradiation of Si with 10 keV Xe ions deviate quite strongly from the SDTrimSP data for all but high angles of incidence. So far as sputtering from a surface with zero transverse curvature K₂₂ is concerned, the simple linear relation between Ψ and θ as well as the assumed proportionality of the sputtering rate and the ion energy both seem to be fairly good approximations. However, the MSMS does rather poorly when applied to sputtering of surfaces with nonzero K₂₂ for higher ion energies.

Both the ion mass and energy differ in our two model examples. To see the effect of independently varying these parameters, we carried out simulations of the sputtering of Si with Ar and Xe ions for a range of ion energies. Figure 12 shows the resulting erosive curvature coefficients $S11erosion$⁠, T₁₁, and $C11erosion$ as well as the MSMS predictions for $C11erosion$⁠. In all cases, T₁₁ becomes negative for sufficiently high angles of incidence, as predicted by the MSMS. For 5 and 10 keV, the agreement between the $C11erosion$ coefficients from the MSMS and SDTrimSP is quite good. There is better agreement for Xe, the more massive of the two ions. For lower ion energies, the agreement between the MSMS and SDTrimSP values of $C11erosion$ is poorer. We attribute this deficiency of the MSMS to an underestimate of the effect of backscattering and the value of α at low energies and higher angles of incidence.

To ensure that our results are not peculiar to Si targets, we carried out SDTrimSP simulations of Xe ion irradiation of amorphous carbon (a-C) for three different ion energies (see Fig. 13). As in the case of Xe bombardment of Si, there is a trend toward improved agreement between the SDTrimSP and the MSMS values of $C11erosion$ as the ion energy increases. The curvature coefficient T₁₁ is once again negative for sufficiently high values of θ, as predicted by the MSMS.

We have carried out Monte Carlo simulations of ion erosion using the program SDTrimSP to determine the erosive moments of the crater function. We not only took into account the curvature dependence of the sputter yield but also the curvature dependence of the erosion crater function itself. Furthermore, we used Monte Carlo simulations to determine the parameters a, α, and β of Sigmund's ellipsoid for given ion-target combinations. With these data, we were able to test the MSMS with Monte Carlo simulations.

An analysis of the fraction of energy carried away by back-scattered ions justified the assumption made in the MSMS of an apparent defection of the ions towards the surface. With this assumption, the MSMS fits the angular dependence of the sputter yield obtained from the Monte Carlo simulations reasonably well, provided that the ion energy is sufficiently high (typically above 1 keV). For the model cases we studied, we found that the simple relation Ψ = ρθ, as introduced in the MSMS, and a linear scaling between the sputter yield and the incident ion energy is sufficient to obtain good quantitative agreement with SDTrimSP simulations for the angular dependence of the sputter yield. Only for very low ion energies is the agreement poor, but it can be improved empirically by adjusting the parameters a, α, β, and Λ.

The sputter yield for ions incident at the origin of a surface with a longitudinal curvature K₁₁ shows the angular dependence predicted by the MSMS for all ion species, ion energies, and target materials we studied. Below a certain critical incidence angle, the sputter yield for a convex surface (K₁₁ < 0) is larger than for a surface with a concave curvature (K₁₁ > 0). Above this critical angle, the situation is reversed. Therefore, the curvature coefficient $C11erosion$ obtained from our simulations includes a destabilizing contribution T₁₁ < 0 at sufficiently large angles that comes from the curvature dependence of the erosive crater function, in agreement with the MSMS. The SDTrimSP simulations indicate that the MSMS actually underestimates this destabilizing contribution. At very low ion energies, SDTrimSP shows almost no K₁₁ curvature dependence of the sputter yield up to about 60°. Therefore, the T₁₁ coefficient obtained from SDTrimSP is quite small under these conditions. It is, however, quite substantial for angles close to grazing incidence.

According to our SDTrimSP simulations, the curvature coefficient C₂₂ has always a stabilizing contribution T₂₂ > 0, in contrast to the prediction of the MSMS for high angles of incidence. However, both the MSMS and Monte Carlo values for T₂₂ are small when θ is close to 90°, and so the difference between them is small as well.

We can gain insight into the differing behavior of T₁₁ and T₂₂ by noting that the sputtering of a surface by an obliquely incident ion beam is influenced by two key factors:

1. Recoil collisions occur closer to the local surface for a convex surface than for a concave surface.

2. For a surface with zero transverse curvature K₂₂, the fraction of ions that is backscattered is higher for a convex surface than for a concave surface (see Figs. 5 and 6). In contrast, for a surface with zero longitudinal curvature K₁₁, the fraction of ions that is backscattered is almost independent of K₂₂ (see Fig. 7).

The first of these two factors is predominant for a surface with K₁₁ = 0. As a consequence, for such a surface, the sputter yield is higher for a convex surface with K₂₂ < 0 than for a concave surface with K₂₂ > 0. Now consider the sputtering of a surface with K₂₂ = 0. Very few ions are backscattered for moderate angles of incidence. Therefore, in this regime, Factor 1 is most important and thus the sputter yield of a convex surface is higher than the yield of a concave surface. On the other hand, for θ close to 90°, the fraction of backscattered ions grows large and Factor 2 plays a crucial role. In this range of angles, the sputter yield of a concave surface is higher than that of a convex surface.

In this paper, we have only discussed the erosive contributions to the curvature coefficients C₁₁ and C₂₂. In order to compare with experiment, the contributions to C₁₁ and C₂₂ that come from mass redistribution must be included. These contributions can also be determined by the crater function formalism.^3,7,9–13

The experimentally relevant coefficients C_jj are the superpositions $Cjjerosion+Cjjredistribution$⁠. Recently, one of us (Hofsäss) introduced an additional stabilizing contribution to $C11redistribution$⁠, which takes into account the angular dependence of the thickness of the irradiated film and the dynamic adjustment between the interface of the unirradiated solid and the surface.^13,26 SDTrimSP simulations for a variety of different ion target combinations and ion energies indicate that typically both erosive and redistributive effects contribute significantly to pattern formation. A comparison of experimentally observed critical angles for ripple formation with predictions from Monte Carlo simulations, taking into account the erosive contributions discussed in this paper as well as all aspects of the mass redistribution, will be presented in a forthcoming publication.

It is natural to ask whether improvements can be made to the MSMS that would lead to better agreement with the simulation results. In fact, in this paper, we only tested the simplest version of the MSMS in which Ψ is taken to be proportional to θ. There is a more general version of the MSMS⁶ in which it is simply assumed that Ψ is a function of θ, i.e., Ψ = f(θ). It is likely that another choice of the function f would result in better agreement between the MSMS and the simulations. For example, consider a function f(θ) = Aθ^p, where A and p are constants. Preliminary work suggests that a choice of p that is greater than 1 yields improved agreement. In addition, in the MSMS, it is assumed that the form of Sigmund's ellipsoidal distribution of deposited energy does not depend on the angle of incidence θ, just as in the SMS. Our simulation results suggest that the parameters α and β that describe this distribution may actually change as θ is varied, however.

Monte Carlo programs that permit the user to simulate sputtering from curved surfaces are not readily available at the present time. Moreover, to obtain the erosive curvature coefficients using our modified version of SDTrimSP required calculations for at least 18 angular steps of 5°, up to 5 different curvatures and at least 10⁵ ions per angle and curvature. Although the computing time of a few hours was short, the data from 90 individual simulations had to be combined for data reduction. This was achieved using an extended version of the program which is not yet generally available. In contrast, it is only necessary to simulate ion impact on flat surfaces in order to obtain the small number (a, α, β, ρ, Λ) of ion-target specific input parameters that appear in the MSMS. These parameters can be easily obtained from commonly used Monte Carlo simulation programs such as SRIM²⁷ for almost any ion-target combination and for a very broad range of ion energies. Experimentally obtained sputter yield data can also be used to obtain the parameters ρ and Λ. Once values for the parameters a, α, β, ρ, and Λ have been acquired for a particular choice of ion species, ion energy, and target material, the MSMS gives reasonably accurate values of the erosive curvature coefficients $S11erosion$⁠, T₁₁, and $C11erosion$ for all angles of incidence.

We have benefitted from discussions with M. P. Harrison, D. A. Pearson, and P. D. Shipman. R.M.B. would like to thank the National Science Foundation for its support through Grant No. DMR-1305449. H.H. would like to thank the German Research Foundation (DFG) for its support through Grant No. HO 1125/20-2.