Sputter yields of surfaces with nanoscale textures: Analytical results and Monte Carlo simulations

When an energetic ion strikes a solid surface, one or more atoms can be ejected from the surface in a process called sputtering. The average number of atoms that are ejected is referred to as the sputter yield and is denoted by $Y$⁠. The sputter yield depends on the local angle of incidence, i.e., the angle that the incident ion’s trajectory makes with the local surface normal. It also depends on the surface curvature at the point of impact $P$⁠.^1,2 In fact, $Y$ depends on all of the spatial derivatives of the surface height at $P$⁠, at least in principle.³ 

In many applications, solid surfaces are sputtered with a broad ion beam. Examples include depth profiling using secondary ion mass spectrometry (SIMS), ion induced Auger electron spectroscopy (IIAES), ion milling, and sputter deposition of thin films. The sputtered surfaces are virtually never planar at the nanoscale, and so there has been longstanding interest in the spatially averaged sputter yield $Y¯$ that results when a broad ion beam is incident on a surface with a nontrivial morphology. This topic is also relevant to the erosion of the vessel walls in nuclear fusion reactors and to the weathering of surfaces in space by the solar wind.

Early analytical studies of the spatially averaged sputter yield of non-planar surfaces were based on Sigmund’s theory of spatially resolved sputtering,^4–7 a theory that is known to be a poor approximation in many circumstances.⁸ Shadowing of the incident ion flux, redeposition of sputtered material, and secondary sputtering were all neglected in these studies. Simulations have been more frequently used to compute the $Y¯$ for non-planar surfaces,^9–19 and attempts have been made to relate the resulting estimates of $Y¯$ to various aspects of the surface texture, including the average local angle of incidence and the root-mean-square surface roughness. In some simulations¹⁹ and in recent analytical work,²⁰ it was assumed that when an ion strikes the surface, the sputter yield depends only on the local angle of incidence. In theories of this type, the curvature dependence of the sputter yield is taken to be negligible. This dependence can, however, have an important effect if there are surface features with sizes approaching the size of a collision cascade.^1,2,21,22

In this paper, we find the spatially averaged sputter yield $Y¯$ analytically for non-planar surfaces that have slowly varying surface heights $h=h(x,y)$⁠. Nonlocal effects are initially neglected. Our theory is quite general and is not based on Sigmund’s theory. Our principal finding is that the leading order corrections to the average sputter yield are proportional to the spatial averages of $(∂h/∂x)2$ and $(∂h/∂y)2$⁠. The constants of proportionality are related to the sputter yield of a flat surface, which means that curvature-dependent sputtering can be neglected when computing the leading order correction to the average sputter yield. We also apply our result to surfaces with sinusoidally varying surface height and to statistically rough surfaces. Finally, we show that the contribution of redeposition to $Y¯$ is proportional to $(A/l)2$⁠, where $A$ and $l$ are the amplitude and characteristic lateral length scale of the surface texture, respectively.

This paper is organized as follows: In Sec. II, we derive an approximate expression for the spatial average of the sputter yield $Y¯$ neglecting nonlocal effects. We also determine the scaling behavior of the contribution of redeposition to $Y¯$ for normal-incidence bombardment. In Sec. III, we carry out Monte Carlo simulations of the sputtering of surfaces with sinusoidally varying surface height. In addition to testing the predictions of our theory on primary sputtering and redeposition, we find the contributions of recoil sputtering and secondary ion sputtering. We discuss our results in Sec. IV and conclude in Sec. V.

Consider the surface of a solid elemental material that is bombarded with a broad beam of noble gas ions. Our goal is to develop an approximate expression for the spatial average of the sputter yield $Y¯$ at time $t=0$ that is valid for a surface that has a small slope everywhere. For simplicity, we will take the target material to be amorphous and neglect any influence of the implanted ions on the sputter yield.

We will adopt a continuum description of the solid surface. Let $h=h(x,y,t)$ denote the height of the surface above the point $(x,y)$ in the $x−y$ plane at time $t$⁠. We assume that $h$ is a single valued, smooth function of $x$ and $y$ at time $t=0$⁠. The partial derivative of $h$ with respect to $x$ will be denoted by $hx$⁠, and $hy$ and $ht$ are defined analogously. We let the incident ion current above the solid be $J=−Je^$⁠, where $J$ is a positive constant and $e^$ is a unit vector. We orient the $x$ and $y$ axes so that $e^=z^cos⁡θ+x^sin⁡θ$⁠, where $θ$ is the angle of ion incidence.

Let $a$ denote the characteristic length scale of a collision cascade. In the Sigmund model of ion sputtering, for example, $a$ can be taken to be the mean depth of energy deposition. We assume that $h$ varies slowly with $x$ and $y$ at time $t=0$⁠. Stated more precisely, we assume that $h(x,y,0)$ only varies appreciably over a characteristic lateral length scale $l$ that is much larger than $a$⁠. Additionally, we take the amplitude $A$ of the surface disturbance $h(x,y,0)$ to be of order $a$ for all $x$ and $y$⁠. Thus,

where $ϵ≡a/l$ and $H$ is dimensionless and of order 1. We will study the limit in which $ϵ$ is small.

$ht(x,y,0)$ is the rate that the surface recedes above the point $(x,y)$ in the $x−y$ plane at time $t=0$⁠. Clearly, $ht(x,y,0)$ depends on the angle of incidence $θ$⁠. It also depends in principle on the shape of the entire surface at time $t=0$⁠, or, equivalently, on all of the spatial derivatives of $h(x,y,0)$⁠. We will write

All of the partial derivatives of $h$ that appear on the right-hand side of Eq. (2) are to be evaluated at $x$⁠, $y$ and time $t=0$⁠.

There are nonlocal effects that can affect the value of $ht(x,y,0)$⁠, including shadowing of the incident ion flux, redeposition of sputtered material, and secondary sputtering. In the limit that $ϵ$ tends to zero while $θ<90°$ is held fixed, the surface slope $hx$ is of order $a/l=ϵ$⁠, which is much smaller than the slope of the trajectory made by an incident ion, $cot⁡θ$⁠. Shadowing, therefore, cannot occur. Secondary sputtering occurs when an ion is reflected from the solid surface, strikes the surface again, and then sputters material at the new point of impact. It also occurs when a previously sputtered atom strikes the surface and then causes sputtering itself. Both secondary sputtering and redeposition of sputtered atoms can contribute to the value of $ht(x,y,0)$⁠. In principle, double reflection of ions and surface reflection of sputtered atoms could also occur, but these are negligible higher-order effects.

For the time being, we will neglect secondary sputtering and redeposition, but later we will discuss what effect this has on our final results. We will also make the customary assumptions that the effect of ion implantation is negligible and that the density of the solid is spatially uniform and independent of time.

Inserting Eq. (1) in Eq. (2), we obtain

where $X≡ϵx$ and $Y≡ϵy$⁠. In the language of the method of multiple scales, $X$ and $Y$ are the “slow” variables, whereas $x$ and $y$ are the corresponding “fast” variables.

Because we are neglecting redeposition and secondary sputtering, $f$ is an analytic function of its arguments. We will expand Eq. (3) out to second order in $ϵ$⁠. So that the result of this expansion can be written down succinctly, we define

and so forth. Similarly, for positive integers $i$ and $j$⁠, $fi,j(θ)$ will denote the partial derivative of $f(θ;hx,hy,hxx,hxy,…)$ with respect to the $i$th and $j$th arguments that appear after the semicolon, evaluated for all the arguments after the semicolon set equal to zero. For example,

We obtain

In writing Eq. (9), we have omitted all of the terms that violate the requirement that $ht$ be invariant under the transformation $y→−y$ since their coefficients must be zero by symmetry. When recast in terms of the original variables, Eq. (9) becomes

where $Λ1≡12f1,1$ and $Λ2≡12f2,2$⁠. All of the partial derivatives of $h$ that appear in Eq. (10) have the arguments $x$⁠, $y$ and $t=0$⁠.

Equation (10) is valid to order $ϵ2$⁠. Clearly, $f0$ is the rate that a flat surface recedes. The terms $f1hx$⁠, $Λ1hx2$⁠, and $Λ2hy2$ come from the slope dependence of the sputter yield. Finally, the terms $f3hxx$ and $f5hyy$ include contributions from curvature-dependent sputtering and mass redistribution. Equation (10) is just the usual anisotropic Kuramoto–Sivashinsky equation that describes the surface dynamics,^3,23,24 except that the term proportional to $∇2∇2h$ does not appear because it is of order $ϵ4$⁠.

Suppose that the surface occupies the domain with $0≤x≤L$ and $0≤y≤L$⁠. We require $h$ to satisfy periodic boundary conditions. At time $t=0$⁠, the spatial average of the sputter yield $Y¯$ is given by

where $Ω$ is the atomic volume. We now insert Eq. (10) into the right-hand side of Eq. (11) and make use of the periodic boundary conditions. We obtain

where $⟨hx2⟩$ denotes the spatial average of $hx2$ at time $t=0$ and $⟨hy2⟩$ is defined analogously. Explicitly,

where $i$ is $x$ or $y$⁠.

The term $f1hx$ in Eq. (10) gives no contribution to $Y¯$⁠. This is to be expected, since this term simply makes a disturbance propagate over the surface without changing its form.²⁵ What is perhaps more surprising is that the terms $f3hxx$ and $f5hyy$ that stem in part from curvature-dependent sputtering also give no contribution to $Y¯$ since their effect averages to zero. Intuitively, this is because the increased sputtering from the surface depressions is offset by reduced sputtering from protrusions.

Let $Y(θ)$ denote the sputter yield if $h(x,y,0)=0$ for all $x$ and $y$⁠, i.e., $Y(θ)$ is the sputter yield if the surface is perfectly flat. $f0$⁠, $Λ1$⁠, and $Λ2$ may be written in terms of $Y(θ)$ and its derivatives.^26,27 The results are

and

where the primes denote derivatives with respect to $θ$⁠. Equation (12) now becomes

where

and

Note that $μ1(0)=μ2(0)=Y′′(0)/2$⁠, and so $μ1$ and $μ2$ coincide for normal-incidence bombardment, as they must by symmetry.

On the right-hand side of Eq. (17), $Y(θ)$ is the sputter yield of a flat surface. The second and third terms on the right-hand side of Eq. (17) are the order $ϵ2$ corrections to this zeroth order result. To this order, the spatially averaged sputter yield $Y¯$ depends on two descriptors of the surface roughness: $⟨hx2⟩$ and $⟨hy2⟩$⁠. As we shall see, typically $μ1(θ)$ is positive for $θ$ near zero and for sufficiently high angles of incidence but is negative for a range of intermediate values of $θ$⁠. $μ2(θ)$⁠, on the other hand, is positive for $θ$ less than a critical value of $θ$ and is negative for larger values of $θ$⁠. If $μ1$ is positive, we have the behavior that we would intuitively expect: if $⟨hx2⟩$ is increased, then $Y¯$ increases. Conversely, if $μ1$ is negative, we come to the surprising conclusion that $Y¯$ decreases if $⟨hx2⟩$ increases. Similar conclusions can, of course, be drawn about the dependence of $Y¯$ on $⟨hy2⟩$⁠.

It is instructive to rewrite our main result [Eq. (17)] in Fourier space. We have

where $x≡xx^+yy^$⁠. The sum in Eq. (20) runs over the wave vectors $k$ consistent with the periodic boundary conditions, i.e., $k=2π(nxx^+nyy^)/L$⁠, where $nx$ and $ny$ are integers. The Fourier coefficient $hk$ is given by

We find that

Equation (22) shows that the contributions of the different Fourier modes simply add to give $Y¯$ to second order in $ϵ$⁠.

Important special cases of Eq. (22) are obtained if the surface height is simply a sinusoidal function of $x$ or $y$⁠. The characteristic lateral length scale $l$ is then the wavelength $λ$⁠. If

then

where $k≡2π/λ$ is the wave number. Similarly, if

then

For Eqs. (24) and (26) to be valid, $A$ must be of order $a$ and $λ$ must be much larger than $a$⁠, so that the correction terms are small. In addition, $A/λ$ must be smaller than $cot⁡θ/(2π)$⁠, because otherwise shadowing will occur and Eq. (24) will, therefore, not apply. In contrast, shadowing cannot occur for a surface described by Eq. (26), no matter how large $A/λ$ is.

As we saw in the Introduction, there has been a great deal of interest in the sputter yield of rough surfaces.^5,6,9–20 Our result Eq. (17) can also be applied to surfaces of this kind. Suppose that $L=NL0$⁠, where $N$ and $L0$ are both large. We divide the $L×L$ domain in the $x−y$ plane into $N2$ square cells, each having side length $L0$⁠. Clearly,

where $⟨hx2⟩n$ denotes the spatial average of $hx2$ over the $n$th cell. Because $N2$ is large, the right-hand side of Eq. (27) is the statistical average of $⟨hx2⟩$⁠. Let $⟨Q⟩stat$ denote the statistical average of the random variable $Q$⁠. The height–height correlation function at time $t=0$ is given by

We assume that

where $g$ is a scaling function, $A$ is a measure of the amplitude of the surface roughness, and $ξ∥$ (⁠$ξ⊥$⁠) is the longitudinal (transverse) correlation length. Equation (29) holds for the surface of a thin film that has kinetically roughened during deposition, for example.²⁸ We differentiate $C(x−x′)$ with respect to $x$ and $x′$ and then take the limit as $x′→x$⁠. This gives

where $gii$ is the second derivative of $g$ with respect to its $i$th argument. We conclude that

In exactly the same way,

Equation (17) now yields our final result,

For a rough surface that is statistically isotropic, $ξ∥=ξ⊥≡ξ$⁠, $g11(0,0)=g22(0,0)$ and

Thus, the deviation of the spatially averaged sputter yield $Y¯$ from the sputter yield for a flat surface $Y(θ)$ is proportional to $A2$ and inversely proportional to the square of the correlation length $ξ$⁠.

It is possible to make analytical progress on the contribution that redeposition makes to the average sputter yield $Y¯$ for the case of normal-incidence ion bombardment. Our starting point for this is the work on redeposition that appears in Ref. 29.

Let $Q(x)$ be the contribution of redeposition to $ht(x,0)$⁠. Expanding $Q(x)$ in powers of $ϵ$ and retaining only the leading order term, one finds that²⁹ 

Here, $h≡h(x,0)$⁠, $h′≡h(x′,0)$⁠, $d2x′≡dx′dy′$⁠, and $γ$ is a positive constant that characterizes the angular distribution of sputtered material (see Ref. 29 for the definition of $γ$⁠). In addition, $s(x+hz^,x′+h′z^)$ is equal to 1 if there is a line of sight between the surface points $x+hz^$ and $x′+h′z^$ and is otherwise equal to zero.

Let $Y¯redep$ denote the contribution of redeposition to $Y¯$⁠. We have

in analogy with Eq. (11). Making use of Eqs. (1) and (35), we obtain

where

Here, $X≡Xx^+Yy^$⁠, $X′≡X′x^+Y′y^$⁠, $d2X≡dXdY$⁠, $d2X′≡dX′dY′$⁠, $∇X≡x^∂X+y^∂Y$⁠, and $∇X′≡x^∂X′+y^∂Y′$⁠. Additionally, the limits on the integrals over $X$⁠, $Y$⁠, $X′$⁠, and $Y′$ in Eq. (38) are all zero and $aL/l$⁠. $F$ is, therefore, a function of $aL/l$⁠. In the limit that $L$ becomes large, $Y¯redep$ must tend to a limit that is independent of $L$⁠. It follows that for large $L$⁠, $F=F(aL/l)$ must be approximately equal to $α(aL/l)2$⁠, where the constant $α$ does not depend on $a$⁠, $l$⁠, or $L$⁠, but it does depend on the detailed form of the function $H=H(X,Y)$⁠. We conclude that

We note that $Y¯redep$ is of order $ϵ2$ and so it is of the same order as the corrections to $Y¯$ coming from direct sputtering that are given by Eq. (17). $Y¯redep$ must also be negative, since redeposition reduces the net amount of material sputtered from the solid surface.

An alert reader may wonder why we did not include the effect of redeposition in the function $f$ that appears on the right-hand side of Eq. (2). The reason for this is that redeposition makes a non-analytic contribution to $f$⁠, and so the expansion that yielded Eq. (9) would not have been possible had we done so.

Let us see why the contribution of redeposition to $f$ is, in general, non-analytic. The contribution of redeposition to $ht(x,y,0)$⁠, or, equivalently, to $f$ will be denoted by $fredep$⁠. $fredep$ depends on the angle of incidence $θ$⁠. It also depends in principle on the shape of the entire surface at time $t=0$⁠, or, equivalently, on all of the spatial derivatives of $h(x,y,0)$⁠. Thus,

All of the partial derivatives of $h$ that appear on the right-hand side of Eq. (40) are to be evaluated at $x$⁠, $y$ and time $t=0$⁠. Consider the value of $fredep$ for $x=y=0$⁠. For the sake of specificity, we consider the normal-incidence bombardment of a surface with $h(x,y,0)=12Kx2$⁠. In this case, $fredep=fredep(0;0,0,K,0,0,…)$ for $x=y=0$⁠. Sputtered material is deposited at the origin for any $K>0$⁠. In contrast, no material is deposited at the origin for $K<0$ because there is no line of sight between any surface point and the origin. Thus, $fredep(0;0,0,K,0,0,…)$ is positive for any $K>0$ and is equal to zero for any $K<0$⁠. Such a function is not an analytic function of $K=hxx(0,0,0)$⁠, and so the assertion has been established.

To test the analytical results presented in Sec. II, we performed Monte Carlo simulations of targets with sinusoidal surfaces, computed the average sputter yield $Y¯$⁠, and then compared the results with the predictions of Eqs. (24) and (26). The coefficients $μ1$ and $μ2$ given by Eqs. (18) and (19) were obtained by Monte Carlo simulations of the angular dependent sputter yields $Y(θ)$ of flat targets. The numerical procedure for obtaining the derivatives $Y′$ and $Y′′$ required for evaluating Eqs. (18) and (19) is described in the  Appendix. Since the calculation of the derivatives requires accurate values of the sputter yield, $108$ ion impacts were simulated for each condition.

For the Monte Carlo simulations, we used the IMplant and Sputter sImuLator (IMSIL).^8,30–33 IMSIL is capable of simulating 1D and 2D targets, the latter with arbitrary geometries defined by polygons. Simulations of irradiation of amorphous Si targets with 1 keV Ar ions were performed, resulting in a mean energy deposition depth of $a=28.25$ Å. We note that for Si bombarded by 1 keV Ar ions at an angle of incidence of $60°$⁠, the maximum concentration of Ar has been shown to be below 4%,³⁴ and so the influence of implanted Ar should be negligible. The 2D targets were defined by one period of a cosine function in the $x$ direction as described by Eq. (23) with wavelength $λ=100a$⁠, varying amplitudes and periodic boundary conditions. A bottom surface, which is required by IMSIL, was placed at $h=−10a$⁠, which was deep enough that forward sputtering was negligible. The cosine function was approximated by a polygon with 1000 segments. Internally, the geometry is represented by a distance function defined on an auxiliary grid³⁵ whose spacing was chosen to be $0.1a$⁠.

The projection of the incidence direction onto the $x$–$y$ plane was chosen to be either parallel or perpendicular to the $x$ axis. We call these parallel and perpendicular modes, respectively. The latter case is equivalent to bombarding a sample with surface height given by Eq. (25) with the projected ion direction along the $x$ axis. The models for nuclear and electronic stopping (the ZBL interatomic potential³⁶ and the mixed Lindhard³⁷ and Oen–Robinson³⁸ model, respectively) as well as the values of cutoff and surface binding energies were the same as in our recent publication Ref. 33. Details of the selection of collision partners and the treatment of the surface are described in Refs. 8 and 32.

The results for the sputter yield $Y$ as a function of incidence angle $θ$ are shown in Fig. 1. In this and all of the figures which follow in which error bars are not shown, the statistical errors are smaller than the symbol sizes. The functional form of $Y(θ)$ is that typical of amorphous targets since $Y′(0)=0$ and there is a single maximum. The behavior near $θ=90°$ is uncertain since it depends on the surface roughness,¹⁷ which cannot be predicted by Monte Carlo simulations. In our simulations, we neglected surface roughness, which leads to $Y′′>0$ near $θ=90°$⁠, while with sufficient roughness $Y′′≈0$⁠.¹⁷ Because of this uncertainty, we will not discuss incidence angles larger than $80°$⁠.

The results for $μ1$ and $μ2$ are shown as a function of the incidence angle $θ$ in Fig. 2. Both coefficients are positive for moderate incidence angles, which means that the average sputter yield of a sinusoidal surface increases with $A/λ$⁠. For the parallel mode, $μ1$ becomes negative and large in magnitude at larger angles, i.e., the sputter yield decreases rapidly with $A/λ$⁠. For the perpendicular mode, on the other hand, $μ2$ is negative but small for higher incidence angles.

According to Eqs. (24) and (26), the change in the average sputter yield $Y¯$ of a sinusoidal surface is quadratic in $A/λ$⁠. We therefore plotted $Y¯$ as a function of $(A/λ)2$ (see Fig. 3). Since redeposition and secondary sputtering are neglected in the theory, we terminated ion and recoil trajectories when they left the target in the Monte Carlo simulations. Excellent agreement between theory and simulation is observed except for the incidence angle of $75°$ in parallel mode. Closer inspection of this case shows that the quadratic approximation given by Eq. (24) only holds for very small $A/λ$⁠, while a quartic term can be identified for $A/λ≳0.01$⁠.

In Fig. 4, the influence of various sputtering mechanisms on the average sputter yield is illustrated for the parallel mode. The blue bullets and lines represent the same data as in Fig. 3, except that additional data points for larger $A/λ$ values are shown. For $θ=0°$⁠, $15°$⁠, $30°$⁠, and $60°$⁠, the excellent agreement between the theory and the simulations extends to $A/λ=0.07$⁠, above which higher-order terms apparently become increasingly important. In the plots for $θ=60°$ and $75°$⁠, the shadowing limits are indicated by dashed vertical lines. Above this limit, the deviation of the theoretical prediction from the simulation results increases significantly. This is to be expected since shadowing of the incident ion flux is not taken into account in the theory.

The orange crosses and lines show simulation results in which recoils that re-entered the target were not considered to have been sputtered, and their trajectories were not followed after re-entering the sample. The difference between orange and blue lines, therefore, represents the contribution of redeposition. As expected, the contribution of redeposition to $Y¯$ is negative and increases in magnitude with $A/λ$⁠. It is largest for moderate ion incidence angles $θ$⁠.

For the data shown with green squares and lines, recoils were followed after re-entering the target, so that they could cause additional sputtering events. Thus, the difference between green and orange lines represents the effect of sputtering by recoils that strike the solid after being ejected from it. This effect is not negligible for moderate incidence angles. In fact, it counteracts roughly 20%–30% of the effect of redeposition.

Finally, for the data shown with red diamonds and lines, reflected ions were also allowed to re-enter the target so that they too could cause secondary sputtering. Thus, the difference between red and green lines represents the effect of secondary sputtering by ions, which becomes more significant for larger incidence angles because the ion reflection coefficient increases with the incidence angle.

In Fig. 5, data analogous to Fig. 4 are shown for the perpendicular mode. As expected from the coefficients $μ1$ and $μ2$ (Fig. 2), the dependence of the average sputter yield on the incidence angle is weaker for the perpendicular mode than for the parallel mode. The excellent agreement between theory and simulation without ion and recoil re-entry (blue lines) extends to $A/λ=0.07$ for all incidence angles shown. Moreover, no shadowing occurs, and so the $A/λ$ dependence is smoother for large incidence angles. For large $θ$⁠, the contribution of redeposition (the difference between orange and blue lines) becomes very pronounced, although it is again partially offset by the contribution of recoil sputtering and secondary ion sputtering: see, in particular, the plot for $θ=75°$⁠.

In Sec. II B, it was shown that for normal-incidence bombardment, the contribution of redeposition to the average sputter yield is proportional to $(A/λ)2$ for small $A/λ$⁠. Comparing the dashed blue line with error bars and the upper black line in Fig. 6, we see that this prediction falls within the statistical error of the Monte Carlo simulations. The other data (which are shown in different colors) indicate that the quadratic dependence on $A/λ$ also holds approximately for oblique ion incidence.

We will begin by discussing our findings for the case in which nonlocal effects are neglected. Later in this section, we will touch on the influence nonlocal effects have on the spatially averaged sputter yield $Y¯$⁠.

As we have seen, the spatial average of the sputter yield $Y¯$ is an increasing function of $⟨hx2⟩$ if $μ1(θ)$ is positive. Conversely, it is a decreasing function of $⟨hx2⟩$ if $μ1(θ)$ is negative. Similar statements apply to the dependence of $Y¯$ on $⟨hy2⟩$⁠, except that in that case, the role of $μ1$ is played by $μ2$⁠.

Figure 2 shows that for the bombardment of an amorphous silicon target with 1 keV Ar$+$ ions, $μ1(θ)$ is positive for $θ$ near zero and that it is negative for a range of intermediate values of $θ$⁠. This figure also shows that $μ2(θ)$ is positive for $θ$ less than a critical value of $θ$ and is negative for larger values of $θ$⁠.

We will show that these conclusions are true for any choice of amorphous target material and noble gas ion species. Only a single, natural assumption is needed: we assume that the sputter yield for a flat surface $Y(θ)$ has a single local maximum. To our knowledge, this assumption is valid for all amorphous target materials and noble gas ions. (For crystalline target materials—which are not considered here—$Y(θ)$ can have multiple local minima along channeling directions.^39,40 Local maxima appear between these local minima.)

Suppose that $Y(θ)$ attains its maximum value at $θ=θc$⁠. Thus, $Y′(θ)$ is positive for $θ<θc$ and is negative for $θ>θc$⁠. It follows from Eq. (19) that $μ2(θ)>0$ for $θ<θc$ and that $μ2(θ)<0$ for $θ>θc$⁠.

Next, we consider the dependence of $μ1$ on $θ$⁠. The sputter yield for a flat surface $Y(θ)$ has a local minimum at $θ=0$ and so $Y′(0)=0$ and $Y′′(0)>0$⁠. It follows from Eq. (18) that $μ1(0)$ is positive. We can also determine the sign of $μ1(θ)$ for $θ=θc$⁠. Because $Y(θ)$ attains a maximum at $θ=θc$⁠, we have $Y′(θc)=0$ and $Y′′(θc)<0$⁠. Equation (18) then shows that $μ1(θc)=Y′′(θc)/2<0$⁠. We conclude that $μ1(θ)$ is positive for $θ$ less than a critical value $θc′$ that lies between zero and $θc$ and that $μ1(θ)$ is negative for a range of $θ$ values that includes the angle $θ=θc$⁠.

We can also gain some physical insight into these conclusions. For simplicity, we consider the cases in which the surface height is given by Eq. (23) or (25), i.e., it is a sinusoidal function of either $x$ or $y$⁠.

$Y¯$ is an increasing function of $A/λ$ for normal-incidence bombardment, i.e., for $θ=0$⁠. The reason for this is that the sputter yield for a flat surface $Y(θ)$ is an increasing function of $θ$ for $θ<θc$⁠. As $A/λ$ is increased from zero, the local angle of incidence $ϕ$ increases everywhere on the surface except at the maxima and minima of the sinusoid. As a consequence, $Y¯$ increases and $μ1(0)=μ2(0)>0$⁠.

To see why $μ1(θ)$ is negative for intermediate values of $θ$⁠, consider bombarding the surface given by Eq. (23) at the angle $θc$ where $Y(θ)$ attain its maximum value. If we start to increase $A/λ$ from zero, the local angle of incidence $ϕ$ will change at all surface points except those where the slope $hx$ is zero. Therefore, the sputter yield will decrease at all surface points except at the maxima and minima of the sinusoid. As a consequence, $Y¯$ decreases and $μ1(θc)<0$⁠.

Now, consider the situation when the surface given by Eq. (25) is bombarded by an ion beam with $θ>θc$⁠. In this regime, $Y(θ)$ is a decreasing function of $θ$⁠. If $A/λ$ is small, then the local angle of incidence $ϕ≅θ+12hy2cot⁡θ$⁠.²⁷ Thus, $ϕ$ is an increasing function of $hy2$⁠. It follows that $ϕ$ is an increasing function of $A/λ$ at all points on the surface except at the maxima and minima of the sinusoid. We conclude that $Y¯$ decreases with $A/λ$⁠, in accord with the result $μ2(θ)<0$ for $θ>θc$⁠.

Makeev and Barabási assumed that the Sigmund model is a good approximation in their work on the sputtering of sinusoidal surfaces.^4,7 All of our analytical results can be immediately applied if the Sigmund model is assumed to be valid since $Y(θ)$ is known for this model.² For the sinusoidal surface given by Eq. (23), for example, Eq. (24) gives $Y¯$ once $μ1(θ)$ has been computed for the sputter yield $Y(θ)$ for the Sigmund model. In particular, the sputter yield correction $Y¯−Y(θ)$ is proportional to $(A/λ)2$ for small $A/λ$⁠, as Makeev and Barabási concluded. We emphasize, though, that we have demonstrated that this result holds even if the Sigmund theory is a poor approximation.

Equation (34) applies for a rough surface that is statistically isotropic. As we have seen, $μ1(0)=μ2(0)>0$⁠, $μ1(θc)<0$ and $μ2(θc)=0$⁠. As a consequence, $μ1(0)+μ2(0)>0$ and $μ1(θc)+μ2(θc)<0$⁠. This means that $Y¯$ decreases as the correlation length $ξ$ grows for normal and near-normal-incidence bombardment. In contrast, $Y¯$ is an increasing function of $ξ$ for $θ=θc$ and for a range of angles about $θc$⁠.

Curvature-dependent sputtering does not contribute to the leading order correction to the spatially averaged sputter yield $Y¯$ for a surface with nanoscale texture. Physically, the reason for this is that the reduced sputtering from a surface protrusion is offset by the increased sputtering from a depression. However, it can be shown that curvature-dependent sputtering does affect $Y¯$ at the next-to-leading order.

Among the nonlocal effects that can only be studied in detail by MC simulations (i.e., redeposition and secondary sputtering by recoils and reflected ions), the significance of sputtering by recoils is the most remarkable. It is well known that redeposition reduces the effective sputter yield and that it leads to a reduction of sidewall angles and achievable depths of holes and trenches milled by focused ion beams.^41,42 Sputtering by reflected ions is known to cause microtrenches at the bottom of broader trenches or near steep sidewalls.^42–44 Sputtering by re-entered recoils, however, has been largely overlooked so far. The one exception is a study in which simulations showed its influence on the depth of trenches milled by a 30 keV Ga focused ion beam.⁴⁵ Our results show that sputtering by re-entered recoils is also significant on surfaces with significantly smaller slopes and at much lower ion energies.

The role of sputtering by re-entered recoils may come as a surprise since the energy spectrum of sputtered atoms is known to be strongly peaked at low energies.⁴⁶ One reason for the recoils’ significant contribution to sputtering from the moderately sloped surfaces studied in this paper is that they strike the surface at a large angle when they re-enter the target. As shown in Fig. 7, for higher incidence angles, the sputter yield is larger than at normal incidence, and it decreases less rapidly as the energy is decreased than at lower incidence angles. Moreover, analysis of the energy and angular spectra of the recoils shows that the fraction of energetic recoils is quite sizable for large exit angles in the forward direction, i.e., for those directions that may lead to re-entry into the target. This is especially true for large ion incidence angles. For example, the fraction of recoils sputtered with energies greater than 50 eV at exit angles around $75°$ is almost 50% for 1 keV Ar ions incident at $75°$⁠. A more detailed analysis will be given elsewhere.

In the analytical part of this paper, we found the spatially averaged sputter yield $Y¯$ for non-planar surfaces that have slowly varying surface height $h=h(x,y)$⁠. Nonlocal effects like redeposition of sputtered material and secondary sputtering were neglected initially. The leading order corrections to the average sputter yield are proportional to the spatial averages of $(∂h/∂x)2$ and $(∂h/∂y)2$⁠, and the constants of proportionality are related to the sputter yield of a flat surface. For a surface with sinusoidally varying surface height, the leading order correction to $Y¯$ is proportional to $(A/λ)2$⁠, where $A$ and $λ$ are the amplitude and wavelength of the sinusoid, respectively. We demonstrated that $Y¯$ decreases with the amplitude $A$ for a range of angles of incidence. We also applied our result to statistically rough surfaces and determined how $Y¯$ depends on the correlation length of the surface. Finally, for normal-incidence ion bombardment of a low amplitude sinusoidal surface, the contribution of redeposition to $Y¯$ was shown to be proportional to $(A/λ)2$⁠.

Our theory agrees remarkably well with our Monte Carlo simulations of sputtering from surfaces with a sinusoidally varying surface height. The simulations confirm that the leading order correction to $Y¯$ due to redeposition is approximately proportional to $(A/λ)2$ for normal-incidence bombardment and show that this relationship also holds for oblique ion incidence. Our simulations also showed that atoms that are sputtered and then strike the surface can themselves cause significant sputtering. This effect is usually neglected in continuum simulations of the surface evolution of ion-bombarded surfaces. In the case studied here (1 keV Ar$+$ incident on a silicon target), it amounts to 20%–30% of the effect of redeposition. Therefore, the combined effect of redeposition and sputtering by recoils is a net reduction in the sputter yield.

One of us (R.M.B.) thanks Professors Patrick and Stephen Shipman for helpful discussions. This work was supported by Grant No. DMR-2116753 awarded by the U.S. National Science Foundation. G.H. is a member of the COST Action FIT4NANO CA19140, http://www.fit4nano.eu/. The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Programme.

The authors have no conflicts to disclose.

R. Mark Bradley: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Writing – original draft (equal). Gerhard Hobler: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Software (equal); Visualization (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In evaluating Eqs. (18) and (19), the first and second derivatives of the sputter yield $Y$ with respect to the incidence angle $θ$ have to be determined for flat targets. Differentiation is a standard problem of numerical mathematics. Approximations have been given by Bickley⁴⁷ for first and higher derivatives given function values at equidistant abscissae. Using three- and five-point stencils of equidistant $θ$ values, $θj=θi+(j−i)h$⁠, we can write

and

respectively. Given exact values of $Y(θi)$⁠, the quality of the approximations increases with order $h2$ and $h4$⁠, respectively, when the spacing $h$ between $θ$ values is decreased. However, since we are using Monte Carlo simulations to determine $Y(θi)$⁠, each $Y$ value has a statistical uncertainty with standard deviation $δY(θi)$⁠, which depends on the number of simulated impacts. This results in statistical uncertainties of the first and second derivatives, which increase as $h−1$ and $h−2$⁠, respectively, as $h$ is decreased.

We define the overall uncertainty of the sputter yield derivatives $Y′$ and $Y′′$ as the root mean square of the statistical error stemming from the Monte Carlo simulations and the systematic error by the finite order of the approximation. The former follows from the laws of error propagation and the $δY(θi)$ values. The latter we estimate by the difference between three-point and five-point approximation. Generally, we find the minimum errors of $Y′$ for $h=5°$⁠, $2.5°$⁠, and $1.25°$ for incidence angles in the ranges $θ<30°$⁠, $30°≤θ<60°$⁠, and $θ≥60°$⁠, respectively, and twice these spacings for calculating $Y′′$⁠. We, therefore, take these $θ$-dependent $h$ values and the five-point approximation for calculating the derivatives.