Establishing a Best Practice for SDTrimSP Simulations of Solar Wind Ion Sputtering

Solar wind (SW) ion irradiation on airless bodies can play an important role in altering their surface properties and surrounding exosphere. Much of the ion sputtering data needed for exosphere studies come from binary collision approximation sputtering models such as TRansport of Ions in Matter and its more recent extension, SDTrimSP. These models predict the yield and energy distribution of sputtered atoms, along with the depth of deposition and damage of the substrate, all as a function of the incoming ion type, impact energy, and impact angle. Within SDTrimSP there are several user-specific inputs that have been applied differently in previous SW ion sputtering simulations. These parameters can influence the simulated behavior of both the target and sputtered atoms. Here, we have conducted a sensitivity study into the SDTrimSP parameters in order to determine a best practice for simulating SW ion impacts onto planetary surfaces. We demonstrate that ion sputtering behavior is highly sensitive to several important input parameters including the ion impact angle and energy distribution and the ejected atom surface binding energy. Furthermore, different parameters can still result in similarities in the total sputtering yield, potentially masking large differences in other sputtering-induced behaviors such as the elemental yield, surface concentration, and damage production. Therefore, it is important to consider more than just the overall sputtering behavior when quantifying the importance of different parameters. This study serves to establish a more consistent methodology for simulations of SW-induced ion sputtering on bodies such as Mercury and the Moon, allowing for more accurate comparisons between studies.


Overview
Surface sputtering by solar wind (SW) ion irradiation is an important process for understanding the surface composition and exosphere formation of bodies such as the Moon, Mercury, and asteroids. Bodies with neither a significant atmosphere nor an intrinsic magnetic field (e.g., the Moon) are directly exposed to SW ions at fluxes that depend on their orbital position and solar activity (Poppe et al. 2018;Tucker et al. 2021). On Mercury, which has a weak magnetic field, SW ions are deflected toward the magnetic poles and impact the surface at high-latitude dayside cusps (Raines et al. 2013(Raines et al. , 2015(Raines et al. , 2022, while at mid latitudes ions impact on the nightside through magnetotail reconnection (Fatemi et al. 2020). In addition to SW ion sputtering, micrometeorite impact vaporization as well as photon-and electron-stimulated desorption can contribute to exosphere formation. The relative contributions of these various processes remains contested for many observed exospheric species (Bida et al. 2000;Burger et al. 2010;Wurz et al. 2010;Schmidt et al. 2012;Tenishev et al. 2013;Gamborino et al. 2019;Killen et al. 2022). Better constraints on the relative contributions of each are needed to interpret ground-and space-based observations of exospheric compositions and densities. Our focus here is on the SW ion sputtering needs of the community. Due to the cost and complexity of sputtering experiments relevant for planetary science, laboratory studies can only provide a small subset of the needed ion sputtering yields, compositional changes, and desorption energy measurements. Hence, theoretical sputtering models are widely used to study the behavior of the incoming ions, impacted surfaces, and sputtered atoms.

Computational Approaches to Modeling Sputtering
Most computational sputtering models use the binary collision approximation (BCA; Eckstein & Urbassek 2007). BCA models predict the yield and energy distribution of the sputtered atoms, as well as the depth of penetration by the primary ion and the substrate damage caused by the collision cascade. Commonly used BCA sputtering models include TRansport of Ions in Matter (TRIM; Ziegler & Biersack 1985) and its more recent extension, SDTrimSP (Mutzke et al. 2019). This program can be run in static (S) or dynamical (D) modes (the latter tracking compositional changes as a function of fluence) using either serial (S) or parallel (P) processing. SDTrimSP allows for variation of the incoming ion type, impact energy, and impact angle. Several previous studies have demonstrated SDTrimSPʼs increased accuracy relative to TRIM in predicting the yield and angular distribution for sputtered atoms in the keV and sub-keV ranges for different ion and single-component target combinations (Hofsäss et al. 2014;Mutzke et al. 2019). In addition, SDTrimSP has been used to predict the role of SW-induced ion sputtering from different planetary-science-relevant multicomponent silicates and minerals (Schaible et al. 2017;Szabo et al. 2020a;.
Within SDTrimSP there are several user-specific inputs that have been applied differently in previous SW ion sputtering simulations (Schaible et al. 2017;Szabo et al. 2020a). For example, simulations can be run using either static or dynamic modes. Static simulations consider impacts into an identical fresh surface for every incident ion and are often used for simplicity or for monoatomic solids. In contrast, dynamic simulations track the sputtered atoms and adjust the target composition accordingly (including as a function of depth). Despite these differences, previous simulations of the total ion sputtering yield did not observe a significant difference in the total sputtering yield between static and dynamic ion sputtering simulations (Szabo et al. 2020a(Szabo et al. , 2020b. However, these studies did not consider the effect of the simulation mode on the elemental composition of the yield or of the target surface. Another fundamental parameter for all BCA models is the surface binding energy (SBE) of atoms in the impacted substrate (Kelly 1986;Behrisch & Eckstein 2007). The SBE has been shown to significantly affect the yield, surface composition, and energy distribution of sputtered atoms (Thompson 1968;Sigmund 1981;Yamamura & Tawara 1996;Behrisch & Eckstein 2007;). This energy distribution is particularly important for exosphere formation models as it influences the proportion of ejecta exceeding the escape energy of the body, the altitude distribution of the ejected atoms, and their corresponding observability . For monoelemental substrates, good agreement with experimental sputtering yields can be achieved by approximating the SBE as the elemental sublimation energy, equal to the cohesive energy for the ground state of the individual atoms in the substrate (Gschneidner 1964;Behrisch & Eckstein 2007). However, SBE values are not well constrained for multielemental substrates, with the SBE of each element often still assumed to be its monoelemental cohesive energy (Mutzke et al. 2019). This implies that the SBE is independent of any bonds formed with other elements, an approximation that has been shown to be unphysical . Therefore, the use of computational parameters measured in monoelemental substrates are likely not correct for multielemental compounds. When experimental sputtering yield data are available for a compound, the SBE can also be fit to match total experimental mass loss. Laboratory data from some simple oxide targets (e.g., SiO 2 and Al 2 O 3 ) are available for comparison with ion sputtering simulations (Eckstein 2007). However, comparatively little work has been done to experimentally measure ion sputtering yields from complex materials relevant to planetary surfaces.

Approximations in Previous Binary Collision Approximation Models
Previous BCA SW simulations of ion sputtering from silicates and minerals relevant to planetary science have often focused on quantifying the SBE of oxygen (O), as it has the largest atomic abundance in many compounds and is therefore expected to make up the largest percentage of the sputtering yield (Schaible et al. 2017;Szabo et al. 2018Szabo et al. , 2020aSzabo et al. , 2020b. The default O SBE recommended by SDTrimSP is 1 eV, based on comparison to Ta 2 O 5 experimental measurements (Hunt & Seah 1983). However, this value is considerably lower than the O SBE values used for other minerals more relevant for planetary surfaces. Schaible et al. (2017) used dynamic BCA simulations to study the sputtering of SiO 2 and Al 2 O 3 and tested O SBE values of 1-3 eV, based on the recommended value within SDTrimSP and on SBE values reported for other monoelemental species. Their findings indicated good agreement with experiment for total ion sputtering and that this total yield did not vary significantly with SBE in this range. However, they did not report the elemental sputtering yield of O, which one would expect to vary by a factor of ∼3 based on the theories of Yamamura & Tawara (1996) and Eckstein (2007). Szabo et al. (2018) conducted dynamic simulations of Ar + and H 2 + impacts onto wollastonite (CaSiO 3 ) for comparison to experimental data. The compound SBE was computed by averaging the tabulated monoelemental cohesive energies for each component (2.6 eV for O). Good agreement with experiment was observed for the total sputtering yield when the SBE of each component observed in the X-ray photoelectron spectroscopy analysis, including the observed carbon (C) contaminant, was included. However, a more recent analysis suggested that C was only at the surface and not in the bulk (Szabo et al. 2020a). When C is not included in the compound SBE, an O SBE of 6.5 eV must be used to obtain good agreement for simulations of He + on wollastonite and iron (Fe)-rich pyroxenes to experimentally measured total yields (Biber et al. 2020;Szabo et al. 2020aSzabo et al. , 2020b. These studies for O indicate that the field has yet to converge on either the importance of the SBE or the appropriate SBE values to use for silicates relevant to planetary science. More research is needed to understand mineral-specific SBEs along with their potential effect on the entire sputtering behavior beyond just total ion sputtering yields. Several additional approximations are often made when simulating the composition, energy, and impact angle of the incoming SW ions. While the SW composition is typically simulated as consisting of about 96% protons (H + ) and 4% helium ions (He ++ ), it also contains minor heavier components including C, N, O, Ne, Si, Ar, and Fe (Bochsler 2007). These secondary components are often ignored for the sake of simplicity, but recent laboratory experiments and TRIM simulations have suggested that this approximation may underestimate the true sputtering yield and damage production in the substrate (Barghouty et al. 2011;Hijazi et al. 2014Hijazi et al. , 2017Nénon & Poppe 2020;Killen et al. 2022). The SW ion energy is typically modeled assuming an energy of 1 keV amu −1 . However, the measured energy distributions of the slow SW (Poppe et al. 2018) has a mean of ∼940 eV, a FWHM of ∼900 eV, and an extended high-energy tail, as is shown in Figure 1. Finally, the impact angle of the SW ions is poorly known and is instead likely a distribution of impact angles. Previous research has shown that ion sputtering yields are highly dependent on the incidence angle (Behrisch & Eckstein 2007;Szabo et al. 2020b), and the SW impact angle is often approximated as either normal (Nénon & Poppe 2020;Killen et al. 2022) or 45°angle relative to a flat surface (Jäggi et al. 2021). In addition, because the surfaces of airless bodies are not flat but composed of individual grains with varied size distributions, the incoming SW ions impact the surface at many different relative angles, further affecting the yield (Wehner et al. 1963;Szabo et al. 2018Szabo et al. , 2022aBiber et al. 2020Biber et al. , 2022. This has been modeled as a parallel flow of incident ions, which is equivalent to a cosine distribution of impact angles onto a flat surface (see, e.g., Cassidy & Johnson 2005). Further research is needed to quantify the effects and validity of these different incoming SW ion approximations on the sputtering behavior.

Purpose
Here, we have investigated the impact of the different approximations for each simulation parameter discussed above on the subsequent predicted sputtering behavior to determine how to best simulate SW ion sputtering. The parameters varied and the values used are shown in Table 1. Based on these options, we have conducted a sensitivity study into the SDTrimSP parameters in order to determine a best-practice for simulating SW ion impacts onto planetary surfaces. The results can be used to establish a more consistent methodology for simulations of SW-induced ion sputtering on bodies such as Mercury and the Moon. The remainder of this paper is organized as follows. Section 2 presents the simulation approach and parameters considered to quantify sensitivity in SW-induced ion sputtering using SDTrimSP. Section 3 presents the results from these simulations. Using these results, Section 4 concludes with a summary of best-practice recommendations to help standardize BCA SW simulations.

Methodology
The BCA model SDTrimSP was used to conduct the parameter sensitivity study, summarized in Table 1, and produce best-practice recommendations for the space-weathering community. We note that SDTrimSP accounts for collisional (i.e., nuclear) and electronic losses during the collision process. However, SDTrimSP does not account for "potential sputtering," which occurs in insulating targets due to the potential energy of multiply charged ions. Previous research has developed models to account for ion sputtering yield increases due to the combination of collisional and potential sputtering (Barghouty et al. 2011;Meyer et al. 2011;Hijazi et al. 2014Hijazi et al. , 2017Szabo et al. 2020a). We also do not include electron sputtering (i.e., sputtering driven by electrons; Assmann et al. 2007) or electronic sputtering (i.e., sputtering driven by ion-generated repulsive electronic states in the atoms in the system; Johnson 2013). Future work is needed to incorporate potential and electronic sputtering effects into BCA models directly.
For the target, we selected the mineral anorthite (CaAl 2 Si 2 O 8 ), a plagioclase feldspar endmember that is considered to be abundant on the Moon (Rickman & Street 2008) and Mercury (Domingue et al. 2014;Nittler et al. 2018). The elemental composition of each component in the target was set to the initial bulk elemental ratio. Following Möller & Posselt (2001) and Szabo et al. (2020a), the elemental mass density of O in the anorthite was then modified to match the experimentally measured bulk mass density of anorthite of 2.73 g cm −3 . 7 The anorthite substrate was then impacted with either hydrogen (H) or helium (He) at a prescribed energy, angle of incidence, and fluence. Diffusion and outgassing from the target was not considered in SDTrimSP due to a lack of experimentally quantified diffusion coefficients (Szabo et al. 2018(Szabo et al. , 2020a. We first quantified the role of the O SBE and static versus dynamic simulations for SW ions of 1 keV H + impacting the surface at normal incidence (a polar angle of 0°). For the O, we  simulated SBEs of 1, 2.6, and 6.5 eV, capturing the range of previously used values. As recommended by the SDTrimSP manual, SBEs for all other elements in the compound were left at their monoelemental cohesive energies of E coh = 2.39, 3.42, and 4.72 eV for Ca, Al, and Si, respectively. We also considered two settings for calculating compound SBEs, referred to as ISBV 1 and ISBV 2, within SDTrimSP. ISBV 1 uses the SBEs of each element in the compound and does not calculate an overall compound SBE. ISBV 2 uses the SBEs of each element of the compound to calculate an averaged overall compound SBE using Here, q x are the number fractions of O, Ca, Al, and Si, respectively, in anorthite. Simulations were conducted up to a fluence of 1-2 × 10 18 atoms cm −2 , ensuring steady state had been reached. We then quantified the effect of the impact angle on the ion sputtering behavior. In our work, three different impact-angle scenarios were considered: impacts normal to the surface, impacts at 45°to the surface, and a cosine distribution of incidence angles.
Next, we quantified the importance of the energy and composition of the incoming SW on the ion sputtering behavior. For the SW impact energy, we compared the 1 keV amu −1 approximation to the slow SW energy distribution (Poppe et al. 2018), from which the incident ion energies were sampled directly in SDTrimSP. For the SW composition, we compared the approximation of 100% H + versus 96% H + + 4% He ++ . It is important to note that SW compositions can also be dependent on the origin of the SW. For example, fast SW, slow SW, coronal mass ejections, and solar energetic particles can all have different compositions and relative concentrations. Such level of detail will be considered in a future study.
For all parameters considered we quantified their effect on (i) the ion sputtering yield, defined as the number of atoms removed from the substrate per impacting ion; (ii) the surface composition, defined as the ratio of each element type in the first 0.5 nm of the substrate; (iii) the damage, defined as the number of vacancies produced at the end of the simulation as a function of depth; and (iv) the energy distribution of the sputtered atoms, defined as the normalized count per emission energy of sputtered atoms. Twenty-two different simulations were conducted in total.

Effect of the Oxygen Surface Binding Energy
The O SBE has a significant effect on the total sputtering yield. This can be seen in Table 2, which shows the total yields for 1 keV H + impacting an anorthite surface at normal incidence as a function of O SBE, the ISBV method, and static versus dynamic simulations. For a static simulation using ISBV 1, the total yield decreases by a factor of ∼5.6 for an increase in the O SBE from 1 to 6.5 eV. For a dynamic simulation using ISBV 1, the decrease is only a factor of ∼2.2. There are also differences between the static and dynamic simulations conducted using the same O SBE. At lower-O SBEs there is a larger difference between the static and dynamic simulations than at higher SBEs. This is because O is more readily sputtered at low SBE and becomes depleted in the surface region in the dynamic simulations, while the higher SBE elements remain. No such depletion is observed in the composition of the static simulations, where the composition is reset before every incident ion and the yields for low-SBE simulations are correspondingly larger.
There are also differences observed between the ISBV 1 and ISBV 2 approaches to the compound SBE. Based on Equation (1) for ISBV 2 and taking default SBE values for Ca, Al, and Si, the anorthite compound SBEs are 2.1, 3.0, and 5.4 eV for O SBEs of 1, 2.6, and 6.5 eV, respectively. Therefore, the compound SBE for an O SBE of 1 eV is higher for the ISBV 2 case compared to the ISBV 1 case. This results in the largest difference in the total yields between the two ISBVs for both static and dynamic cases occurring at an SBE of 1 eV. Additionally, in the dynamic case the bulk concentration of each element changes with irradiation up to steady state for ISBV 1, whereas in the dynamic case for ISBV 2 the elemental bulk concentrations remain constant with fluence. This is a direct result of using unique elemental SBEs for ISBV 1 compared to the same weighted-average SBE for all elements for ISBV 2. Moreover, research has suggested that the SBE for a given element in a compound is an element-specific value and should not be equal for all atom types being sputtered from the compound (Kubart et al. 2010;). Therefore, approximations using ISBV 2 are expected to oversimplify this sputtering process, and for the rest of this work the remaining simulations will consider dynamic substrates using only the ISBV 1 case.
The influence of the O SBE can also be seen in the composition of the sputtering yield during dynamic simulations, i.e., the individual elemental yields of the substrate and in the changing surface composition as a function of incident ion fluence (Table 3). For example, at an O SBE of 1 eV the total sputtering yield is decreased by a factor of 2.6 for dynamic simulations as compared to static simulations. However, the yield of O is decreased by a factor of 4.2. Therefore, when assessing the effect of different simulation parameters, it is necessary to investigate the composition of the sputtering yield and not just the total yield. This is distinctly different from what has been done in previous comparisons of modeling to experiment, which often track the total mass change in the target when assessing agreement (Schaible et al. 2017;Szabo et al. 2018Szabo et al. , 2020a. As shown in Figure 2, the change in surface composition for all elements is dependent on the O SBE. For an O SBE of 1 eV, initially strong preferential sputtering of O is The large depletion in surface O predicted for an SBE of 1 eV has not been observed in SW irradiation experiments and surface characterization of space-weathered samples. In irradiation experiments using 1 and 100 keV H, Burke et al. (2011) observed 60% more O loss relative to Al for anorthite, but did not directly quantify surface concentrations. In irradiation experiments using 4 keV He + , Dukes and colleagues observed only a minor reduction in O surface concentration (∼2%-5%) for Apollo soils (Dukes & Baragiola 2015) and a slight increase in O surface concentration for olivine (Dukes et al. 1999). In irradiation experiments using 1 keV amu −1 H and He, Laczniak et al. (2021) observed only minimal changes in O surface concentration for Murchison chondrite material, but did not directly quantify surface concentrations. More recently, Chaves & Thompson (2022) quantified the elemental composition as a function of depth for olivine from the spaceweathered asteroid Itokawa. The composition of O was reduced from ∼0.65 at depth to ∼0.5 at the surface, agreeing better with our simulations using an O SBE of 2.6 eV. Based on the lack of a significant O reduction in these previous surface concentration measurements for several different silicates,it appears unlikely that the O SBE could be 1 eV. However, it is unknown if a process such has radiation-enhanced diffusion may replenish the sputtered surface O in those experiments, and thus making these low SBEs possible (Quadery et al. 2015;Chaves & Thompson 2022;Christoph et al. 2022). Here, we consider the 6.5 eV SBE to be the most likely, but for the sake of completeness we include O SBE values of 2.6 subsequent simulations.
Our simulations demonstrate the degree to which ion sputtering yields depend on the choice of the SBE values. Research directly quantifying elemental SBEs from compounds and silicates is limited. Kubart et al. (2010) used a combination of BCA simulations and experimental results to estimate the O SBE from a series of metal oxides. The O SBEs for the oxides were significantly higher than those for the monoelemental cohesive energy, ranging from 6-9 eV depending on the oxide studied. More recent research on SBEs has directly quantified mineral-specific SBEs using molecular dynamics (MD). For example,  used MD to demonstrate that the Na SBE was mineral specific and ranged between 2.6-8.4 eV. Future research is needed to further quantify mineral-specific SBEs for each element type in substrates relevant to planetary science and incorporate them into SDTrimSP to increase simulation accuracy in predicting sputtering yields, preferential sputtering, and ejecta energy distribution.

Effect of Impact Angle
Varying the ion incidence angle also significantly affected the ion sputtering behavior, as indicated in cases 1-3 of Table 4 for dynamic simulations. For irradiation at oblique incidence angles (either 45°or a cosine distribution), both the elemental and total sputtering yields increased compared to normal incidence. For all SBEs, the elemental and total yields are increased for a cosine distribution as compared to 45°. For example, for impacts of 1 keV amu −1 H + and an O SBE of 6.5 eV, the elemental (total) sputtering yields increase by a factor 2.6 and 4.4 (2.5 and 3.5) as compared to normal incidence for the 45°and cosine distribution simulations, respectively.
In addition to the yields, it is also important to consider the role of incidence angle on the elemental surface composition and damage as a function of depth. Figure 3 shows the effect of incidence angle can also be seen in the fluence-dependent surface compositions (Figures 3(A) and (B)) and damage (as measured by the vacancy distribution) as a function of depth (Figure 4(A)). The O surface concentration is reduced with increasing fluence for all incidence angles at an O SBE of 2.6 eV. There is similar behavior between the normal and 45°i ncidence cases, whereas for the cosine distribution case the O surface concentration is reduced faster but plateaus sooner. For an O SBE of 6.5 eV, for the two oblique cases there is a relative increase in the O surface concentration compared to normal incidence. At an O SBE of 6.5 eV, the fraction of O sputtering to the total yield at low fluences is lower than the proportion of O in the anorthite mineral bulk. Therefore, at an O SBE of 6.5 eV the fraction of O at the surface is enhanced relative to the bulk. Moving on to the damage (Figure 4), there was also a decrease in the peak depth and an increase in the peak of the damage distribution. For example, at normal incidence the damage peaks at a depth of ∼90 Å and a vacancy count of ∼1.4/ion. In contrast, at 45°incidence the damage peaks at a depth of ∼30 Å and a vacancy count of ∼1.7/ion. For normal incidence, more particles penetrate deeper into the target, leading to a deeper depth for the peak in damage. For oblique incidence, more projectiles are deposited at a shallower depth, thus concentrating the damage in the near-surface region. There is better agreement in the average vacancy per ion count over the entire depth (15.1, 13.3, 14.0 atoms per ion for normal, cosine, and 45 degree impacts, respectively). The damage was not dependent on the O SBE, which is in line with the SBE only coming into play when energy from the collision cascade reaches a surface atom. Note. The simulations were conducted using 1 keV H + impacting at normal incidence. For the dynamic simulations, the ion sputtering yields were sampled at steady state.
Because the cosine distribution of incidence angles approximates the variation in impact angles and corresponding yields expected for a spherical grain surface, it is considered the best approximation currently available within SDTrimSP for SW interactions with regolith. More recent research has also shown that an upcoming 3D version of SDTrimSP could provide another option to capture the granular surface profile (Szabo et al. 2022b). Previous research has suggested that ion sputtering of regolith due to normally incident SW can be distinctly different than ion sputtering from a   Note. The simulations were conducted using either H + or H + + He ++ impacting at a cosine distribution of impact angles and ISBV 1. The ion sputtering yields were sampled at steady state.
perfectly flat surface or an ensemble of uniform spherical grains (Cassidy & Johnson 2005;Kulchitsky et al. 2018;Sarantos & Tsavachidis 2021). As shown in the present study, when considering sputtering from a spherical-shaped grain there is an increase in the total yield as compared to a flat surface due to the influence of oblique impacts in the cosine distribution. This is due to the different incidence angles considered in the cosine distribution. However, the close packing of these grains can also potentially lead to sputtered atoms being adsorbed onto neighboring grains during emission. Cassidy & Johnson (2005) developed a model to account for the effect of a granular regolith structure on the total ion sputtering yield. They demonstrated that for some sizes 74% of sputtered products from an ensemble of uniform grains are retained by neighboring grains when 100% sticking is assumed. However, this reduction factor is dependent on other factors such as sticking coefficient, roughness, grain size, and substrate porosity. Future research is needed to better understand how grain packing can influence the ion sputtering behavior.

Effect of H + Impact Energy
The effect of using a 1 keV amu −1 impact energy approximation versus incorporating the slow SW distribution of impact energies is shown for a cosine distribution of impact angles in cases 3 and 4 of Table 4, respectively. We only consider a cosine distribution of impact angles as it is the most realistic for regolith and dynamics simulations. For the ion sputtering yield, there are only minor differences between the two incidence energy approaches, with the slow SW impact energy distribution yield being consistently lower (3%-5%) than the 1 keV amu −1 approximation. For each case, the proportion of O in the total sputtering yield is consistent. The effect of incidence energy approach can also be seen in the surface composition (Figure 3(C)) and the damage production as a function of depth (Figure 4(B)). There is good agreement in the surface composition between the 1 keV amu −1 approximation and the slow SW energy distribution. Concerning damage, at intermediate depths (110-210 Å) there is 40% less damage when the slow SW energy . Simulations for vacancies/ion as a function of depth for (A) different incidence angle cases using 1 keV H + impacts, (B) different H + impact energy cases using a cosine distribution of incidence angles, and (C) different SW compositions using 1 keV amu −1 impacts with a cosine distribution of incidence angles. Vacancy counts were independent of O SBE (see text). distribution is used. This is because the slow SW energy distribution peaks at an energy lower than 1 keV ( Figure 1). However, the damage extends deeper (>250 Å) due to the highenergy tail in the distribution. More recently, damage has been shown to be a potentially important indicator for the exposure duration and age of exposed samples (Poppe et al. 2018;Keller et al. 2021). Therefore, when damage is being simulated, it is important to use the appropriate ion impact energy distribution of the incoming SW. Future laboratory simulations that approximate the SW as a single ion energy should be aware of the discrepancies this approximation can introduce.

Effect of Solar Wind Composition
The effect of the SW composition on the elemental and total ion sputtering yield is shown in case 5 of Table 4 for a cosine distribution of impact angles. Accounting for the He ++ component in the SW increases both the total and elemental ion sputtering yields by a factor of 1.2-1.3 as compared to H + impacts only (case 3). Therefore, while He ++ makes up only 4% of the SW, it accounts for over 20% of the ion sputtering yield. When comparing the H + with the H + + He ++ options, the proportion of O in the yield stays consistent. The effects of SW composition on the impacted surface composition as a function of incident radiation fluence and on the damage production as a function of depth are shown in Figures 3(D) and 4(C), respectively. For the two compositions there is excellent agreement in the surface composition of the substrate as a function of fluence, indicating that the presence of He ++ does not significantly change the preferential sputtering of the species present. However, when the SW is 96% H + and 4% He ++ the peak in the damage as a function of depth increases by a factor of 1.25. Therefore, approximating the SW as 100% H + during BCA simulations underestimates both the ion sputtering yields and damage production. In addition to H + + He ++ , the SW contains other heavier ion types in much lower concentrations. Simulations using SRIM/TRIM have estimated that these minor ions can increase the collisional sputtering yield by approximately 7% (Nénon & Poppe 2020). It is important to note that SW compositions can also be dependent on the origin of the SW (i.e., fast SW, slow SW, coronal mass ejections, solar energetic particles, etc.; Gloeckler et al. 1999;Von Steiger et al. 2000;Killen et al. 2012).

Energy Distribution of Sputtered Oxygen
It is also important to consider the energy distribution of the sputtered atoms. This energy distribution can be used to help determine the proportion of ejecta exceeding the escape energy of the body and the altitude distribution of the ejected atoms, and thus quantify the contribution of ion sputtering to the exosphere . Figure 5 shows the energy distribution of sputtered O for different SBEs and impact-angle cases. As predicted by Thompson theory, the energy distribution of sputtered O is highly dependent on the O SBE. In all cases the peak in the distribution occurs at ∼E b /2. For larger SBEs, the peak shifts to higher energies and the width of the distribution increases. Increasing the SBE by a factor of 2.5 (from 2.6 to 6.5 eV) increases the peak of the energy distribution by a factor of 2.5 and the FWHM by a factor of ∼6. In addition, there is good agreement in the energy distributions for static versus dynamic simulations along with normal versus a cosine distribution of impact angles. However, because SDTrimSP is developed from binary collision theory it can only predict Thompson-like energy distributions. Therefore, further research using MD and experimental measurements are needed to investigate departures from BCA.

Summary and Best-practice Recommendation
We have conducted a series of BCA simulations of SWinduced ion sputtering using SDTrimSP to quantify the sensitivity of sputtering behavior to the various inputs. Previous simulations of SW ion sputtering have used a range of simulation parameters and often consider only the total ion sputtering yield when modeling ion sputtering behavior (Schaible et al. 2017;Szabo et al. 2018Szabo et al. , 2020aSzabo et al. , 2020b. We have demonstrated that ion sputtering behavior can be highly sensitive to many important input parameters including SBE and impact angle. Furthermore, it is important to consider the elemental ion sputtering yield, surface concentration, and damage production, in addition to the total yield, when quantifying the importance of different parameters. Based on these findings, we conclude with the following best-practice recommendation for SDTrimSP simulations of SW-induced ion sputtering: Figure 5. Energy distribution of sputtered O from anorthite as a function of emission energy for different simulation setups. Simulations were conducted using 1 keV H + impacts. The jaggedness of the curves is due to the Monte Carlo nature of the simulations.