Second order corrections to the sputter yield of a curved surface

The sputter yield of a solid surface depends not only on the angle of incidence of the impinging ion but also on the surface curvature.^1,2 This effect becomes especially important when surfaces with radii of curvature comparable to the ion penetration depth are sputtered.^3,4 The sputtering of elemental spheres has been studied in great detail^5–14 but the sputtering of arbitrary surface morphologies¹⁵ and of alloy spheres and wires^16–18 has also been investigated. Curvature-dependent sputtering is relevant to ion irradiation of dust grains in space, the sputtering of nanoparticles supported on flat surfaces,^6,18 and to ion implantation of nanowires.^17,18 In addition, when an initially flat solid surface is bombarded with a broad ion beam, curvature-dependent sputtering can lead to a surface instability and to the emergence of self-organized nanoscale patterns.^2,19

Consider an ion that is normally incident on a solid sphere of radius $R$⁠. For $R$ much larger than the mean depth of nuclear energy deposition (NED) $a$⁠, the average sputter yield $Y$ is, to a good approximation, a linear function of $a/R$⁠,

where $Y∞$ is the yield from a flat surface and $η1$ is a positive, dimensionless constant.¹² Molecular dynamics and Monte Carlo (MC) simulations show that as $a/R$ is increased, $Y$ at first increases linearly but then goes through a maximum and finally tends monotonically to zero as $a/R$ becomes large.¹² For small $a/R$⁠, there are corrections to Eq. (1) of order $(a/R)2$⁠.

The generalization of Eq. (1) to the normal-incidence bombardment of a curved surface of arbitrary shape is

where $H$ is the mean curvature of the surface at the point of impact.¹⁵ (For a sphere, $H=−1/R$⁠.) Equation (2) shows that the sputter yield is increased (reduced) for impact on convex (concave) surfaces relative to impact on a flat surface. Based on our knowledge of the sputtering of spherical particles, we would expect there to be corrections to Eq. (2) of order $(aH)2$⁠.

In this paper, we will show that there is indeed a correction to Eq. (2) that is proportional to the square of the mean curvature. However, we will also demonstrate that there is a correction of the same order that is proportional to the Gaussian curvature $K$⁠,

where $η2$ and $η3$ are dimensionless constants. The correction terms proportional to $H2$ and $K$ become increasingly important as the radius of curvature is reduced. We will show that, in general, the coefficients $η1$⁠, $η2$⁠, and $η3$ can be written in terms of moments of the crater function and then explicitly evaluate them for the special case in which the Sigmund model of ion sputtering¹ is valid. We will also use MC simulations to obtain estimates of $η1$⁠, $η2$⁠, and $η3$ for silicon targets and for three different noble gas ion species and three ion energies and then compare their values to the analytical expressions obtained using the Sigmund theory. Finally, the implications of our results for pattern formation on solid surfaces bombarded by a broad ion beam will be discussed.

Consider the impact of a single ion with energy $E$ at a point on the surface of a solid object composed of an elemental material. The angle of incidence of the ion will be denoted by $θ$⁠. It will be helpful to consider general values of $θ$ before specializing to the case of normal incidence, $θ=0$⁠. We place the origin $O$ at the point of impact and put the $z$ axis along the outward-pointing normal to the surface. We orient the $x$- and $y$-axes so that the ion’s direction of incidence $−e^$ lies in the $x−z$ plane and is given by $−e^=−x^sin⁡θ−z^cos⁡θ$⁠.

We wish to compute the sputter yield $Y$⁠, which is roughly speaking defined to be the average number of sputtered atoms per incident ion. However, we must touch on two technical points in order to define $Y$ precisely. First, an atom that is sputtered but then redeposits on the surface of the solid will be counted as a sputtered atom and so it contributes to the sputter yield $Y$⁠. Second, the incident ion could be reflected from the surface or exit from it and then strike the solid surface once again at another point, possibly causing additional sputtering. This additional sputtering will not be included in the sputter yield $Y$ but is expected to be small for $θ=0$⁠. If the solid surface is convex, then neither of these two possibilities can occur, i.e., there can be no redeposition of sputtered atoms and ions cannot reimpinge on the solid surface.

Let $u=u(x,y)$ be the height of the solid surface above the point $(x,y)$ in the $x−y$ plane. To begin, we will assume that $u(x,y)$ is a smooth, single-valued function that is defined for all $x$ and $y$⁠. Later, we will discuss sputtering from objects like spheres that have a double-valued surface height that is not defined for all $x$ and $y$⁠. Note that since the $z$ axis is normal to the surface at the origin, $∂u/∂x=∂u/∂y=0$ for $x=y=0$⁠.

We will develop an approximate expression for the sputter yield for a surface of the form

where $R>0$ has dimensions of length and $U$ is dimensionless. $R$ is the characteristic length scale of the surface. The paraboloid of revolution $u=−(x2+y2)/(2R)$⁠, for example, satisfies Eq. (4). The paraboloid’s characteristic length scale $R$ is just the radius of curvature at its apex.

The value of the crater function $F$ at the point $(x,y)$ is defined to be minus the average change in the surface height $u$ above the point $(x,y)$ in the $x−y$ plane as a result of a single ion impact at $x=y=0$⁠.^20,21 The crater function $F$ depends on $x$⁠, $y$⁠, and the angle of incidence $θ$⁠. It also depends in principle on the shape of the entire surface, or, equivalently, on all of the spatial derivatives of $u(x,y)$ evaluated at $x=y=0$⁠. We will write

where the subscripts on $u$ denote partial derivatives. The partial derivatives of $u$ that appear on the right-hand side of Eq. (5) are all to be evaluated at $x=y=0$⁠. We assume that $F$ is known a priori from another theory or from atomistic simulations.

The sputter yield $Y$ is the average number of sputtered atoms and is given by

where $Ω$ is the atomic volume. For convenience, we let $M≡ΩY$ and write

We also introduce the crater function moments

and so on. Similarly, for positive integers $i$ and $j$⁠, $Mi,j(θ)$ will denote the partial derivative of $M(θ;uxx,uxy,uyy,uxxx,…)$ 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,

Let $a$ be a length scale that characterizes the size of the collision cascade. For the Sigmund model of ion sputtering,¹ for example, $a$ may be taken to be the average depth of NED. In  Appendix A, we find an approximate expression for $Y$ that is valid for a surface that has a characteristic length scale $R$ that is much larger than $a$⁠. More precisely, we let $ϵ≡a/R$ and work to second order in $ϵ$⁠. We also specialize to the case of normal incidence. We obtain Eq. (3) and find that

and

Equations (13)–(15) relate the coefficients $η1$⁠, $η2$⁠, and $η3$ in Eq. (3) to moments of the crater function. Note that because $ux=uy=0$ at the origin, the mean and Gaussian curvatures are $H=(uxx+uyy)/2$ and $K=uxxuyy−uxy2$⁠, respectively.

Equation (3) is valid to order $ϵ2$⁠, no matter what the nature of the crater function. For a flat surface, $Y=Y∞$⁠. The second term on the right-hand side of Eq. (3) is the correction term of order $ϵ$ for a curved surface. The third and fourth terms on the right-hand side are the correction terms of order $ϵ2$⁠.

In the derivation of Eq. (3) given in  Appendix A, $u(x,y)$ is taken to be single-valued function of $x$ and $y$ for all $x$ and $y$⁠. In  Appendix B, we show that Eq. (3) is also valid for sputtering from a solid object like a sphere that has a double-valued height function $u(x,y)$ provided that the crater function falls off with sufficient rapidity away from the point of impact. Naturally, the correction terms in Eq. (3) will only be appreciable if the dimensions of the object do not exceed a few micrometers for ion energies commonly used in sputtering experiments.

In this subsection, we will obtain an approximate expression for the sputter yield $Y$ for $θ=0$ using the Sigmund model of ion sputtering and the results given by Eqs. (3) and (13)–(15). Thus, we will need to find $M0$⁠, $M1$⁠, $M1,1$⁠, and $M1,3$ for this model.

Consider an ion of energy $E$ in an infinite medium composed of the target material. Initially, the ion is at the origin and is moving in the $−z$ direction. Let $ED(x,y,z)$ be the average energy density deposited at the point $(x,y,z)$ by nuclear collisions. The integral of $ED(x,y,z)$ over all space is the total energy $ED$ deposited in the material by nuclear collisions. In the Sigmund model, it is assumed that the NED distribution $ED(x,y,z)$ is a Gaussian,

Here, $a$ is the average depth of NED and $α$ and $β$ are the longitudinal and transverse straggling lengths, respectively. The contours of equal NED are ellipsoids of revolution centered at the point $−az^$ with the $z$ axis as their axis of symmetry.

We now return to the problem of computing the sputter yield $Y$ for the impact of a normally incident ion on a curved solid surface. In the Sigmund model, the average number of atoms sputtered from a small surface element $dA$ centered on a point $r=(x,y,z)$ is assumed to be proportional to the NED density, i.e.,

where the sputtering efficiency $λ$ is a constant. The sputter yield $Y$ is $∫Sn(r)dA$⁠, where the integral is taken over the solid surface $S$⁠.

This discussion shows that Sigmund’s model is based on three assumptions:²² 

1. the sputter yield from a surface element is proportional to the local NED density and the proportionality factor (the sputtering efficiency $λ$⁠) is a constant;

2. the NED distribution is independent of the orientation and shape of the surface and is identical to the one in an infinite medium; and

3. the NED distribution in an infinite medium is Gaussian.

In addition to finding the coefficients $η1$⁠, $η2$⁠, and $η3$ for the Sigmund model, we will also find them for a generalized version of the Sigmund model. In this version of the Sigmund model, assumptions (1) and (2) are made but assumption (3) is not. Instead, we will only assume that $ED(x,y,z)$ is rotationally symmetric about the $z$ axis.

The sputter yield for the generalized Sigmund model is

where $d2x=dxdy$⁠. It follows that

Note the partial derivatives of $u$ that appear in Eq. (19) are to be evaluated at the origin.

The required crater function moments are now readily evaluated. Inserting them in Eqs. (13)–(15), we find that

and

where

and

Here, $E~2D(x~,z~)$ denotes the 2D NED density, which is obtained from $E~D(x~,y~,z~)$ by integration along the $y$ direction,

The tildes indicate scaling with respect to the mean NED depth $a$⁠: $x~≡x/a$⁠, $y~≡y/a$⁠, $z~≡z/a$⁠, $E~2D≡a2E2D$⁠, and $E~D≡a3ED$⁠.

For the Sigmund model, Eq. (16) applies and a simple calculation yields

and

Equations (20) and (21) then give

and

where $aα≡a/α$ and $aβ≡a/β$⁠. Inserting Eq. (31) into Eq. (3) and remaining only terms up to order $ϵ$⁠, we obtain

a result that has already been obtained elsewhere.¹⁵ 

In previous work, sputtering from spheres^11,12 and cylinders¹⁵ was considered. Here, we consider in addition convex parabolic cylinders, $u(x,y)=−y2/2R$⁠, convex paraboloids of revolution, $u(x,y)=−(x2+y2)/2R$⁠, and catenoids. The catenoids we study have their axis of symmetry parallel to the $x$ axis, along the line with $z=−R$ and $y=0$⁠,

This axis was chosen so that the origin lies on the surface of the catenoid. The parabolic cylinders and paraboloids are similar to the (circular) cylinder and the sphere in that they have the same curvatures at the origin. In contrast to cylinder and sphere, they have constant second derivatives $uxx$ and $uyy$ rather than constant curvatures. Moreover, they do not allow forward sputtering, and so avoid the associated contribution to the sputter yield. The catenoid is particularly interesting as it is a minimal surface, i.e., the mean curvature $H$ vanishes everywhere on the surface. For the spheres, convex paraboloids of revolution, and catenoids, $R$ is the radius of curvature along the $x$ and $y$ directions at the origin. For circular and convex parabolic cylinders, $R$ is the radius of curvature along the $y$ direction, and the radius of curvature along the $x$ direction is infinite.

For circular or convex parabolic cylinders, $H=−1/2R$ and $K=0$ at $x=y=0$⁠. Hence, according to Eq. (3),

For spheres or convex paraboloids, $H=−1/R$ and $K=1/R2$ at $x=y=0$⁠, and thus

For the catenoid, $H=0$ and $K=−1/R2$ at $x=y=0$⁠. The only nonvanishing correction term on the right-hand side of Eq. (3) is, therefore, the one containing $η3$⁠,

We will also consider the sputter yields of five concave surfaces. The sputter yield of a concave paraboloid of revolution described by $u(x,y)=(x2+y2)/2R$ is obtained by replacing $R$ by $−R$ in Eq. (37). This is also the sputter yield for the normal-incidence impact of an ion with the wall of a spherical cavity of radius $R$ in an infinite medium. Similarly, the sputter yield of a concave parabolic cylinder with $u(x,y)=y2/2R$ is given by Eq. (36) with $R$ replaced by $−R$⁠. This is also the sputter yield for the normal-incidence impact of an ion with the wall of a cylindrical cavity of radius $R$ in an infinite medium. Finally, we will also consider ions impinging in the negative $z$ direction on the wall of a catenoidal cavity in an infinite medium. The wall of the cavity is given by Eq. (35). The sputter yield is given by Eq. (38) in this case, just as it is for a solid catenoid.

In this section, we present MC simulations to demonstrate that the second-order approximation given by Eq. (3) significantly extends the range of curvatures for which the theory provides a good approximation to the sputter yield. Moreover, we quantitatively compare the predictions of the Sigmund theory and the generalized Sigmund theory (Sec. II B) with the MC results.

We used the MC simulator IMSIL^23,24 to determine sputter yields. IMSIL is able to consider targets with translational symmetry in the $x−y$ plane [one-dimensional (1D) layered targets] and targets with translational symmetry along the $y$ axis [two-dimensional (2D) geometries]. For 1D targets, the surfaces and any interfaces between materials are specified by $z$ values, while the target is assumed to be homogeneous in the $x$ and $y$ directions. For 2D targets, the surfaces and interfaces are specified by closed polygons in the $x−z$ plane, while the target is assumed to be homogeneous in the $y$ direction. As amorphous Si is the only material considered in this paper, there are no interfaces.

For the purpose of this work, IMSIL was enhanced by the implementation of a rotational symmetry mode with the symmetry axis parallel to one of the coordinate axes. In this mode, the distance $r$ from the symmetry axis takes the role of $z$ when the position of an ion or recoil is checked with respect to the target geometry. Thus, 1D geometry definitions specify surfaces with constant $r$⁠, and so cylinders and cylindrical cavities may be defined. 2D geometry definitions describe surfaces by specifying polygons in the $x−z$ plane and then rotating them about an axis parallel to the $x$ axis, and so paraboloids of revolution, catenoids, and spheres may be approximated. Other orientations of the symmetry axis are realized by permutations of $x$⁠, $y$⁠, and $z$⁠. For spheres, a spherical symmetry mode combined with a 1D geometry¹⁶ is also available. In this mode, the distance from the symmetry center takes the role of $z$ in geometry checks.

For the simulations presented in this work, we used 1D geometry definitions together with rotational and spherical symmetry to describe cylinders and spheres, respectively. Parabolas and catenaries were approximated by polygons with a resolution of about $a/200$ to describe parabolic cylinders, paraboloids, and catenoids. Since the polygons must be closed in IMSIL, we cut off the bodies at $z=±10a$⁠. While it would have been easy to hard-code these geometries, we used the approximation by polygons, since this is an existing feature of IMSIL. The results for the parabolic cylinder and paraboloid converge very nicely to those for the cylinder and sphere, respectively, at small curvatures. This provides evidence that not only was the resolution of the polygons sufficient, but also that the sputtering from a polygonal surface is consistently implemented in IMSIL.²⁵ 

To describe nuclear scattering, we used the pair-specific Ziegler–Biersack–Littmark (ZBL) potential²⁶ for the Si–Si interactions and the universal ZBL potential²⁶ for the ion-target atom interactions. The Lindhard model²⁷ was used for the electronic stopping power and the Oen–Robinson model²⁸ to describe the impact parameter dependence of electronic stopping, assuming an equipartition rule. NED distributions were obtained by recording the energy transferred to target atoms that are not displaced from their sites (subthreshold collisions) and the remaining energies of the ion and the recoils when their trajectories are terminated.²² Before being sputtered, the Si atoms had to surmount a potential well that followed the surface at a distance equal to the maximum impact parameter.²⁵ The height of the well was taken to be equal to the heat of sublimation (4.7 eV/atom). Because of the abruptness of the assumed well, refraction and reflection may be treated as if the well were planar^29,30 with a normal vector in the direction of the gradient of the distance function.²⁵ The number of reflections from the well was limited to one in order to avoid having atoms reflect back and forth between the surface potential well and the repulsive potentials of the target atoms. These oscillations would probably not occur in a real system due to the inevitable surface roughness that develops on an atomic scale during sputtering. The effect of suppressing these oscillations in the simulations was found to be negligibly small.

We set the cutoff energy for ion and recoil trajectories equal to 4.7 eV, since we were mainly interested in the sputter yield. A recoil with lower energy can neither leave the target nor can it generate another recoil that is able to leave the material. Consideration of lower energies might cause a slight broadening of the NED distribution. However, treatment of such low-energy recoils in the binary collision approximation is questionable and would have increased computation times. For details on the selection of collision partners and the treatment of the surface, we refer the reader to Refs. 22 and 31.

To find the NED distribution, its mean depth $a$⁠, longitudinal standard deviation $α$⁠, and transverse standard deviation $β$⁠, we performed simulations in an infinite medium in which the ions were started at the origin. To determine the parameters $η1$⁠, $η2$⁠, and $η3$⁠, we simulated the flat surface plus convex parabolic cylinders, convex paraboloids, and catenoids with $a/R=2−k$⁠, where $k=3,4,…,8$⁠. In each simulation, the target was bombarded with $108$ ions directed toward the origin along the negative $z$ axis. The sputter yields obtained have an accuracy of about $10−4$⁠. We then simultaneously fit the simulated sputter yields for a given ion species and energy with the functions $Y(a/R)$ given by Eqs. (36)–(38) using the method of least squares with $Y∞$⁠, $η1Y∞$⁠, $η2Y∞$⁠, and $η3Y∞$ as the fitting parameters. It is important to note that the same values of the fitting parameters were used in the fits of the data for the convex parabolic cylinders, convex paraboloids, and catenoids. This yielded our estimates of $η1$⁠, $η2$⁠, and $η3$⁠. In order to determine error bars, we repeated the fit 100 times with Gaussian deviations added to the mean value of each sputter yield, using standard deviations of $Y$ as determined by IMSIL.

Additional simulations, including bombardments of cylinders and spheres, were performed to investigate the range of validity of the approximation given by Eq. (3), which is valid to second order in $ϵ=a/R$⁠. To provide additional tests of the theory, simulations were also carried out for concave parabolic cylinders, concave paraboloids, and for spherical, cylindrical and catenoidal cavities in an infinite medium.

MC simulations and fits were performed for three ion species (Ne, Ar, and Xe) and three ion energies (1, 5, and 20 keV). The energies 1 and 5 keV are lower than the energies that have typically been studied in the past^12,15,16,31 and were included because pattern formation experiments are generally performed with ions that have energies lower than 20 keV. As an example of the accuracy of the fits, the results for 1 keV Ar ions normally incident on convex surfaces are shown in Fig. 1. Note the almost perfect fit for all surfaces up to curvatures of $a/R=0.25$ (dashed lines) even though the same values of $η1$⁠, $η2$⁠, and $η3$ were used in the fits of the data for the convex parabolic cylinders, convex paraboloids, and catenoids. This shows the accuracy of the second-order approximation given by Eq. (3). Beyond $a/R=0.25$⁠, the results for spheres, cylinders, and catenoids start to deviate from the predictions of Eq. (3), which is most likely primarily due to forward sputtering.^12,15 In contrast, the second-order approximation performs much better at large $a/R$ for the convex parabolic cylinder and the convex paraboloid, from which forward sputtering is not possible. In the inset in Fig. 1, one can also see that the first-order results Eq. (2) (dotted lines) start to deviate noticeably from the data around $a/R=0.05$⁠, and the second-order correction terms become increasingly important for larger values of $a/R$⁠. We obtained very similar results for six other choices of ion species and energy, as shown in the supplementary material.

For each choice of ion species and energy, we computed estimates of the coefficients $η1$⁠, $η2$⁠, and $η3$ using three different methods:

1. The MC results for the sputter yields of parabolic cylinders, paraboloids, and catenoids were simultaneously fit using the same values of $η1$⁠, $η2$⁠, and $η3$⁠, as in Fig. 1.

2. The 2D NED distribution in an infinite medium was determined by MC simulations. $η1$⁠, $η2$⁠, and $η3$ were then computed using Eqs. (20)–(27), which are valid for the generalized Sigmund model.

3. The values of $a$⁠, $α$⁠, and $β$ were determined from the 2D NED distributions and then input into Eqs. (31)–(33), the equations that apply for the Sigmund model.

Of these methods, Method 1 is expected to be the most accurate and Method 3 is expected to be the least.

To implement Method 2, we saved 2D NED distributions with a resolution of $a/25$ in our MC simulations in the infinite medium, and applied Eqs. (20)–(26). The first and second derivatives of $E~2D$ with respect to $z~$ were approximated by finite differences using the first two and three layers below $z~=0$⁠, respectively, of the NED distribution.

The estimates of $η1$⁠, $η2$⁠, and $η3$ are shown in Fig. 2. The agreement between the values of $η1$⁠, $η2$⁠, and $η3$ computed using the three different methods is quite good for 20 keV Ar ions. This is consistent with previous work in which the first-order theory based on the Sigmund model showed excellent agreement with MC simulations for 20 keV Ar ion impacts on silicon spheres¹² and cylinders.¹⁵ [In Fig. 4(a) of Ref. 15, the theoretical lines were drawn with incorrect slopes; the correct slopes would have given even better agreement at small curvatures.] The agreement between the results obtained using the three methods is only fair for 5 keV Ar ions. The results computed for 1 keV Ar ions using the Sigmund model differ widely from the values determined using Method 1, however. The same is true of the results for 1 keV Ne and Xe ions. This is not unexpected, because the Sigmund model is an asymptotic approximation in terms of the inverse energy.²⁹ As a consequence, as the ion’s energy is reduced, the predictions of the Sigmund model are expected to deviate to an increasing degree from the MC results. For the 1 keV ions, the values of $η1$⁠, $η2$⁠, and $η3$ computed using Methods 2 and 3 are relatively close to one another. This suggests that the inaccuracy of the Sigmund model at this energy does not stem primarily from the Gaussian approximation to the NED distribution in an infinite medium. Instead, the problem must lie with either Assumption 1 or Assumption 2, or both.

The values of $η1$⁠, $η2$⁠, and $η3$ for the 20 keV ions computed using Method 1 increase with the ion mass. The estimates determined using Method 3 do not reproduce this trend, but the estimates obtained using Method 2 do. This suggests that the Gaussian approximation to the NED distribution is responsible for this failing of the Sigmund model.

In Fig. 3, we show the results for 1 keV Ar ions impinging on the surfaces of the concave parabolic cylinders, the concave paraboloids, and the catenoidal cavities. Once again, the ions’ incidence is perpendicular at the origin. The dotted and dashed lines represent the model predictions. The same parameters $η1$⁠, $η2$⁠, and $η3$ were used as in Fig. 1. Excellent agreement is observed at moderate curvatures with slightly larger deviations than in Fig. 1 at the largest curvatures. The differences between the concave paraboloid and the spherical cavity, and between the concave parabolic cylinder and the cylindrical cavity, are much smaller than for the corresponding convex surfaces, since for all four of these concave surfaces, there is no forward sputtering.

Equations (1) and (2) were obtained analytically in past work. In the derivations of these results, however, it was assumed that the Sigmund model is valid.^11,12,15 This is problematic, since the Sigmund model is known to be a rather poor approximation, especially for light ions or high ion energies.²² In addition, the derivation of Eq. (1) given in Refs. 11 and 12 is not rigorous because a quadratic approximation to the shape of the surface was employed.^11,12 Similarly, the derivation of the more general result (2) in Ref. 15 is not rigorous because the dependence of the crater function $F$ on third and higher order spatial derivatives of $u$ was neglected.

In the present work, Eq. (2) was rigorously shown to be valid by carrying out a systematic expansion in the small parameter $ϵ≡a/R$⁠. The result was obtained without assuming that the Sigmund model is valid, and the dependence of $F$ on spatial derivatives of $u$ of arbitrarily high order was taken into account. In addition, the correction terms to the sputter yield $Y$ of order $ϵ2$ were found. We demonstrated that one of these corrections is proportional to the square of the mean curvature $H$ and the other is proportional to the Gaussian curvature $K$⁠. The correction terms become increasingly important as the characteristic length scale of the surface $R$ is reduced.

Our principal result, Eq. (3), also has noteworthy implications for the nanoscale patterns that emerge when an initially flat solid surface is bombarded with a broad ion beam. For example, suppose that the solid is bombarded with a broad ion beam with angle of incidence $θ$ and that the solid is rotated with angular velocity $ω$ about its normal. For simplicity, we assume that the target material is elemental and that a layer at the surface of the solid is amorphized by the ion bombardment. An instability results in the formation of a nontrivial surface morphology if $θ$ exceeds a critical value.³² If $ω$ is sufficiently large, the equation of motion (EOM) for the surface is rotationally invariant about the $z$ axis. In that regard, the EOM resembles the one that applies for normal-incidence bombardment, although there is no surface instability in the latter case. The EOM for the problem with oblique-incidence bombardment and sample rotation that is usually adopted is

where $u(x,y,t)$ denotes the height of the solid surface above the point $(x,y)$ in the $x−y$ plane at time $t$⁠, and $v0>0$⁠, $A>0$⁠, $B>0$⁠, $λ$⁠, and $c$ are constants.^32,33 Equation (39) is expected to apply only for low amplitude surface disturbances with small slopes.

The term $−A∇2u$ in EOM (39) stems from curvature-dependent sputtering² and ion-induced mass redistribution.^34–36 The physical origin of the term $−B∇2∇2u$ is thermally activated surface diffusion or ion-induced viscous flow,³⁷ while the term $λ(∇u)2$ describes the effect of the slope dependence of the sputter yield to lowest nontrivial order.^38,39 The term $c∇2(∇u)2$ is the so-called conserved Kuramoto–Sivashinsky (CKS) nonlinearity; it is second order in $u$ and fourth order in the gradient $∇$⁠. This term leads to coarsening of the surface pattern with time, i.e., to an increase in the characteristic lateral length scale and the surface width.^19,32

To lowest nontrivial order, the mean curvature of the surface $H$ is $12∇2u$⁠. Equation (39), therefore, includes the effect of curvature-dependent sputtering to the lowest nontrivial order in $H$⁠. Since the sputter yield $Y$ depends on $H$⁠, however, terms proportional to $H2$ and higher powers of $H$ should, in general, appear on the right-hand side of Eq. (39). Note that $H2≅14(∇2u)2$ is second order in $u$ and fourth order in $∇$⁠, i.e., it is of the same order as the CKS nonlinearity $c∇2(∇u)2$ that appears in Eq. (39). This means that Eq. (39) is actually inconsistent as it stands, and a term proportional to $H2$ must, in general, be included. In addition, the Gaussian curvature $K$ is $uxxuyy−uxy2$ to lowest order. This too is second order in $u$ and fourth order in $∇$⁠, and so, for the sake of consistency, a term proportional to $uxxuyy−uxy2$ must also appear in the equation of motion. To simplify the resulting EOM, we set $u~(x,y,t)=u(x,y,t)+v0t$⁠, where $v0$ is the rate of recession of a flat surface. After dropping the tilde, we obtain

where $μ$ and $ν$ are constants.

The terms proportional to $H2$ and $K$ that appear in Eq. (40) are analogs of the second-order correction terms that appear in Eq. (3), and they can lead to dramatic effects on the surface morphology. Bernoff and Bertozzi have studied Eq. (40) for the special case in which $c=0$ as a model of a moving, nearly planar interface between two stable phases of a material.⁴⁰ Through both rigorous analytical and numerical work, they showed that for a range of parameter choices and initial conditions, spike singularities develop in finite time. Naturally, singularities would not actually occur when a rotating surface is bombarded with a broad ion beam. Instead, either the CKS nonlinearity or an additional term not included in Eq. (40) would come into play as a spike grows and prevent the formation of a true singularity. However, sharp spikes could develop on the surface. Spike formation has not yet been observed in experiments in which a sample was concurrently rotated and bombarded with an obliquely incident broad ion beam,⁴¹ but it could be in the future.

In this paper, we studied the sputter yield $Y$ of a curved surface that is struck by a normally incident ion. We focused on the limit in which the characteristic length scale of the surface $R$ is large compared to the size of the collision cascade $a$⁠. The leading order correction to $Y$ is proportional to the mean curvature $H$ at the point of impact and is of order $ϵ=a/R$⁠. In addition, the correction terms of order $ϵ2$ were found. We demonstrated that one of these corrections is proportional $H2$ and that the other is proportional to the Gaussian curvature $K$ at the point of impact. The coefficients $η1$⁠, $η2$⁠, and $η3$ of the correction terms were shown to be related to moments of the crater function and were explicitly determined for the Sigmund model of ion sputtering and for a generalized version of the Sigmund model.

To test these predictions, we carried out Monte Carlo simulations of the sputtering of silicon targets with a variety of shapes and surface morphologies for three different noble gas ion species and three ion energies. The approximation to the sputter yield that is valid to order $ϵ2$ [Eq. (3)] was found to be in good agreement with our MC results for $a/R$ values up to 0.25 and is accurate for a much broader range of surface curvatures than the result that is valid to order $ϵ$⁠. The Sigmund values for $η1$⁠, $η2$⁠, and $η3$ are in reasonable accord with the corresponding MC estimates for 20 keV Ar ions. The Sigmund model does not reproduce the dependence of $η1$⁠, $η2$⁠, and $η3$ on the ion mass, however. The generalized Sigmund model yields better agreement with the ion mass dependence of the MC results for 20 keV ions. Finally, both the Sigmund model and the generalized Sigmund model give values for $η1$⁠, $η2$⁠, and $η3$ that differ markedly from the MC results for 1 keV ions. MC simulations rather than the Sigmund model should, therefore, be used to estimate the coefficients in Eq. (3) whenever possible.

We also discussed our results’ implications for the nanoscale patterns that emerge when an initially flat solid surface is bombarded with a broad ion beam. Our work suggests that more systematic experiments with oblique-incidence ion bombardment and concurrent sample rotation would be worthwhile and could lead to the observation of sharp spikes that protrude from the solid surface.

See the supplementary material for the values of the parameters $a$⁠, $α$⁠, and $β$ as determined with IMSIL, and for the sputter yield as a function of curvature for Ne, Ar, and Xe ions with energies of 1, 5, and 20 keV.

We are grateful to Peter Sigmund for illuminating correspondence. The work of R.M.B. was supported by Grant No. DMS-1814941 from the U.S. National Science Foundation.

The data that support the findings of this study are available within the article and its supplementary material.

In this Appendix, we derive the approximate expression for the sputter yield (3) and the formulas (13)–(15) that relate the coefficients $η1$⁠, $η2$⁠, and $η3$ in Eq. (3) to crater function moments. We begin by introducing the dimensionless quantities $u~≡u/a$⁠, $x~≡x/a$ and $y~≡y/a$⁠. Equation (4) becomes

where $ϵ≡a/R$ is small. We also let $X~≡ϵx~$ and $Y~≡ϵy~$⁠, so that

Equation (7) becomes

We will work consistently to second order in $ϵ$⁠. To that order,

or, in terms of the original, unscaled variables

If the solid is reflected about the $x−z$ plane, the sputter yield does not change. Therefore, $M$ must be invariant under the transformation $y→−y$⁠. It follows that $M2=M5=M7=M1,2=M2,3=0$⁠. Equation (A5) reduces to

For later use, we consider the sputter yield $Y′=M′/Ω$ obtained by changing the azimuthal angle of the ion beam from zero to $ϕ$⁠. Let the coordinate axes obtained by rotating the axes $x$⁠, $y$⁠, and $z$ through the angle $ϕ$ about the $z$ axis be the $x′$⁠, $y′$⁠, and $z′$ axes. Then,

where

and

So far, the results we have obtained apply for all $θ$⁠. We will now specialize in the case of normal incidence, $θ=0$⁠. In this case, if the solid is reflected about the $y−z$ plane, the sputter yield does not change, and so $M$ must be invariant under the transformation $x→−x$⁠. As a consequence, $M4=M6=0$⁠. Moreover, by symmetry, $M3=M1$ and $M1,1=M3,3$⁠. Equation (A6) reduces to

As $θ→0$⁠, $M−M′$ tends to zero and $M$ and $M′$ coincide for $θ=0$⁠. Therefore,

We now take the partial derivative of Eq. (A11) with respect to $uxx$ twice and set $uxx=uxy=uyy=uxxx=uxxy=uxyy=uyyy=…=0$⁠. Making use of Eqs. (A8) and (A9), we find that

This shows that

This identity allows us to recast Eq. (A10) as

Let $H$ and $K$ denote the mean and Gaussian curvature of the surface at the origin, respectively. Because $ux=uy=0$ there, we have $2H=uxx+uyy$ and $K=uxxuyy−uxy2$⁠, where the partial derivatives are evaluated at $x=y=0$⁠. Equation (A14) now yields

Dividing Eq. (A15) by $M0=ΩY∞$⁠, we obtain Eq. (3) and the expressions (13)–(15), which relate the $η1$⁠, $η2$⁠, and $η3$ to moments of the crater function.

In deriving Eq. (3), we took $u(x,y)$ to be single-valued and to be defined for all $x$ and $y$⁠. Let us now consider sputtering from a solid object like a sphere that has a double-valued height function $u(x,y)$ that is defined only on a domain $D$ that is a subset of the $x−y$ plane. We will compute the yield from backward sputtering; the forward sputtering yield is very small for $R≫a$ in any event. The height function $u(x,y)$ that describes the upper surface of the object is single-valued and is defined only on the domain $D$⁠. As before, we assume that Eq. (4) is satisfied. We will also take the boundary of $D$ to be given by an equation of the form $s(x/R,y/R)=0$⁠. For the upper hemisphere of a solid sphere of radius $R$⁠, for example,

and $D={(x,y)|x2+y2≤R2}$⁠. Equation (B1) shows that Eq. (4) is satisfied. In addition, the boundary of $D$ is given by $(x/R)2+(y/R)2−1=0$ and so has the required form.

Equation (A15) remains valid when the domain $D$ is finite, but in evaluating the crater function moments, we must take into account the fact that the integration is over $D$ rather than over all $x$ and $y$⁠. We will estimate the size of the error we incur by taking the integration in the crater function moments to be over all $x$ and $y$ rather than $D$⁠. This error is negligibly small if the crater function falls off with sufficient rapidity away from the point of impact. The estimate will be made for $M0$ for the sake of definiteness. We assume that the disk $D0={(x,y)|x2+y2≤R02}$ lies within $D$⁠, where $R0=cR$ and $c$ is a positive constant of order 1. Let $r≡x2+y2$⁠. We will also assume that

for all $r≥R0$⁠, where $n>4$ and $γ$ is dimensionless and positive. The error is

where $D′$ is the complement of $D$ and $d2x≡dxdy$⁠. Using the inequality (B2), we have

Because $n$ is by assumption greater than 4 and we are only working to order $ϵ2$⁠, the error $E0$ is negligible, as claimed.

For the Sigmund model,

$F(x,y,0;0,0,…)$ tends to zero faster than any power of $r$⁠. This means that a tiny error is incurred by replacing the integration over $D$ in a crater function moment by integration over the entire $x−y$ plane.

We asserted that the forward sputtering yield is very small when $a≪R$⁠, i.e., for $ϵ≪1$⁠. It is possible to be more precise if the Sigmund model of ion sputtering is adopted. We will show that for the Sigmund model, the contribution of forward sputtering to the total sputter yield is smaller than order $ϵ2$ for a spherical object of radius $R$⁠. To see this, note that the yield from forward sputtering $Yf$ is given by Eq. (18) with the NED distribution given by Eq. (16) and with $u=−(R+R2−x2−y2)$ and the domain of integration $r≤R$⁠. Since $u≤−R$ on the lower hemisphere and the area of the hemisphere is $2πR2$⁠,

This shows that $Yf$ is smaller than $ϵn$ for any finite $n$⁠. In particular, its contribution to the total sputter yield is negligible because we discarded contributions of order $ϵ3$ and higher.