Unraveling energy loss processes of low energy heavy ions in 2D materials

Structuring of 2D materials and their heterostructures with ion beams is a challenging task, because typically low ion energies are needed to avoid damage to a substrate. In addition, at the very first monolayers of a material, ions are not yet in charge equilibrium, i.e. they may either charge up or neutralize depending on their velocity. The change in electronic structure of the ion during scattering affects the energy, which can be transferred to the recoil and therefore the energy available for defect formation. In order to make reliable use of ion beams for defect engineering of 2D materials, we present here a model for charge state and charge exchange dependent kinetic energy transfer. Our model can be applied to all ion species, ion charge states, and energies. It is especially powerful for predicting charge state dependent stopping of slow highly charged ions.Defects engineering plays a central role in manipulating the electronic, optical and structure properties of materials and for device fabrication. The authors propose and experimentally verify a method to use ion beams for defect engineering of 2D materials based on a time dependent interatomic potential.

T ailoring electronic, optical, and structural properties of 2D materials by defect engineering is a promising strategy for device fabrication 1 . Point defects with a missing lattice atom 2 , substitutional defects with foreign atoms 3 , rotational defects 4 , line defects 5 , and complex defects as a combination of the former ones 6 , constitute a rich toolbox for material property manipulation. The importance of defects was shown in context of electronic, optical, vibrational, magnetic, elastic, and chemical properties of graphene 7-10 and transition metal dichalcogenides 11 . To introduce defects deliberately and reliably, a variety of processes can be used: direct host atom displacement by electron knock-on in a transmission electron microscope 12 , wet or dry chemical etching 13,14 , as well as focused and broad ion beam irradiation 15 . The latter approach is especially efficient in introducing specific structural defects and substitutional defects with foreign atoms [16][17][18][19] . The efficiency for structural defect formation (defect per ion) increases both with ion mass and decreasing kinetic energy 20,21 , because the ion scattering cross section increases. Thus, slow heavy ions may be ideally suitable in application.
The energy transfer to the target material, governing the defect formation, is determined by the sum of three terms: elastic energy transfer to the target nucleus (nuclear stopping), inelastic energy transfer to the target electrons (electronic stopping), and potential energy deposition into the target electronic system. The latter is a result of the neutralization of the ion, while the former two are traditionally treated as independent contributions to the stopping force 22 .
Recent experimental work showed that heavy ions, especially in high initial charge states, undergo ultrafast neutralization and de-excitation within less than 10 fs when transmitted through a freestanding monolayer of graphene 23 or 1 nm thick polymeric carbon nanomembranes 24 . The neutralization is accompanied by a strong enhancement of the ion stopping 25 , which depends on both, the ion incident charge state Q in and the amount of charge exchange ΔQ = Q in − Q out . Established stopping force models [26][27][28] cannot describe this charge state enhanced stopping and only recently attempts to improve this were made 29 . However, the largely enhanced kinetic energy losses can not be predicted by any existing reliable theory, yet.
Here we show that charge exchange between the target material and the ion, which drives the neutralization, also affects the scattering potential of the ion-target system and consequently influences elastic and inelastic energy transfer. In fact, it leads to a close coupling between nuclear and electronic losses such that also nuclear energy transfer becomes inelastic and both terms cannot be treated independently anymore. The effect of the timedependent potential alteration due to charge exchange is especially important in the very first monolayer at a surface and thus in 2D materials. For deeper surface layers the ion accommodates a constant (apart from small fluctuations) mean charge state, i.e. charge equilibrium is reached quickly. In a 2D material the ion is typically in the pre-equilibrium regime of charge exchange, which received only little attention so far. In this work we propose a Time-Dependent Potential (TDPot) model, which gives the total and the inelastic nuclear energy losses of slow ions impacting on 2D materials. We distinguish both as: the energy lost by the projectile (total energy loss) and the energy received by the target recoil (inelastic nuclear energy loss). The difference is dissipated into the electronic subsystem. Nuclear losses are typically the ig. 1 Maps of some quantities extracted from Time Dependent Potential (TDPot) simulations for 40 keV Xe 40+ ions on a graphene single layer. a Values of the distance of closest approach R min . b Final charge of the projectile Q out . c Total energy loss. Possible excitation of C atoms is not considered here. a-c are for scattering into Ω = 4π and d-f for δ = 1.6°detector acceptance angle dominating contribution to atom displacement and therefore structural defect creation. The key point of our model is that the scattering of ion and target is determined by the screened Coulomb potential, which becomes time-dependent due to charge exchange. Also the details of level population and their dynamical change due to de-excitation processes in the ion and excitations of the target are included. Our de-excitation model is based on Interatomic Coulombic Decay (ICD) 30,31 quenching an excited atom/ion upon collision with a surface. The distance-dependent de-excitation rate γ is taken from literature where-ever possible to reduce the number of free parameters of our model. In fact, our model depends only on one unknown free parameter (the target excitation parameter ξ). If we do not consider target atom excitation during collision, the only uncertainty of the model is given by the extrapolation of γ from literature values. On the one hand, ongoing research into ICD in a collisional system will reduce the uncertainty of our model results. On the other hand, our model can be used to gauge ICD rates in these systems by comparing with experimental data from ion scattering difficult to obtain with other methods. The TDPot model predicts ion stopping for all types of 2D materials, ion species, charge states, and kinetic energies. It also gives an atomistic picture of impact-parameter dependent charge exchange, i.e. it predicts the exit charge state of an ion transmitted through a 2D material. Since only slow ions contribute significantly to atom displacement, we focus on keV heavy ions in the following.

Results
40 keV Xe 40+ ions on graphene. We use our TDPot model (see Methods section) to investigate the extreme case of highly charged Xe ions transmitted through a freestanding single layer of graphene, because it nicely illustrates the importance of charge exchange for ion stopping. As an example, we consider 40 keV Xe 40+ ions, because experimental data exists for this case. Figure 1 shows the values of the distance of closest approach R min between the ion and a carbon atoms in the graphene lattice, exit projectile charge Q out and total energy loss for each impact position of 40 keV Xe 40+ ions on a graphene layer. All parameters used in the calculation are depicted in Table 1. The positions of the C atoms are in the middle of each triangular unit cell (18 unit cells are visible in the plot). As expected, ions impinging with small impact parameter on a C atom approach it to smaller minimal distances (see Fig. 1a). Also charge exchange is larger, i.e. Q out is smaller at close encounters. Additionally, the energy loss is impact-parameter dependent and also highest for central collisions. To compare with experiments, one should scan over all impact parameters within a unit cell in Fig. 1b and c to determine the histograms of outgoing charge state and energy loss. The histograms yield the mean energy loss and mean outgoing charge state. Since each impact parameter is also associated with a scattering angle, the detector acceptance angle δ needs to be taken into account. The white parts of the plots in Fig. 1d-f correspond to cases where the Xe ions are not detected since their angular deflections are larger than the angular acceptance, here used as δ = 1.6°. Then, the maximum electronic energy loss values correspond to regions between two C neighbors.
Our models is based on the time-dependent change of the number of electrons bound to the ion, which is expressed by three values: N core , the number of frozen core electrons; N cap (t), the number of electrons captured into highly excited Rydberg states; and N stab (t), the number of electrons stabilized at the ion after deexcitation. The values are determined by the rate Eqs. (6)- (8). Figure 2a shows the results for the solution of the rate equations using the values of λ and γ given by Eqs. (9) and (10), respectively, for 40 keV Xe 40+ on graphene (see Methods section). N core = 14 is constant. N cap has a fast rise at about −1.2 nm (graphene layer at 0 nm), increasing by 40 electrons. From this point the Xe atom is neutralized as a hollow atom. Close to the membrane, the stabilization phase starts. Depending on the impact parameter p, a different number of electrons captured in high-Rydberg states decays to low lying empty states of the projectile (stabilization). As observed in Fig. 2a for p = 0.07 nm, 22 electrons are stabilized maintaining 8 electrons in high-Rydberg states. These electrons are assumed to leave the atom on the way towards the detector via intra-atomic Auger processes (μs time of flight) and therefore the final charge of the projectile Q out will be close to Q in − N stab (∞) = N cap (∞). The changed number of electrons after interaction alters the scattering potential accordingly, see Fig. 2b. Consequently, the kinetic energy of the ion is lower on the outgoing trajectory and the kinetic energy difference is inelastic (Fig. 2b is for the center-of-mass).
The results for the total and nuclear energy losses are displayed in Fig. 3a for 40 keV Xe 40+ ions impinging on a single layer of graphene as a function of the final ion charge state Q out . They are calculated from the velocities of the Xe and C atoms after a distance of 50 atomic units from the graphene membrane. Our model is based on the time-dependent scattering potential and charge exchange, which is a highly dynamic process, where level population of the projectile is important. We assume classicalover-barrier charge transport above the surface forming a hollow atom 32 . Subsequently, captured electrons de-excite into atomic states with higher binding energy and thus form a compact lowly charged ion (i.e. ion with Q out < Q in in the ground state) 33 (see Fig. 2). On the one hand, the result of the calculated energy loss is insensitive to the chosen classical-over-barrier transport rate, because as a resonant process its rate is much larger than the inverse of the time the ion spends above the surface (we consider slow ions). On the other hand, the electron stabilization phase, i.e. the time scale of the hollow-atom collapse is of fundamental importance. This time scale is taken by the ICD rate γ as described in the Methods section. The ICD rate γ is strongly atom-atom distance dependent (mainly γ ∝ 1/R 6 ) 30,34-36 and consequently highest for smallest R min in our case. A change of this quantity from Eq. (10) by a factor of 2 affects essentially the value of the exit charge-state Q out by 8-10 charge units. We take the ICD rate γ from literature, where time-resolved experiments of ICD in Ne 2 and Ar 2 dimers and theoretical calculations for a variety of diatomic molecules were performed (see Methods section). Thus, it is not a free parameter in our model. The sensitivity of our results to the ICD rate is an advantage, because ion transmission experiments together with the TDPot model can be used to extract the ICD rate in collisional systems not easily accessible with other methods.
To evaluate the influence of the ICD rate γ on our results quantitatively, we used our extrapolation function (see Eq. (10)) and varied it by a factor 1/2 and 2, respectively (see Fig. 3a). The max , μ, R c for the formation of the hollow atom, γ for the stabilization process, α for the target excitation and δ for the opening angle of the detector calculated energy loss varies by a factor of 2-3. Due to lack of better knowledge, we use γ as obtained from an extrapolation from literature values without correction factor in the following. In Fig. 3a, we can also observe that both energy loss terms decrease for increasing exit charge Q out . For the electronic energy loss (difference between total and nuclear loss) this is a consequence of the sudden change in the potential energy of the system when the ions cross the membrane and therefore the electronic energy loss should be (mainly) proportional to where Q in is the initial projectile charge, Z 2 is the atomic number of the target, R min the distance of closest approach and ϕ is a screening function (see Methods section). However, we would expect a different behavior for the nuclear energy loss since the momentum transfer to the target atoms is proportional to the interaction strength and thus it should increase with Q out (Q on average is higher). Nevertheless, there is a strong correlation between Q out and the distance of closest approach R min as can be observed in Fig. 1a, b. Because of the stabilization phase, smaller values of Q out arise from smaller values of R min and therefore larger nuclear energy loss. Hence, for larger exit charges the nuclear stopping and corresponding damage will be smaller. To shed more light on the energy loss processes, it should be noted that the kinetic energy is not conserved in a scattering event where the potential is explicitly time dependent. However, the momentum is conserved. Consequently, we solve Newton's equations of motion and define the nuclear energy loss as the kinetic energy of all target atoms after scattering. The electronic energy loss is therefore the difference between projectile kinetic energy after scattering and the nuclear energy loss. This energy is contained in the change of the internal energy of the system (projectile and target material) and is reflected by charge redistribution and excitation. The latter is also explicitly taken into account by a dynamical change of the target atom screening in our model (see below, c.f. stretching parameter ξ). In a physical picture, which goes beyond results of our present model, target excitation may result in electron excitation and emission (ionization) and subsequently in material damage due to Coulomb explosion or more complex damage formation scenarios (electronic sputtering, defects induced by electronic transitions (DIET) 37 ). It also changes optical and electronic properties transiently due to a large excitation of electrons from the valence to the conduction band of the material.
Charge-state dependent energy-loss in graphene. The energy loss of slow Xe ions in graphene was measured recently 23 as a function of the incident and exit Xe charge states. The Ion position (nm) experimental energy loss results increase with increasing charge exchange ΔQ = Q in − Q out . This trend is corroborated by the present calculations and was not predicted by any other energy loss model before.
Our TDPot calculations underestimate the experimental values by about a factor of 2 if we neglect target atom excitation. In fact, due to charge exchange, a large number of electrons is drawn locally from graphene followed by an ultrafast neutralization of the C atoms nearby (in-plane current) 23 . During the collision time, these C atoms will not be in the ground-state anymore as required by the typical statistical models used to describe the interatomic potential. The easiest way to describe an excited potential in the framework of statistical models is the introduction of a stretching parameter ξ, which stretches the screening length of the (excited) scattering partners. Details are described in the Methods section.
The effect of excited C atoms is shown in Fig. 3b for the 40 keV Xe 40+ ions on graphene as a function of Q out . The results from Fig. 3a for non-excited C are also shown in this figure for comparison. The excitation of the C atoms strongly increases both energy losses and maintains the general trend of decreasing energy loss for increasing exit charge. The excitation of the C atoms before the impact on graphene and fast quenching afterwards lead to a strong change of the time-dependent potential and thus to a larger inelastic energy loss. The same holds true for the nuclear energy loss, where the reduced screening before the collision enlarges the momentum transfer to the C atoms. The process of ICD is responsible for the rapid quenching of the hollow atom and subsequent excitation of the target material. Thus, ICD causes directly (i) target electron excitation by transfer of the projectile potential energy and indirectly (ii) charge redistribution due to a stretching of the target screening, which acts on the scattering potential and subsequently adds to nuclear energy loss (to a minor extent) and electronic energy loss. Certainly, other de-excitation channels are present and in competition with ICD. Intra-atomic (within the projectile alone) non-radiative and radiative de-excitation processes would not directly cause target excitation. This would only be facilitated by secondary processes, like electron-electron scattering. The importance of target excitation during scattering is shown in Fig. 3. Further, the importance of other competing processes was already ruled out based on comparison of the process lifetimes/rates 31 . Other inter-atomic processes also contribute only to a minor extend for a highly asymmetric scattering system like Xe-C, where (quasi-)molecular orbital formation is suppressed. As a consequence, side-feeding of electrons directly into low-lying projectile states is absent, because states of C and highly excited Xe are not in resonance. Similarly penning ionization and exchange transition mediated decay (ETMD 38 ), direct Auger neutralization are not important. The picture may change for more symmetric projectile-target systems. Target excitation and charge redistribution from these processes may be similar to the ones from ICD. However, we can conclude from our results here, that the dominant inter-atomic deexcitation mechanism must follow an inter-atomic distance (R) dependence similar to Eq. (10) from ICD. Note, that at R ≪ 1 Å for similar mass scattering partners, the distinction between ICD, ETMD and penning ionization may blur out.
The results of the energy loss and exit charge abundance for an excited graphene layer are present in Fig. 4a, b. The magnitude of the experimental energy loss values is well reproduced by the calculations, but charge exchange is somewhat underestimated for small Q in and energy loss is underestimated for large Q in ≳ 30. A larger discrepancy occurs for Q in = 20, which may be due to an uncertainty in the measurement. A scatter in the measured data in similar magnitude can also be seen for the Q in = 30 data. Note, that data acquisition in experiment takes about 24-72 h per incident charge state, where small fluctuations in beam energy and focus may occur and influence the results. Since all parameters but one (α or ξ, respectively) used in the calculations are either taken from literature (γ(R)) or show only little influence on the results within a reasonable range of variation (λ, R c , ϕ, μ, β), we see this as evidence for (1) a strong excitation of the graphene layer due to charge extraction such that target excitation always need to be considered, (2) charge extraction appears even when the projectile is still above the surface, and thus before the actual impact, and (3) ultrafast interatomic electron excitation energy release from the projectile without recharging, i.e. ICD during the collision. Discrepancies between our model results and experimental values could in the future be reduced by treating α and γ(R) as entirely free parameters.

Discussions
Our results show that the kinetic energy loss of ions is significantly enhanced by the charge state. Electron capture of the ion reduces the charge state during collision and therefore reduces the charge state effect of the total energy loss. However, the energy sharing between nuclear losses (energy transfer to target atoms) and electronic losses (energy transfer to target electrons) is altered by charge exchange. In case of strong charge exchange, electronic losses are increased more than nuclear losses (see Fig. 3b). Note, that in cases where charge exchange is important, also the potential energy (internal energy) of the ion is released contributing additionally to target damage. Heavy ions in higher charge states ≥1 should therefore be used in experiments in order to increase defect creation efficiency (defect per ion). It is pointed out that the electronic excitations due to kinetic energy of the projectile also take place but are of minor importance for the example depicted in Fig. 4a. They can be estimated by using the transport cross-section approach (TCS) 26,28 and amount to few tens of eV for neutral Xe projectiles or few hundreds of eV for a hollow Xe ion.
The most sensitive parameter in TDPot is the stabilization rate γ taken from ICD models. The sensitivity to the ICD rate can be used to determine the interatomic deexcitation rate of neutralized ions in a collisional system from transmission experiments, not directly accessible with other methods. It should be emphasized that the rate γ is a deexcitation rate of the ion, where it was shown that deexcitation is mainly driven by non-radiative processes (such as ICD). One should be aware, that at close collisions, the term "interatomic" in ICD may loose its meaning, even though multi-electron dynamics is still active leading to the quenching of the hollow atom.
The microscopic picture given by our TDPot model allows us also to interpret measured data from ion transmission through 2D materials, where highly differential data is acquired. Ions are detected under a certain scattering angle with a charge state Q and a kinetic energy E, significantly different from their initial values before transmission. Our model can predict further angleresolved charge exchange and energy loss measurements (charge exchange dependent scattering pattern). When comparing measured data to models, it is imperative to consider the measurement geometry, i.e. detector acceptance angle, which is already included in TDPot. Based on TDPot, we can also calculate differential cross sections, and therefore probabilities for charge exchange, which we do not discuss in this paper.
It should be emphasized here, that for 2D materials multiple scattering is of minor importance and therefore it exists a strong correlation between R min , Q out and ΔE. When thicker targets are considered, this correlation gets lost. Still, TDPot could be implemented in a Binary Collision Approximation simulation for few-layer 2D materials where the charge state of the ion is still higher than its equilibrium value. TDPot in its present form could further be useful for thick materials if swift ions are used, where charge exchange at the surface proceeds in the opposite direction, i.e. ions charge up. Due to the high ion velocities in this case, a charge state far from equilibrium will be maintained for several nanometers.
TDPot is an easy to handle, low-computational cost model to predict energy deposition of ions in 2D materials. It goes well beyond standard software packages like SRIM and its results show clearly that charge exchange of keV ions typically used in sputtering experiments of 2D materials is significantly enhanced. With TDPot we finally understand the origin of this enhancement and have a microscopic picture of ion charge exchange in the 2D lattice. TDPot further enables a comprehensive understanding of ion beam spectroscopy measurements with wide applicability in techniques like Low and Medium Energy Ion Scattering, High-Resolution Rutherford Backscattering Spectroscopy and Elastic Recoil Detection where quantification of charge fractionization is important.

Methods
Time-dependent interatomic potential. Here we describe the model for the electronic and nuclear energy losses of slow heavy ions on graphene using the concept of the time-dependent interatomic potential. The input of this model is the time-dependent number of electrons bound to the projectile N 1 (t) = N core + N cap (t) + N stab (t), which is separated into 3 parts. The first one is the number of unchanged or frozen electrons during the collision (N core ), which populate the inner-shells of the ion. The second part is the number of captured electrons (N cap ) from the graphene into highly excited Rydberg states in the ion. Part of these electrons (N stab ) are stabilized via Auger decay or via extremely fast ICD 31 .
The time-dependent interaction potential proposed here is given by (atomic units are used throughout this section) which is the straightforward generalization of the potential used to calculate the charge-dependent energy loss in elastic collisions 29 . It is based on the statistical description of an ion with a positive nucleus (charge Z 1 ) and surrounding electron cloud (with N 1 electrons) interacting with a neutral target atom (charge Z 2 ). The distance between the two nuclei is R; ϕ and ϕ cap hollow are screening functions discussed below. The first and second terms of Eq. (1) describe the interaction of the screened part of the projectile's nuclear charge (Z 1 − Q(t) = N 1 (t)) with the neutral target atom. The last term is the remaining interaction with the total (unscreened) ion charge Q(t) = Z 1 − N 1 (t) with the (initially) neutral target atom. Different statistical models (e.g. Moliere 39,40 , Ziegler-Biersack-Littmark (ZBL) 41 , Krypton-Carbon KrC 42 ) can be utilized for the screening function ϕ(x) in connection with corresponding screening lengths 41 . For energies larger than few tens of keV, the difference between them is of minor importance and therefore we adopted here the KrC potential to describe the interaction between a heavy ion and the C atoms from graphene. The corresponding screening lengths are given by 29 where N 2 = Z 2 for a a neutral target atom. We advance the statistical description of atoms in order to take electrons properly into account which occupy levels with large principal quantum numbers n. A neutral atom with many empty inner shells, and consequently many electrons in high-n Rydberg states is called a hollow atom and forms due to classical-overbarrier charge transport in front of a solid surface 32 . The weakly bound electrons N cap lead to a different screening described by the function ϕ cap hollow ðR; r 0 Þ. It is obtained by calculating the interaction of the target atom with the projectile captured electrons, which are distributed over a spherical shell of radius r 0 and density ρ cap ðRÞ ¼ N cap =ð4πr 2 0 ÞδðR À r 0 Þ and reads ϕ cap hollow ðRÞ ¼ X i a i e Àb i R=a 2 ðZ 2 Þ À sinhðb i r < =a 2 ðZ 2 ÞÞ b i r < =a 2 ðZ 2 Þ e Àb i r > =a 2 ðZ 2 Þ R r > ; ð4Þ where r > = max(R, r 0 ), r < = min(R, r 0 ). The coefficients a i and b i are taken from the KrC potential 42 ( P i a i e Àb i x ¼ ϕ KrC ðxÞ). The value of r 0 is taken to be the critical capture distance for first electron transfer above the surface, R c . It is a function of Q in , and could be determined from the classical-over-barrier model directly 32 . To describe the experiments from Gruber et al. 23 we used the following expression determined from TDDFT calculations 23 and differing only little from the result of Burgdörfer et al. 32 .
Projectile electrons. The time-dependent potential Eq. (1) depends on the fixed number of core electrons N core and the number of captured (N cap (t)) as well as stabilized electrons N stab (t), which explicitly depend on the time. These numbers are obtained from the solution of the following coupled rate equations, which include the formation and collapse of the hollow atom as dN cap ðtÞ dt ¼ λðtÞ Z À N core À N cap ðtÞ À N stab ðtÞ À γðtÞN cap ðtÞ ð6Þ and very close to the membrane they are stabilized by ICD 31 . This rapid deexcitation of the hollow atom takes place with the rate γ.
According to the over-barrier-model 32 , the capture rate λ should be constant for distances smaller than the critical capture distance R c and therefore λ is taken as the rectangular profile where λ max is the maximum value of λ, erf(x) is the error function, v p is the projectile velocity and μ is related to the rise and fall time (steepness of the onset of the capture rate). R c is given by Eq. (5). The results do not depend very much on the parameters λ max and μ, and therefore details of the formation of the hollow atom before the impact with the membrane are of minor importance. Consequently we choose λ max = μ = 1 (see Table 1). However, the ICD rate γ is important for the results of our model. As the ICD line width is unknown for interatomic distances of less than about 2 Å, we estimate its value from extrapolating literature values rather than treating it as a free parameter. Figure 5 shows the line width of ICD for various dimers found in literature. Most of the calculations are based on ab initio calculations solving Schrödinger's equation for the multi-electron system. The CaNe and MgNe system was calculated using the most reliable Green's function method treating the manybody system. Some experimental results are also available and agree with the calculations within a factor of 10. So, in general the trends are correct, but the absolute agreement is not satisfactory. Further, for our extrapolation to small R, we get a consistent rate of about 500-4000 meV (for R = 0−1.4 Å) for CaNe, MgNe, and Ar 2 . In case of HeNe a different functional form is found in literature and therefore we used a slightly different extrapolation. All extrapolations are based on exponential R-dependencies as this fits literature values best and it can clearly be seen in Fig. 5 that none of the literature results follow the simple 1/R 6 scaling over the entire range. However, we decided to use a 1/R n (n = 8) functional form, to keep a simple expression for the ICD rate, which may in future be compared to ab initio calculations. As it turns out, ICD in a collisional process is dominated by the rate at small R anyway (where it is roughly constant) and therefore the exact R dependence is of minor importance. Our model here is very sensitive to the ICD rate (see Fig. 3a) and extrapolation of literature values points in the right direction. We want to emphasize, that the ICD rate γ and the stretching parameter ξ related to the target excitation (see below) are the only uncertain ingredients to our model. We see this as an opportunity to determine the ICD rate in a collisional system for the first time, by varying γ with a scaling factor leaving the shape constant. Current models for ICD cannot easily treat interatomic distances below a few Å, also because electronic states may loose their atomic character. Experimental determination would demand a truly time-resolved ion scattering experiment with fs timing resolution not feasible today.
We used the following expression for the ICD rate with R % ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðv p tÞ 2 þ p 2 q ), where p is the impact parameter relative to the C atom at the center of the unit cell. Note, that a R −6 dependence, as typically assumed, does not fit literature data well and it was further shown, that the ICD rate may increase even stronger than R −6 for small R if we go beyond the virtual photon model 34 .
Monte Carlo calculations. The nuclear and electronic energy losses were calculated by simulating the ion and target classical trajectories. Random numbers are used to select the projectile impact parameters inside a Wigner-Seitz cell, where each C atom is located at the center and the nearest first neighbors are outside. The effect of thermal vibrations can also be simulated by randomizing the initial target positions. After selecting the position of the ion and target C atoms, the corresponding classical trajectories are determined deterministically by solving the equations of motion through the Runge-Kutta method.
The ions impinge a single layer of graphene under normal incidence with an initial velocity v p and charge Q in . The Wigner-Seitz cell is an equilateral triangle, whose apothem corresponds to the half of the smallest distance d between two C atoms in graphene (d = 1.42 Å). The area of the unit cell is A ¼ n À1 A , where n A = 3.87 × 10 15 C cm −2 . Besides the main interaction between the ion and the C atom at the center of the triangle, the interactions between the ion and the first 3 C neighbor atoms were also taken into account. Therefore, 5 equations of motion were solved (ion + 4 C atoms) numerically using the timedependent potential Eq. (1) coupled with the rate Eqs. (6)- (8). Specifically, the resulting force on the ion is calculated considering the interaction with the first four nearest C atoms. The force on each of these C atoms, which are initially at rest, is exerted from the projectile. Of course the sum of all forces is zero and hence the total momentum is conserved. The sum of all C kinetic energies is ascribed as the nuclear energy loss. The difference of kinetic energies of the projectile before and after the impact is the total energy loss. For a time independent potential (the case of a Xe with a frozen charge) the total energy loss of the projectile is equal to the nuclear energy loss and therefore the inelastic energy loss is zero. Our model calculates the inelastic energy loss upon a variation of the projectile charge state.
For each ion trajectory, we obtain the total energy loss ΔE and the nuclear energy loss ΔE n , where ΔE was determined from the final energy of the projectile and ΔE n from the kinetic energy given to C atoms. Since the potential is timedependent, the total mechanical energy is not conserved and therefore ΔE − ΔE n corresponds to the electronic energy loss. The final projectile charge state Q out was determined from the simulation as shown in Fig. 2b according to Only projectile trajectories having angular deviations smaller than the detector acceptance of 1.6°were used to generate the data of energy loss as a function of Q out . The results from Eq. (11) are used as mean values for a Gaussian distribution of FWHM = 3 unit charges.
Excited potential. The excitation of the target atoms is modelled by changing the screening length between the projectile and the target atom. This screening length depends on the size of the electron cloud around the target atom and hence it increases for an excited electron density. Therefore, the screening length should be stretched when the target electrons are excited. Thus, the effect of an excited C atom can be easily included in Eq. (1) by changing the screening lengths a 1 and a 2 according to a 1 ðN; Z 2 Þ ! a 1 ðN; Z 2 =ξÞ a 2 ðZ 2 Þ ! a 2 ðZ 2 =ξÞ ð12Þ with a time dependent stretching parameter ξ given by where α is the parameter to describe the maximum stretching effect and β is the parameter responsible for the switching time, i.e. the rate of excitation. The excitation of the target atoms is asymmetric in time (before and after the interaction) because of the stabilization (quenching) process. In the beginning the excitation is small. When ICD sets in, target electrons are excited and at the late stage of the interaction, excitation is still large but the ion is (almost) in its ground state and the potential is consequently screened efficiently irrespective of target electron excitation. The "switching time" is taken from the stabilization time, which depends on the rates λ and γ. The maximum stretching parameter was obtained from a best fitting of the experimental data. Here we used α = 1-25, and β = 2v p , which corresponds to a switching time of about 0.4 fs for 40 keV Xe ions.

Data availability
The data sets generated during and/or analysed during the current study are available from the authors on reasonable request.

Code availability
The custom computer code used in this study is based on standard procedures to solve Eqs. (6)- (8) including the steps outlined in the Methods section. A standard Monte Carlo approach is used to vary the impact parameter. The full code is available from the authors on reasonable request.