Sputtering of surfaces by ion irradiation: A comparison of molecular dynamics and binary collision approximation models to laboratory measurements

We compare various sputtering simulation methods to experimental results in both the low energy (<1 keV) and high energy ( ≥ 1 keV) impact regimes for argon ions impacting a pure copper substrate at normal incidence. Our results indicate that for high energy impacts, both binary collision approximation (BCA) and molecular dynamics methods can be used to generate reasonable predictions for the yield and energy distribution of the sputtered atoms. We also find reasonable agreement between the theoretical and experimental results down to impact energies of 600 eV. However, at 200 eV impact energies, significant discrepancies appear between the experimental and theoretical ejecta energy distributions in the peak position, the width of the energy distribution, and the magnitude of the high energy tail. These discrepancies appear to arise from the experimental results being only for atoms sputtered normal to the substrate surface, whereas the theoretical results are integrated over all 2 π solid angles above the surface. Using the BCA code SDTrimSP and limiting the results to only atoms emitted within ±15° of the surface normal brings theory and experiment into reasonable agreement. These results suggest that for low energy impacts, the energy distribution of sputtered atoms is highly dependent on the emission angle of the ejecta.


I. INTRODUCTION
Sputtering due to ion irradiation is a process important to a diverse range of fields, such as surface and exosphere evolution of airless planetary bodies [1][2][3] and the production of thin films for optical coatings and electronic devices. 4,5Advances in these fields build on our experimental and theoretical understanding of sputtering.However, laboratory simulations of surface sputtering are expensive and complex to conduct and can provide only a small fraction of the required data.Theoretical calculations are needed to provide the bulk of the data required to move these fields forward.
The most physically realistic simulations use molecular dynamics (MD) methods, 4 which consider all atomic interactions in the system.However, MD calculations produce significant computational loads, thus limiting the spatial size of the simulated volume and the duration of the simulated event.In contrast, binary collision approximation (BCA) models consider sputtering to be a result of binary collision cascades. 5These models can be used to predict the energy distribution and yield of sputtered atoms as a function of incoming ion type, energy, and angle, with only modest computational requirements.For this reason, the majority of theoretical sputtering studies have been conducted using BCA methods. 5However, it is important to benchmark these theoretical methods against more physically realistic MD simulations and experimental results.
BCA models rely on traditional linear collision cascade theory and assume that the impact event can be approximated as a series of independent collisions between pairs of atoms. 5For this reason, BCA models work best at high energies (≥1 keV) where the motion of the incident ion is faster than the vibrational motions of the substrate. 6One of the simplest BCA models, the Thompson distribution, 7,8 was derived to provide an estimate for the energy distribution of the sputtered particles and has been found to be independent of the incoming ion type or energy.Another BCA model is the Transport of Ions in Matter (TRIM) software, 9 which can be used to predict the sputtering yield and energy distribution for a wide range of ion/target combinations.
As impact energy moves below the keV regime, the motion of the incident particle becomes comparable to the motion of the atoms in the substrate and many-body effects must be included. 5herefore, the approximation of a moving atoms colliding with only a single target atom (binary collisions) begins to break down.In this regime, processes become increasingly dependent on the orientation of the substrate relative to the incoming particle.This is commonly referred to as the anisotropic regime.In an effort to account for these low energy effects, SDTrimSP extends the more commonly used TRIM software by incorporating a krypton-carbon type interaction potential, all within the BCA. 5,10The underlying physics of SDTrimSP is summarized by Eckstein. 11 Previous comparisons of SDTrimSP results to experimental data has shown that SDTrimSP can be used to predict the yield and angular distribution of sputtered atoms above 1 keV. 12However, we are unaware of any studies that compare both the predicted yield and energy distribution from SDTrimSP to experimental results for both high and low impact energies.
A key parameter in predicting accurate energy distributions and sputtering yields for BCA models is the surface binding energy (SBE) of the substrate. 10,13The SBE is a user defined value that describes the binding strength of surface atoms to a target.As an example of the importance of the SBE, the sputtered atom energy distribution given by Thompson theory is Here, E is the energy of the sputtered ion and E b is the SBE of the substrate.The peak in this distribution corresponds to half the SBE.Despite the clear importance of SBE in simulating sputtering, its actual meaning and value are not well known for many substrates.For single component substrates, the SBE is often approximated as the cohesive (i.e., sublimation) energy for the atoms in the substrate. 14,15However, different authors have argued, using both MD simulations and experiment, that this approach can underestimate the true SBE by 30%-40%. 13,14,16A 40% underestimation of the SBE translates to a factor of ∼1.7 reduction in the peak energy for the Thompson energy distribution.Furthermore, Stepanova and Dew 10 reviewed published experimental energy distributions for various monoatomic substrates.They found that for a given substrate, the experimental data could be fit to a Thompson distribution, but the different measurements exhibited a range of values for E b .These values averaged approximately to the cohesive energy of the substrate.Stepanova and Dew concluded that E b is better considered a fitting parameter.Therefore, it is unknown whether E b represents a true SBE, a cohesive energy, or is instead fitting parameter.Such uncertainties can have a significant impact in planetary science on the predicted structure and composition for exospheres formed by solar wind ion irradiation (e.g., Mercury or the Moon). 17,18here is no universal approach to estimating the SBE for multicomponent substrates.Given that BCA methods rely on a user defined SBE, this can be a significant source of error for more complex substrates.BCA methods are also limited to amorphous substrates and thus do not take into account the effects of crystallinity or temperature.As such, while we expect the vast majority of studies to employ BCA methods, situations arise where they are less suitable.It is, therefore, important to benchmark BCA methods with more physically realistic methods and experimental results.
MD, in contrast to the assumptions of BCA models, uses an interatomic potential to define all atomic interactions within the system during the collision cascade.A typical MD simulation begins with the known initial positions and velocities for all atoms within the system.The atomic accelerations are then calculated using an interatomic potential that defines the forces between the interacting particles.With the atomic acceleration known, Newton's equations of motion are then used to predict subsequent positions and velocities.However, care must be taken to ensure that the selected potential is adequately parameterized for the ion/target interactions and free surfaces used in the simulations.Because all interactions are captured in the simulation, MD removes the assumption of binary collisions and can be used to simulate both amorphous and crystalline samples.Furthermore, unlike BCA approaches, SBE is not a user input but is instead a function of the atomistic interactions on the surface of the substrate. 19This allows for the simulation of more complex substrates that can account for the effects of pre-existing damage and surface roughness.This is particularly important when simulating multicomponent planetary surfaces where impacts with solar wind ions can lead to sputtering and surface modifications. 20owever, MD simulations require several user specific choices to allow for a more realistic and accurate simulation, such as the interatomic potential, equilibration methods, and boundary conditions.
While significant research exists on MD simulations of sputtering, studies that consider both yield and energy distribution are more limited.The study by Yang et al. 13 is such an example, where they simulated the impact of argon ions with energies of 200-500 eV on a tungsten substrate.Their findings suggested that MD could be used to accurately predict both the yield and energy distribution of sputtered particles for low energy impacts (<1 keV).However, higher energies were not considered and the results were also not compared to experimental values to benchmark the method.As recommended by Urbassek, 19 before MD can be used as an additional sputtering simulation tool, research is needed on its accuracy relative to experimental values.
Building on this need, the purpose of this study is to benchmark MD and BCA simulations with published experimental data to provide a better understanding of the applicability and accuracy of current sputtering simulation methods.Specifically, we consider a copper substrate impacted by argon ions with kinetic energies between 200 and 3000 eV and compare the MD and BCA results to the experimental values for the resulting sputtering yield and energy distributions.This energy range allows for results in the high energy range as well as in the low energy range, where BCA methods may be less suitable.
The rest of this paper is organized as follows.Section II presents the methodologies used to simulate sputtering using the MD and BCA approaches.Section III presents the results for sputtering yields and energy distributions compared to published experimental results and then considers potential low energy dependencies.Section IV concludes the paper with a summary of the findings.

II. METHODOLOGY A. Molecular dynamics simulations
We began by constructing a suitable substrate for our simulations.In specific, a single crystal face centered cubic (FCC) copper substrate was formed from 10 × 10 × 55 unit cells in x, y, and z axes, respectively, with a lattice unit of 3.61 Å, 21 and was oriented, using Miller indices, with the [100], [010], and [001] directions along the x, y, and z axes, respectively.Boundary conditions were periodic in the x and y directions and fixed in the z direction to simulate an infinite slab with a constant thickness and a free surface.The simulation domain extended 250 Å above the free surface to capture any sputtered atoms before they leave the simulation volume.Next, a Berendsen number, volume, and temperature (NVT) thermostat 22 was used to equilibrate the slab to room temperature by keeping the volume and number of atoms fixed while allowing the temperature to reach the desired value.The bottom three layers of the substrate were frozen in time and space to keep the substrate set (Fig. 1).To simulate the temperature dissipation found in a macroscale sample and to eliminate shock wave effects, two layers of atoms on the sidewalls and above the frozen bottom of the substrate were kept at room temperature using the Berendsen thermostat. 13,23,24The remaining atoms were free to interact and increase in temperature via argon induced collisions.
Impacts were simulated one at a time by adding an argon ion above the surface.For each iteration, an ion was initiated from a random x and y location with a velocity corresponding to the desired energy (Fig. 1).All simulated impacts occurred normal to the {100} surface on a 7 × 7 unit cell area.Previous research has shown that channeling is insignificant for {100} FCC copper. 23uring each simulated impact, the region above the slab was scanned for emitted atoms.If sputtered atoms were observed, their atom types, locations, and velocities were recorded.If the impacting argon atom was reflected off the substrate, the atom was not counted as part of the sputtering yield.The scanning continued until no further sputtering was observed and the argon atom was either embedded or reflected off the substrate and out of the simulation box.The MD simulations followed every sputtered atom for the duration of the simulated time.Velocity changes at every time increment were recorded.The selected simulation time per impact was determined by increasing simulation time until the results had plateaued.The total simulation time for each impact was 20 ps corresponding to 20 000 time steps and computation time of ∼3 min per impact.Based on the simulation duration and plateauing of the results, we believe that all ejected atoms were accounted for.Argon impact energies of 200, 600, and 1000 eV were tested.These energies were selected based on the availability of experimental data, with the goal of capturing the sputtered atom energy distribution for both low and high impact energies.To ensure adequate sampling and statistics, 6000 argon impacts were simulated for 200 eV and 2000 argon impacts for all higher impact energies (see Table I).The substrate was kept in rest before each new simulated argon ion impact.For each simulated energy, the total number of sputtered atoms was divided by the number of impacts modeled to obtain a total sputtering yield, i.e., ejected atoms/incident ions.Similarly, the energy distribution of sputtered atoms was obtained by binning the sputtered atom energies into 1 eV bins.
We also simulated a 3000 eV impact energy.Due to the temperature spike caused by the higher impact energy, the simulations were conducted on a larger substrate.The properties of the substrate were identical to those used earlier except that now we used a single crystal FCC copper substrate formed from 34 × 34 × 55 unit cells in x, y, and z axes, respectively.Given the importance of the SBE in the predicted energy distribution and yield of sputtered atoms, MD was also used to calculate the SBE for pure copper.The SBE was defined as the minimum energy needed to completely remove a surface atom from the substrate.Using an identical substrate to the above simulations, a random surface copper atom was given a specific amount of energy and its subsequent position and remaining energy were tracked as a function of time.An iterative method was used to determine the minimum energy needed to remove one copper atom completely from the substrate.Five random surface copper atoms were tested.
In all simulations, the interactions between atoms in the system were calculated using a hybrid interatomic potential.Cu-Cu interactions were simulated using an embedded atom method (EAM) interatomic potential originally developed by Mishin et al. 25 This potential type has been previously used to describe sputtering yields of copper substrates at impact energies between 0.1 and 5 keV 23 and was originally parameterized to describe lattice, defect, and surface properties.Short distance interactions in EAM potentials are often obtained by extrapolation from equilibrium, but this renders them unreliable for the high energy repulsive forces observed during short distance interatomic separation.The potential of Mishin et al. was specifically developed to replicate short approach energy spikes and has been shown to have an improved reliability for shockwave and other high energy simulations.As such, we deemed this potential suitable for capturing the close atomic approaches potentially observed during atomic impacts.Furthermore, Mishin et al. 25 also demonstrated that this potential predicted a cohesive energy of pure copper within the range of previously reported values as determined via ab initio and experimental methods (3.49-3.547][28] All of our simulations were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulation (LAMMPS) package. 29 BCA model simulations BCA simulations were conducted using SDTrimSP 30 and the Thompson equation.For SDTrimSP, a copper substrate was again normally impacted with argon ions at energies of 200, 600, 1000, and 3000 eV.We simulated 500 000 argon ion impacts on a 100-layer thick amorphous copper slab, resetting the substrate for each new impacted simulation (see Table I).Unlike MD, the user does not actually need to consider the substrate crystallinity or decide upon dimensions or boundary conditions.BCA simulations are less computationally intensive than MD simulations, allowing for significantly more iterations to be conducted at each impact energy.A value of 3.49 eV was used for the SBE, based on the pure copper cohesive energy, as recommended by SDTrimSP.31 The Thompson distribution [Eq.( 1)] was also used to predict the energy distribution of the sputtered copper particles.Thompson theory predicts that the energy distribution is solely dependent on the SBE of the substrate.Again, an SBE of 3.49 eV, corresponding to the pure copper cohesive energy, was used.31

III. RESULTS
A. Surface binding energy MD simulations were conducted to determine the SBE of pure copper based on the minimum energy needed to completely remove an atom from the surface.For each surface copper atom tested, the SBE was exactly 4.60 eV.This value is ∼30% higher than the cohesive energy of copper reported by various simulations and experiments (3.49-3.547][28] Given that the interatomic potential used in the MD simulation also predicted a 3.49 eV cohesive energy, the higher SBE is unlikely to be an artifact of an inaccurate potential.A similar finding was reported for the MD simulations by Yang et al. 13 who also observed that the minimum energy needed to remove surface atoms from pure tungsten and beryllium substrates was 30%-40% higher than the individual atom's cohesive energy.

B. Sputtering yield
Argon impacts on the substrate surface led to the sputtering of copper atoms and a damage region below the surface, as can be seen in the MD simulations shown in Fig. 2. The resulting simulated sputtering yields between 200 and 3000 eV using MD and SDTrimSP are compared to experimental values for argon impacts on copper at a normal angle of incidence in Fig. 3.The experimental values represent a compilation of several studies reported by Eckstein. 32for normally incident argon atoms on a polycrystalline copper substrate.The MD and SDTrimSP results at all energies considered agree well with experimental results, though the MD simulations tend to be on the high end of the experimental results and the SDTrimSP on the low end.The MD sputtering yields increase from 1.2 to 4.8 atoms/ion at 200 and 3000 eV, respectively, while for SDTrimSP, the predicted sputtering yields increase from 0.88 to 4.8 atoms/atoms.
In addition, it is interesting to note that if we use the MD SBE for the SDTrimSP calculations, the resulting SDTrimSP sputtering yield is reduced by about 30%.This would make the SDTrimSP results discrepant with the measurements and suggests that for SDTrimSP the cohesive energy is the more appropriate value to use for the simulations.
Sputtering yield results can also be used to investigate the importance of crystallinity in ion sputtering.The experimental studies were conducted on polycrystalline substrates, representing yields for a random distribution of crystal orientations.In contrast, the MD simulations were conducted on the {100} surface of a single crystal and the SDTrimSP simulations were with amorphous copper.Given the good agreement between the sputtering yields for the polycrystalline, single crystal, and amorphous substrates, it appears that the crystal orientation has an insignificant influence on the sputtering yield for normal angles of incidence.This finding is supported by previous experimental 32 and MD studies 23 on both copper and other metals.In a future study, we plan to further investigate the important of crystallinity on sputtered atom energy and angular distribution.

C. Sputtered atom energy distribution
The energy distributions of the sputtered copper atoms at impact energies of 1000 and 600 eV are shown in Fig. 4 and at 200 eV in Fig. 5.In these figures, we compare our MD, SDTrimSP, and Thompson results to the experimental results by Brizzolara et al. for polycrystalline copper sputtered by normally incident argon. 31They measured the energy distribution for particles ejected normal to the surface, and so their results are most readily directly comparable to the Thompson predictions, which assume that the energy distribution is independent of the emission angle.This is to be contrasted with the MD and SDTrimSP results, which yield the energy distribution as a function of emission angle.Comparisons between the experimental results and the MD and SDTrimSP results thus potentially provide information about the importance of the emission angle on the resulting energy distribution.In order to compare our MD and SDTrimSP results at 1000 and 600 eV to the measurements, we have first summed over all 2π sr of the solid angle above the surface.For the 200 eV results, in addition to this, we have also performed simulations where we limited the SDTrimSP results to only those particles ejected within ±15°of the surface normal.For all impact energies, the sputtered particle energy distributions were normalized to one by dividing the counts in each energy bin to the counts in the peak energy bin.
Figure 4 shows that at 1000 and 600 eV, there is good agreement between the BCA simulations and the experimental results.This good agreement suggests that at these impact energies, there is no significant angular dependence for the energy of the ejecta.For the MD results, the general agreement with experiment and BCA simulations is reasonable, though the MD results do peak at an energy of ∼0.5 eV higher.It is difficult to interpret the meaning of this difference as the low numbers in the MD results required the data be binned into 1 eV bins to improve the statistical accuracy of the results.Our findings suggest that for certain systems, both the BCA method and MD may still be reliable for impact energies somewhat below 1000 eV.In addition, the experimental, BCA, and MD results assume three different atomic arrangements: polycrystalline, amorphous, and single crystal copper, respectively.Given the close overlap in energy distributions for all three arrangements, these results suggest that above 600 eV, the substrate crystallinity does not significantly affect the energy distribution.
The close agreement between MD and BCA methods in Fig. 4, despite their differing underlying assumptions, can also be used to better understand the meaning of the E b term in BCA methods.Although MD simulations require a reliable interatomic potential to model many-body effects, BCA methods are limited to independent collisions and approximate the SBE with a cohesive energy.We find that with the BCA methods using an E b value that is 30% lower than the SBE predicted from MD simulations, there is still a close agreement between the two methods.The close agreement observed between BCA and MD simulations despite the different SBEs suggests that the E b term in the Thompson distribution and SDTrimSP may instead be better considered a fitting parameter.For single component systems, the E b term appears well approximated as the cohesive energy.However, we note that according to BCA theory, increasing the SBE from 3.49 to 4.9 eV results in a 30% decrease in the sputtering yield and a broadened energy distribution with a peak that is 30% higher in energy.Potential reasons for this close agreement despite the different SBEs will be explored in a future work.In addition, further research is needed on estimating E b for multicomponent substrates.
Moving to 200 eV, the lowest impact energy studied here, Fig. 5 shows significant differences between experiment and theory in the peak position, the width of the energy distribution, and the magnitude of the high energy tail.Starting with the peak position, there is good agreement between the experimental results and the Thompson distribution.The SDTrimSP peak is about 0.5 eV higher than these, while the MD peak is about 2 eV higher than the experimental and Thompson results.The experimental distribution also narrows as compared to the higher impact energy experimental results and to the Thompson and SDTrimSP results.The full width at half maximum (FWHM) for the 200 eV experimental results is 5.8 eV.For the higher impact energy experimental results and Thompson theory, it is 6.7 eV, while for SDTrimSP, it is 6.9 eV.The FWHM of the MD results at 200 eV is 6 eV, but the binning results in an uncertainty of ∼±0.5 eV, making it hard to infer the presence or absence of a discrepancy.Finally, starting from slightly above the peak energy, all three of the theoretical results are higher than the experimental data.
The discrepancies between the experimental results and the Thompson distribution are potentially due to the theory not accounting for the angular dependence of the ejecta with impact energy.A behavior similar to what we observe was been seen by Goehlich et al. 33,34 They measured the sputtered particle energy distribution for atoms emitted normal to the surface due to argon ions incident at normal angles on tungsten and aluminum substrates for both high and low impact energies.Their experimental results found good agreement with Thompson distributions at high impact energies, but at low impact energies the distributions were narrower and had a lower magnitude for the high energy ejecta.
Moving on to the SDTrimSP and MD results, the source of the observed discrepancies with the measurements potentially lies in the experimental conditions and in a dependence of the ejecta FIG. 3. Sputtering yield for normally incident argon impact on a copper substrate using MD (blue circles) and SDTrimSP (gray triangles), as compared to the experimental results sampled from Ref. 31 (red squares).The blue and gray curves are fits to the theoretical results.
energy distribution with the emission angle of the sputtered atom.To further investigate this dependency, we used SDTrimSP to calculate ejecta energy distributions for emission angles ±15°from normal.We compare these results to the experimental data in Fig. 5(b).The fluctuations in the SDTrimSP results are due to lower sample counts during binning and make it difficult to comment on whether the emission angle dependency is the source of the peak position discrepancies.However, the emission angle does appear to explain the remaining discrepancies.When only those atoms emitted close to the surface normal are considered, the SDTrimSP width and magnitude of the energy distribution are now in good agreement with the measurements.These results suggest that for low energy impacts, the energy distribution of sputtered atoms is dependent both on the impact energy and on the angle of emission.As such, at low impact energies, it is important to account for possible dependencies in the energy distribution with the impact energy angle of emission of the ejecta.In a future study, we will elaborate on these dependencies and attempt to further explain the physical reasons behind them.

IV. SUMMARY
We have performed new theoretical calculations for the sputtering yield and energy distribution resulting from argon impacting a pure copper substrate at normal incidence with impact energies between 200 and 3000 eV.For energies of 600 eV and higher, the sputtering yields and energy distribution as predicted by the MD and BCA methods agree well with the experimental results.Furthermore, despite BCA methods using an E b value that is 30% lower than the SBE predicted from MD simulations, there is still close agreement between the two methods.This close agreement suggests that the E b term in BCA methods may instead be considered a fitting parameter.Further research is needed on estimating E b for multicomponent substrates.At an impact energy of 200 eV, there were significant discrepancies between the experimental and theoretical energy distributions in the peak position, the width of the energy distribution, and the magnitude of the high energy tail.The sources of these observed discrepancies appear to be due to a combination of factors.The experiments only measured atoms sputtered normal to the surface, and there is a dependence of the ejecta energy distribution with the emission angle of the sputtered atom.When only those atoms emitted close to the surface normal were considered, the SDTrimSP width and magnitude of the energy distribution were in good agreement with the measurements.More theoretical and experimental studies are needed to understand the origin of these low energy dependencies.

FIG. 1 .
FIG.1.Illustration of the MD setup for argon impacts (black) on a copper substrate showing the impact, thermostat, and fixed regions, in green, orange, and red, respectively.

FIG. 2 .
FIG. 2. MD calculations showing a section of the copper substrate surface (A) before impact with an argon ion at 1000 eV and (B) after.Copper atoms are shown in orange and the argon ion in blue.The x axis is into the page.

FIG. 4 . 31 FIG. 5 .
FIG. 4. Normalized sputtered atom energy distributions at (a) 1000 and (b) 600 eV.The blue bars are the MD results; the black tops show the ±one-sigma counting statistics uncertainty; the red lines show the SDTrimSP results; the purple dashed lines show the Thompson distribution; and the green lines show the experimental results.31

TABLE I .
Total numbers of impacts and sputtered atoms for our MD and SDTrimSP simulations.