††thanks: currently working at the European Research Council Executive Agency
(ERCEA)

#  On the effect of electronic stopping in molecular dynamics simulations of
collision cascades in Gallium Arsenide

Johannes L. Teunissen  Univ. Lille, CNRS UMR 8520 - IEMN - Institute of
Electronics, Microelectronics and Nanotechnology, F-59000 Lille, France  Royal
Belgian Institute for Space Aeronomy, Av Circulaire 3, 1180 Brussels, Belgium
<EMAIL_ADDRESS>Thomas Jarrin  CEA, DAM, DIF, F-91297, Arpajon, France
Nicolas Richard  CEA, DAM, DIF, F-91297, Arpajon, France  Université Paris-
Saclay, CEA, Laboratoire Matière en Conditions Extrêmes, F-91680 Bruyères-le-
Châtel, France  Natalia E. Koval  CIC Nanogune BRTA, 20018 Donostia-San
Sebastián, Spain  Daniel Muñoz Santiburcio  CIC Nanogune BRTA, 20018 Donostia-
San Sebastián, Spain  Instituto de Fusión Nuclear “Guillermo Velarde”,
Universidad Politécnica de Madrid, c/José Gutierrez Abascal 2, 28006 Madrid,
Spain  Jorge Kohanoff  Instituto de Fusión Nuclear “Guillermo Velarde”,
Universidad Politécnica de Madrid, c/José Gutierrez Abascal 2, 28006 Madrid,
Spain  Atomistic Simulation Centre, Queen’s University Belfast, Belfast
BT71NN, Northern Ireland, United Kingdom  Emilio Artacho  CIC Nanogune BRTA,
20018 Donostia-San Sebastián, Spain  Donostia International Physics Center
DIPC, 20018 Donostia-San Sebastián, Spain  Theory of Condensed Matter,
Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United
Kingdom  Ikerbasque, Basque Foundation for Science, 48011 Bilbao, Spain
Fabrizio Cleri  Univ. Lille, CNRS UMR 8520 - IEMN - Institute of Electronics,
Microelectronics and Nanotechnology, F-59000 Lille, France  Fabiana Da Pieve
Royal Belgian Institute for Space Aeronomy, Av Circulaire 3, 1180 Brussels,
Belgium

###### Abstract

Understanding the generation and evolution of defects induced in matter by ion
irradiation is of fundamental importance to estimate the degradation of
functional properties of materials. Computational approaches used in different
communities, from space radiation effects to nuclear energy experiments, are
based on a number of approximations that, among others, traditionally neglect
the coupling between electronic and ionic degrees of freedom in the
description of displacements. In this work, we study collision cascades in
GaAs, including the electronic stopping power for self-projectiles in
different directions obtained via real time Time Dependent Density Functional
Theory in Molecular Dynamics simulations of collision cascades, using the
recent electron-phonon model and the previously developed two-temperature
model. We show that the former can be well applied to describe the effects of
electronic stopping in molecular dynamics simulations of collision cascades in
a multielement semiconductor and that the number of defects is considerably
affected by electronic stopping effects. The results are also discussed in the
wider context of the commonly used non-ionizing energy loss model to estimate
degradation of materials by cumulative displacements.

##  I Introduction

The description of radiation-induced effects in materials is a challenging
multiscale problem at both time and length scales. At very short time scales
($\sim$100 as), a particle passing through a material loses energy by
transferring it to electronic degrees of freedom of the target (electronic
stopping).

On longer time scales, predominantly when the particle has been slowed down by
the target’s electrons, the impacting particle undergoes nuclear elastic
collisions, thus displacing atoms in the target. These first displaced atoms
are called primary knock-on atoms, (PKAs). Such PKAs act as additional self-
projectiles and collide with other atoms creating a collision cascade [1] (in
the ps regime). Different types of defects are created as a result, such as
vacancies and interstitials (Frenkel pairs), and defect clusters followed by a
slow partial recombination on a longer time scale up to nanoseconds and
longer.

Several physically realistic metrics for the quantification of radiation
damage exist [2, 3, 4]. In the high-energy particle physics and space
radiation effects communities, the Non Ionizing Energy Loss (NIEL) concept is
commonly used to describe the damage due to cumulative atomic displacements.
The NIEL represents the portion of the energy of an incident particle that
goes into non-ionizing phenomena in the target, being displacement damage
effects. Apart from the high energy range, where non-ionizing phenomena are
due to inelastic nuclear scattering, in the medium-low energy range atomic
displacements nuclear stopping represents the main contribution to the NIEL,
with phonons not directly leading to a displacement also eventually
contributing. On the basis of a considerable set of experimental observations,
it is generally assumed that the degradation of output parameters of a
semiconductor device under ion-irradiation is linearly correlated with the
NIEL [5, 6, 7, 8, 9]. According to the commonly used Norgett-Robinson-Torrens
model (NRT) [10] which is based on the binary collision approximation (BCA),
there is a linear relationship between the NIEL and the number of
displacements being subsequently a linear function of the PKA energy. The NRT
model is implemented in the widely used SRIM/TRIM [11] software [12] and in
wider-purpose Monte Carlo (MC) particle transport codes (at low energies) [13,
14].

The linearity of the NIEL means that only the number of defects should give a
measure of the damage, irrespective of their nature, e.g. whether defects are
clustered in small regions (as for neutron-induced damage) or homogeneously
scattered over a relatively large volume (as for low-energy proton or
$\gamma$-ray-induced damage) [15]. This independence allows to obtain the
damage produced by different particles and with different energies via a NIEL
scaling [5, 6, 16, 7, 17, 9]. Deviations from the linear dependence of the
number of defects on the PKA energy were found via MD simulations on Si [18,
19, 20, 21] for low energy PKAs. These works [18, 19, 20, 21] also highlighted
that atomic displacements occur even for energies lower than the commonly used
displacement threshold energy, implying that not all the energy transfer is to
phonons below such energy but that also collective motion can generate some
damage [22]. The non-linearity of the NIEL was also observed in experiments
[7] and MC simulations [23]. Deviations from the linearity of the degradation
parameters and the NIEL were also found for different energy ranges [24, 25].
Efforts have been made [22, 26, 27, 28, 29, 30], to propose an effective or
adjusted NIEL model to correct the deviations from linear dependence, either
via MD studies [22, 31, 32, 33] or experimental works [17, 34].

The manner in which the energy dissipation by the PKAs along the incident ion
trajectory is determined and partitioned into energy loss to electrons and
into atomic motion is fundamental for the counting and further evolution of
defects. NIEL calculations in SRIM and wider-purposes MC codes are done
without considering any coupling between ionic and electronic degrees of
freedom, using amorphous targets and using a constant threshold displacement
energy for each element that is independent on the material in which the
element is inserted. For the nuclear stopping, the universal screening
potential (Ziegler-Biersack-Littmark, ZBL) is used in SRIM, while in wider-
purpose MC codes such as Geant4 [13] either the ZBL [35] or a modified
screening function given by the Moliere’s approximation of the Thomas-Fermi
model are used.

The electronic stopping of heavy ions is calculated via a charge scaling of
the proton stopping. The latter is calculated in SRIM on the basis of fitting
(and extrapolation) of experimental data [36, 12] using the (linear-response)
Lindhard-Sharff model. In Geant4 this model is used for low energy protons
(below $\sim$10keV/amu) only, and a parametrization based on ICRU tables and
then the Bethe Bloch formula are used for increasing energies [37]. All these
approaches, however, inherit the main approximations mentioned above.

The accuracy of such stopping power values, in particular for low energy, for
self-irradiated mono- and multi-element systems, is actually unknown and still
needs to be largely verified, given the paucity of available experimental
data. Several recent works highlight the limitations of SRIM in the prediction
of the electronic and nuclear stopping power [12, 38, 36, 39].

Real Time-Time Dependent Density Functional Theory (RT-TDDFT) represents a
promising parameter-free approach to obtain the electronic stopping in any
target for any projectile, at energies where the assumptions of SRIM can no
longer be justified. Recent studies on monoelemental systems [40, 41, 42, 43,
44, 45] showed that the inclusion of electronic effects obtained via RT-TDDFT
in MD cascade simulations affected both the number of defects and the cascade
morphology.

In this context, the recently proposed electron-phonon (EPH) model for non-
adiabatic dynamics [46, 47, 48, 41, 49, 40, 42] improves upon the previously
developed Two-Temperature model (TTM) (by Duffy et al. [50]). The EPH model
treats the electron-phonon coupling and electronic stopping as two
manifestations of the same physical process. In contrast with the TTM, the
friction term now depends on the local structural environment allowing for
better agreement with ab-initio electronic stopping data. Recent work on
cascade simulations in Si showed that the EPH model correctly reproduces the
density-dependent ESP and that the final radiation damage differs considerably
from the TTM predictions [45].

In this work, we perform an ab-initio study of collision cascades induced by
self-recoils of 100 eV, 1 keV and 10 keV, using both the EPH model and the TTM
in the paradigmatic case of GaAs, whose damage is often studied in the space
engineering community being the most active layer of state-of-the-art solar
cells used in space missions. We devise a more robust procedure as put forward
by Caro et al. [41] to fit the EPH model to the electronic stopping determined
via RT-TDDFT. A different two-tier method was used that immediately prevents
unbounded stopping and that relies on two global optimization steps to reduce
the problem of local optima. We show for the first time that the EPH model
provides a good description of the electronic stopping in the semiconductor
material GaAs. Moreover, we show that it is possible to fit the EPH model for
multiple-element materials. We additionally show that the description of the
electronic stopping and electron-phonon coupling influences the number of
defects, and is thus crucial to correctly model radiation damage. MD collision
cascade simulations are performed for PKAs with different energies in GaAs and
are analyzed in terms of number of defects and clustering. How the number of
defects depends on the counting criteria is also discussed. Lastly, the number
of defects is compared to the NRT model and previous works on GaAs, also
discussing eventual deviations of the linearity of the number of defects with
PKA energy.

This work is part of a larger scheme as envisioned by M. Raine et al.[51]
spanning multiple time and length scales to model the radiation impact. Within
this scheme the distribution of PKAs is initially evaluated from high energy
Monte Carlo simulations, resulting in PKA energy distributions for different
species. In this work the impact of the individual Ga and As PKAs is evaluated
leading to damage configurations only a few picoseconds after the PKA caused
the collisional cascade. After our work, the obtained damage configurations
are evolved up to a timescale of seconds using another Monte Carlo simulation
(k-ART) [52].

##  II Methods

This work is divided in three sections. Section A discusses the determination
of the ESP via RT-TDDFT calculations. Section B, the obtained ESP results are
used to train the EPH model. An iterative procedure is applied to optimize the
dissipation functions that are part of the EPH model to best reproduce the ab-
initio ESP. In the last section, the EPH model and TTM are used to run MD
simulations of collision cascades.

###  II.1 Calculation of the electronic stopping power

For projectiles with high velocity, the main kinetic energy loss mechanism
occurs via non-adiabatic excitations of the electronic system. This loss is
usually quantified as the ESP, $S_{e}$ that is defined as:

$S_{e}(E)=-\frac{dE_{e}(x)}{dx}$ (1)

where $dE_{e}$ is the energy loss to the electronic system, $dx$ is the
distance travelled by the projectile, i.e., the penetration depth. We compute
the ESP by simulating the passing of an ion with a constant velocity through
the target material while keeping the rest of the atoms frozen, which allows
to readily obtain $dE_{e}$ as the change of the total energy of the system
(thus neglecting the nuclear stopping component). This approximation is
completely valid for the projectile velocities investigated in this work [53,
54, 55], and it is widely employed in most of the computational works
determining the ESP by RT-TDDFT methods [41, 44, 56].

The RT-TDDFT calculations are performed with the plane-wave Qb@ll code [57,
58] using the PBE functional [59]. Previous experience of calculating
electronic stopping power using LDA and GGA (PBE) functionals show quite
quantitative agreement with experiments [44, 56, 47, 60]. Due to the high
velocities of the PKAs, leading to relatively deep-level electronic
excitations, pseudopotentials with 13 and 15 valence electrons are used for
gallium and arsenide respectively. The core electrons are especially relevant
for highly energetic projectiles [55]. The aim of our RT-TDDFT calculations is
to parametrize the EPH model for collision cascades with PKA energies well
below 100 keV.

Hamann-type norm-conserving pseudo-potentials[61] as obtained from the
quantum-simulation pseudopotential library have been used [62, 63]. The
pseudopotential cut-offs are 0.6 Å for gallium and 0.56 Å for arsenic. The
wavefunction basis cutoff was set at 100 Rydberg and a smearing was added
using 4 empty states at 1000 K. A 2x2x3 supercell of GaAs is used for all
calculations containing 96 atoms and its Brillouin zone is sampled only at the
$\Gamma$-point.

In our simulations, the wavefunctions are propagated using the enforced time
reversal symmetry propagator. In this work, the projectile, unless stated
otherwise, is always neutral, i.e., it starts in an unionized form since we
assume that the electronic environment of the PKA at the start of its
trajectory is equivalent to that of a neutral atom in its ground state. See
Supplemental Material S1 [64] for ESP evaluations of highly ionized As PKAs.
Below 1000 keV we used a time step of 0.4 attoseconds. At higher energies the
time step is chosen such that the projectile displaces 0.005 Å per simulation
step.

There is an initial transient period before the projectile charge is
stabilized. For this reason, only the results after this transient phase are
used to obtain the average ESP. Every time the projectile passes close to an
atom, the potential energy of the projectile shows a peak. See Supplemental
Material S2 [64] for detailed information about how the average ESP is
obtained.

###  II.2 Including electronic stopping in MD

In MD simulations, the nuclear stopping can be calculated relatively
straightforwardly as it corresponds to the energy loss of the projectile due
to all the interatomic collisions. However, especially at higher projectile
energies ($E>$ 100 keV), the most important energy loss mechanism is the
electronic stopping. Several methods are constructed to account for the
electronic stopping in MD simulations without taking the electrons into
account explicitly [65, 50]. One of the earlier models applied only an extra
friction force to the atoms, but the energy lost by the projectile is removed
from the system and not returned as thermal energy via electron-phonon
coupling.

The TTM allows both, the inclusion of the electronic stopping and the return
of the energy adsorbed by the electronic system back to the nuclear system
[50, 66]. TTM adds an effective electronic system to the classical atomic MD
system with its own temperature. The normal MD equation of motion is exchanged
by a Langevin-type equation by adding a friction term and a stochastic force
term. The friction term scales linearly with the atomic velocity (in agreement
with Lindhard theory [67]) and is scaled with a $\gamma_{s}$ parameter related
to the electronic stopping. The value of $\gamma_{s}$ is normally obtained
from SRIM. The energy transport in the electronic system is covered by an
electronic heat equation using the electronic specific heat and electronic
thermal conductivity.

The advantage of the TTM is that it allows to include electronic stopping and
electron-phonon interactions into MD simulations, however, there are a few
disadvantages. Firstly, it applies the electronic stopping such as a
projectile would experience in a uniform electron gas and thus effects such as
channeling are not taken into account [68]. Secondly, many parameters are
included for which there is little experimental data, making it difficult to
obtain reliable results and predictability [69].

Application of the TTM proved that the inclusion of the electronic temperature
elongates the duration of the thermal spike observed during collision
cascades, and that this enables more thermal annealing leading to decreased
defect numbers. Additionally, the TTM allows to investigate radiation
resistant materials that have large electron-phonon couplings [66]. The TTM as
implemented in LAMMPS [70] only allows to give one $\gamma_{s}$ parameter, so
for GaAs, the friction term will be the same for Ga and As projectiles.

The recently developed EPH model [48, 40, 41] improves upon the TTM by
modelling the electronic system by spherical electronic densities placed on
each atom and thus the ESP can depend on the local electron density
experienced by the projectile. The electron-phonon coupling is treated as an
electronic stopping process at low energies (meV scale). Thirdly, it was
proven that TDDFT can give accurate results for the ESP and the electron-
phonon coupling [55, 46] making it possible to parametrize the EPH model using
TDDFT results. In [41], it is shown that ESP profiles from a well fitted EPH
model closely match the RT-TDDFT results.

In the EPH model, the friction and random force are described by many-body
tensors, contrary to the TTM in which they are scalar [40]:

$m_{I}\frac{\partial\mathbf{v}_{I}}{\partial
t}=\mathbf{F}_{I}-\sum_{J}B_{IJ}\mathbf{v}_{J}+\sum_{J}W_{IJ}\bm{\xi}_{J}$ (2)

The second and third terms on the right-hand side correspond to the friction
force and to the random force acting on the $I^{th}$ atom. The forces acting
on the $I^{th}$ atom now depend on the velocity of all other atoms $J$, and
thus all forces are correlated. The set of $\bm{\xi}_{J}$ are white-noise
mutually-uncorrelated Gaussian random variables normalized to $2k_{B}T_{e}$
[41]. $B_{IJ}$ and $W_{IJ}$ are two tensors that describe the spatial
correlations and they are related via the fluctuation-dissipation theorem such
that:

$B_{IJ}=\sum_{K}W_{IK}W_{JK}^{\text{T}}$ (3)

Both the ESP and the electron-phonon coupling are determined by the terms in
the $W$ matrix.

The terms in $W$ are constructed such that the random forces do not produce
any local net translational nor rotational momentum. This condition is
motivated by the preservation of the correct phonon lifetimes and
polarization. $W$ is then defined as:

$W_{IJ}=\left\\{\begin{aligned}
&-\alpha_{J}(\bar{\rho_{J}})\frac{\rho_{I}(r_{IJ})}{\bar{\rho_{J}}}\vec{e}_{IJ}\times\vec{e}_{IJ}&\;\;(I\neq
J)\\\ &\alpha_{I}(\bar{\rho_{I}})\sum\limits_{K\neq
I}\frac{\rho_{K}(r_{IK})}{\bar{\rho_{I}}}\vec{e}_{IK}\times\vec{e}_{IK}&\;\;(I=J)\\\
\end{aligned}\right.$ (4)

where $\rho_{I}(r)$ is the value of the spherical density $\rho$ (centered at
atom $I$) at a distance $r$. $\bar{\rho_{I}}=\sum_{J\neq I}\rho_{J}(r_{IJ})$
is the electron density experienced by atom $I$ and is the sum of
contributions from all its neighbors. $\alpha_{I}(\bar{\rho_{I}})$ is a
function defining the coupling between the electrons and ions [42].

The shape of the $\alpha$ function is free and depends on the atomic species.
Thus for every element a different $\alpha$ function is used. The $\alpha$
function has to be defined over the whole relevant range, from low densities
occurring at crystal vacancies and in the center of wide channel directions,
to high densities encountered at close collisions. This function is optimized
using the RT-TDDFT data. The ESP has to be defined in a wide range of
electronic densities in order to obtain an accurate $\alpha$ function along
the whole density range.

The electronic heat equation is very similar to the TTM and the same values
for $C_{e}$ and $\kappa_{e}$ are used [42]. The EPH model is available as an
external plug-in for LAMMPS.

###  II.3 Collision Cascade MD simulations

Following the approach of Gao et al. [32], we use the bond order potential as
developed for GaAs by Albe et al. [71] combined with the ZBL functional [72]
at small distances to describe the short range repulsive interactions that
occur during strong collision impacts. (See Supplemental Material S3 [64] for
the potential visualization.)

For all MD simulations, periodic boundary conditions are applied. The total
simulation box is a large cubic supercell. All the outer cells of the
simulation supercell form a thermostat keeping the atomic temperature at 300 K
via a Nose-Hoover thermostat [73], thus absorbing the thermal wave. The EPH or
TTM equations only apply within the interior region.

For each simulation type, the system is first initiated at 300 K and
equilibrated for a sufficient time of 100 ps. Afterwards, the equilibrated
system is used as a starting point for a set of collision cascades. At the
beginning of the cascade, the time step is set to 0.001 fs. Subsequently, the
time step is updated such that no atom moves more than 0.05 Å per time step.
All simulations are run for several picoseconds until the temperature returns
back to 301 K, so the thermal spike is sufficiently dissipated. This allows us
to compare the direct effect of using different simulation types but for
experimental comparison longer time scale methods would be required to include
longer timescale effects.

To get a good statistical average, for each MD setting, 78 different cascades
are run with different initial PKA directions chosen according to the
procedure outlined in reference [69]. The PKA is chosen such that its initial
velocity vector points as much as possible towards the center of the interior
region to keep the collision damage into the box and preventing atoms from
leaving the interior during the cascade. The clustering analysis is done by
using the Louvain algorithm for community detection [74].

##  III Results and discussion

###  III.1 Calculation of the electronic stopping power

####  III.1.1 Choosing the PKA directions

Different crystal directions are chosen for the PKA trajectories such that a
variety of crystal environments and electron densities are probed to enable
accurate subsequent fitting of the EPH model. We included three full
channeling directions where the projectile travels in the middle of the
channel, thus traversing a lowest density path through the crystal: $\langle
001\rangle$, $\langle 011\rangle$ and $\langle 111\rangle$.

Off-center channeling directions, i.e., directions parallel to the channel but
closer to the atoms, are added in the $\langle 001\rangle$ channel. We adopt
the following naming scheme: $\langle 001\rangle_{X}^{N}$ indicates a
trajectory moved away from the $\langle 001\rangle$ center towards the atoms
of type $X$, located at $\frac{1}{N}$ of the atom-center distance. For the
included off-center directions see Figure 1.

Figure 1: Schematic representation of the $\langle 001\rangle$ channel. Blue
dots signify arsenic and red dots signify gallium atoms. The green dot
indicates the center channeling direction $\langle 001\rangle$. The turquoise
dots indicate the position of the $\langle 001\rangle_{Ga}^{7}$ and $\langle
001\rangle_{As}^{7}$ off-center channeling directions, located at one seventh
of the atom-center distance. The orange dots indicate the position of the
$\langle 001\rangle_{Ga}^{4}$ and $\langle 001\rangle_{As}^{4}$ off-center
channeling directions, located at a quarter of the atom-center distance.

Random trajectories that are incommensurate with any crystallographic symmetry
vector are also included. Since these directions are not periodic, they are
selected to probe each possible density environment of the crystal. These so-
called off-channel directions are chosen using several conditions:

  * •

The starting point of the PKA has to be not too close to another atom to get a
good starting wavefunction.

  * •

Due to the choice of a $2\times 2\times 3$ supercell, the PKA direction has to
have the largest velocity component in the $z$ direction.

  * •

Since the off-channel directions are sampled for a longer trajectory, the
projectile has to re-enter the unit-cell as far as possible from where it has
already been before.

  * •

The trajectory should not contain a too close encounter (frontal collision)
with another atom.

  * •

The trajectory should be as incommensurate with a crystallographic direction
as possible. To test this, we use the auto-correlation function, which tests
if there is any periodicity in a signal. We record the closest distance to
another atom as a function of the trajectory along the first 60 Å. If there
are no large peaks (except the self correlation) in the auto-correlation
function, the chosen direction is sufficiently incommensurate with any low-
index crystallographic direction.

In the EPH model, the electron-phonon interactions are treated as an
electronic stopping process that happens around the atomic reference
positions. The electronic stopping at these (low atomic density) reference
positions is sampled by including an incommensurate PKA trajectory moving
through a vacancy.

####  III.1.2 Electronic stopping power calculations

The ESP for several trajectories was evaluated for a range of different
energies. Since the results for gallium and arsenic PKAs are qualitatively
similar, in Figure 2 we show only the ESP for arsenic PKAs for simplicity.
Also the stopping power as obtained by SRIM is included.

Figure 2: Electronic stopping of Arsenic PKAs in Gallium Arsenide. (Some
points are omitted to cut the computational cost.)

The ESP values for all channeling directions are smaller than the SRIM
results, as is also observed for many other materials [56, 75, 76, 55]. We
observe that the ESPs along the centers of the three channels are very
similar. For example, for 1 MeV As PKAs, the ESP values are 26.0, 26.5 and
27.4 eV/Å for the $\langle 001\rangle$, $\langle 011\rangle$ and $\langle
111\rangle$ channels respectively. This deviates from other works on materials
with a similar diamond-like structure [75, 76] where the wider $\langle
011\rangle$ channel resulted in slightly lower ESP values than the $\langle
001\rangle$ and $\langle 111\rangle$ results.

The ESP results for the off-center channeling close to the atoms, have a very
high ESP due to the high electron densities encountered during the projectile
trajectory. The ESP for $\langle 001\rangle_{\text{As}}^{7}$ is slightly
higher than for the $\langle 001\rangle_{\text{Ga}}^{7}$, so the stopping on
the projectile passing an As atom is higher, in agreement with its higher
electron density.

The off-channeling results are in between the channeling and off-center
channeling results. This is expected as the off-channeling trajectory probes
the electron densities most evenly (time averaged) while channeling and off-
center channeling mainly probe either only low or mainly high electron
densities.

Since SRIM does not take into account channeling, i.e., it considers only a
series of binary collisions, the off-channeling result should match best with
the SRIM result. We confirm this assumption at lower energies ($\leq$ 10 MeV).
At higher projectile energies, there is a larger deviation from SRIM. This is
also observed in several other works [55, 44] and can be attributed to the
increasing importance of core electrons at higher projectile velocities. Our
calculations only explicitly treat 13/15 valence electrons of the 31/33
electrons of Ga and As respectively so the effect of core electrons is not
included.

To fit the EPH model, we focus on the results at the PKA energy of 100 keV,
since it offers the best balance between accuracy and computational time, as
going to lower energies (such as 20 keV and lower) drastically increases the
computational costs while at much higher energies, the linear relationship
between the ESP and the PKA velocity breaks down [50].

Figure 3 shows the ESP values at multiple trajectories both for Ga and As
PKAs, all with the kinetic energy of 100 keV. Some additional trajectories are
considered such as $\langle 001\rangle_{Ga}^{4}$ and $\langle
001\rangle_{As}^{4}$. Another off-channeling trajectory is included that does
not exhibit a close collision as does the offchannel1 direction around 38 Å
(See Supplemental Material S2 [64]). Last, an off-channeling trajectory moving
through a vacancy is included (vacancy1).

Figure 3: Electronic stopping values of Ga and As PKAs in GaAs along different
trajectories with PKA energy of 100 keV.

The ESPs along the off-center channeling directions $\langle
001\rangle_{X}^{4}$ are lower than along the $\langle 001\rangle_{X}^{7}$
trajectories (with $X$ being Ga or As). The ESP along the second off-channel
direction without a close collision is clearly lower than the on on the first
off-channeling direction and is very close to the SRIM result. The trajectory
moving through a vacancy yields a lower ESP than along the other off-
channeling directions.

In general, the differences in ESP between Ga and As are small and most
pronounced for the channeling directions. The ESP for As is typically larger
than for Ga. Off-center channeling trajectories close to the As atoms have
also higher ESPs than those close to the Ga atoms, in accordance with the
higher electron density around As.

###  III.2 Fit of the EPH model to the ESP results

Using the set of ESP results along the different trajectories, we fit the EPH
model to include the ESP into the MD simulations. Since it would be too costly
to run a self-consistent quantum-mechanical calculation within each MD step,
in the EPH model the electron density is approximated as the sum of the
spherical densities centered on each atom. Thus, one needs to provide the
spherical densities of Ga and As. These atom-like densities are obtained using
the OPIUM program on the isolated atoms [77] (See Supplemental Material S4
[64] for the atomic density visualisations).

####  III.2.1 Optimizing the $\beta$ function

The $B$ matrix in Equation 2 represents the friction coefficient and is the
Hermitian square of $W$ (Equation 3). For this reason, the friction scales
with $\alpha^{2}$. In the following section, we optimize this $\alpha^{2}$
function, which we call $\beta$, following the same notation as in reference
[41].

In principle, $\beta$ depends not only on the local density, but also on the
directions of the moving particles. Especially at higher densities, it can be
observed (see especially at the peaks of the middle plot of Figure 4) that at
the beginning of a close collision, when two atoms get closer to each other,
the ESP is high, while further on, when the particles move away from each
other, the ESP is lower. Nevertheless, the time-independent density profile is
symmetric. This means that the ESP cannot be expressed exactly as a function
that only depends on the time-independent density [41]. It is thus an
approximation of the EPH model to let $\beta$ be only a function of the time-
independent density, approximated by the sum of the spherical densities of
each atom called $\rho_{0}$.

In the following section, we try to find a $\beta(\rho_{0})$ function. We
define that the best $\beta$ function is the one that, when the same
trajectories are run with EPH-MD and RT-TDDFT, gives the smallest difference
in energy dissipation due to the ESP.

From the EPH-MD simulations, we directly obtain the energy dissipated by the
system due to the ESP at each point along the trajectory. To obtain the
instantaneous dissipation, or instantaneous ESP, we have to subtract the
effect of the crystal structure from the RT-TDDFT energy curves. The structure
effect can be removed by computing the energy corresponding to the ground
state (BOA) of each atomic configuration, and then subtracting it from the
energy obtained during the RT-TDDFT run at the same configuration: Now, the
energy increase caused by the non-adiabatic effects only, $E_{elec}(r)$, which
equals the dissipated energy to electronic stopping, is obtained as:

$E_{elec}(r)=E_{RT-TDDFT}(r)-E_{BOA}(r)$ (5)

In Figure 4, we show for three trajectories how the effect of the structure is
removed to get $E_{elec}$. For the $\langle 001\rangle$, the $E_{elec}$ curve
is very smooth, so the ESP is very stable along the channel. From the other
two trajectories, it can be observed that most energy is dissipated during
close collisions.

Figure 4: Three examples of how the instantaneous dissipation energy (orange)
due to electronic stopping is obtained by subtracting the BOA energy (black)
from the RT-TDDFT (red). Here, they are shown for Arsenic PKAs of 100 keV,
along the $\langle 001\rangle$, $\langle 001\rangle_{Ga}^{4}$ and offchannel1
trajectories.

We choose to define the error for a given $\beta$ function as the difference
between the EPH-MD and RT-TDDFT dissipation. For a single trajectory the mean
absolute error is defined as:

$MAE=\frac{1}{r_{12}}\int_{r_{1}}^{r_{2}}\left|E_{MD}(r)-E_{elec}(r)\right|dr$
(6)

And the total error to be minimized is the sum of all $MAE$s from all the
trajectories taken into account during the minimization.

There are different strategies to solve the challenge of finding a $\beta$
function that minimizes the total $MAE$. Generally, we assume that when there
are no electrons, there is no electronic stopping so the function has to pass
through the origin. We also expect the ESP to increase when the density
increases, but only until a certain density is reached after which the ESP
either saturates, or decreases. Different strategies can be applied. Caro et
al. [41] constructed a $\beta$ function by using a cubic spline with 6 knots
and afterwards added an exponential decay to prevent unbounded stopping at
high electron densities. Jarrin et al. [45] used a function of the form:
$\beta(\rho)=c_{1}\rho\cdot e^{c_{2}(\rho-c_{3})}$ where the $c$’s are
optimization parameters.

Both approaches have their disadvantages. The optimization space using the
splines is non-linear so the optimization problem is prone to find a local
optimum while the functional form is too restricted. We decided to combine
both approaches. First, a simple functional form is optimized to get a general
shape of the $\beta$ function. Subsequently, a second optimization is run that
adds a cubic spline function onto the earlier function. This second
optimization essentially fine-tunes the first function by making small
adjustments to the function but keeps the general function shape. The
advantage of the two-tier optimization procedure is the reduced dimensionality
of the problem. In both optimization steps, the DIRECT global derivative-free
optimization algorithm [78] is used as implemented in the non-linear
optimization software NLOPT [79]. During the optimization procedures in each
iteration, NLOPT proposes a trial $\beta$ function. Next, short LAMMPS
calculations are run for all the relevant trajectories using the trial $\beta$
function. Subsequently, the total error is evaluated.

Note that in our case, for GaAs, we have two elements and so two $\beta$
functions are optimized simultaneously. As we noted that at high densities the
ESP is smaller but does not drop to zero, instead of the exponential decay we
decided to saturate to a constant $\beta$ function at higher densities. The
use of a constant $\beta$ function is also used in [45] and [40]. (The
trajectories $\langle 001\rangle_{\text{Ga}}^{7}$ and $\langle
001\rangle_{\text{As}}^{7}$ were not used during the EPH model training as the
high density regions sampled along these trajectories are rarely encountered
during the collisional cascades and the dissipation functions will be constant
at these densities. In the first optimization step, we use the simple
functional form:

$\beta(\rho)=\left\\{\begin{aligned} &c_{1}\cdot\rho&\;\;(\rho<c_{2})\\\
&c_{1}\cdot c_{2}&\;\;(\rho>c_{2})\\\ \end{aligned}\right.$ (7)

In the second optimization step, a 4-point cubic spline is added in the range
[0, $c_{2}$] that has to be zero at both ends of the range. In Figure 5, the
optimized $\beta$ functions are displayed.

Figure 5: The optimized $\beta(\rho)$ functions for gallium and arsenic with
the density in $e^{-}/$Å3. The $\beta_{1}$ curves are the result of the first
optimization step, following Equation 7. The $\beta_{2}$ curves are the result
of the sum of $\beta_{1}$ and the splines correction curves. $\beta$ values
are in eV$\cdot$ps/Å2.

In Figure 6, the comparison of the EPH-MD and TTM-MD runs with the ab-initio
results are shown for 4 trajectories for Ga (top) and As (bottom) PKAs. The
TTM result has significant errors along most trajectories, as it does not take
into account any local crystal information. The EPH model is clearly more in
sync with the RT-TDDFT data, although the agreement could certainly be
improved by using more advanced fitting methods.

Figure 6:  Comparison of the dissipation energy due to electronic stopping
evaluated via RT-TDDFT ($E_{elec}$), EPH-MD and TTM-MD vs. distance along the
trajectory. Upper row for gallium PKAs and lower row for arsenic PKAs along
the indicated trajectories.

We decided to optimise (using NLOPT) the $\gamma_{s}$ parameters of the TTM
using the same error measure as for the EPH model to get a fair comparison
with the EPH model. The optimized $\gamma_{s}$ parameters are 45.82 and 43.34
g/mol/ps for Ga and As respectively.  The SRIM obtained $\gamma_{s}$ values
are 54.54 and 51.51 g/mol/ps for Ga and As, respectively. These values are
clearly higher than our RT-TDDFT based values and thus would result in a too
high ESP.

The total MAE for EPH is 1421 eV while for TTM it is 4680 eV. Thus with the
EPH as fitted here, a 70% reduction in the dissipation error with respect to
the TTM is obtained.

###  III.3 Collision Cascade MD simulations

In this section, we present the collision cascades obtained using the TTM and
EPH model as well as cascades without electronic stopping. In both models we
used specific heat, $C_{e}$ and thermal conductivity, $\kappa_{e}$ values of
$6.45\cdot 10^{-6}$ eV/K and $1.0\cdot 10^{-9}$ eV/K/Å/ps respectively. As
shown by Jarrin et al [45], the EPH-MD results do not strongly depend on these
electronic parameters, while for the TTM-MD results, only changing $C_{e}$ had
a large impact. In the TTM model, the critical velocity was set to 96.28 Å/ps

We ran MD simulations of collision cascades with PKA energies of 100 eV, 1 keV
and 10 keV. The MD systems were $N*N*N$ supercells with $N$ being 20, 30 or 50
conventional GaAs unit cells respectively. Above 10 keV, the cascade
structures become self-similar due to subcascade formation [80].

In Figure 7, we show how the energies evolve over time during the collision
cascade. All curves start with a peak of kinetic energy due to the
projectile’s direct transfer of kinetic energy followed by a peak of the
potential energy, the so-called heat spike. The dissipation of the heat is
much faster with TTM-MD than with EPH-MD. In the evolution of the kinetic
energy in the EPH-MD simulations, we can observe a second slight increase
around 2-3 ps, caused by the the electron-phonon coupling, where the hot
electrons dissipate heat to the ions, thus increasing again the kinetic energy
of the system. As we ran the cascade until the temperature was back at 301
Kelvin after the heat spike, running times for EPH-MD were significantly
longer than for TTM-MD. The average run time for the MD calculations with PKA
energies of 10 keV were 4 picoseconds for TTM-MD vs. 51 picoseconds for EPH-
MD.

Figure 7: Average evolution of the energies indicated in the figure during the
MD simulations. Each curve is the average of 78 MD runs and the shaded area
the standard deviation.

There are two main methods to evaluate the number of defects.

  * •

The Wigner-Seitz method [81]: The cascade is initiated on a perfect lattice.
After the cascade, every atom is assigned to its closest reference crystal
position. Each reference position that has two or more atoms for which it is
the closest position is counted as an interstitial. Every reference position
that has no atom for which it is the closest lattice position is a vacancy.
Antisites are defined as reference positions with only one affiliated atom but
of another type as the original type at that reference position.

  * •

Cut-off methods [82]: A sphere is centered on each original crystal position.
After the cascade, each sphere that is empty is labelled a vacancy while each
atom that is located outside any sphere is an interstitial. Antisites are
defined as spheres with the cut-off radius centered on the original crystal
atoms that have one atom within the sphere that is of a different element than
the original element.

Here, we choose to use the second method with a cutoff of 1.0 Å, which is
within the stable regime where the number of defects w.r.t. the cutoff radius
is most stable. See Supplemental Material S5 [64] for how the number of
defects changes w.r.t. the defect definition and cut-off value. The number of
defects and antisites determined via the cut-off method as well as the cluster
sizes and the PKA penetration depths are shown in Table 1. Three typical final
damage configurations obtained with EPH-MD for the Ga PKAs are presented in
Figure 8. As expected, higher PKA energies lead to higher numbers of defects,
antisites clusters and PKA depths.

Figure 8: Three examples of final damage configurations shown in the
$xy$-plane, ran with EPH-MD with different Ga PKA energies of 100 eV, 1 keV
and 10 keV respectively. Insterstitials are red, vacancies green and antisites
are purple. The axes dimensions are in Å. The original, non-defect crystal
atoms are only shown in proximity of the defects and depicted as small dots:
gallium (blue) and arsenic (orange).

The number of defects and antisites determined via the cut-off method as well
as the cluster sizes and the PKA penetration depths are shown in Table 1. The
results of the EPH-MD simulations for As and Ga PKAs are similar. The average
PKA penetration depths and number of clusters formed is only slightly higher
for Ga than As, being always within the respective error bars of each other.
In the EPH-MD runs without electronic stopping, the number of defects is
around 10% higher (see Table 1), so the electronic stopping effectively
removes energy from the nuclei, resulting in lower final radiation damage.

A higher number of defects, antisites, and clusters is obtained for MD with
the EPH model as compared to TTM. This indicates that with the TTM, more
energy is dissipated thanks to the ESP and thus as heat, while with EPH model,
since there are regions with lower ESP (channeling-like conditions), more
energy goes into nuclear stopping, leading to more defects. The higher
dissipation to electronic stopping in the TTM-MD case goes in the same
direction as in SRIM, which also overestimates the ESP [83, 38].  The SRIM
obtained values for $\gamma_{s}$ are higher. They would also be higher when we
would have optimized $\gamma_{s}$ using only a random off-channeling
trajectory $\gamma_{s}$. Running TTM-MD with higher $\gamma_{s}$ values would
result in a stronger ESP and the number of defects would be even lower.  Also,
our $C_{e}$ and $\kappa_{e}$ values are relatively low, but using higher
values would presumably lead to even lower number of defects, in analogy with
reference [45]. The number of clusters from the TTM-MD is clearly lower, even
when taking into account the lower amount of defects; an indication that TTM-
MD leads to larger defect clusters. This behaviour is also observed for Si
[45].

Table 1: Number of defects and antisites as determined via the cut-off method
($1$ Å) as well as cluster sizes

and PKA depths for the MD calculations for specified PKA and initial PKA
energies. All values are mean values of 78 MD simulations. Values in between
brackets are standard error of the mean (SEM) values.  EPH-MD NVE TTM-MD
$E_{0}^{PKA}$ As Ga Ga Ga N defects 100 eV 6 (2) 7 (2) - 6 (2) 1 keV 59 (11)
54 (8) 61 (8) 52 (8) 10 keV 542 (86) 542 (93) 592 (111) 464 (52) N antisites
100 eV 0 (0.5) 0 (0.5) - 0 (0.4) 1 keV 5 (2) 5 (2) 5 (2) 2 (2) 10 keV 64 (15)
64 (13) - 26 (6) N clusters 100 eV 2.9 (0.7) 3.6 (1.0) - 3.4 (0.8) 1 keV 14.8
(3.0) 16.3 (2.8) 16.7 (3.1) 13.0 (2.6) 10 keV 96.3 (13.0) 100.1 (13.3) 106.3
(14.3) 67.3 (9.2) PKA depth (Å) 100 eV 4 (3) 8 (4) - 7 (3) 1 keV 21 (12) 29
(16) 31 (17) 24 (12) 10 keV 81 (59) 89 (55) 112 (69) 98 (64)

We ran additional EPH-MD calculations on 0.5, 2, 5 and 20 keV PKAs. The
evolution of the number of defects is shown in Figure 9. A clear linear
relationship is observed for the number of interstitials and antisites. We
also included the number of defects expected from the NRT formula:
$N_{\mathrm{defects}}=\frac{0.8E_{PKA}}{2*E_{d}}$ where $E_{d}$ is the
threshold displacement energy. $E_{d}$ is taken as 13 eV, the average of the
reported values for Ga and As [84]. As one can see, the number of defects, as
obtained with the cut-off method is considerably higher than what predicted
via the NRT model, but has the same linearity.

The number of defects at 1 keV and 10 keV increases with approximately a
factor of 9-10, in agreement to what is observed in MD simulations by Gao et
al. [32] and Nordlund et al. [80]. However, contrary to our work, in the work
by Gao et al [32], the number of defects was much higher than what predicted
by the NRT model as shown in Figure 9. The authors associated this results
with the direct-impact amorphization that occurs in GaAs, which prevents the
fast annihilation of interstitials with vacancies during the relaxation phase.
The defect production efficiency, defined as the ratio of $N_{F}$ to
$N_{NRT}$, ranged from 3.0 to 5.0 and increased with increasing PKA energy,
rather than decreasing with increasing PKA energy as revealed in SiC and
metals. This suggested that the defect production efficiency in GaAs is higher
than those found in metals or alloys because of the absence of a thermal spike
in GaAs. When determined via the WS definition, our number of defects is
slightly smaller than predicted via the NRT model. In both cases, there is a
high linearity in the number of defects vs. PKA energy. The number of defects
reported in Nordlund et al. [80] for GaAs is higher than ours but much closer
to our reported values than those from Gao et al [32]. As in their case, the
number of antisites is an order of magnitude lower than the number of standard
defects (vacancies and interstitials). We note that great care should be taken
to the comparison of absolute defects numbers as there are still large
uncertainties also indicated by Nordlund et al. [80] not only due to different
electronic stopping models but also due to the use of different force field
potentials and defect detection algorithms.

As can be seen from the linear fit to the lower energy defect numbers given in
Figure 9, at higher energies, the agreement with the linear NRT model breaks
down to some extent. This can be caused by the larger amount of energy lost by
electronic stopping and also because of non-linear effects due to the start of
overlapping defect cascades.

Figure 9:  Number of defects for multiple PKA energies. Number of frenkel pair
defects as determined via the Wigner-Seitz (WS) method are compared with the
numbers reported in Gao et al. [32] and Nordlund et al. [80]. The number of
defects as predicted via the NRT model
($N_{\mathrm{defects}}=\frac{0.8E_{PKA}}{2E_{d}}$) is evaluated with
$E_{d}=13$ eV. A fit to the number of EPH-MD defects at energies up to 5 keV
is shown in grey, which has an effective $E_{d}$ of 17.1 eV.

##  IV Conclusions

In this work, we have shown that it is possible to get proper ab-initio
electronic stopping behaviour of PKAs in GaAs via RT-TDDFT calculations. It
has to be noted that while RT-TDDFT simulations are becoming more widely used,
simulations of heavy ions are anything but routine. We devised a more robust
procedure as put forward by Caro et al. [41] to fit the EPH model to the ESP
determined via RT-TDDFT. A different two-tier method was used that immediately
prevents unbounded stopping and that relies on two global optimization steps
to reduce the problem of local optima. Although theoretically based on metals,
we have shown that the EPH model provides a good description of the electronic
stopping in the semiconductor material GaAs. Moreover, we have demonstrated
the application of the EPH model for multiple elements materials.

We conclude that the inclusion of the electronic stopping and electron-phonon
coupling influences the number of defects, and is thus crucial to correctly
model the radiation damage. The temperature dissociation in the TTM-MD was
proven to be much faster, while the EPH model dissipated the heat spike much
more slowly.

In comparison with the commonly used NRT model, the cut-off method gives a
higher number of defects, while the Wigner-Seitz method gives a lower one. Our
work shows that the number of defects obtained with the EPH model differs from
the results obtained with other models and within the Wigner-Seitz definition
of defects, this number is ”attenuated”. Although our work suggests that the
realistic, non-adiabatic cascades, would not strongly affect the linear
dependence of $N_{non-adiab}$ vs. $E_{PKA}$, deviations from linearity start
to increase beyond PKA energies of 10 keV. Such non linear effects could be
caused by the larger amount of energy lost by electronic stopping and also
because of non-linear effects due to the start of overlapping defect cascades.
Non linearity effects observed in experiments could also be induced by
intrinsic defects in the target, an aspect that should be further
investigated. It has to be stressed that an overestimation of the number of
defects at the nanometer scale can affect the lifetime estimation of the
functions that the material serves to.

##  V acknowledgements

This work received funding from the Research Executive Agency under the EU’s
Horizon 2020 Research and Innovation program ESC2RAD (grant ID 776410). We
acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France
(PRACE Project: 2019215186). The author thankfully acknowledges the computer
resources at MareNostrum and the technical support provided by Barcelona
Supercomputing Center (RES-FI-2020-1-0012, RES-FI-2020-2-0021, RES-
FI-2020-3-0014, RES-FI-2021-1-0035).

Jorge Kohanoff was supported by the Beatriz Galindo Program (BEAGAL18/00130)
from the Ministerio de Educación y Formación Profesional of Spain, and by the
Comunidad de Madrid through the Convenio Plurianual with Universidad
Politécnica de Madrid in its line of action Apoyo a la realización de
proyectos de I+D para investigadores Beatriz Galindo, within the framework of
V PRICIT (V Plan Regional de Investigación Científica e Innovación
Tecnológica).

We thank Alfredo Correa from the Lawrence Livermore National Laboratory for
fruitful discussions. Funding from Spanish MICIN is also acknowledged through
grant PID2019-107338RB-C61/AEI/DOI:10.13039/501100011033 as well as a María de
Maeztu award to Nanogune, Grant CEX2020-001038-M funded by
MCIN/AEI/10.13039/501100011033.

For Fabiana Da Pieve, the information and views set out in this article are
those of the author and do not necessarily reflect the official opinion of the
ERCEA.

## References

  * Was [2017] G. S. Was, Fundamentals of Radiation Materials Science: Metals and Alloys (Springer New York, New York, NY, 2017) Chap. The Damage Cascade, pp. 131–165.
  * Nordlund  _et al._ [2018a] K. Nordlund, S. J. Zinkle, A. E. Sand, F. Granberg, R. S. Averback, R. E. Stoller, T. Suzudo, L. Malerba, F. Banhart, W. J. Weber, F. Willaime, S. L. Dudarev, and D. Simeone, Primary radiation damage: A review of current understanding and models, Journal of Nuclear Materials 512, 450 (2018a).
  * Nordlund  _et al._ [2018b] K. Nordlund, S. J. Zinkle, A. E. Sand, F. Granberg, R. S. Averback, R. Stoller, T. Suzudo, L. Malerba, F. Banhart, W. J. Weber, F. Willaime, S. L. Dudarev, and D. Simeone, Improving atomic displacement and replacement calculations with physically realistic damage models, Nature Communications 9, 1084 (2018b).
  * Griffin [2014] P. J. Griffin, Relationship between metrics used to represent displacement damage in materials 10.2172/1617630 (2014).
  * Moll [2018] M. Moll, Displacement Damage in Silicon Detectors for High Energy Physics, IEEE Transactions on Nuclear Science 65, 1561 (2018).
  * Srour [2003] J. R. Srour, Review of displacement damage effects in silicon devices, IEEE Trans. Nucl. Sci. 50 (2003).
  * Summers  _et al._ [1993] G. P. Summers, E. A. Burke, P. Shapiro, S. R. Messenger, and R. Walters, Damage correlation in semiconductors exposed to gamma, electron and proton radiations, IEEE Trans. Nucl. Sci. 40, 1372 (1993).
  * SR  _et al._ [2001] M. SR, S. GP, B. EA, W. RJ, and X. MA, Modeling solar cell degradation in space: a comparison of the NRL displacement damage dose and the JPL equivalent fluence approaches, Progress in Photovoltaics: Research and Applications 9, 103 (2001).
  * GP  _et al._ [1994] S. GP, W. RJ, X. MA, B. EA, M. SR, S. P, and S. RL., A new approach to damage prediction for solar cells exposed to different radiations, In IEEE Proceedings of 1st World Conference on Photovoltaic Energy Conversion,Waikoloa, HI , 2068 (1994).
  * Norgett  _et al._ [1975] M. Norgett, M. Robinson, and I. Torrens, A proposed method of calculating displacement dose rates, Nuclear Engineering and Design 33, 50 (1975).
  * Ziegler  _et al._ [2010a] J. F. Ziegler, M. Ziegler, and J. Biersack, SRIM – the stopping and range of ions in matter (2010), Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 268, 1818 (2010a), 19th International Conference on Ion Beam Analysis.
  * W.J.Weber and Zhang [2019] W.J.Weber and Y. Zhang, Predicting damage production in monoatomic and multi-elemental targets using stopping and range of ions in matter code: Challenges and recommendations, Current Opinion in Solid State and Materials Science 23, 100757 (2019).
  * Allison  _et al._ [2016] J. Allison, K. Amako, J. Apostolakis, P. Arce, M. Asai, T. Aso, E. Bagli, A. Bagulya, S. Banerjee, G. Barrand, B. Beck, A. Bogdanov, D. Brandt, J. Brown, H. Burkhardt, P. Canal, D. Cano-Ott, S. Chauvie, K. Cho, G. Cirrone, G. Cooperman, M. Cortés-Giraldo, G. Cosmo, G. Cuttone, G. Depaola, L. Desorgher, X. Dong, A. Dotti, V. Elvira, G. Folger, Z. Francis, A. Galoyan, L. Garnier, M. Gayer, K. Genser, V. Grichine, S. Guatelli, P. Guèye, P. Gumplinger, A. Howard, I. Hřivnáčová, S. Hwang, S. Incerti, A. Ivanchenko, V. Ivanchenko, F. Jones, S. Jun, P. Kaitaniemi, N. Karakatsanis, M. Karamitros, M. Kelsey, A. Kimura, T. Koi, H. Kurashige, A. Lechner, S. Lee, F. Longo, M. Maire, D. Mancusi, A. Mantero, E. Mendoza, B. Morgan, K. Murakami, T. Nikitina, L. Pandola, P. Paprocki, J. Perl, I. Petrović, M. Pia, W. Pokorski, J. Quesada, M. Raine, M. Reis, A. Ribon, A. R. Fira, F. Romano, G. Russo, G. Santin, T. Sasaki, D. Sawkey, J. Shin, I. Strakovsky, A. Taborda, S. Tanaka, B. Tomé, T. Toshito, H. Tran, P. Truscott, L. Urban, V. Uzhinsky, J. Verbeke, M. Verderi, B. Wendt, H. Wenzel, D. Wright, D. Wright, T. Yamashita, J. Yarba, and H. Yoshida, Recent developments in Geant4, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 835, 186 (2016).
  * Fasso  _et al._ [2011] A. Fasso, A. Ferrari, G. Smirnov, V. Vlachoudis, and F. Sommerer, Fluka realistic modeling of radiation induced damage, in  _Progress in Nuclear Science and Technology; Journal Volume: 2; Conference: SNA+MC 2010: Joint international conference of the 7th supercomputing in nuclear application and the 3rd Monte Carlo, Tokyo (Japan), 17-21 Oct 2010_ (2011) pp. 769–775.
  * Dawson  _et al._ [2021] I. Dawson, F. Faccio, M. Moll, and A. Weidberg,  _Radiation effects in the LHC experiments: Impact on detector performance and operation_, Tech. Rep. (CERN, Geneva, 2021) report CERN-2021-001.
  * Srour and Lo [2000] J. Srour and D. Lo, Universal damage factor for radiation-induced dark current in silicon devices, IEEE Transactions on Nuclear Science 47, 2451 (2000).
  * Messenger  _et al._ [2001a] S. R. Messenger, G. P. Summers, E. A. Burke, R. J. Walters, and M. A. Xapsos, Modeling solar cell degradation in space: A comparison of the NRL displacement damage dose and the JPL equivalent fluence approaches, Progress in Photovoltaics: Research and Applications 9, 103 (2001a), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pip.357 .
  * Nordlund  _et al._ [1998] K. Nordlund, M. Ghaly, R. S. Averback, M. Caturla, T. Diaz de la Rubia, and J. Tarus, Defect production in collision cascades in elemental semiconductors and fcc metals, Phys. Rev. B 57, 7556 (1998).
  * Hobler and Otto [2003] G. Hobler and G. Otto, Amorphous pocket model for silicon based on molecular dynamics simulations, Nucl. Instrum. Meth. Phys. Res. B 206, 81 (2003).
  * Santos  _et al._ [2006] I. Santos, L. A. Marqués, and L. Pelaz, Modeling of damage generation mechanisms in silicon at energies below the displacement threshold, Phys. Rev. B 74, 174115 (2006).
  * Santos  _et al._ [2007] I. Santos, L. A. Marques, L. Pelaz, and P. Lopez, Molecular dynamics study of damage generation mechanisms in silicon at the low energy regime, in  _2007 Spanish Conference on Electron Devices_ (2007) pp. 37–40.
  * Inguimbert  _et al._ [2010] C. Inguimbert, P. Arnolda, T. Nuns, and G. Rolland, “effective niel” in silicon: Calculation using molecular dynamics simulation results, IEEE Transactions on Nuclear Science 57, 1915 (2010).
  * Huhtinen [2002] M. Huhtinen, Simulation of non-ionising energy loss and defect formation in silicon (2002).
  * Moll  _et al._ [2002] M. Moll, E. Fretwurst, M. Kuhnke, and G. Lindstroem, Relation between microscopic defects and macroscopic changes in silicon detector properties after hadron irradiation, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 186, 100 (2002).
  * Sumita  _et al._ [2003] T. Sumita, M. Imaizumi, S. Matsuda, T. Ohshima, A. Ohi, and H. Itoh, Proton radiation analysis of multi-junction space solar cells, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 206, 448 (2003), 13th International Conference on Ion Beam Modification of Materials.
  * Akkerman  _et al._ [2001] A. Akkerman, J. Barak, M. Chadwick, J. Levinson, M. Murat, and Y. Lifshitz, Updated niel calculations for estimating the damage induced by particles and $\gamma$-rays in si and gaas, Radiation Physics and Chemistry 62, 301 (2001).
  * Akkerman and Barak [2006] A. Akkerman and J. Barak, New partition factor calculations for evaluating the damage of low energy ions in silicon, IEEE transactions on nuclear science 53, 3667 (2006).
  * Messenger  _et al._ [2003] S. Messenger, E. Burke, M. Xapsos, G. Summers, R. Walters, I. Jun, and T. Jordan, Niel for heavy ions: an analytical approach, IEEE Transactions on Nuclear Science 50, 1919 (2003).
  * Weller  _et al._ [2004] R. Weller, M. Mendenhall, and D. Fleetwood, A screened coulomb scattering module for displacement damage computations in geant4, IEEE Transactions on Nuclear Science 51, 3669 (2004).
  * Inguimbert and Gigante [2006] C. Inguimbert and R. Gigante, Nemo: A code to compute niel of protons, neutrons, electrons, and heavy ions, IEEE Transactions on Nuclear Science 53, 1967 (2006).
  * Inguimbert and Messenger [2012] C. Inguimbert and S. Messenger, Equivalent Displacement Damage Dose for On-Orbit Space Applications, IEEE Transactions on Nuclear Science 59, 3117 (2012).
  * Gao  _et al._ [2017] F. Gao, N. Chen, E. Hernandez-Rivera, D. Huang, and P. D. LeVan, Displacement damage and predicted non-ionizing energy loss in gaas, Journal of Applied Physics 121, 095104 (2017).
  * Salzberger  _et al._ [2017] M. Salzberger, C. Nömayr, P. Lugli, S. R. Messenger, and C. G. Zimmermann, Degradation fitting of irradiated solar cells using variable threshold energy for atomic displacement, Prog. Photovolt: Res. Appl. (2017).
  * Lu  _et al._ [2011] M. X. Lu, R. C. Wang, Y. H. Liu, W. T. Hu, Z. Feng, and Z. Han, Adjusted NIEL calculations for estimating proton-induced degradation of GaInP/GaAs/Ge space solar cells, Nuclear Instruments & Methods in Physics Research Section B-beam Interactions With Materials and Atoms 269, 1884 (2011).
  * Ziegler  _et al._ [2010b] J. F. Ziegler, M. D. Ziegler, and J. Biersack, Rim - the stopping and range of ions in matter (2010), Nucl. Instr. Meth. in Phys. Res. B 268, 1818 (2010b).
  * Wittmaack [2016] K. Wittmaack, Misconceptions impairing the validity of the stopping power tables in the srim library and suggestions for doing better in the future., Nucl. Instr. and Methods in Phys. Res. B 380, 57–70 (2016).
  * Allison  _et al._ [2006] J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce Dubois, M. Asai, G. Barrand, R. Capra, S. Chauvie, R. Chytracek, G. A. P. Cirrone, G. Cooperman, G. Cosmo, G. Cuttone, G. G. Daquino, M. Donszelmann, M. Dressel, G. Folger, F. Foppiano, J. Generowicz, V. Grichine, S. Guatelli, P. Gumplinger, A. Heikkinen, I. Hrivnacova, A. Howard, S. Incerti, V. Ivanchenko, T. Johnson, F. Jones, T. Koi, R. Kokoulin, M. Kossov, H. Kurashige, V. Lara, S. Larsson, F. Lei, O. Link, F. Longo, M. Maire, A. Mantero, B. Mascialino, I. McLaren, P. Mendez Lorenzo, K. Minamimoto, K. Murakami, P. Nieminen, L. Pandola, S. Parlati, L. Peralta, J. Perl, A. Pfeiffer, M. G. Pia, A. Ribon, P. Rodrigues, G. Russo, S. Sadilov, G. Santin, T. Sasaki, D. Smith, N. Starkov, S. Tanaka, E. Tcherniaev, B. Tome, A. Trindade, P. Truscott, L. Urban, M. Verderi, A. Walkden, J. P. Wellisch, D. C. Williams, D. Wright, and H. Yoshida, Geant4 developments and applications, IEEE Transactions on Nuclear Science 53, 270 (2006).
  * Crocombette, Jean-Paul and Van Wambeke, Christian [2019] Crocombette, Jean-Paul and Van Wambeke, Christian, Quick calculation of damage for ion irradiation: implementation in iradina and comparisons to srim, EPJ Nuclear Sci. Technol. 5, 7 (2019).
  * Stoller [2000] R. E. Stoller, The role of cascade energy and temperature in primary defect formation in iron, Journal of Nuclear Materials 276, 22 (2000).
  * Tamm  _et al._ [2018] A. Tamm, M. Caro, A. Caro, G. Samolyuk, M. Klintenberg, and A. A. Correa, Langevin dynamics with spatial correlations as a model for electron-phonon coupling, Phys. Rev. Lett. 120, 185501 (2018).
  * Caro  _et al._ [2019] M. Caro, A. Tamm, A. A. Correa, and A. Caro, Role of electrons in collision cascades in solids. i. dissipative model, Phys. Rev. B 99, 174301 (2019).
  * Tamm  _et al._ [2019] A. Tamm, M. Caro, A. Caro, and A. A. Correa, Role of electrons in collision cascades in solids. ii. molecular dynamics, Phys. Rev. B 99, 174302 (2019).
  * Sand  _et al._ [2019] A. Sand, R. Ullah, and A. Correa, Heavy ion ranges from first-principles electron dynamics, npj Comput Mater 5, 43 (2019).
  * Lee  _et al._ [2020] C.-W. Lee, J. A. Stewart, R. Dingreville, S. M. Foiles, and A. Schleife, Multiscale simulations of electron and ion dynamics in self-irradiated silicon, Phys. Rev. B 102, 024107 (2020).
  * Jarrin  _et al._ [2021] T. Jarrin, J. Teunissen, F. Da Pieve, A. Hémeryck, and N. Richard, Integration of electronic effects into molecular dynamics simulations of collision cascades in silicon from first-principles calculations, Phys. Rev. B 104, 195203 (2021).
  * Caro  _et al._ [2015] A. Caro, A. A. Correa, A. Tamm, G. D. Samolyuk, and G. M. Stocks, Adequacy of damped dynamics to represent the electron-phonon interaction in solids, Phys. Rev. B 92, 144309 (2015).
  * Caro  _et al._ [2017] M. Caro, A. A. Correa, E. Artacho, and A. Caro, Stopping power beyond the adiabatic approximation, Scientific Reports 7, 2618 (2017).
  * Caro  _et al._ [2018] M. Caro, A. Tamm, A. Correa, and A. Caro, On the local density dependence of electronic stopping of ions in solids, Journal of Nuclear Materials 507, 258 (2018).
  * Correa [2018] A. A. Correa, Calculating electronic stopping power in materials from first principles, Computational Materials Science 150, 291 (2018).
  * Duffy and Rutherford [2006] D. M. Duffy and A. M. Rutherford, Including the effects of electronic stopping and electron–ion interactions in radiation damage simulations, Journal of Physics: Condensed Matter 19, 016207 (2006).
  * Raine  _et al._ [2016] M. Raine, A. Jay, N. Richard, V. Goiffon, S. Girard, M. Gaillardin, and P. Paillet, Simulation of single particle displacement damage in silicon –part i: Global approach and primary interaction simulation, IEEE Transactions on Nuclear Science 64, 133 (2016).
  * Béland  _et al._ [2011] L. K. Béland, P. Brommer, F. El-Mellouhi, J.-F. m. c. Joly, and N. Mousseau, Kinetic activation-relaxation technique, Phys. Rev. E 84, 046704 (2011).
  * Pruneda  _et al._ [2007] J. M. Pruneda, D. Sánchez-Portal, A. Arnau, J. I. Juaristi, and E. Artacho, Electronic stopping power in lif from first principles, Phys. Rev. Lett. 99, 235501 (2007).
  * Zeb  _et al._ [2012] M. A. Zeb, J. Kohanoff, D. Sánchez-Portal, A. Arnau, J. I. Juaristi, and E. Artacho, Electronic stopping power in gold: The role of $d$ electrons and the $\mathbf{H}/\mathrm{He}$ anomaly, Phys. Rev. Lett. 108, 225504 (2012).
  * Ullah  _et al._ [2018] R. Ullah, E. Artacho, and A. A. Correa, Core electrons in the electronic stopping of heavy ions, Phys. Rev. Lett. 121, 116401 (2018).
  * Schleife  _et al._ [2015] A. Schleife, Y. Kanai, and A. A. Correa, Accurate atomistic first-principles calculations of electronic stopping, PHYSICAL REVIEW B 91, 014306 (2015).
  * Gygi [2008] F. Gygi, Architecture of qbox: A scalable first-principles molecular dynamics code, IBM Journal of Research and Development 52, 137 (2008).
  * Andrade and Draeger [2012] X. Andrade and E. Draeger, Qb@ll, https://github.com/LLNL/qball (2012).
  * Perdew  _et al._ [1996] J. P. Perdew, M. Ernzerhof, and K. Burke, Rationale for mixing exact exchange with density functional approximations, The Journal of Chemical Physics 105, 9982 (1996).
  * Reeves  _et al._ [2016] K. G. Reeves, Y. Yao, and Y. Kanai, Electronic stopping power in liquid water for protons and $\alpha$ particles from first principles, Phys. Rev. B 94, 041108 (2016).
  * Hamann  _et al._ [1979] D. R. Hamann, M. Schlüter, and C. Chiang, Norm-conserving pseudopotentials, Phys. Rev. Lett. 43, 1494 (1979).
  * pse [2020] Pseudopotential libary, quantum-simulation.org (2020), accessed: June 2020.
  * Fuchs and Scheffler [1999] M. Fuchs and M. Scheffler, Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory, Computer Physics Communications 119, 67 (1999).
  * [64] See Supplemental Material at link for ESP evaluations of highly ionized As PKAs, information about how the average ESP is obtained and the close collision peak occurring during the off-channeling trajectory, the force field potential visualisation, atomic density visualisations and how the number of defects changes w.r.t. the defect definition and cut-off value.
  * Zhong  _et al._ [1998] Y. Zhong, K. Nordlund, M. Ghaly, and R. S. Averback, Defect production in tungsten: A comparison between field-ion microscopy and molecular-dynamics simulations, Phys. Rev. B 58, 2361 (1998).
  * Rutherford and Duffy [2007] A. M. Rutherford and D. M. Duffy, The effect of electron–ion interactions on radiation damage simulations, Journal of Physics: Condensed Matter 19, 496201 (2007).
  * J. Lindhard and Thomsen [1963] M. S. J. Lindhard, V. Nielsen and P. V. Thomsen, Integral equations governing radiation effects, Mat. Fys. Medd. Dan. Vid. Selsk. 33 (1963).
  * Lohmann  _et al._ [2020] S. Lohmann, R. Holeňák, and D. Primetzhofer, Trajectory-dependent electronic excitations by light and heavy ions around and below the bohr velocity, Phys. Rev. A 102, 062803 (2020).
  * Jarrin  _et al._ [2020] T. Jarrin, A. Jay, A. Hémeryck, and N. Richard, Parametric study of the two-temperature model for molecular dynamics simulations of collisions cascades in si and ge, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 485, 1 (2020).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117, 1 (1995).
  * Albe  _et al._ [2002] K. Albe, K. Nordlund, J. Nord, and A. Kuronen, Modeling of compound semiconductors: Analytical bond-order potential for ga, as, and gaas, Phys. Rev. B 66, 035205 (2002).
  * Ziegler  _et al._ [1985] J. Ziegler, J. Biersack, and U. Littmark, The stopping and ranges of ions in solids, vol. 1, The Stopping and Ranges of Ions in Matter (1985).
  * Shinoda  _et al._ [2004] W. Shinoda, M. Shiga, and M. Mikami, Rapid estimation of elastic constants by molecular dynamics simulation under constant stress, Phys. Rev. B 69, 134103 (2004).
  * Blondel  _et al._ [2008] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, Fast unfolding of communities in large networks, Journal of Statistical Mechanics: Theory and Experiment 2008, P10008 (2008).
  * Koval  _et al._ [2020] N. E. Koval, F. Da Pieve, and E. Artacho, Ab initio electronic stopping power for protons in ga0.5in0.5p/gaas/ge triple-junction solar cells for space applications, R. Soc. Open Sci. 7, 200925 (2020).
  * Ullah  _et al._ [2015] R. Ullah, F. Corsetti, D. Sánchez-Portal, and E. Artacho, Electronic stopping power in a narrow band gap semiconductor from first principles, Phys. Rev. B 91, 125203 (2015).
  * Grinberg  _et al._ [2000] I. Grinberg, N. J. Ramer, and A. M. Rappe, Transferable relativistic dirac-slater pseudopotentials, Phys. Rev. B 62, 2311 (2000).
  * Jones  _et al._ [1993] D. R. Jones, C. D. Perttunen, and B. E. Stuckman, Lipschitzian optimization without the lipschitz constant, Journal of Optimization Theory and Applications 79, 157 (1993).
  * [79] S. Johnson, The nlopt nonlinear-optimization package, github.com/stevengj/nlopt.
  * Nordlund  _et al._ [2001] K. Nordlund, J. Peltola, J. Nord, J. Keinonen, and R. S. Averback, Defect clustering during ion irradiation of gaas: Insight from molecular dynamics simulations, Journal of Applied Physics 90, 1710 (2001).
  * Wigner and Seitz [1933] E. Wigner and F. Seitz, On the constitution of metallic sodium, Phys. Rev. 43, 804 (1933).
  * Hensel and Urbassek [1998] H. Hensel and H. M. Urbassek, Implantation and damage under low-energy si self-bombardment, Phys. Rev. B 57, 4756 (1998).
  * Messenger  _et al._ [2001b] S. Messenger, R. Walters, E. Burke, G. Summers, and M. Xapsos, Niel and damage correlations for high-energy protons in gallium arsenide devices, IEEE Transactions on Nuclear Science 48, 2121 (2001b).
  * Chen  _et al._ [2017] N. Chen, S. Gray, E. Hernandez-Rivera, D. Huang, P. D. LeVan, and F. Gao, Computational simulation of threshold displacement energies of gaas, Journal of Materials Research 32, 1555–1562 (2017).