We find the spatially averaged sputter yield Y¯ analytically for non-planar surfaces that have slowly varying surface heights h=h(x,y). To begin, nonlocal effects like redeposition of sputtered material and secondary sputtering are neglected. We show that the leading order corrections to Y¯ are proportional to the spatial averages of (h/x)2 and (h/y)2. The constants of proportionality can be written in terms of the first and second derivatives of the sputter yield of a flat surface with respect to the ion incidence angle θ. For a range of θ values, Y¯ is a decreasing function of the amplitude of the surface texture. We also determine how the contribution of redeposition to Y¯ depends on the amplitude and characteristic lateral length scale of the surface morphology. As a test of our theory and to quantify the roles of redeposition and secondary sputtering, we performed Monte Carlo simulations of sputtering from Si targets with sinusoidal surfaces by 1 keV Ar+ ions. The theory agrees remarkably well with our Monte Carlo simulations. Our simulations also lead to the notable result that atoms that are sputtered and then strike the surface can themselves cause significant sputtering.

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.3 

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.8 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 simulations19 and in recent analytical work,20 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 xy 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,

h(x,y,0)=AH(ϵx,ϵy),
(1)

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 xy 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

ht(x,y,0)=f(θ;hx,hy,hxx,hxy,hyy,hxxx,hxxy,hxyy,hyyy,hxxxx,).
(2)

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

ht=f(θ;ϵAHX,ϵAHY,ϵ2AHXX,ϵ2AHXY,ϵ2AHYY,),
(3)

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

f0(θ)f(θ;0,0,),
(4)
f1(θ)hxf(θ;hx,0,0,)|hx=0,
(5)
f2(θ)hyf(θ;0,hy,0,0,)|hy=0,
(6)
f3(x,y,θ)hxxf(x,y,θ;0,0,hxx,0,0,)|hxx=0,
(7)

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 ith and jth arguments that appear after the semicolon, evaluated for all the arguments after the semicolon set equal to zero. For example,

f1,3(x,y,θ)hxhxxf(x,y,θ;hx,0,hxx,0,0,)|hx=hxx=0.
(8)

We obtain

ht=f0+ϵAf1HX+ϵ2Af3HXX+ϵ2Af5HYY+12ϵ2A2f1,1HX2+12ϵ2A2f2,2HY2.
(9)

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

ht=f0+f1hx+f3hxx+f5hyy+Λ1hx2+Λ2hy2,
(10)

where Λ112f1,1 and Λ212f2,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 22h does not appear because it is of order ϵ4.

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

ΩY¯JL2cosθ=0L0Lht(x,y,0)dxdy,
(11)

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

Y¯=1ΩJcosθ(f0+Λ1hx2+Λ2hy2),
(12)

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

hi21L20L0Lhi2(x,y,0)dxdy,
(13)

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.25 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

f0=ΩJY(θ)cosθ,
(14)
Λ1=12ΩJ[Y(θ)cosθ2Y(θ)sinθ],
(15)

and

Λ2=12ΩJY(θ)cosθcotθ,
(16)

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

Y¯=Y(θ)+μ1(θ)hx2+μ2(θ)hy2,
(17)

where

μ1(θ)12[Y(θ)2Y(θ)tanθ]
(18)

and

μ2(θ)12Y(θ)cotθ.
(19)

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

h(x,0)=1Lkhkeikx,
(20)

where xxx^+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

hk=1L0L0Leikxh(x,0)dxdy.
(21)

We find that

Y¯=Y(θ)+1L2k[μ1(θ)kx2+μ2(θ)ky2]|hk|2.
(22)

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

h=h(x)=Acos(kx),
(23)

then

Y¯=Y(θ)+2π2μ1(θ)(Aλ)2,
(24)

where k2π/λ is the wave number. Similarly, if

h=h(y)=Acos(ky),
(25)

then

Y¯=Y(θ)+2π2μ2(θ)(Aλ)2.
(26)

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 xy plane into N2 square cells, each having side length L0. Clearly,

hx2=1L20L0Lhx2(x,y,0)dxdy=1N2n=1N2hx2n,
(27)

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

C(xx)=[h(x,0)h(x,0)]2stat.
(28)

We assume that

C(xx)=A2g(xxξ,yyξ),
(29)

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.28 We differentiate C(xx) with respect to x and x and then take the limit as xx. This gives

hx2(x,0)stat=A2ξ2g11(0,0),
(30)

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

hx2=A2ξ2|g11(0,0)|.
(31)

In exactly the same way,

hy2=A2ξ2|g22(0,0)|.
(32)

Equation (17) now yields our final result,

Y¯=Y(θ)+μ1(θ)A2ξ2|g11(0,0)|+μ2(θ)A2ξ2|g22(0,0)|.
(33)

For a rough surface that is statistically isotropic, ξ=ξξ, g11(0,0)=g22(0,0) and

Y¯Y(θ)[μ1(θ)+μ2(θ)]A2ξ2.
(34)

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 that29 

Q(x)=ΩJY(0)γd2x[hh(xx)h][hh(xx)h]|xx|4×s(x+hz^,x+hz^).
(35)

Here, hh(x,0), hh(x,0), d2xdxdy, 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+hz^) is equal to 1 if there is a line of sight between the surface points x+hz^ and x+hz^ and is otherwise equal to zero.

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

ΩJL2Y¯redep=0Ldx0LdyQ(x),
(36)

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

Y¯redep=Y(0)γ(aL)2F,
(37)

where

F[HH(XX)XH][HH(XX)XH]|XX|4×s(X+Hz^,X+Hz^)d2Xd2X.
(38)

Here, XXx^+Yy^, XXx^+Yy^, d2XdXdY, d2XdXdY, Xx^X+y^Y, and Xx^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

Y¯redep=JY(0)αγ(Al)2.
(39)

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,

fredep=fredep(θ;hx,hy,hxx,hxy,hyy,hxxx,hxxy,hxyy,hyyy,hxxxx,).
(40)

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%,34 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 grid35 whose spacing was chosen to be 0.1a.

The projection of the incidence direction onto the xy 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 potential36 and the mixed Lindhard37 and Oen–Robinson38 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,17 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 Y0.17 Because of this uncertainty, we will not discuss incidence angles larger than 80°.

FIG. 1.

Sputter yield Y of a flat target as a function of incidence angle θ as obtained from the Monte Carlo simulations.

FIG. 1.

Sputter yield Y of a flat target as a function of incidence angle θ as obtained from the Monte Carlo simulations.

Close modal

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.

FIG. 2.

Coefficients μ1 and μ2 given by Eqs. (18) and (19) as a function of incidence angle θ.

FIG. 2.

Coefficients μ1 and μ2 given by Eqs. (18) and (19) as a function of incidence angle θ.

Close modal

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.

FIG. 3.

Average sputter yield of a sinusoidal surface divided by the sputter yield of the flat target as a function of (A/λ)2 for incidence angles of 0°, 15°, 30°, 45°, 60°, and 75°. (a) Parallel mode and (b) perpendicular mode. In the simulations, re-entry of ions and recoils into the target was suppressed.

FIG. 3.

Average sputter yield of a sinusoidal surface divided by the sputter yield of the flat target as a function of (A/λ)2 for incidence angles of 0°, 15°, 30°, 45°, 60°, and 75°. (a) Parallel mode and (b) perpendicular mode. In the simulations, re-entry of ions and recoils into the target was suppressed.

Close modal

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.

FIG. 4.

Average sputter yield of a sinusoidal surface for the parallel mode as a function of A/λ for the incidence angles 0°, 15°, 30°, 45°, 60°, and 75° (going from the top left to the bottom right panel). For the meaning of the different simulation modes, see the main text. The vertical dashed lines in the plots for θ=60° and θ=75° indicate the shadowing limit. The shaded areas around the theory lines indicate the statistical error that comes from the uncertainty in the Monte Carlo values for Y(θ).

FIG. 4.

Average sputter yield of a sinusoidal surface for the parallel mode as a function of A/λ for the incidence angles 0°, 15°, 30°, 45°, 60°, and 75° (going from the top left to the bottom right panel). For the meaning of the different simulation modes, see the main text. The vertical dashed lines in the plots for θ=60° and θ=75° indicate the shadowing limit. The shaded areas around the theory lines indicate the statistical error that comes from the uncertainty in the Monte Carlo values for Y(θ).

Close modal

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°.

FIG. 5.

Same as Fig. 4, but for the perpendicular mode.

FIG. 5.

Same as Fig. 4, but for the perpendicular mode.

Close modal

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.

FIG. 6.

Contribution of redeposition to the average sputter yield ratio Y¯/Y of a sinusoidal surface as a function of A/λ, shown in a double-logarithmic plot (colored lines with error bars). The solid lines represent proportionality to (A/λ)2. (a) Parallel mode and (b) perpendicular mode.

FIG. 6.

Contribution of redeposition to the average sputter yield ratio Y¯/Y of a sinusoidal surface as a function of A/λ, shown in a double-logarithmic plot (colored lines with error bars). The solid lines represent proportionality to (A/λ)2. (a) Parallel mode and (b) perpendicular mode.

Close modal

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θ.27 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.2 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.45 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.46 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.

FIG. 7.

Dependence of the sputter yield of a flat Si target for Si ions between 10 eV and 1 keV and incidence angles as given in the legend. The simulations were performed with 105 incident ions.

FIG. 7.

Dependence of the sputter yield of a flat Si target for Si ions between 10 eV and 1 keV and incidence angles as given in the legend. The simulations were performed with 105 incident ions.

Close modal

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 Bickley47 for first and higher derivatives given function values at equidistant abscissae. Using three- and five-point stencils of equidistant θ values, θj=θi+(ji)h, we can write

Y(θ0)=Y(θ1)Y(θ1)2h+O(h2),
(A1)
Y(θ0)Y(θ1)2Y(θ0)+Y(θ1)h2+O(h2),
(A2)

and

Y(θ0)Y(θ2)8Y(θ1)+8Y(θ1)Y(θ2)12h+O(h4),
(A3)
Y(θ0)Y(θ2)+16Y(θ1)30Y(θ0)+16Y(θ1)Y(θ2)12h2+O(h4),
(A4)

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 h1 and h2, 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.

1.
P.
Sigmund
,
J. Mater. Sci.
8
,
1545
(
1973
).
2.
R. M.
Bradley
and
J. M. E.
Harper
,
J. Vac. Sci. Technol. A
6
,
2390
(
1988
).
3.
4.
M. A.
Makeev
and
A.-L.
Barabasi
,
Appl. Phys. Lett.
72
,
906
(
1998
).
5.
M. A.
Makeev
and
A.-L.
Barabási
,
Appl. Phys. Lett.
73
,
2209
(
1998
).
6.
M. A.
Makeev
and
A.-L.
Barabási
,
Nucl. Instrum. Meth. B
222
,
316
(
2004
).
7.
M. A.
Makeev
and
A.-L.
Barabási
,
Nucl. Instrum. Meth. B
222
,
335
(
2004
).
8.
G.
Hobler
,
R. M.
Bradley
, and
H. M.
Urbassek
,
Phys. Rev. B
93
,
205443
(
2016
).
9.
M.
Küstner
,
W.
Eckstein
,
V.
Dose
, and
J.
Roth
,
Nucl. Instrum. Meth. B
145
,
320
(
1998
).
10.
M.
Küstner
,
W.
Eckstein
,
E.
Hechtl
, and
J.
Roth
,
J. Nucl. Mater.
265
,
22
(
1999
).
11.
D. N.
Ruzic
,
Nucl. Instrum. Meth. B
47
,
118
(
1990
).
12.
T.
Kenmotsu
,
Y.
Yamamura
,
T.
Muramoto
, and
N.
Hirotani
,
Nucl. Instrum. Meth. B
228
,
369
(
2005
).
13.
A.
Hu
and
A.
Hassanein
,
Nucl. Instrum. Meth. B
281
,
15
(
2012
).
14.
J.
Drobny
,
A.
Hayes
,
D.
Curreli
, and
D. N.
Ruzic
,
J. Nucl. Mater.
494
,
278
(
2017
).
15.
V. I.
Shulga
,
Nucl. Instrum. Meth. B
339
,
8
(
2014
).
16.
V.
Shulga
,
Radiat. Eff. Def. Sol.
172
,
610
(
2017
).
17.
S. J.
Lindsey
,
G.
Hobler
,
D.
Maciążek
, and
Z.
Postawa
,
Nucl. Instrum. Meth. B
393
,
17
(
2017
).
18.
19.
C.
Cupak
,
P.
Szabo
,
H.
Biber
,
R.
Stadlmayr
,
C.
Grave
,
M.
Fellinger
,
J.
Brötzner
,
R.
Wilhelm
,
W.
Möller
,
A.
Mutzke
, and
M. V.
Moro
,
Appl. Surf. Sci.
570
,
151204
(
2021
).
20.
P. S.
Szabo
,
C.
Cupak
,
H.
Biber
,
N.
Jäggi
,
A.
Galli
,
P.
Wurz
, and
F.
Aumayr
,
Surf. Interfaces
30
,
101924
(
2022
).
21.
M. L.
Nietiadi
and
H. M.
Urbassek
,
Appl. Phys. Lett.
103
,
113108
(
2013
).
22.
H. M.
Urbassek
,
R. M.
Bradley
,
M. L.
Nietiadi
, and
W.
Möller
,
Phys. Rev. B
91
,
165418
(
2015
).
23.
R.
Cuerno
and
A.-L.
Barabási
,
Phys. Rev. Lett.
74
,
4746
(
1995
).
24.
M.
Makeev
,
R.
Cuerno
, and
A.
Barabási
,
Nucl. Instrum. Meth. B
197
,
185
(
2002
).
25.
If Dirichlet boundary conditions are adopted and instead of periodic boundary conditions, the contribution of f1hx to Y¯ is of order L1. This contribution tends to zero as the domain size L becomes large.
26.
D. A.
Pearson
and
R. M.
Bradley
,
J. Phys.: Condens. Matter
27
,
015010
(
2015
).
27.
H.
Hofsäss
and
O.
Bobes
,
Appl. Phys. Rev.
6
,
021307
(
2019
).
28.
A.-L.
Barabási
and
H. E.
Stanley
,
Fractal Concepts in Surface Growth
(
Cambridge University Press
,
1995
).
29.
R. M.
Bradley
,
Phys. Rev. B
83
,
75404
(
2011
).
30.
G.
Hobler
,
Nucl. Instrum. Meth. B
96
,
155
(
1995
).
31.
C.
Ebm
and
G.
Hobler
,
Nucl. Instrum. Meth. B
267
,
2987
(
2009
).
32.
G.
Hobler
,
M. L.
Nietiadi
,
R. M.
Bradley
, and
H. M.
Urbassek
,
J. Appl. Phys.
119
,
245105
(
2016
).
33.
R. M.
Bradley
and
G.
Hobler
,
J. Appl. Phys.
129
,
194301
(
2021
).
34.
D. W.
Oh
,
S. K.
Oh
,
H. J.
Kang
,
H. I.
Lee
, and
D. W.
Moon
,
Nucl. Instrum. Meth. Phys. Res. B
190
,
598
(
2002
).
35.
G.
Hobler
and
S.
Selberherr
,
IEEE Trans. Comp.-Aided Des.
8
,
450
(
1989
).
36.
J. F.
Ziegler
,
J. P.
Biersack
, and
U.
Littmark
,
The Stopping and Range of Ions in Matter
(
Pergamon
,
New York
,
1985
).
37.
J.
Lindhard
and
M.
Scharff
,
Phys. Rev.
124
,
128
(
1961
).
38.
O. S.
Oen
and
M. T.
Robinson
,
Nucl. Instr. Meth.
132
,
647
(
1976
).
39.
J. M.
Fluit
,
P. K.
Rol
, and
J.
Kistemaker
,
J. Appl. Phys.
34
,
690
(
1963
).
40.
41.
J. D.
Casey
,
A. F.
Doyle
,
R. G.
Lee
,
D. K.
Stewart
, and
H.
Zimmermann
,
Microelectron. Eng.
24
,
43
(
1994
).
42.
S.
Lindsey
and
G.
Hobler
,
Nucl. Instrum. Meth. B
282
,
12
(
2012
).
43.
S. M.
Rossnagel
and
R. S.
Robinson
,
J. Vac. Sci. Technol. A
1
,
426
(
1983
).
44.
E.
Platzgummer
,
A.
Biedermann
,
H.
Langfischer
,
S.
Eder-Kapl
,
M.
Kuemmel
,
S.
Cernusca
,
H.
Loeschner
,
C.
Lehrer
,
L.
Frey
,
A.
Lugstein
, and
E.
Bertagnolli
,
Microelectron. Eng.
83
,
936
(
2006
).
45.
S. J.
Lindsey
, “Computer simulation of sputtering,” Ph.D. thesis (TU Wien, Vienna, 2015).
46.
47.
W. G.
Bickley
,
Math. Gaz.
25
(
263
),
19
(
1941
).