Why are physical sputtering yields similar for incident ions with different masses?—physical sputtering yields of the Lennard–Jones system

Plasma etching of nano-meter-scale complex structures for semiconductor device manufacturing requires a deeper understanding of etching mechanisms. For example, it is known experimentally that the sputtering yield of a material tends to have weak dependence on the mass of incident ions except for extremely light ions such as helium. To understand this property, the sputtering yield of a system of atoms interacting with Lennard–Jones (LJ) potentials was evaluated with molecular dynamics simulation. As the simplest possible case involving two atomic species, a single-element face-centered-cubit (fcc) LJ solid surface interacting with purely repulsive atoms was examined, which emulates a solid surface sputtered by noble-gas ions. The sputtering of such a system at specific incident ion energy depends only on two parameters, i.e. the mass ratio and a parameter representing the relative interaction range between the surface atom and the incident ion. For real materials of our concern used in plasma etching, the range of these two parameters was found to be relatively limited. It was also found that the physical sputtering yield of the LJ system weakly depends on the mass ratio in this relatively narrow parameter range. Because the simple model predicts the weak yield dependence on the incident ion mass, it is considered as a generic property of physical sputtering, independent of the detailed atomic interactions of the surface material and incident ion species.


Introduction
Plasma etching has significantly contributed to the miniaturization of semiconductor devices as it is one of the main driving * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
technologies for the fabrication of micro/nano-scale structures on wafer surfaces [1][2][3][4][5]. Plasma etching refers to the process of removing atoms from a surface by the impact of energetic and/or chemically reactive incident species (ions and radicals, typically) generated in a plasma. Chemically reactive species adsorbed on the surface modify the surface chemical compositions and the impact of energetic incident ions causes collision cascades of atoms inside the material, whose billiard-like multiple collisions can lead to the ejection of some atom from the surface. Synergetic effects between the chemical modification of the surface and the ion impact determine the sputtering yield (i.e. the average number of atoms removed from the surface per ion impact) and, therefore, such etching is called reactive ion etching (RIE) [6][7][8]. If no modification of the surface chemical composition takes place, the process is called physical sputtering or physical etching [9,10].
When deep and narrow structures such as high-aspect-ratio holes and slits used, for example, in three-dimensional (3D) NAND flash memory devices [11] are formed with RIE processes, ion incident energies are typically set high. It is known that, even in RIE processes, the higher the ion incident energy becomes, the more physical the nature of etching becomes. For example, with high ion incident energies, the sputtered (or desorbed) species from the surface are typically single atoms or fractions of molecules with small numbers of atoms rather than large molecules [12]. Such sputtered species typically have high sticking probabilities with material surfaces, resulting in the redeposition of sputtered materials on the sidewalls of etched deep structures [13][14][15]. Sputtering deposition processes are another example where physical sputtering plays an important role in material processing. Sputtered species from the target surface are used to form thin films over micro/nano-scale structures on a wafer facing the target surface. For example, in the case of ionized magnetron sputtering, the deposited thin film on the wafer is subject to simultaneous etching and deposition processes, depending on the fluxes and energies of incident ion and neutral species [16][17][18][19].
The experimental and simulation data obtained so far indicate that, for a given surface material and ion incident energy, the sputtering yields tend to depend weakly on the mass of incident ion species, except for extremely small ions such as helium (He), if the sputtering is dominated by physical sputtering. For example, figure 1 shows the comparison of the experimentally obtained sputtering yields of silicon (Si), iron (Fe), nickel (Ni) and copper (Cu) for self-sputtering and sputtering by noble gas ions, i.e. helium (He + ), neon (Ne + ), argon (Ar + ), krypton (Kr + ) and xenon (Xe + ). Silicon (Si) and iron (Fe) have diamond and body-centered cubic (bcc) crystalline structures whereas nickel (Ni) and copper (Cu) have facecentered cubic (fcc) crystalline structures. The data are taken from [21] and the solid curves are the empirical fitting curves proposed by Yamamura and Tawara [21]. It is seen that, except for He + sputtering, the sputtering yield hardly depends on the incident species in an energy range from about 100 to several keV. At larger or lower incident energies, a dependence on the incident species is observed. On the other hand, the sputtering yields by He + ions are lower in most of the energy range. Similar weak mass dependence is also observed for the sputtering yields of composite materials [30,31].
The weak dependence of sputtering yields on the incident ion mass in the intermediate ion-incident-energy range seems a generic property of sputtering. The goal of this study is to understand why it is so and clarify how the sputtering yield depends on the nature of the material and incident species as well as the incident ion energy. To achieve this goal, we performed MD simulations of atoms interacting with simple interatomic potential functions. The reason we use simple interatomic potential functions, rather than realistic ones, is that we conjecture that this generic property does not depend on the details of specific atomic interactions of the system we consider.
One of the simplest and most widely used interatomic potential functions to describe materials and fluids is the Lennard-Jones (LJ) potential. The simplicity of the LJ potential function makes it a popular choice to model approximate covalent and metallic bonds as well as van der Waals interactions. The LJ potential has been widely used for the study of thermodynamics properties [75], transport properties, and equations of state [76] for solids, fluids, and gases, and the phase transitions among them [77][78][79][80][81].
A recent study [82] has shown that the LJ system, i.e. the system of particles interacting with the LJ potential functions, can reproduce the self-sputtering yields of real materials reasonably well except for a high ion-incident-energy range, where the discrepancy between the repulsive potential of the LJ system and those of real atoms becomes significant [82,83]. Therefore, in this study, we extend the LJ-system of a single element used in [82] to that for general sputtering, i.e. the system where the surface material and incident particles are of different species, and compare the sputtering yields of the model system with those of real materials in a relatively low ion-incident-energy range. To keep the system as simple as possible, we construct a single-element solid consisting of LJ particles (i.e. particles interacting with the LJ potential) and inject energetic particles having an LJ-like repulsive potential and a different mass to the solid surface. The system is intended to emulate a solid interacting with noble gas ions such as Ar + . As we shall discuss in the following section as well as in appendix A, in this way, we can minimize the number of system parameters while keeping the essential physical mechanisms that govern the sputtering phenomena. Such a system may be considered as a minimal extension of the LJ system of a single element used in [82]. Because the LJ system is one of the most extensively studied models for materials and fluids, determining its sputtering properties may also be considered valuable in itself. Therefore, in this study, we shall examine sputtering properties of the LJ system of a single element by repulsive particles of a different species.
The rest of this paper is structured as follow: section 2 presents the LJ system used in this study. Section 3 outlines the simulation method. In section 4, the sputtering yields of the LJ system obtained from MD simulations are presented and compared with the corresponding experimental values. The conclusions are presented in section 5.

LJ potentials and equations of motion
In this study, the interactions among atoms are represented by the LJ potential function [84]. No Coulomb interaction was taken into account in this study because, in real systems, incident ions are expected to be charge-neutralized by an Auger emission process right before they impact the surface [85]. The LJ potential function depends on two parameters: the 'length' parameter σ, which is the interatomic distance at which the LJ potential is null, and the 'binding' parameter ε, which represents the strength of the interaction. The simplicity of the LJ potential function makes it a popular choice to model approximate covalent and metallic bonds as well as van der Waals interactions, as discussed above. It is, however, known that the short-range part of the LJ potential is too repulsive compared with the interatomic potential of real atoms [82,83].
In this study, we consider two atomic species. The interaction among atoms that constitute the surface material is modelled with a full LJ potential function V αα (r). The interaction between an atom of the surface material and an injected atom is modelled with a repulsive potential function V αβ (r). The interaction among injected atoms is modelled with another repulsive potential function V ββ (r). The expressions of these functions are given by Here r denotes the interatomic distance. The subscript α or β indicates that atomic species. It should be noted that the parameters ε αβ σ αβ 12 and ε ββ σ ββ 12 always appear as these combinations and, for example, ε αβ and σ αβ are not individually defined in the framework of this study. However, we keep these parameters as they are because, as we shall discuss later, the corresponding values of ε αβ and σ αβ can be defined separately for most atoms, including noble gases, and can be used for the comparison between the simulation results and experimental observations. Weakly bound atoms such as noble gas atoms tend to have small ε values while σ values are similar to those of other atoms. If the interactions among all atoms are represented by full LJ potential functions, the potential energy system of equation (1) is valid in the limit of ε αβ /ε αα → 0 and ε ββ /ε αα → 0 (see appendix A). The dynamics of the LJ system of N atoms of two species is governed by the following equations of motion; It is understood that the summation over j is taken in such a way that j ̸ = i. Here m α and m β are the atomic masses of species α and β, v i is the velocity of the ith atom, r i is the position of the ith atom, r ij = |r i − r j | is the interatomic distance between the ith and jth atoms. i ∈ α means the ith atom is of species α. ∑ j∈α and ∑ j∈β denote the summation over all atoms j ̸ = i belonging to species α only and β only.
By setting the following normalised variables and system parameters The normalised equations of motion can then be written as: Because incident species have only repulsive interaction with other atoms, they are most likely to leave the surface together with the surface atoms when sputtering takes place. Therefore, we expect incident species hardly interact among themselves and the second term of the LHS of equation (3b) hardly affects the sputtering yield. Therefore, for the sake of simplicity, we assume in this study With these assumptions, the system now depends only on two parameters µ and ρ and the equations of motion (3) can be written as with the normalised interatomic functions defined bỹ Now we examine the interatomic potential functions of this system. The functionṼ αα takes its minimum value at ξ = ξ M = 2 1/6 , i.e.Ṽ αα (ξ M ) = −1/4 and becomes null at ξ = 1. The repulsive potential functionṼ αβ becomes unity at ξ = ρ, which indicates that the parameter ρ represents a typical interaction range between a surface atom α and an incident atom β relative to the interaction range between two surface atoms. Therefore, in this article, we call ρ the relative interaction range of incident ions. Figure 2 shows the normalised LJ potentialṼ αα and three different repulsive potentialṼ αβ with three different values of ρ = 1.0, 0.5, and 2.0.
With the normalization of equation (2), any energy quantity is normalised with 4ε αα , as we have seenṼ αα = V αα /4ε αα . The normalised kinetic energy of an atom is also defined as where E = mv 2 /2 is the kinetic energy in dimensional units, with v being the magnitude of the incident particle velocity v = |v| andv = vT 0 /σ αα .

Outline of MD simulation
MD simulation was performed to examine sputtering phenomena of the LJ system. All MD simulations in this study were performed with the classical MD simulations code LAMMPS [86,87] distributed under the GNU General Public License (GPL) license by Sandia National Laboratories, USA.
The material used in this study was made of 27 000 atoms arranged in a face-centered cubic (fcc) crystalline structure [88] of a lattice constant of l a = 1.54σ. (Here and in what follows, the subscript αα of equation (1a) is omitted for conciseness unless needed.) It had a rectangular shape with a dimension of 15, 15, and 30 times the lattice constant l a in the x, y and z directions, respectively. In all cases, the Miller index of the top surface was (100). The material was placed at the bottom of a simulation box having the same size in the horizontal directions, i.e. x and y directions, while having a size of 35 times the lattice constant l a in the z direction to maintain some space above the material for particle injection and surface dynamics. Periodic boundary conditions were applied in the horizontal directions to model an infinitely large surface. The position of the atoms in a few layers at the bottom of the material were fixed to prevent the material from moving downward due to the particle injection. All atoms leaving the simulation box were considered lost and removed from the system. Since the material is thick enough, all particles leaving the simulation box did so through the top boundary, meaning that they are either sputtered atoms from the surface or reflected incident particles.
MD simulations were performed with a material having a surface four times larger (108 000 atoms) or twice deeper (54 000 atoms) to confirm that the use of a larger material gives essentially the same sputtering yields. The material was initially thermalised at a normalised temperature of T = k B T/4ε = 0.006, with k B being the Boltzmann constant and T being the material temperature. For a typical material with ε = 1 eV, this normalised temperature corresponds to a temperature of T = 300 K.
At the beginning of each injection cycle, a single atom was injected with normalised incident energy E at normal incidence, i.e. along the z axis. The simulation was performed with a normalised timestepδt = 1.7 × 10 −3 , which corresponds to a timestep in physical units δt = 0.1 fs for a material with m = 58.693 g·mol −1 (corresponding to the mass of Ni), σ = 1.5 Å and ε = 1 eV. Up to τ = 12.0, MD simulation was performed with constant-energy conditions (NVE). From τ = 12.0 to τ = 30.8, a Langevin cooling [89] was applied, followed by the application of constant-temperature conditions (NVT) between τ = 30.8 and τ = 34.2 to force the system to reach a thermalised state atT = 0.006. The relaxation of the system after a single ion bombardment would take a much longer time in real materials. The damping parameter (i.e. inverse of the normalised friction coefficient) τ damps and the target temperatureT of the Langevin cooling were τ damp = 17.1 andT = 0.006. At the end of the cycle, i.e. at τ = 34.2, all atoms out of interaction with the substrate bulk are considered as sputtered atoms and removed from the simulation to avoid the accumulation of floating atoms above the substrate.
This injection cycle was repeated around 1000 times (i.e. around 1000 injections or a normalised ion dose of 1.87) until the sputtering yield reached steady state, i.e. became essentially independent of the ion dose. A summary of the simulation parameters with the corresponding dimensional values is given in table 1. The substrate used in this study is shown in figure 3(a).

Sputtering yield
As mentioned earlier, the sputtering yield Y of surface atoms of species α is defined as the average number of atoms of species α removed from the material per ion impact, i.e.
with N gas α and N in β denoting the total number of atoms of species α removed from the material surface (and therefore in the gas phase) and the number of incident ions of species β. The sputtering yield is usually measured in steady state, i.e. when Y averaged over a certain period of ion incidence hardly depends on the ion dose.
If the incident ions are of the same species as surface atoms, the experimentally measured sputtering yield (i.e. selfsputtering yield) is defined differently because the incident and surface species cannot be distinguished experimentally. The self-sputtering yield Y self is defined [82] as , which correspond to the sputtering of Ni by Ne + , Ni + , and Xe + ions, respectively. The blue, red, and black dashed curves are the empirical fitting functions by Yamamura and Tawara [21] for experimentally obtained sputtering yields of Ni by Ne + , Ni + , and Xe + ions. It should be noted that the red curve represents the self-sputtering yield Y self of Ni, rather than the sputtering yield of surface atoms Y.
with N gas all and N in denoting the total number atoms removed or reflected from the material surface and the number of incident ions. It should be noted that, if Y self < 1, the net deposition takes place even if surface atoms are sputtered by ion impact.
With α = β, the relation between Y and Y self may be given formally as with N gas β denoting the total number atoms of species β desorbed or reflected from the material surface and s = 1 − N gas β /N in β being the sticking probability of the incident species. Most self-sputtering cases that we are interested in are those of stable solid materials such as metal. In such cases, energetic ions with normal incidence typically penetrate deeply into the material and get deposited, so their sticking probability is close to unity, i.e. s ≃ 1 and sputtered species are mostly surface atoms. Therefore, the sputtering yield of surface atoms Y is essentially equal to the self-sputtering yield Y self . However, when the incident energy is low or the angle of incidence is large from the surface normal, the sticking probability s becomes small and the difference between Y self and Y can become significant.
Furthermore, in our system of equation (4), the incident atoms only have repulsive potential and therefore cannot be the same as the material atoms even with parameters of µ = 1 and ρ = 1. However, as discussed above, if the incident energy is sufficiently high and the injected atoms are embedded in a deeper region of the material, we expect the yield Y of sputtered species from the surface is comparable with the selfsputtering yield Y self of the full LJ system [82] because the repulsive potential is similar. A more detailed discussion is found in appendix B.

Energy dependence
MD simulations were performed to determine the dependence of the sputtering yield Y on the normalised kinetic energy E. The angle of ion incidence was normal to the surface. Figure 4 shows the obtained physical sputtering yield for a couple of combinations of system parameters µ and ρ. We did not plot the yield for the energy higher than E > 10 3 because, at such high energies, the sputtering yields obtained from the LJ system significantly deviate from those of the real materials. For the model to emulate a real system better, the repulsive potential needs to be modified in general [82].
It is seen that the sputtering yield increases with the normalised energy E and is similar to those with different parameters in an energy range of around E > 200. The dotted lines are the empirical fitting curves to the experimentally observed sputtering yields of Ni by Ne, Ni, and Xe ions. The comparison with such experimental values will be discussed in section 4.4. The parameters µ and ρ selected here also correspond to those of the actual systems, as will be discussed in section 4.4.

Dependence on parameters µ and ρ
The dependence of sputtering yield on the two parameters µ and ρ at E = 250 is shown in figure 5. The angle of ion incidence was normal to the surface. The filled circle, filled up-pointing triangle, and filled down-pointing triangle of this figure are the cases corresponding to the data plotted in figure 4. The corresponding experimental values are plotted with crosses, which will be discussed in section 4.4.
It is seen that the yield value Y varies up to an order of magnitude as ρ increases from 0.33 to 2. More specifically, for each value of µ, as ρ increases, the sputtering yield of surface atoms Y increases until it reaches its maximum and then decreases. In other words, if the incident ions are too small compared with the interatomic distance of the material surface (e.g. ρ < 0.5), the sputtering yield is low probably because the incident ions penetrate the material too deeply and release their energies and momenta along longer collision paths, making the energy and momentum transfer to each surface atom less efficient. On the other hand, if the incident ions are too large (e.g. ρ > 1.2), the sputtering yield is again low probably because they cannot penetrate the material deeply, thus mostly affecting only the top surface atoms they collide with.
It is also seen in figure 5 that the yield Y is the largest around µ = 1 and ρ = 1 possibly because the momentum transfer among atoms with the same mass and size is the most efficient.
However, this observation is not universal and depends on the incident energy, as we see in figure 4, where the sputtering yield for µ = 0.34 and ρ = 0.91 is the highest among the three cases in an energy range of E < 100. In general, if the value of µ is extremely small (e.g. µ < 0.1), the sputtering yield is also extremely small as light ions cannot transfer much energy or mass in a single impact.

Comparison with experiments
So far we have evaluated the sputtering yield for a range of µ and ρ without reference to real materials. For real materials with specific incident ions, these two parameters are not independent. For example, for smaller incident ions, both mass ratio and relative interaction range are typically small. In this subsection, we discuss how these parameters can be set for real materials.
One can determine ε and σ of the LJ potential from the basic properties of a solid that the corresponding atoms constitute, especially if the solid is a face-centered cubic (fcc) crystal that the LJ particles are known to form under low-temperature and low-pressure conditions.
The parameter σ can be calculated from the closest neighbour distance d, which is known for most materials [90]. For an fcc crystal, the closest neighbour distance d, the lattice constant l a , and σ are related by  [92] The parameter ε can be obtained from the cohesive energy E C , which is defined as the energy per unit amount of substance needed to separate all the constituent atoms into noninteracting neutral atoms. The cohesive energy is also known for most materials [90]. The cohesive energy E C of an fcc crystal and ε are related by Both relations (6) and (7) are in good agreement with previous calculations [91]. For simplicity, we extend these definitions of σ ε for any material (even if it does not form an fcc lattice) for which the closest neighbour distance d and cohesive energy E C are defined. It is understood that, with such extended definitions of σ and ε, the actual atomic interaction for some materials may not be well represented by the LJ potential. The parameters σ and ε determined in this way for various gaseous and material species are listed in table 2. (These parameters of gaseous species at standard temperature and pressure were determined from their solid states at sufficiently low temperature.) For interactions between atoms of different species, the values of σ and ε may be estimated from those parameters of each species. In this study, we use the following mixing rules, the Lorentz rule [93] for σ and the Berthelot rule [94] for ε; Tables 3 and 4 show the values of the mass ratio µ and relative interaction range ρ obtained from the Lorentz-Berthelot combining rules above with the parameters listed in table 2 for the interaction of materials C, Al, Si, Fe, Ni, Cu, and Au with incident ionic species of He, Ne, Ar, Kr, and Xe.
The correlation between the parameters µ and ρ is seen in figure 6, where the values of ρ and µ in tables 3 and 4 are plotted. It is seen, while µ varies by nearly three orders of magnitude, the variation of ρ is rather limited. The correlation is nearly linear as indicated by the dotted curve. By comparing figure 6 with the sputtering yields summarized in figure 5, we observe that the sputtering yield indeed depends strongly on these parameters in general. However, if we confine our interest to materials of relatively large atomic masses (such as Al or larger) and exclude extremely light or heavy incident species relative to the mass of a surface atom (such as He for most surface materials or Xe for some light surface materials), the values of the parameters fall into the ranges of 0.3 < µ < 3 and 0.8 < ρ < 1.2, as we see in tables 3 and 4. Figure 5 shows that the sputtering yields at E = 250 are relatively similar within this parameter range; although the µ varies by an order of magnitude in this range, the yield varies only up to a factor of 2.
In figure 4, the sputtering yields for (µ, ρ) = (0.34, 0.91), (1.00, 1.00), and (2.24, 1.19) correspond to those of Ni by Ne + , Ni + , and Xe + ions, according to tables 3 and 4. The blue, red, and black dotted curves are the empirical fitting functions proposed by Yamamura and Tawara [21] for experimentally obtained sputtering yields of Ni by Ne + , Ni + , and Xe + ions. It should be noted that the red dotted curve represents experimentally obtained self-sputtering yields of Ni. The three curves take similar values in an energy range of 100 < E < 1000 although their values are lower than the yield values of the LJ system for larger E values. As discussed in [82], this is expected because the repulsive potential of the LJ system is too strong compared with those of actual atoms. Figure 7 compares the sputtering yields of the LJ system obtained in this study with experimentally obtained sputtering yields of various materials. The LJ-based yield values obtained from MD simulations are the same as those plotted in figure 5 and their symbols are also the same. The experimental data were obtained from the empirical fitting functions by Yamamura and Tawara presented in [21]. Especially, the experimental yield values of Ni by Ne + , Ni + , and Xe + ions are plotted by blue, red, and black crosses, whose values of µ are 0.34, 1.00, and 2.24, respectively. It is seen that these values are lower than the corresponding yield values of the LJ system.
Except for the experimental yield values of Ni by Ne + , Ni + , and Xe + ions mentioned above, no simulation data is presented with the same µ and ρ values as those of the corresponding experimental data. However, the mass ratios of the experimental data presented here mostly fall in the range of 0.3 < µ < 3 with a few exceptions. In most cases, the experimental values seem close to the corresponding simulation values within a factor of 2 to 3, except for Al and Si, whose yields are much lower than those of the corresponding LJ system. The mass ratio of Al to Kr + or Xe + is larger than 3 and Si has a diamond crystalline structure, which may partially account for the discrepancies. For each material, however, we observe the weak dependence of the yield on the mass ratio µ, which varies more than a factor of 7 in this figure. Figure 4 also shows that, in a lower energy range of E < 100, the sputtering yields of these three cases obtained from MD simulations deviate from each other and so do the three dotted curves. It is seen that the circles follow the corresponding experimental curves. The predicted sputtering yields of Ni by Ni + and Xe + ions deviate significantly from the experimental curves. As discussed earlier, the sputtering yield Y plotted here is different from the self-sputtering yield Y self , which the red dotted curve represents. The disagreement between the red dotted curve and red up-pointing triangles may arise from this difference in definition. For more discussion on this, the reader is referred to appendix B. The reason for the discrepancy between the black dotted curve and black down-pointing triangles, referring to the sputtering yield of Ni  figure 5 and their symbols are the same. The experimentally obtained yields of Al, Si, Fe, Ni, Cu, and Au are denoted by the symbols listed in the legends. For each material, from the left to right, the experimentally obtained sputtering yields by Ne + , Ar + , Kr + and Xe + are plotted, except for the self sputtering at ρ = 1. Among them, the experimentally obtained yields of Ni by Ne + , Ni + , and Xe + ions are shown with blue, red, and black crosses, whose values of µ are 0.34, 1.00, and 2.24, respectively. The corresponding yields obtained from the LJ-based simulations are plotted by the filled blue circle, filled red up-pointing triangle, and filled black down-pointing triangle, as in figure 5. It should be noted that the red cross represents the self-sputtering yield Y self . All experimental data were obtained from the empirical fitting functions by Yamamura and Tawara [21].
by Xe + ions, especially at low energies, is not clear. When the incident ion energy is low, the sputtering yield may depend more sensitively on the details of actual atomic interactions.

Conclusions
We evaluated sputtering yields of the (100) surface of an fcc solid whose atoms interact with the LJ potential subject to the bombardment of ions that interact with surface atoms with LJ-type repulsive potential, using MD simulation. No Coulomb interaction was taken into account in this model because, in real systems, incident ions are expected to be charge-neutralized by an Auger emission process right before they impact the surface [85]. The equations of motion solved in this study are given by equation (4) in the normalised form and the system depends only on two parameters, i.e. the ratio of the incident ion mass to the mass of surface atom µ and the relative interaction range of the incident ion ρ.
This study is essentially an extension of an earlier study of [82] on the self-sputtering of the LJ system. In [82], the incident species are assumed to be the same as those of the surface atoms and obey the full LJ potential whereas, in the present study, the incident species have only repulsive interactions with surface atoms. Therefore, strictly speaking, the system of this study does not include that of [82]. However, as the repulsive potential of energetic incident species, rather than their attractive interactions with surface atoms, is expected to affect the collision cascade dominantly, the average number of surface atoms ejected from the surface by knock-on collisions should be similar between the two systems. Using the LJ system for two species defined above, we evaluated the sputtering yields by MD simulations. In the energy and parameter range we examined, i.e. 10 < E < 1000, 0.07 ⩽ µ ⩽ 2.24, and 0.33 ⩽ ρ ⩽ 2, the sputtering yield vary by orders of magnitude. As well known from earlier studies, the sputtering yield is an increasing function of the incident energy and varies up to 10 or more. In general, if the value of µ is extremely small (e.g. µ < 0.1), the sputtering yield is also extremely small as light ions cannot transfer much energy or momentum to the surface in a single impact. For a given incident energy E and a mass ratio µ, as the incident ion becomes smaller or larger than certain values (e.g. ρ < 0.5 or ρ > 1.2), the sputtering yield becomes lower. This is probably because, if the radius of incident ions is too small, they penetrate the material too deeply and release their energy and momentum along longer collision paths, thus making the energy and momentum transfer to each surface atom less efficient. If the radius is too large, they cannot penetrate the material deeply, thus mostly affecting only the top surface atoms they collide with.
The goal of this study is to explain why physical sputtering yields tend to have weak dependence of the mass of incident ions except for extremely light ions such as He + . This is generally the case when the incident energy is less than a couple of thousands eV, i.e. the energy range of interest for semiconductor device manufacturing processes. If the incident ion energy is low and the sputtering yield is less than 0.1, the sputtering yield seems to depend on the mass of incident species more strongly, but the yield is so small and the difference among small yields hardly matters in most practical applications. To determine the system parameters µ and ρ of the LJ system most relevant for real materials, we used equations (6) and (7) as well as the Lorentz-Berthelot combining rules. For real physical systems, the values of µ and ρ have a strong correlation. It was found that, for the combination of materials and ions of our interest (excluding extremely light ions), the set of parameters µ and ρ falls into a narrow range of around 0.3 < µ < 3 and 0.8 < ρ < 1.2. The MD simulations of this study have shown that the sputtering yields are relatively similar in this parameter range. Although the µ varies by an order of magnitude in this range, the yield varies only up to a factor of 2. This accounts for the similarity of sputtering yields among different incident ion species. Admittedly, how similar or different the sputtering yields among different incident ionic species are rather subjective. As we see in figure 1 and discussed above for the LJ system, the sputtering yield does depend rather strongly on the mass of incident energy in the high or low ion-incident-energy ranges. Probably because, for most processing applications, we use the ion incident energy ranging from around 100 to a few thousand eV and care less about extremely low sputtering yields, we might have formed a rather biased impression that the sputtering yield is less dependent on the mass of incident ions. What we have clarified in this study using the LJ system is that, as emphasized above, in the limited parameter range of around 0.3 < µ < 3 and 0.8 < ρ < 1.2, the sputtering yield varies within a factor of 2 or so.
Concerning the discussion above, another general trend that we may find in figure 1 is that, for each surface material, the sputtering yield by heavier ions is higher at high ion incident energy and lower at low ion incident energy. Even for extremely light He + ions, its sputtering yields are close to those by heavier ions at low ion incident energy. Except for He + ions, the sputtering yields by heavier ions overtake those by lighter ions in the intermediate energy range as the ion incident energy increases. This means that the sputtering yields by most ions are nearly equal in the intermediate ionincident-energy range, which we have attempted to explain in this study. A similar trend is seen in figure 4, where the sputtering yields by heavier ions are shown to be much lower than those by lighter ions at E < 100 for the LJ system. The reason for this trend is not clear and left as a subject of future work.
How representative the LJ system is of real systems is another question. It has been known that the repulsive potential of the LJ system is much stronger than those of real atoms [82,83], which tends to give unrealistically high sputtering yield at high ion incident energies. This is the reason that, in this study, we focused our study on a normalised energy range of E < 1000, which corresponds to an energy range of less than a couple of thousands eV in physical units. Even in this range, our LJ model tends to overestimate the sputtering yield of real systems. However, our goal of this study is to understand the qualitative nature of sputtering and therefore discussion on the validity of the LJ's repulsive potential is deferred to a future study.
A recent study [95] showed that, based on analysis of a large amount of experimentally obtained sputtering yield data, the sputtering yield and related properties depend essentially on a relatively small number of physical parameters. In other words, the sputtering yield may not depend sensitively on the details of interatomic potentials. Based on this premise, we surmise that the present study based on the LJ system has captured some essence of the sputtering properties of real materials.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Instead of the repulsive potentials of equations (1b) and (1c), if we employ the full LJ system for incident particles as  Figure B1. Comparison between the sputtering yield of surface atoms Y of the LJ system of this study for µ = ρ = 1, denoted by triangles (as given in figure 4), and the self-sputtering yield Y self of the full LJ system from [82], denoted by circles. The dashed curve is the empirical fitting functions proposed by Yamamura and Tawara for experimentally obtained self-sputtering yields of Ni given in [21].
in addition to equation (1a), the normalised equations of motion corresponding to equations (3) become (i ∈ β). (A.1b) Integrating these equations by numerical simulation is not much more time-consuming than integrating equations (3) or (4). However, equation (A.1) have more system parameters than those of equations (3) or (4). The reason we employ a simpler set of equations (4) in this study is to limit the range of our parameter survey to that of a smaller number of important parameters. Equations (A.1) above approach equation (3) in the limit of ε αβ /ε αα → 0 and ε ββ /ε αα → 0. Figure B1 shows the sputtering yield of surface atoms Y for µ = ρ = 1 (denoted by filled triangles; data are the same as those in figure 4) and the self-sputtering yield Y self of the full LJ system of single species given in [82] (denoted by the empty circles). The solid curve, following the filled triangles, is a guide to the eye. The definitions of Y and Y self are not the same and they are related via equation (5), as discussed in section 4.1. It is seen that yield data of filled triangles and empty circles agree well in a high energy range of E > 100, where we expect the sticking probability of the incident species is high, i.e. s ≃ 1. Under such conditions, Y ≃ Y self . On the other hand, when the incident energy is low, the incident species cannot penetrate the material deeply and may leave the surface after impact, even if they have attractive interaction with surface atoms, resulting in s < 1. Therefore the difference in definition between Y self and Y, at least partially, accounts for the disagreement between filled triangles and empty circles at low incident energies.

Appendix B. Difference between Y self and Y
This difference can also arise from the difference in interatomic potential. In the LJ system of this study, incident species react with surface atoms only repulsively whereas, in the full LJ system of [82], both incident species and surface atoms have the same full LJ potential. In the former system, incident species are unlikely to be deposited. In the latter system, if incident species are deposited and the surface structure changes, the sputtering yield may change accordingly. Discussion on such effects on the sputtering yield is out of the scope of this work and will be analysed elsewhere.
The experimental self-sputtering yields are also plotted by the dotted curve. [21]. It is seen that the empty circles show good agreement with the experimental values for E < 10, as both represents Y self . In a higher energy range of around E > 200, the yields obtained from the LJ system are consistently higher than the experimental data. As discussed in [82], this is due to the inaccurate representation of the LJ's repulsive potential for real systems.