Low energy ion-solid interactions: a quantitative experimental verification of binary collision approximation simulations

Ultra-low energy ion implantation has become an attractive method for doping of two-dimensional materials and ultra-thin films. The new dynamic Monte Carlo program IMINTDYN based on the binary collision approximation allows a reliable prediction of low energy implantation profiles and target compositional changes, as well as efficient simulation of high energy light ion scattering. To demonstrate the quality of these predictions and simulations, we present a model case experiment where we implanted W ions into tetrahedral amorphous carbon with low (10 keV) and ultra-low (20 eV) ion energies and analyzed the W implantation profiles with high resolution Rutherford backscattering spectrometry (HR-RBS). This experiment is compared with a complete simulation of all aspects of ion-solid-interactions of the experiment using the new IMINTDYN program. A unique novel simulation option, also relevant for implantation into 2D materials, is the inclusion of the vacancy as target species with dynamic vacancy generation and annihilation. Whereas simulations neglecting vacancy formation cannot reproduce the measured implantation profiles, we find excellent agreement between simulated and measured HR-RBS spectra. We also demonstrate the important role of simultaneous weak collisions in the binary collision approximation at low projectile energies.


Introduction
Low energy ion irradiation and ion implantation of solids becomes increasingly important due to the advances in ultra-low energy ion implantation of two-dimensional (2D) materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16], implantation of shallow junctions [17][18][19][20] and the broad applied field of ion erosion [21][22][23][24][25]. Ion energies required for doping of 2D materials are as low a 20 eV. Implantation into 2D materials like graphene or metal dichalcogenides (e.g. MoS 2 or MoSe 2 ) at ultra-low energies results in ion ranges of only few atomic layers, accompanied by backscattering of ions and recoils from the underlying substrate and significant intermixing of the atomic layers of a 2D material. Typical implantation fluences between 10 14 and 10 15 ions cm −2 cause significant intermixing of the 2D layers as well as damage formation. Simulation of implantation profiles in 2D materials requires an accurate simulation of low energy atomic collisions, as well as a realistic representation of the complete target structure, (2D materials on a substrate) and its dynamical evolution. A well-established, versatile and fast computational approach for simulating ion-solid interactions are Monte Carlo (MC) simulation codes based on the binary collision approximation (BCA) [26][27][28][29]. Some popular BCA simulation programs [26,27] were designed for simulation of higher energy ion irradiation, the typical case for implantation doping of bulk semiconductors, and fail to correctly predict ultra-low energy ion implantation or ion erosion [30]. BCA programs usually assume an amorphous homogeneous or layered target structure, where collision details are selected by random numbers. This is well justified if the implantation fluences lead to significant damage and intermixing. Furthermore, graphene or metal dichalcogenides are rather hollow structures with empty spaces in the center of the hexagon structures. Multilayer graphene has a large layer spacing of 3.35 Å but the carbon layer thickness itself is only about 1-1.5 Å. It is therefore necessary to take the hollow volumes into account, for example by introducing filler vacancies. Vacancies as target components have not yet been implemented in BCA simulation code, except for our novel IMINTDYN program [31].
For a successful and reliable application of MC-BCA simulation codes to low and very low energy ion irradiation one has to consider three important aspects, which will be addressed in this publication and are provided by the novel BCA code IMINTDYN [31]. These aspects are (i) the dynamic vacancy generation and annihilation with a vacancy appearing as a real target component rather than a 'counter', (ii) the implementation of simultaneous weak collisions with next-nearest target atoms in binary collisions, and (iii) the correct choice of a cutoff energy to terminate the elastic collision cascade.
The aim of this work is to evaluate the accuracy of BCA simulations of low energy ion implantation by comparison with well characterized experimental systems. Therefore, we have chosen tetrahedral amorphous carbon (ta-C) as chemically stable large area uniform target with an atomically smooth surface. As ion we have chosen W, which forms stable carbide bonds upon implantation and does not show clustering or out-diffusion. The implantation profiles of W into ta-C, implanted at low (10 keV) and ultra-low (20 eV) energies, can be analyzed with high precision and depth resolution by high resolution Rutherford backscattering spectrometry (HR-RBS). The IMINTDYN BCA code also allows to simulate high energy ion scattering spectra such as highresolution Rutherford-backscattering spectra for the simulated low energy implantation profiles. We can therefore compare HR-RBS experiments with simulations of HR-RBS. The novelty of this work is, that experiment and simulations are carried out in exactly the same sequential way (10 keV implantation followed by 20 eV implantation, followed by HR-RBS analysis). The dynamic formation of vacancies should result in local fluctuations of the atom density and should have an influence on the implantation profiles. Therefore, the experimental model case system is ideal suited to be compared with BCA simulations, in particular to evaluate the role of vacancy formation during implantation.

Comparison of different BCA simulation programs
One of the first popular MC-BCA computer codes was TRIM (Transport of ions in matter) [26]. TRIM was further developed and later renamed to SRIM (Stopping and ranges in matter) [27]. The popularity of SRIM is mainly due to its easy-to-use graphical user interface and the extensive database of electronic stopping power data. SRIM computes sequences of single binary collisions using the Ziegler-Biersack-Littmark (ZBL) interaction potential [32] and the so-called magic formula for solving the scattering integral [see [33] pages 22 ff.]. TRIM and SRIM were essentially developed to predict ion ranges in materials for medium energy ions (several 10 keV or higher). The use of BCA programs to reliably predict low energy collisions, for example occurring during sputter erosion, requires certain extensions, which will be briefly discussed in the following sub-sections.

Simultaneous weak collisions
It was recognized by Eckstein and Biersack in 1984, that at low projectile energies single binary collisions (as it is used in SRIM and TRIM) are not sufficient, because a significant fraction of weak collisions with larger impact parameter (small deflection angles and small energy transfer to target atoms) are neglected. As a consequence, the collision cascade simulation remains incomplete and influences in particular sputter yield calculations. Eckstein and Biersack developed the software code TRIM.SP [34], which takes simultaneous weak collisions with next nearest and next-to-next nearest neighboring target atoms with respect to that target atom selected for the binary collision with smallest impact parameter into account. Furthermore, TRIM.SP uses the Kr-C interaction potential [35], which turned out to be much better suited to describe low energy collision cascades and sputtering processes. The main focus of the study by Biersack and Eckstein was an adequate simulation of sputtering processes, which is dominated by very low energy collisions with energies comparable to typical surface binding energies as threshold energy for sputtering, e.g. 7.4 eV for carbon.
A closer look at range profiles calculated with different BCA software codes indicates that SRIM uses no simultaneous weak collisions, and it was shown that SRIM is not suited to reliably predict sputter yields [30]. In the paper on TRIM.SP [34] it is mentioned that one additional weak collision was used in a precursor program called TRSPCR1. More recent software codes such as TRIDYN [28,36,37] and SDTrimSP [29], which include dynamic stoichiometry changes of the target with main focus on low energy ion-solid interactions, also offer the feature of simultaneous weak collisions and offer the choice of different interaction potentials and integration methods. IMINTDYN [31], which is based on the SDTrimSP program, and also takes simultaneous weak collisions into account.

Cutoff energies
In general, BCA simulation programs treat the collision cascade approximately as a sequence of elastic collisions. One may subtract a bulk binding energy from the recoil energy after an elastic collision. Furthermore, an atom leaving the surface has to overcome a surface binding energy, typically taken from sublimation enthalpies. If a projectile or recoil atom of a collision cascade reaches an energy lower than a certain cutoff energy, it is assumed to be stopped and is incorporated at its actual position. The cutoff energy is an important simulation parameter which influences the simulated ion range profiles. The TRIDYN version prior to 2017 uses a default value for the cutoff energy (often labelled as final energy E F ) of only 0.1 eV. The recent TRIDYN version uses the minimum of the involved surface binding energies as cutoff energy E F . However, if noble gas atoms with zero surface binding energy are involved, E F will be set to a default minimum value, most likely 0.1 eV. In SDTrimSP and IMINTDYN the cutoff energy can be specified for each projectile and recoil species and is 1 eV as default or even smaller if taken from the element property tables. For the novel bulk binding energy model (BBE-model) in IMINTDYN, the cutoff energy is no longer a free parameter but is determined as 1/3 of the atomic sublimation enthalpy, or 1eV in case of a vanishing sublimation enthalpy. In a dynamic TRIM version developed by Biersack, a cutoff energy E F = 3 eV was chosen [38]. SRIM uses a built-in cutoff energy of undocumented value, but most likely in the range E F ≈ 1-3 eV. Chapter 6.2.1 in the book of Eckstein is devoted to the cutoff energy [33]. For sputter yield simulations a cutoff energy equal to the surface binding energies (up to 8 eV) can be chosen. For damage production a cutoff energy equal to a certain displacement energy (25-40 eV) is recommended. In these cases, simulation results related to the collision cascade will be wrong. The recommendations by Eckstein are essentially made to optimize the computing time to obtain a specific simulation result.

Vacancy generation and annihilation.
Up to now, vacancies appear in BCA simulation program only as counter of damage production based on the Kinchin-Pease model [39]. If a collision with a target atom takes place and a recoil is generated with recoil energy E recoil > E displ (displacement energy) or E recoil > E subl (sublimation energy), then the program has to decide if a vacancy is generated. In the simple Kinchin-Pease model such a collision simply increases a vacancy counter. This leads to output values such as generated vacancies per ion. However, a recoil is generated without creating an actual vacancy. This only leads to an undercoordinated local structure with a reduced atomic density. In the case of a dynamic relaxation, the atomic density relaxes towards the initial density and leads to a relaxation of the layer thickness and a modification or adjustment of the layer composition. IMINTDYN, however, considers vacancies as real target components which can be generated and annihilated depending on specific generation and annihilation probabilities. A vacancy occupies a specified atomic volume and increases the free flight path of a projectile by an additional atomic spacing but does not change the overall atomic density N of the material (as vacancies are considered a special atomic species) and therefore the maximum impact parameter in the BCA collision remains unchanged. The effect of vacancies in the collision cascade is therefore different than a reduced density, which only slightly increases the free flight path l = -( ) N 2 1 3 / and the maximum impact parameter. If we consider the hollow structures in 2D materials as filler vacancies, we find a filler vacancy concentration of 33% for graphene and 67% for metal dichalcogenides (each hexagon of the triple layer structure only contains three atoms). Furthermore, multilayer graphene has a layer spacing of 3.35 Å but the carbon layer thickness itself is only about 1Å. Therefore, a multilayer graphene target should be represented by a carbon layer, separated by to layers of filler vacancies. For simulation of 2D materials, we can now introduce filler vacancies to represent the hollow volumes in the hexagon atomic rings as well as vacancy spacer layers in between individual graphene layers.
The concept of vacancies as target component can also be applied to our model-case target amorphous carbon. The structural and electronic properties of amorphous carbon materials have been extensively studied theoretically as well as experimentally [40][41][42][43]. In the case of tetrahedral amorphous carbon with an atomic density of about 1.5·10 23 at cm −3 and a high fraction of fourfold coordinated sp 3 bonds, irradiation does not create vacancies similar to a crystalline material. However, the bond structure of atoms around a recoiled missing C atom will rearrange and most likely form sp 2 bonds, generating a local vacancy-like region with lower density. Therefore, the concept of dynamic vacancy or void formation can also be applied to amorphous carbon materials.

Experimental
We have implanted ta-C substrates (ta-C films on Si) with W ions at 10 keV ion energy and normal ion incidence and subsequently 20 eV ions at normal ion incidence using our low energy ion beam deposition system [4,5,44]. The nominal W fluence, based on the beam current measurement, was 1·10 15 /cm 2 for each energy. Due to deceleration of the ions towards the target at very low energies, the secondary electron suppression is very efficient and the implantation fluence for the 20 eV W implantation can be measured rather accurate. For the 10 keV W implantation, the implanted fluence is about 25% lower compared to the nominal fluence due to secondary electron emission. The ta-C films on Si were grown by the Fraunhofer Institute for 'Werkstoffund Strahltechnik' (IWS) using the Laser-Arc ® deposition and exhibit a high sp 3 bond fraction and an atomic density of about 1.5·10 22 at cm −3 . The films are transparent with a purple interference color corresponding to a film thickness of 250-260 nm, which was verified by Rutherford backscattering. Atomic force microscopy measurements reveal an atomically smooth surface with rms roughness of typically ω rms ≈ 0.5 nm.
As reference sample we also prepared a sputter-deposited W film on a ta-C film on Si with a thickness of approximately 10 nm. This film reveals the position of the backscattering energy of He + in the HR-RBS spectra, backscattered from W at the surface of the film.
HR-RBS analyses were done at our 500 keV ion accelerator [45] using He + ions at an energy of 446 keV, a beam size of 1 mm diameter and a cylindrical electrostatic analyzer (ESA) with 300 mm central radius, 6 mm plate spacing and 90°sector angle [46]. The experimental chamber is shown in figure 1. The backscattering angle of the system is 127°and the nominal energy resolution is ΔE/E ≈ 0.5% or 2 keV at 400 keV backscattering energy. The total path length of about 55 cm through the cylindrical capacitor between target and ESA detector was evacuated to 2·10 -5 Pa to make energy loss and energy straggling of the He ions passing through the detector negligible. A line shaped 50 mm × 1mm Si surface barrier detector serves as an ion counter behind the electrostatic analyzer, which can detect He ions with energies above about 100 keV. The analyzer was operated with triangular oscillating plate voltages with oscillation period of 10 s. The two plate voltages U + = −0.98·U − were chosen slightly asymmetric, so that zero potential is in the center of the capacitor plates. The circuit diagram of the ESA system is shown in figure 2. A signal in the semiconductor detector behind the analyzer exit slits is amplified and converted to a trigger signal using a single channel analyzer. This is used to trigger the gate of a gated amplifier which selects the momentary control voltage of the analyzer voltage supply and converts it to a signal for pulse height analysis digitization using a Canberra Multiport II analog-to-digital converter into 2048 energy bins covering an energy range of about 560 keV. In this way a HR-RBS spectrum is recorded sequentially by scanning periodically a given analyzer voltage range (corresponding to a range of backscattering energies).
Fluctuations of the beam current are sufficiently averaged in this way. A characteristic feature of an ESA energy spectrum is that each energy bin covers an energy range given by the energy resolution of the cylindrical analyzer (as mentioned in our case ΔE/E ≈ 0.5%) . In our case the bin width is 257 eV for 2048 bins, but the analyzer resolution is 2 keV at 400 keV or only 1 keV at 200 keV. Therefore, a given backscattering event at 400 keV may be recorded in one out of 8 adjacent energy bins. This will be considered in the simulation of ESA HR-RBS spectra.
The energy calibration of the accelerator was also done using the electrostatic analyzer He + ions and a Au film on Si as substrate. First, the perfect linear relation between function generator voltage and the gate output amplitude was verified by manually varying the function generator voltage. For a fixed function generator voltage events are recorded in just one single bin of the 2048 energy bins. Then the position of the Au surface edge in the RBS spectra was determined as a function of different control voltages of the high voltage power supply of the 500 kV accelerator. In this way a precise energy calibration of the 500 kV high voltage power supply was achieved and the energy for the HR-RBS experiment was set to 446 keV (corresponding to 7.5000 V control voltage).
HR-RBS spectra of backscattered He + ions were recorded with a beam current of 1-10 μA for typically 10 min of data acquisition time. Higher irradiation fluences are avoided because a slight degradation of the irradiated samples sets in.

Simulation of vacancy generation and annihilation
The vacancy was added in IMINTDYN as a new target element with atomic number Z = 0, Mass A = 0 and an atomic volume based on the local target atomic density. The atomic density of a vacancy can be specified in the input script file, or is calculated based on the local target atomic density. A vacancy cannot act as a projectile and cannot be recoiled. Vacancies can be easily used as additional target components in static simulations. If a certain vacancy concentration c Vac is assumed, the overall mass density of the target is reduced to a value of Vac 0 and the atomic density of massive particles is = - Vac 0 The overall atomic density N including the vacancies remains unchanged. Based on N we calculate the mean free flight path between two collisions (equation (2)). Since a vacancy just extends this mean free flight path, the straight path of projectiles is increased by a factor of two for a fraction c Vac of all collisions. The maximum impact parameter for collisions with massive particles is calculated based on the atomic density of massive particle only, and is given using equation This ensures that collisions with massive particles, in particular simultaneous weak collisions (see next chapter) with more distant target atoms are correctly taken into account. The path of a projectile 'hitting' a vacancy continues with small deflection angle only determined by weak distant collisions. In contrast, just reducing the mass density and/or the atomic density N only leads to a slightly larger value of the maximum impact parameter P max and a slightly increased free flight path, because both depend on N −1/3 . Therefore, inserting vacancies have a different effect than reducing the density, because a vacancy has a much stronger influence on the free flight path of projectiles. The formation of a vacancy usually requires an activation energy for vacancy formation. In addition, vacancies may diffuse out or immediately annihilate with nearby interstitial atoms. The final situation whether a vacancy is generated or not is hard to predict. In the IMINTDYN software, a vacancy can be either defined as a static component of a target (filler vacancy), or vacancies can be created or annihilated dynamically. In a dynamic simulation, we tried to implement the dynamic formation and annihilation of vacancies within the process of a collision cascade. The energy needed to create a vacancy by recoiling a target atom is supplied by the projectile kinetic energy.
In our approach, a vacancy is generated with a certain probability p vac . In the bulk binding energy model in IMINTDYN [31] this leads to a recoil kinetic energy reduced by the sublimation energy of the target atom. The probability for vacancy formation p vac is calculated based on a general formation probability q vac (between 0 and 1) and the local vacancy concentration c vac (between 0 and 1). The probability for vacancy formation is given by selectable functions for each massive target species. Some functions reduce the probability p vac with increasing c vac , others let p vac increase towards a value 1 with increasing c vac . A suitable function for our system ta-C which rapidly stabilized the vacancy concentration towards a steady state turns out to be = ⋅ - If a projectile or recoil is stopped within the target, then the stopped particle may annihilate with a vacancy. This annihilation of course must depend on the local vacancy concentrations c Vac . We introduce a bond coordination number n coord to count the nearest neighbors of an arbitrary atom site. If one of these neighbor atoms of a particle stopping position is a vacancy, then the vacancy is annihilated. Therefore, the probability of vacancy annihilation is given as For our simulation we use q vac = 0.7 and n coord = 4, which results in a rapidly developing steady state vacancy concentration of about 14% for irradiation of ta-C with 10 keV W ions.
Vacancies are not expected to be stable close to the surface of a target. In order to avoid a too high vacancy concentration close to the surface or even a vacancy containing layer at the surface in a dynamic simulation, we have implemented an error function profile, which reduced the maximum allowable vacancy concentration towards the surface. If, in a dynamic simulation, the concentration of vacancies exceeds the value given by the error function profile, the excess vacancies are annihilated. For our simulations the target is split into ≈ 1Å thick layers and the error profile is centered at layer 6 with a width of 3 layers.

Simultaneous weak collisions and cutoff energies
A detailed description of the physics of ion-solid interactions which is used in the BCA computer simulations can be found in the seminal book of Eckstein from 1991 [33]. In 1984, Biersack and Eckstein discuss simulation of sputtering, including simultaneous weak collisions, using the newly developed BCA software TRIM.SP [34]. In all BCA simulation codes, the motion of the incident projectile and also the recoiling target atoms are treated in the following way: Each particle first propagates along a mean free path of length λ (see equation (2)). Next, the particle encounters the next elastic collision, with a collision partner selected by random numbers and an impact parameter between zero and a maximum value P max , (see equation (1)) also determined by random numbers. This process is then repeated, until the particle energy is lowered to the cutoff energy (see below) and the particle comes to rest.
This choice ensures that we have only one target atom in a cylindrical volume V (called the inner ring cylinder) of length λ and cross section area pP max However, for small values P max one neglects collisions with nearby atoms with larger impact parameter and only consider a limited range of scattering and recoil angles. This can lead to an underestimated energy transfer to target atoms, artificially forward directed collisions, larger ion ranges and larger depth of recoil distributions.
To take simultaneous collisions into account, atoms are randomly chosen in ring cylinders of equal volume around the central inner ring cylinder volume given by equation (6). As shown in figure 3, each ring cylinder contains one atom. The volumes of these ring cylinders, including the inner ring cylindrical volume, are then given as lp The number of additional ring cylinders (n max > 0), which have to be applied, has to be checked individually for each ion-target system chosen. In TRIDYN the default value is n max = 1, SDTrimSP and IMINTDYN use the default value n max = 2 for incident ions as projectiles as well as for recoils as projectiles. Eckstein recommended a value n max = 3 (page 93 of [33]). SDTrimSP and IMINTDYN use a threshold energy based on the continuous inelastic loss cross section S LS proposed by Lindhard and Scharff (see equation 5.2.1 of chapter 5.2 of [33]) below which simultaneous weak collisions treatment is activated. The inverse of S LS has a dependence on projectile atomic number and mass (Z 1 , A 1 ), target mass A 2 and projectile energy E 0 and scales with the amount of energy transfer to recoils. Typically, based on the value of S LS , simultaneous weak collision treatment is activated for particle energies below about 10 keV. It becomes relevant for heavy projectiles, light target atoms and dense targets (small value P max ).
For the case of n > 0 a projectile actually experiences (n + 1) collisions with target atoms in one BCA collisions. The total energy loss of a projectile will be therefore larger compared to the case n = 0. As example we consider low energy W (M 1 = 184 amu) projectiles and C (M 2 = 12 amu) target atoms. The maximum scattering angle Θ max is obtained from evaluating the energy and momentum conservation in elastic collisions. Heavy projectiles cannot only be forward scattered with a maximum possible scattering angle of (equation The scattering angle as function of impact parameter for W projectiles and C targets and for the screened Kr-C interaction potential is shown in figure 4. For 100 eV projectiles of mass M 1 , the impact parameters have values between 0 and P max with scattering angles between 0 and up to 3.7°. The kinematic factor K, describing the ratio of projectile energies before and after a collision, is obtained from equation 3. The minus sign holds for M 1 > M 2 . We find K(0°) = 0.77 and K(Θ max ) = 0.88. As average we find a kinematic factor of 〈K〉 ≈ 0.84. For the 1st and 2nd weak collision, the scattering angles are 3.7°−3.0°and 3.0°−2.0°with kinematic factor of about K 1 ≈ 0.83 and K 2 ≈ 0.80. The total energy loss of the BCA collision with n = 2 simultaneous weak collisions is then K eff = 〈K〉·K 1 ·K 2 ≈ 0.55. Without simultaneous weak collisions (n = 0) we have a reduced the energy transfer to target atoms, a stronger forward directed scattering and the projectiles reaches 10% of its initial energy after 13-14 collisions. In contrast, for n = 2 10% of the initial energy is reached after only 4 collisions. Between each strongly forward directed collision is a mean free path of about 1.9 Å in ta-C. Therefore, if simultaneous weak collisions are neglected and in addition a very small cutoff energy is chosen, then the ion range of low energy ions increases drastically.
To demonstrate this effect, we have also carried out simulations for different cutoff energies and different numbers of simultaneous weak collisions for 10 keV W in ta-C at normal incidence. Figure 5 shows the average number of collisions of W projectiles in a collision cascade until the atoms come to rest. Whereas for 1 or 2 simultaneous weak collisions and sufficiently large cutoff energy the average number of collisions is about 44-45, the value increases to 60-63 for low cutoff energies and neglected simultaneous weak collisions. The projectile range is therefore increased by about 2.5-3 nm, which is a relative increase of about 30% IMINTDYN can optionally use the bulk binding energy approach for cutoff energies (1/3 of the sublimation energy) also for simulations with the surface binding energy model. Cutoff energies below 1 eV, which could lead to artificially increased ion ranges should be avoided.

Electronic stopping
In BCA simulations using the SRIM code electronic stopping is based on the SRIM-2013 data base derived from fits to experimental stopping data for H and He ions in various target materials. Electronic stopping data are provided as analytical fit functions. For low energy projectiles (less than about 10 keV) the stopping functions are  extrapolated values. SDTrimSP and TRIDYN use the empirical stopping models introduced by Lindhard and Scharff (LS) [50] and also Oen and Robinson (OR) [51], often described as non-local and local stopping models. These models only apply to low energies since they describe a stopping power proportional to the projectile velocity. Recent versions of SDTrimSP allow to select the Ziegler-Biersack (ZB) stopping (see chapter 5.2 in [33]), which is a precursor of the SRIM-2013 stopping data base, and is described in SRIM documentation book [52]. IMINTDYN provides the choice of all stopping models including the SRIM-2013 electronic stopping. For our implantation simulations we use the SRIM-2013 stopping power data. For the HR-RBS simulation we also use the SRIM-2013 stopping data for He ions. The SRIM-2013 stopping data were extracted using a PYTHON script reading the SRModule application provided with the SRIM software [52] and 92 stopping tables for all projectile elements and all target elements and energies up to 2 GeV were created. IMINTDYN uses these tables and calculates stopping powers for specific energies by interpolation of tabulated values.

Simulation of HR-RBS spectra
IMINTDYN has the capability of simulating various types of ion scattering configurations with high energy MeV light ions (Elastic recoil detection (ERD), elastic backscattering (EBS), even with non-Rutherford cross sections, coincidence-ERD). For this we implemented the features 'enforce-scattering' and 'free_flight_path' and some more related simulation options. In our work we simulated ion implantation of W with first 10 keV and then 20 eV ion energy and the resulting implantation profiles were then used as input for the HR-RBS simulation.
'enforce-scattering' strongly reduces the impact parameter for selected collisions so that scattering at desired large angles is forced to occur. The position within the target where enforced scattering occurs is chosen by random numbers. Each projectile undergoes just one enforced scattering event, so that we automatically follow multiple scattering processes of the in-and out-going projectiles. To keep track of the cross sections and local concentrations, the output data of scattered and recoiled particles contain corresponding weight factors. For our case of 127°backscattering of 446 keV He ions the impact parameter range is reduced from about 0-1 Å (0 -P max ) to about 0-10-3 Å. The individual impact parameter is determined using a random number 0 r 1 with = P P r.
max Since the scattering cross section scales with the square of the impact parameter, the reduced range of impact parameters leads to an amplification of the large angle scattering probability of 10 6 compared to a regular BCA simulation. The program allows to select a range of primary scattering angles as well as a range of exit angles and does book-keeping of small, medium and large angle collisions. In this way contributions due to dual scattering and multiple scattering events can be evaluated. In our simulations we consider a scattering angle > 90°as huge angle scattering, scattering angles between 5°and 90°as medium angle scattering and scattering angles below 5°as small angle scattering. The simulated RBS spectra shown in figures 6, 10 and 11 consist of 1 single huge angle scattering event (single scattering). The events with 1 huge angle and 1 medium scattering angle (dual scattering) are very rare and can be neglected.
The feature 'free_flight_path' is similar to the option 'ion distribution and quick calculation of damage' in the SRIM software. IMINTDYN allows to increase the mean free path between two subsequent collisions based on the maximum allowed deflection angles or the maximum allowed relative energy transfer. Typically, we use a maximum deflection angle of 1°and a maximum relative energy transfer <1%. The actual mean free path between two collisions is then smeared out by ±25% using random numbers to avoid oscillatory interferences of the energy loss increments between subsequent collisions and the energy binning of the recorded spectra. For the present simulation with backscattering close to the surface, we do not need to use the 'free_flight_path' option, but it becomes important for MeV incident light ions and targets with thickness exceeding 1 μm.
IMINTDYN automatically determines the accessible depth range of a scattering experiment, the required range of reduced impact parameters for enforced scattering and calculates the cross-section and concentration dependent weighting factors. As mentioned before, the program also takes energy and angular straggling into account. A post processing program then analyzes the data files of output data of scattered and recoiled emitted particles and generates the scattering spectra for single, dual and multiple scattering. For the simulation of an energy spectrum recorded with an electrostatic analyzer, the analyzer energy resolution is specified, and a simulated backscattering event at a given energy is randomly sorted into the accessible energy bins given by the ESA energy resolution, so that each backscattering event is only counted once.

Results and discussion
As reference sample we first analyzed a W-film sputter-deposited on ta-C with an estimated thickness of 10 nm. The HR-RBS spectrum is shown in figure 6 together with two simulated spectra. The simulations were done with IMINDYN, assuming (a) a 10 nm thick W-film shown with gray symbols and (b) a film with a thin oxide surface layer and a thickness inhomogeneity of 8.1 ± 1.2 nm shown as open symbols. The simulations include energy straggling, multiple scattering and the energy resolution ΔE/E = 0.5% of the electrostatic analyzer. The surface edge for W is at 416 keV, which is expected for backscattering of 446 keV He at W with average atomic mass 184 and kinematic factor K = 0.9326 @127°. We have also simulated the spectra for the different W isotopes, but the relative variation of the masses is very small, so that isotope effects can be neglected (see figure 6). The experimental spectrum has a significant smaller slope at the W surface edge and also a broadening around 400keV. Both can be modeled with a simulation assuming an oxide surface layer consisting of 1.2 nm WO 3 + 1.2 nm WO as well as a thickness inhomogeneity of 8.1 ± 1.2 nm.
The experimental HR-RBS spectra for the W-implanted ta-C samples is shown in figure 7 together with the experimental HR-RBS spectrum of the W-film on ta-C. The two peaks of the implantation profiles, corresponding to implantation with 10 keV and with 20 eV, are clearly resolved. The profile for 20 eV appearing close to 416 keV is rather sharp and indicates that the profile is very close to the surface and no further broadening e.g. due to oxidation has occurred. The profile for 10 keV implantation appears at 400 keV backscattering energy and has a Gaussian shape with a width of σ = 2.3 keV of FWHM = 5.42 keV. The integral counts of the shown profiles are 2780 counts for the 20 eV profile and 2200 counts for the 10 keV profile. Including a correction of 1.05 taking the different He ion energies before backscattering at W close to the surface and at a depth of the 10 keV W profile into account, we find that the implantation fluence of the 10 keV implantation is 75% compared to the 20 eV implantation.
A HR-RBS spectrum of the W implanted ta-C film recorded for a broad range of energies is shown in figure 8. Besides the He + backscattering signal originating from the two W implantation profiles around 400 keV we observe a duplicated spectrum at about 200 keV. This spectrum is due to backscattered He 2+ ions, which pass the analyzer at half the analyzer voltage. The effective charge state of He in C at 400 keV is Z eff = 1.43, so we expect a significant contribution of He 2+ backscattering. We can also identify backscattering of He + from Carbon at about 120-150 keV. The cutoff below 120 keV is mainly due to the dead layer of the Si surface barrier detector used in the Analyzer.
We now compare the experimental HR-RBS spectrum with IMINTDYN simulations for five different simulation conditions and target structures listed in table 1. The simulations were done using the surface binding energy model (SBE-model) as well as the bulk binding energy model (BBE-model), with binding energies taken from elemental sublimation enthalpies. Furthermore, we use the SRIM-2013 electronic stopping data, Gauss-Legendre integration, the Kr-C interaction potential and 2 simultaneous weak collisions. Displacement energies and bulk binding energies were set to zero in SBE-model simulation mode. For each implantation simulation we used 1.2·10 5 incident W ions and a dynamic development of target Figure 6. Experimental HR-RBS spectrum (black data points) for a W film on ta-C measured with 446 keV He + ions @ 127°. The experimental spectrum is less rectangular compared to the simulated spectrum for a 10 nm film (blue data points). A good agreement is obtained for simulation assuming a thin oxide layer at the surface and a thickness variation of ±1.2 nm (red data points). The surface edge in the experimental spectrum agrees well with the calibrated and simulated value of 416 keV. The contributions from different W isotopes to the simulated HR-RBS spectrum of the W-film are also indicated as small dotted lines.
composition. The targets which were implanted with 10 keV W in a first step were then used as input for the subsequent 20 eV W ion implantation. Simulation condition (e) was also done in dynamic mode but otherwise it resembles a SRIM-like simulation with ZBL interaction potential, magic formula integration and no simultaneous weak collisions. Simulation done with the SRIM program are static simulations with displacement energies and bulk binding energies set to zero.  The simulated depth profiles are shown in figure 9 for simulation conditions (a)-(d). The results for the SRIM-like simulation condition (e) are listed in table 2. Also shown in table 2 are results obtained for the BBE model, which do not significantly differ from the ones obtained using the SBE model. The last two rows in table 2 contain the result of simulations with SRIM-2013 using zero values for the displacement energy and the bulk binding energy. In all cases we only find sputter erosion of C atoms for the 10 keV W implantation. Otherwise the sputter yields are zero.
The simulation results exhibit a large variation of the 10 keV ion range between 7.8 nm and 13 nm with a width σ varying between 0.9 and 1.6 nm. The ion range for 20 eV implantation varies between 0.25 nm and 2.4 nm almost by a factor of 10. SRIM simulations and SRIM-like simulations (condition (e)) produce comparable results for the 10 keV implantation. Because of neglected simultaneous weak collisions, the 20 eV profiles are much deeper compared to condition (b) and (c).
For a graphite-like a-C target we find the 10 keV Gaussian profile at position 10.5 nm and width σ = 1.3 nm. The maximum of the W concentration in this profile is about 3 at%. The profile for 20 eV appears close to the surface at position 0.32 nm and width σ = 0.135 nm. The maximum W concentration reaches 24 at% so that we can expect the formation of tungsten carbide, because W atoms come to rest in a carbon environment, the formation enthalpy is slightly exothermic (≈ −40 kJ mol −1 e −1 ), and carbide formation has been observed in similar earlier work [53,54]. For ta-C with higher atomic density the 10 keV W profile shifts towards the surface with position 7.8 nm and width σ = 0.95 nm. The profile for 20 eV is at position 0.25 nm and width σ = 0.13 nm. Here, the maximum W concentration reaches 20 at%.
If we allow the formation of vacancies, e.g. for simulation condition (c), we find the 10 keV W profile much deeper at position 12 nm with a larger width σ = 1.6 nm. The maximum W concentration for this profile reaches about 2 at%. For simulation condition (c) we obtain an almost constant steady state vacancy concentration of 14% up to a depth of 14 nm, except close to the surface where the vacancy concentration drops to 4%-8% due to the implemented error function profile and due to efficient annihilation of vacancies by the subsequently implanted 20 eV W ions. The profile for 20 eV is at position 0.38 nm with width σ = 0.29 nm and a maximum W concentration of 12 at%.
If we additionally suppress simultaneous weak collisions and reduce the cutoff energy for the W projectiles to 0.1 eV, which corresponds to simulation condition (d), we find the 10 keV W profile at 12.9 nm and a rather narrow width of σ = 1 nm. The maximum W concentration is 3 at%. Interestingly, the C sputter yield is  significantly higher compared to simulation condition (c). This is mainly due to the missing energy transfer in weak collisions, resulting in a larger number of created C recoils per W projectile (1938 compared to 1426) and a lower average energy per sputtered recoil (20 eV compared to 25 eV). Compared to the other simulation conditions, the 20 eV W profile appears at a significantly larger depth of 2.37 nm with a width of σ = 0.35 nm. The maximum W concentration reaches 8 at%. The neglected simultaneous weak collisions in combination with a too low cutoff energy leads to an artificially increased ion range of more than 600%. The simulated HR-RBS spectra are now compared with the experimental spectra in figures 10 and 11. For this comparison, W implantation fluences of 7.5·10 14 /cm 2 and 1·10 15 W cm −2 were used in the simulations for the 10 keV and 20 eV implantations, respectively. For the backscattering simulations we used 6·10 6 446 keV He + projectiles incident normal to the surface and detection at 127°backscattering angle. The simulated HR-RBS spectra for a-C and ta-C overlap almost perfectly. Although the ion ranges of 10.5 nm and 7.8 nm are quite different due to the different atomic densities, the electronic energy loss of the He projectiles also scales with the density of the target. Therefore, the backscattering energies become independent of the density of the carbon target material and the HR-RBS spectra look identical. The simulated 10 keV W profiles for a-C (condition (a)) and ta-C (condition (b)) appear at 404 keV with width FWHM = 3.53 keV, significantly higher in energy and narrower than the experimental profile at 400 keV. The 20 eV profiles fit very well to the experimental data.
The 10 keV profile simulated for SRIM conditions (condition (e)) appears at about 402.6 keV and width FWHM = 3.30 keV. The shift to lower backscattering energies (or larger implantation depth) occurs because SRIM neglects simultaneous weak collisions. This is also seen for the 20 eV profile and SRIM-like simulation conditions, which appears at a too low energy of 414 keV. The simulated profiles under SRIM conditions are also in poor agreement with the experimental ones.
The simulated profiles for implantation into ta-C including vacancy generation, but neglecting simultaneous weak collisions and reducing the cutoff energy to 0.1eV (condition (d)) exhibit some remarkable features. The 10 keV profile appears at 399 keV with a width of FWHM = 3.06 keV, which is significantly narrower than the experimental profile. Moreover, the profile for 20 eV implantation appears at a too low energy of 412.9 keV and has little overlap with the experimental 20 eV profile.
Finally, a very good agreement between experiment and simulation is obtained for implantation into ta-C allowing simultaneous weak collisions and dynamic vacancy generation (condition (c)), resulting in a steady state vacancy concentration of 14% ( figure 11). The simulated profiles for 10 keV W ions appear at 400 keV with a width of FWHM = 5.18 keV, which is very close to the experimental profile position and shape (FWHM = 5.42 keV). The simulated 20 eV implantation profile at position 415 keV and FWHM 2.2 keV is also in very good agreement with the experiment, indicating that the simulation correctly reflects the energy resolution of the electrostatic analyzer and that no oxidation of W implanted with 20 eV close to the surface has occurred. Figure 12 shows the almost perfect linear relation between depth of backscattering event versus He + Figure 10. Comparison of the experimental HR-RBS spectrum with simulations (colored lines) for different conditions. All of these the simulations are in poor quantitative agreement with the experimental data.
backscattering energy, in addition to effects of limited ESA energy resolution and energy straggling. Events around 415 keV correspond to backscattering of W implanted with 20 eV, the broad distribution around 400 keV to backscattering at W implanted with 10 keV. Such a plot is not available using common analysis tools such as SIMNRA [55].
These results demonstrate the strong influence of vacancies on the ion range and range distribution, which cannot be modeled by simply reducing the target atomic density. Figure 11. Comparison of the experimental HR-RBS spectrum (black data points) and a simulation (red data points) for ion implantation into ta-C with initial carbon density of 1.5·10 23 at cm −3 , allowing dynamic vacancy generation. The steady state vacancy concentration is rapidly established and is about 14%. The carbon atomic density is 1.3·10 23 at cm −3 , much higher than the atomic density of graphite. The analysis of HR-RBS spectra can be done using popular software codes like SIMNRA, which provide easy and fast analytical approximation and fitting to experimental spectra. However, the interpretation remains difficult and often ambiguous. Figure 13 shows a quantitative perfect fit of the experimental HR-RBS spectrum obtained with SIMNRA by manually defining six different layers with varying thickness and varying W content. Each layer is then fitted to optimize agreement with the experimental spectrum. The result of the analysis is shown in figure 14. The center position of the 6 layers, and the error bar crosses were manually extracted from  SIMNRA. The thickness is given as area density in units of 10 15 at cm −2 . For the 10 keV W profile we find a center position of 1.55·10 17 at cm −2 , corresponding to an ion range of 13.8 nm for a-C and 10.3 nm for ta-C. Here we have to assume a certain target atomic density to derive ion ranges. According to the BCA simulation (conditions (a) and (b)) we expect 10.5 nm or 8 nm, a deviation of about 30%. SIMNRA can only provide a rough estimate of the width of the profile. For the 20 eV W profile we obtain an ion range of about 0.5 nm assuming a-C and 0.34 nm assuming ta-C, which is about 60% larger than the BCA prediction for conditions (a) and (b)). It is obvious that SIMNRA provides a perfect fit, but the interpretation remains ambiguous. In the worst case one would conclude that the BCA simulations have a large error. Moreover, the mainly manual SMINRA analysis of the spectra and the construction of the graph in figure 15 is time consuming and slower than the computation time of all simulations with IMINTDYN.
In table 2 we have also shown the results for SRIM-2013 simulations. It is interesting to note that the 20 eV profiles from SRIM has a segmented structure due to discrete free flight path steps. This is shown in figure 15 for 20 eV W implantation into a-C. Moreover, although the sputter yield from SRIM for 10 keV W implantation has a reasonable value (1.45 is close to our SRIM-like simulation of 1.67 for ta-C), a closer look, however, shows that the angular distribution of sputtered carbon atoms derived with SRIM is strongly peaked in direction of the surface normal, an artefact of SRIM which has been recognized earlier [30].

Conclusions
We have developed the versatile Monte-Carlo binary collision approximation simulation code IMINTDYN capable of (a) predicting the implantation profiles for low and ultra-low energy ion implantation and (b) calculating high energy light ion scattering spectra of the implantation profiles to compare with experimental scattering spectra obtained for experimentally implanted samples. A novel feature of IMINTDYN is the availability of vacancies as target components and the dynamic generation and annihilation of vacancies during evolution of a collision cascade. IMINTDYN provides a dynamic treatment of vacancies as target components for the first time. In particular the use of vacancies as special target components in IMINTDYN allows a much more precise simulation of ultra-low energy ion implantation into 2D materials such as graphene and 2-dimensional metal dichalcogenides. As model case system we have chosen the subsequent implantation of 10 keV and 20 eV W ions into tetrahedral amorphous carbon (ta-C) targets as well as a sputtered W film on ta-C as reference sample. The experimental film profiles and implantation profiles were analyzed with Rutherford backscattering spectrometry of 446 KeV He + ions using a high-resolution electrostatic analyzer (HR-RBS). Ion implantation was simulated Figure 15. Implantation depth profile for 20 eV W ions incident into ta-C, calculated with SRIM-2013. The profile has a segmented structure due to discrete interaction increment of about 1.9 Å, the average free flight path in ta-C. Displacement energy E D and Bulk binding energy E B were set to zero. for different simulation conditions and different carbon target materials. Steady state vacancy concentrations of 10%-20% strongly increase the ion ranges and the width of an implantation profile. This effect is much more pronounced compared to a comparable reduction of the target density. We find excellent quantitative agreement between simulation and experiment for the implantation profiles in a ta-C sample with dynamically generated steady state vacancy concentration of 14%. The simulated HR-RBS profiles for 10 keV W ions appear at 400 keV with a width of FWHM = 5.18 keV, which is very close to the experimental profile at position 400 keV and width of FWHM = 5.42 keV. Simulated profiles for ta-C and lower density graphitic a-C without vacancies are not in agreement with the experiment. We have also shown that simultaneous weak collisions are necessary to correctly predict the implantation profiles. Neglecting simultaneous weak collisions, as it is the case in SRIM simulations, leads to significantly larger ion ranges in particular at ultra-low ion energies. The comparison between experiment and simulation for the W-film on ta-C reveals a thin surface oxide as well a small thickness inhomogeneity both of about 1-2 nm. Besides, IMINTDYN can easily handle target element isotopes. For sputter simulations IMIMTDYN offers the bulk binding energy model as an alternative approach to the commonly used surface binding energy model [31]. For the present simulations, we find excellent agreement between both models for predicting the ion range distributions as well as sputtering yields.
IMINTDYN has novel features such as 'enforce_scattering' and 'free-flight_path' which makes it possible to efficiently simulate various ion scattering spectra using high energy (even MeV) ions. The simulation allows to evaluate dual and multiple scattering contributions, energy and angular straggling, as well as relations between scattering depth and backscattering energy. The simulation can be done for previously simulated ion irradiated targets with complex target structures and requires short a computing time comparable to the time needed to acquire the spectra experimentally. We have simulated the backscattering from the W implanted samples using 446 keV He + ions, a backscattering angle of 127°and energy analysis characteristic for an electrostatic analyzer with given energy resolution ΔE/E = constant.
Besides, IMINTDYN has many more features which are not discussed in this publication. IMINTDYN is capable of simulating coincidence ERDA spectra with coincident detection of forward scattered and recoiled atoms, non-Rutherford scattering spectra and even nuclear reaction analysis spectra using cross section data from the IBANDL data base [56]. Furthermore, IMINTYN is able to calculate crater function moments up to 6th order for erosion, redistribution and implantation crater functions, which are required to predict ion-induced surface pattern formation [57].