Analytical model for the sputtering of rough surfaces

The sputtering yields of solids under ion bombardment are highly sensitive to the roughness of their surfaces. Understanding how sputtering is exactly affected by different surface morphologies is relevant especially for plasma-wall interaction in fusion reactors and space weathering of planetary surfaces. We present an analytical model that allows to calculate sputtering yields of random gaussian rough surfaces under arbitrary angles of incidence, taking into account local incidence angles, shadowing and redeposition of sputtered materials. Sputtering yields of a rough surface can then be calculated with the sputtering yield’s dependence on the ion incidence angle for a flat surface and a single statistical parameter, which characterizes the surface roughness. The model supports previous findings that the mean surface inclination angle δm is a well-suited parameter to describe the sputtering behavior of rough surfaces. Comparisons of the results to previous experiments and numerical simulations for various cases are presented, showing that the model allows to quantitatively reproduce sputtering yields of different samples over a wide range of roughness regimes.


Introduction
Bombardment of solids with energetic ions leads to the emission of surface atoms due to energy transfer from the ions during a large number of collisions in the target [1]. This erosion process is called sputtering and it is quantified by the sputtering yield Y, which denotes the mean number of atoms being sputtered per impinging ion. The sputtering yield depends on various parameters, such as atomic masses of the projectiles and the target atoms, the energy of the ions or the surface binding energy [2]. One particularly pronounced dependence is also observed for the angle of incidence of the ion in relation to the surface normal. Under oblique incidence, more collisions happen closer to the surface, which leads to more atoms being sputtered. These processes are generally well understood for ideally flat samples both from fundamental theory [1], as well as for numerical simulations that even allow precise reproduction of dynamic effects for multi-component targets [3,4].
Sputtering is an important technique for established applied methods such as the creation of thin coatings by sputter deposition [5], sample preparation with a Focused Ion Beam (FIB) [6,7] or sample analysis with Secondary Ion Mass Spectrometry (SIMS) [8,9]. Recent efforts in sputtering research focused on understanding how the surface roughness affects the sputtering yields. This has largely been motivated by nuclear fusion research, where ions from the fusion plasma erode reactor walls [10]. Here, roughness due to the sample preparation is of relevance, but especially plasma-wall-interaction-induced roughening can be expected, e.g., in the form of needle structures in EUROFER steels [11] or complex W fuzz structures [12]. The surface structures of plasma-facing components then strongly affect how material is eroded and redeposited at other positions [13,14]. Another application for sputtering research is found for space weathering of airless planetary objects in the solar system: The sputtering by ions mainly from the solar wind significantly contributes to the alteration of planetary surfaces [15] or the formation of thin exospheres around the Moon or Mercury [16]. The surfaces of these planetary bodies are mostly made up of regolith grains, which also represent a significant difference to an ideally flat surface [17,18]. Both fields therefore strongly profit from an understanding of the effects of surface morphology on sputtering yields to correctly assess the ion-induced erosion.
Sputtering of rough surfaces has mostly been investigated by numerical simulations, building upon established 1D-codes based on the Binary Collision Approximation (BCA). For example, Küstner et al. simulated sputtering of rough surfaces by using surface inclination angles from Scanning Tunneling Microscopy (STM) as inputs for 1D-simulations [19,20]. Several different efforts were also based on including surface roughness in BCA codes via fractal geometries [21][22][23][24] or varying density in the near-surface region [25]. Recent developments focused on full 3D-BCA codes such as SDTrimSP-3D [26], TRI3DYN [27] and IM3D [28] or geometric ray-tracing codes such as SPRAY [29]. For application in fusion devices, ERO2.0 also includes impurity transport to take into account the conditions of plasma exposure [30,31]. These simulation programs are capable of reproducing a large number of experimental results [26,29,[31][32][33][34][35][36]. It was shown how roughness can change the dependence of the sputtering yield on the projectile's impact angle and that the sputtering of very rough samples is significantly reduced due to redeposition of sputtered atoms [26,29,32,36]. However, such numerical approaches can require long computation times and make it challenging to deduce fundamental dependencies of the sputtering yield. The latter aspect is especially of interest due to the ongoing discussion regarding the best-suited parameters for surface characterization other than the often used root-mean-square (RMS) roughness σ.
Similar to Küstner et al. [19,20], the studies by Stadlmayr et al. [37] and Arredondo et al. [38] also used the surface inclination distribution for characterizing their experimentally investigated surfaces. Cupak et al. recently showed how the mean surface inclination angle δ m is much better suited to categorize rough surfaces in regard to their sputtering behavior than the RMS roughness σ [29]. Their observations are in line with the work by Li et al. [32], who used the mean surface inclination angle to derive an approximate analytical formula for the sputtering yield under normal incidence, achieving good agreement with their IM3D simulations. However, their analytical approach does not include a full treatment of surface characteristics and does not describe non-normal incidence. Another analytical approach was developed by Makeev and Barabási by modifying Sigmund's theory for both self-affine and rippled surfaces [39,40]. Their model predicts both sputtering increases and decreases due to surface roughness, but it does not include shadowing or redeposition effects. Stepanova et al. also discussed an analytical model for the sputtering of rough surfaces with a similar approach as the model presented in this manuscript [41], but their method only takes into account periodic, idealized surface structures. Comparisons to their Cu sputtering yield measurements indicate that this is not universally applicable.
This manuscript presents an analytical model for the sputtering of random gaussian rough surfaces, which can be solely characterized by the mean surface inclination angle δ m . Based on previous geometrical models aimed at describing the shadowing of light at rough surfaces, expressions for the sputtering yield are derived that include sputtering due to different local angles of incidence, shadowing of surface features and redeposition of sputtered material. Using surface characteristics and the sputtering yields of flat samples, we show how, despite the assumption of a random gaussian rough surface, various experimentally and numerically derived sputter data can be reproduced very well with this model. Given a valid input for the yield's dependence on the angle of incidence for a flat surface, sputtering yields for any projectile-target combination can be calculated.

Surface roughness characteristics
Following the approaches by He et al. [42] and Smith [43] for calculating the scattering and shadowing of light on a rough surface or the work by Li et al. about sputtering simulations of rough surfaces [32], a restriction to random isotropic gaussian rough surfaces is made. The mean height of such a gaussian surface is 0 and its height z follows the probability distribution (see [43]) with the RMS roughness σ. Even for random gaussian surfaces, σ alone is not sufficient to describe the surface structure. Instead additional information is included in the joint distribution functions of surface slopes p = ∂z/∂x and q = ∂z/∂y, which is given by (see [43]): P pq (p, q) dp dq = 1 2πw 2 exp Due to the isotropy of the surface, the standard deviations of the distributions of p and q are described by the single parameter w. It can be calculated by using the autocorrelation function ρ(r), which describes the correlation of heights at different surface points with distance r (see Appendix A). Evaluating the second derivative of the autocorrelation function ρ(r) at r = 0 then gives w 2 = − ρ ′′ (0) [43]. Technical surfaces often show gaussian autocorrelation functions characterized by an autocorrelation length τ and therefore w 2 = 2σ 2 /τ 2 (see Appendix A). In the following, the ratio of RMS roughness and autocorrelation length σ/τ will be used, with larger values of σ/τ representing larger surface roughness, but this remains interchangeable with the more general parameter w.
As described in detail in Appendix A, the distribution of surface inclination angles P θ (θ) dθ can be calculated from the slope distribution P pq (p, q) dp dq: Figure 1 demonstrates how P θ varies for different ratios σ/τ. The distributions are characterized by a significant peak and shift to larger angles for rougher surfaces. The calculated distribution are similar to those determined by Stadlmayr et al. [37] and Arredondo et al. [38] from Atomic Force Microscopy (AFM) measurements of their investigated samples.
For the mean surface inclination angle δ m in relation to the z-axis, the following expression can be derived: The dependence of δ m on σ/τ is shown in Fig. 2. The same behavior as in the approximative derivation by Li et al. [32] can be observed, when the profile width of their gaussian surfaces W = 100 nm is used as the autocorrelation length τ. δ m is a monotonously increasing, unique function of the roughness parameter σ/τ. δ m thus contains the same information as σ/τ and therefore fully describes a random gaussian rough surface. It also projects the possible values of σ/τ to the finite interval [0 ∘ , 90 ∘ ].

Sputtering effects on rough surfaces
After presenting the mathematical descriptions of the rough surfaces that are of interest, their influence on the sputtering yield Y is discussed in the upcoming sections. In particular, following effects are included (see the sketch in Fig. 3): • local incidence angles • shadowing of surface areas • redeposition of sputtered material Secondary sputtering by reflected ions or increased sputtering on sharp edges can play a role for very rough surfaces, but these phenomena are not implemented in the presented model. Similarly the change of the surface structure due to dynamic erosion and redeposition of material is not discussed here.
Overall, the mean sputtering yield Y(α) over the whole surface as a function of the global incidence angle α is the quantity of interest and should include the three above mentioned effects. The expression for Y(α) is of following form: ×F v (p, q, z, α) F sp (p, q, z) P pq (p, q) P z (z) dz dq dp, where Y(θ loc (p, q, α)) represents the sputtering due to the local incidence angle θ loc at a certain surface segment with the slopes p and q.
F v (p, q, z, α) represents a weight factor for a surface segment to be visible to the incoming ions and thus not being shadowed from another surface segment. During ion irradiation under oblique incidence, some surface areas will be more exposed, while others will not be irradiated at all.
F sp (p, q, z) represents the probability for sputtered atoms to escape the surface instead of being redeposited at another feature. Both F v and F sp are dependent on z due to higher features being more exposed and having a larger free solid angle for both incident ions and sputtered atoms. To calculate the mean sputtering yield Y for a rough surface, the product of these elements is then integrated with the distribution functions P pq (p, q) and P z (z).

Sputtering at local incidence angles
The local angle of incidence θ loc is given as the angle between the direction of the ion beam with angle α to the nominal surface normal and the normal vector of the local surface segment with slopes p and q. Due to the isotropy of a gaussian surface, a coordinate system can be chosen where the ion beam is in the xz-plane and then θ loc can be calculated as (see Appendix B) For the general angular dependence of the sputtering yield, the Eckstein fitting formula [44] is used to describe the sputtering yield Y, when a surface is hit under an angle β: The fitting parameters β 0 , b, c and f allow a precise description of the sputtering yield's initial increase with the incidence angle, its maximum and the subsequent decrease for incidence angles towards 90 ∘ . In the integral over p and q in Eqs. (5), (6) for θ loc is inserted into the Eckstein formula to calculate sputtering yields Y(θ loc ) depending on the local angle of incidence.

Shadowing
The shadowing of an ion trajectory hitting a rough surface under a global angle α shall be described by a weight factor F v (p, q, z, α) indicating how many ions hit a surface segment at height z and the slopes p and q. For a similar problem, Eq. (22) by Smith [43] derived the probability of a surface segment to be illuminated by a light beam. Based on this expression, F v (p, q, z, α) can be constructed by combining Eq. (22) from Smith [43] with a normalization using Eq. (24) from Smith [43] because every ion still hits the surface. This then gives with the Heaviside function H and (see [43]) Λ(cot α) represents a measure for the likelihood that a surface segment is shadowed under an angle α. It is 0 for α = 0 and becomes infinite towards α = π/2. In particular, 1/(1 + Λ(cot α)) is the probability that a surface segment with the local normal vector pointing in the direction of the incoming ion beam (angle α) is not shadowed [43].
The Heaviside function in Eq. (8) guarantees that only slopes facing the ion beam can get irradiated. The fraction ensures the normalization of F v and the expression in squared brackets favors surface areas that are higher and therefore less likely to be shadowed by other features. F v is not dependent on q and the dependence on p is restricted to the Heaviside function. Apart from this aspect, the surface slopes in principle do not influence the probability of shadowing, which is in line with the assumption of the height probabilities P z and the slope probabilities P pq being independent of each other.

Redeposition of sputtered material
Going towards larger roughness, the redeposition of sputtered material plays an important role in reducing the sputtering yield. Redeposition has been identified as a key influence on the sputtering yield of rough surfaces [19,20] and has also been investigated in regard to pattern formation [45,46]. This effect is included in a full binary collision approximation (BCA) treatment in SDTrimSP-3D [26] or TRI3DYN [27], but can also be formulated analytically for the regarded isotropic random rough surfaces.
The probability that a recoil atom sputtered under a global polar angle θ r can leave the surface without being redeposited is the same as the probability that the surface point, where an atom is sputtered from, is visible under the same angle. The expressions discussed in the previous section on shadowing can thus be used in a similar fashion, as they were also applied by Smith to model the thermal emission of the Moon [47]. As described in more detail in Appendix C, the total probability F sp that an atom from a surface with slopes p and q is not redeposited can be calculated to result in following expression: with the global polar angle θ r under which the atom is sputtered: For the angular distribution of sputtered particles f sp (θ ′ , ϕ ′ ), we use (see A polar coordinate system of the local surface segment described by the angles θ ′ and ϕ ′ is chosen for the easiest implementation of the distribution of sputtered particles f sp (θ ′ , ϕ ′ ). The Heaviside function restricts possible sputtering directions to the upper hemisphere in the global coordinate system. The expression in square brackets is the same as for the shadowing in the previous section, describing the probability that a surface segment at height z is visible under the angle θ r . For the distribution of sputtered atoms f sp (θ ′ , ϕ ′ ), the cosine distribution represents a first order approximation for the distribution of sputtered atoms that was also applied by Küstner et al. [19,20]. In principle, a more precise with azimuthal dependence can also be implemented.

Evaluation of the dependence of the sputtering yield on roughness
Following the above presented derivation of mathematical expressions for sputtering on local incidence angles θ loc , shadowing F v and redeposition F sp , the averaged sputtering yield Y(α) can be calculated with Eq. (5). The Eckstein sputtering yield formula from Eq. (7) has to be evaluated at θ loc from Eq. (6), F v has to be inserted from Eq. (8) and F sp has to be taken from Eq. (10). This leads to integrals over p, q, z, θ ′ and ϕ ′ over several functions that depend on these parameters. Following the calculations by Smith [43], the integral over z can be calculated analytically because of Then After the integration over z, this expression only includes dependence on roughness parameters in the form of the ratio σ/τ. This is consistent with the purely geometrical effects that are considered within this model and should therefore be independent of any scaling of the surface. The unique description with σ/τ also means that the sputtering behavior of a random gaussian rough surface is also uniquely described by the mean inclination angle δ m , following Eq. (4). This result therefore supports the usage of δ m as a parameter for characterizing the surface roughness in regard to sputtering behavior, as previously discussed by Cupak et al. [29]. For small surface roughness, the sputtering will mostly be influenced by the effect of local incidence angles. However for high δ m values, Fig. 4. The mean fraction of redeposited material 1 − 〈F v F sp 〉 is plotted. For rougher surfaces, redeposition significantly increases, which reduces the sputtering yields for large δ m . Due to different shadowing, redeposition is less important under oblique incidence. The result of redeposition calculations by Diddens and Linz is included for comparison as the dashed black line [46]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) redeposition of sputtered particles will strongly change the total sputtering yield. The mean fraction of redeposited atoms independent of the projectile-target combination can be determined by calculating 1 − 〈F v F sp 〉, with the angle brackets denoting the integral over the product of F v F sp with the distributions P z (z) and P pq (p, q): F v F sp is the probability that a sputtered atom can leave the surface from a chosen surface area and 〈F v F sp 〉 represents the average over the whole surface. The redeposition fraction 1 − 〈F v F sp 〉 is depicted in Fig. 4 for different angles of incidence. The plot clearly shows, how redeposition becomes more important for very rough surfaces. For the redeposition model of Diddens and Linz [46], the authors also calculated the fraction of redeposited atoms and found a similar dependence on the aspect ratio of their investigated surfaces. Figure 4 includes their result for normal incidence with a dot-like surface and a cosine angular distribution of sputtered material and their surface aspect ratio is taken as equal to σ /τ.
Their calculation predicts the same behavior as our model and underlines the suitability of our approach. The fraction of redeposited atoms calculated by Diddens and Linz is slightly dependent on the distribution of sputtered material [46]. As mentioned previously, future studies should also test model results with other distributions of sputtered atoms f sp . Furthermore, Fig. 4 also makes the role of shadowing in the redeposition process evident: For oblique incidence angles, shadowing favors surface positions at larger heights z to be hit by impinging ions. At these points, the solid angle of free paths away from the surface is larger and therefore sputtered atoms are less likely to be redeposited.

Comparison to experiments and numerical simulations
Equation (14)

Sputtering of CaSiO 3 by 2 keV Ar + ions
The mineral wollastonite (CaSiO 3 ) has been used to investigate the sputtering of planetary bodies in the solar system [48][49][50][51], such as the Moon or Mercury, which are continuously exposed to solar wind ions. In the scope of this study, further sputtering yield measurements with a rough wollastonite sample were performed under irradiation with 2 keV Ar + ions. Sputtering yields were determined with a quartz crystal microbalance technique [52], in experiments analogous to those described in [48][49][50]. From AFM images, a mean surface inclination angle δ m = 15.9 ∘ was determined, which was used as an input for the theoretical calculation. For the sputtering yield's angular dependence, the experimental values for flat samples from Szabo et al. [48] were fitted with the Eckstein formula from Eq. (7). With β 0 = π /2 for noble gas projectiles (see [44]), this gives following fit parameters: Y 0 = 21.06, c = 0.75, f = 4.81, b = 1.85.
The blue dots in Fig. 5a show the angular dependence of the measured sputtering yield of the rough wollastonite sample. For composite materials such as wollastonite, only the sputtered mass per impinging ions can be measured with the QCM technique. The yield is therefore given as mass sputtering yield y with units [u/ion] to differentiate it from atomic sputtering yields Y [atoms/ion]. Compared to a flat surface (black dashed line), the yield is increased at normal incidence and decreased at very oblique incidence angles around 70 ∘ . The model calculation with δ m = 15.9 ∘ is included as the solid red line, which agrees very well with the experimental data. The inset in Fig. 5a shows that the surface inclination angle distribution P θ deduced from AFM measurements (blue) differs from the distribution P θ described by Eq. 3 (red). However, this demonstrates that the model can also be used reasonably well to describe the sputtering of non-gaussian surfaces. Further calculations with other δ m values are depicted as red dashed lines, showing how the yield is predicted for even rougher surfaces. Further development of the observed tendencies occur, before the sputtering yield decreases at all incidence angles at δ m towards 80 ∘ .
The behavior of the sputtering yield as the ratio y rough /y flat is shown in Fig. 5b. This quantity represents the change of the sputtering yield for a rough surface compared to the yield for a flat surface. The x-axis describes different mean inclination angles δ m and the y-axis gives different angles of incidence α. For the case of 2 keV Ar + ions sputtering wollastonite, a significant increase under normal incidence up to about twice the initial value is predicted by our model. In contrast, sputtering yields under oblique incidence immediately decrease for rougher surfaces, especially at α ≈ 70 ∘ , where the yield's angular dependence has its maximum. A similar variation of both sputtering yield increases and decreases for different incidence angles has also been found by Makeev and Barabási for rippled surfaces [40]. The decrease is more prominent in our model due to redeposition becoming the dominant effect for very rough surfaces.

Sputtering of rough Be surfaces
Küstner et al. applied a simulation approach to calculate the angular dependence of the sputtering yield of rough Be surfaces based on local incidence angles and redeposition [19,20]. STM images of a size of several µm were used as an input for the surface topography, from which the distribution of incidence angles and a redeposition factor were determined. Based on the incidence angle distribution, sputtering yields with TRIM.SP were calculated. For rough Be surfaces, sputtering yields were measured for the irradiation with Be + and He + ions, leading to a good agreement with the simulation for all presented cases.
For comparisons of our model with the results by Küstner et al., mean surface inclination angles δ m were determined from the given surface inclination angle distributions. Inclination angle values of δ m = 29 ∘ for the sample roughened by prolonged ion bombardment and δ m = 37 ∘ for the unpolished sample from [20] were calculated, even though the limited resolution of the given surface inclination angle distributions ( Fig. 4c and e in [20]) leads to some uncertainties for these parameters. It also hinders the assessment of how well these surfaces correspond to a gaussian rough surface, but the surface inclination angle distributions derived from STM images in [20] are somewhat broader than those described by Eq. (3). Using the derived δ m values and the angular dependence of the sputtering yields of flat Be samples from TRIM.SP, which are also given in [20], sputtering yields of rough Be samples were calculated with our model under different angles of incidence. Figure 6 shows comparisons of the different approaches for the irradiation of Be samples with 3 keV Be + ions (Eckstein fit parameters for a flat surface: β 0 = π/2, Y 0 = 0.31, c = 0.92, f = 3.17, b = 0.68) and 300 eV He + (β 0 = π/2, Y 0 = 0.11, c = 0.89, f = 5.65, b = 2.04).
Experimental results (solid blue dots) and simulation results (open blue dots) from Küstner et al. [20] are compared to our model (red line) for Be + for unpolished samples (δ m = 37 ∘ , Fig. 6a) and ion-beam-roughened samples (δ m = 29 ∘ , Fig. 6b). The TRIM.SP simulation for a flat surface (black dashed line) is included for comparison. In both cases, the deviation from the flat sample is very well reproduced: Due to roughness effects, the yields are increased for normal incidence and significantly decreased for oblique angles of incidence. The almost linearly increasing yield is also reflected in our calculation. Given the uncertainties of the surface characterization inputs, the quantitative agreement especially for the roughened sample in Fig. 6b is very good. For a different case,  Fig. 6c compares the results of 300 eV He + irradiations of the roughened sample (δ m = 29 ∘ , same inputs as for Fig. 6b). Again, the sputtering yields from the model show a similar angular dependence as the data from Küstner et al. [20].

Sputtering of rough W surfaces by 2 keV Ar + ions
In their study on different ways for characterizing the sputtering behavior of rough surfaces, Cupak et al. [29] investigated the sputtering of rough W surfaces in both experiment and modeling with the geometrical SPRAY code. SPRAY simulations take into account similar effects that are included in the present manuscript, with the addition of precise particle distributions from BCA-codes and secondary sputtering effects. Cupak et al. found that SPRAY is well-suited for simulating experimentally derived sputtering yields and subsequently modeled a large number of artificially generated surfaces. Their results revealed very similar sputtering yields for surfaces with the same mean surface inclination angle δ m , showing that this parameter is much better suited for surface characterization than the RMS surface roughness σ.   (blue) of normal incidence irradiations for different roughness parameters δ m very precisely. Deviations arise only above δ m ≈ 60 ∘ , most likely due to the increasing importance of secondary sputtering, which denotes the additional sputtering by reflected projectiles. Our theoretical model neglects this effect and therefore underestimates the total sputtering yield. Similar results are observed for comparing the sputtering yields of varying rough surface under a nominal angle of incidence of 60 ∘ , which is shown in Fig. 7b. The model gives slightly lower results than the SPRAY simulations, but the general dependence of the yield, which continuously decreases with increasing roughness, is very well reflected. Figure 8a presents model results for a variation of mean inclination angles δ m and angles of incidence α for the case of 2 keV Ar + projectiles and a W sample. Similar to Fig. 5b, the ratio of yields for rough and flat targets is depicted here. In contrast to wollastonite, normal incidence sputtering yields of W only slightly increase for rougher surfaces. This is caused by the less pronounced angular dependence of the sputtering yield of W, which can be seen in the Eckstein fitting parameter f (2.804 for W and 4.81 for wollastonite). The effect of different local incidence angles is therefore not very significant in increasing normal incidence sputtering yields for W, while the decrease at oblique incidence is similarly observed for W and wollastonite. In the same fashion, redeposition of sputtered atoms strongly decreases all sputtering yields for very rough surfaces with large δ m values. The highest relative yield increase is observed for angles of incidence α of around 85 ∘ , but this is also influenced by the yield being almost zero at these angles for a flat surface. Figure 8b compares the measured sputtering yields from Cupak et al. [29] with model results. Measurements (dots) for different samples with low roughness (blue, δ m = 7.6 ∘ ), medium roughness (green, δ m = 19.7 ∘ ) and larger roughness (orange, δ m = 36.5 ∘ ) are compared to respective calculations (lines). Yields at 0 ∘ and 60 ∘ are in very good agreement with all experiments, in accordance to the comparison with SPRAY simulations shown in Fig. 7a and b. However, the exact angular dependence for the medium rough sample shows a different behavior in the experiment. The inset in Fig. 8b indicates some differences for the surface inclination angle distribution P θ from the AFM measurement in [29] and P θ from Eq. (3). Though, deviations at a similar δ m value have also been observed for wollastonite in Fig. 5a with much better agreement between model and experiment. Instead, we assume the observed discrepancies are more likely connected to the larger valley-like structures that are visible in the AFM images presented in [29]. SPRAY simulations that take into account the full geometry from AFM images are able to reproduce the whole angular dependence more precisely [29]. Nevertheless, the model predicts the qualitative developments of the sputtering yields for rougher surfaces also for the cases shown in Fig. 8b well.

Conclusions
In this manuscript, a full analytical model for the sputtering of random gaussian rough surfaces is presented. For this class of surfaces, expressions for the surface inclination angle distribution and the mean surface inclination angle δ m can be derived. We show that δ m includes all information on the roughness of gaussian surfaces and is therefore wellsuited as a characterization parameter, especially for surfaces with random roughness. Building on this description, the geometrical aspects of sputtering such as different local incidence angles, shadowing of surface features and redeposition of sputtered atoms can be combined in integral expressions for the sputtering yield of a rough sample. Sputtering yields can then be calculated by numerical integration and the model is able to reproduce the effect of roughness on the sputtering yield for different cases of samples relevant for space-weathering and nuclear fusion research. The model describes how large the roughness-related increase of the sputtering yield at normal incidence can be for different sample materials. Good agreement is also found, when the surfaces deviate from random gaussian rough surfaces. This indicates that the approximation by a gaussian rough surface with the same mean inclination angle δ m is a valid approach for many cases. Some discrepancies in the exact angular dependence can occur if the differences to a gaussian rough surface are too large, but overall the model accurately describes sputtering yield behaviors for different target-projectile combinations and different roughness regimes.
Further improvements of this theoretical approach could be implemented by including sputtering by secondary particles and more precise simulation results will of course be possible with SPRAY [29], SDTrimSP-3D [26] or TRI3DYN [27]. However, our above presented model has been shown to predict sputtering yields of rough surfaces well without the need for large-scale numerical simulations. The sputtering yields can also be calculated with only the mean surface inclination angle δ m being necessary as a single parameter to describe the surface roughness. This advantage of our model makes its application especially favorable for larger calculations on plasma-wall interaction and exosphere formation, where sputtering represents only one of several processes that have to be simulated.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
P pq (p, q) dp dq = 1 2πw 2 exp with the RMS roughness σ and the standard deviation of the slope distributions w. w is given by (see [43]) with the autocorrelation function (see [43]) r and R are vectors in the horizontal xy-plane. Due to the isotropy of the surface, ρ only depends on the distance r = |r|. For a gaussian autocorrelation function ρ(r) = σ 2 exp w can be evaluated as w 2 = 2σ 2 /τ 2 , with τ as the autocorrelation length. Equation (A2) can be rewritten as: Another common function to describe a surface structure is the distribution P θ of local surface inclination angles θ. Its relation to the slopes p and q is sketched in Fig. A1 and can be described with p = tanθ p and q = tanθ q as well as the azimuthal angle φ: The distribution P pq (p, q) dp dq can be transformed to a distribution P θφ (θ, φ) dθ dφ by calculating the Jacobi determinant: with p 2 + q 2 = tan 2 θ. P θφ is then written as It is not dependent on φ, which agrees with the isotropy of a random gaussian rough surface, and integration over φ yields the distribution P θ (θ) of surface inclination angles θ: P θ is shown in Fig. 1 in the main text for different values of σ/τ. For surfaces with larger roughness, the distribution shifts towards higher angles. Additionally, the width of the distribution significantly varies. For intermediate roughness of σ/τ = 0.5, P θ becomes very broad compared to σ /τ = 0.05 or σ/τ = 5.
As described in the main text, the mean inclination angle δ m can be calculated from Eq. (A10):

Appendix B. The local incidence angle
The local angle of incidence θ loc of an ion hitting a segment of a rough surface, can be calculated from the direction vector of the ion beam α and the normal vector n of the local surface segment sketched in Fig. A1. Due to the assumed isotropy of the rough surface, we can choose a coordinate system where the ion beam direction vector α is confined to the xz-plane. α is then given as: n can be calculated as P.S. Szabo et al.

Appendix C. Redeposition of sputtered material
Redeposition of sputtered materials can be calculated similarly to shadowing, but the distribution of sputtered particles f sp has to be taken into account. The visibility probability S for a surface segment under an angle θ r from Smith [43] is given as: The probability that an atom is sputtered in a certain direction follows within a first approximation the distribution f sp (θ ′ ) = cos(θ ′ ) /π, with θ ′ describing the angle between the sputter direction and the local surface normal n. Due to the slopes of the surface, θ ′ does not equal the global polar angle θ r . The probability that a sputtered atom actually escapes the surface and is not redeposited on another slope, has to be calculated by integration over one of these angles. With the distribution of sputtered particles f sp (θ ′ ) being only dependent on the local polar angle, this is best done by defining such a local coordinate system, where the local surface normal vector n determines the z-axis and the global ẑ vector lies in the xz-plane. In this coordinate system, ẑ can be written as . The direction vector of a sputtered atom v sp is simply given by polar and azimuthal angles θ ′ and ϕ ′ : The angle θ r represents the global polar angle of the direction vector and is thus given as the angle between v sp and ẑ: The probability that an atom is sputtered in a direction described by (θ ′ , ϕ ′ ) from a surface with slopes p and q and not redeposited on another surface segment is thus given as the product of the visibility probability S and the angular distribution of sputtered particles f sp : The Heaviside function in S was here replaced. The original Heaviside function H(cot θ r − p) describes surface segments that cannot be illuminated, which is not needed here as all sputtered particles start from a surface segment. Instead H[π/2 − θ r (θ ′ , ϕ ′ )] constrains the angles θ ′ and ϕ ′ so that the direction vector v sp points away from the surface. Particles sputtered in the direction of θ r > π/2 will be redeposited every time independent of surface roughness or slope orientation. The total probability F sp that an atom from a surface with slopes p and q is not redeposited can be calculated with integration over the angles θ ′ and ϕ ′ , where the Jacobi determinant sinθ ′ has to be included: (C6)

Supplementary material
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.surfin.2022.101924. P.S. Szabo et al.