# Examining Different Regimes of Ionization-Induced Damage in GaN Through Atomistic Simulations

Miguel Sequeira, Flyura Djurabekova, Kai Nordlund, Jean-Gabriel Mattei, Isabelle Monnet, Clara Grygiel, Eduardo Alves, Katharina Lorenz

► To cite this version:

Miguel Sequeira, Flyura Djurabekova, Kai Nordlund, Jean-Gabriel Mattei, Isabelle Monnet, et al. Examining Different Regimes of Ionization-Induced Damage in GaN Through Atomistic Simulations. Small, 2022, 18 (49), pp.2102235. 10.1002/smll.202102235. hal-03922054

HAL Id: hal-03922054

https://hal.science/hal-03922054v1

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Examining Different Regimes of Ionization-Induced Damage in GaN Through Atomistic Simulations

Miguel C. Sequeira,* Flyura Djurabekova, Kai Nordlund, Jean-Gabriel Mattei, Isabelle Monnet, Clara Grygiel, Eduardo Alves, and Katharina Lorenz

The widespread adoption of gGaN in radiation-hard semiconductor devices relies on a comprehensive understanding of its response to strongly ionizing radiation. Despite being widely acclaimed for its high radiation resistance, the exact effects induced by ionization are still hard to predict due to the complex phase-transition diagrams and defect creation-annihilation dynamics associated with group-III nitrides. Here, the Two-Temperature Model, Molecular Dynamics simulations and Transmission Electron Microscopy, are employed to study the interaction of Swift Heavy Ions with GaN at the atomic level. The simulations reveal a high propensity of GaN to recrystallize the region melted by the impinging ion leading to high thresholds for permanent track formation. Although the effect exists in all studied electronic energy loss regimes, its efficiency is reduced with increasing electronic energy loss, in particular when there is dissociation of the material and subsequent formation of N₂ bubbles. The recrystallization is also hampered near the surface where voids and pits are prominent. The exceptional agreement between the simulated and experimental results establishes the applicability of the model to examine the entire electronic energy loss spectrum. Furthermore, the model supports an empirical relation between the interaction cross sections (namely for melting and amorphization) and the electronic energy loss.

## 1. Introduction

Group-III nitride semiconductors have contributed to some of the most significant technological achievements of the last decades, setting new standards and often changing the landscape of entire industries. Light Emitting Diodes and Laser Diodes are perhaps the two most accomplished examples.[1] More recently, they have seen another surge of interest due to their use in High Electron Mobility Transistors (HEMT). When

comparing GaN HEMT to Si transistors, the former not only offer higher efficiency, work at higher voltages, and have faster-switching transitions but also enjoy better thermal stability and are more resistant to radiation.[2,3] When compared to SiC, GaN-based HEMT offer much higher electronic mobility.[4,5] Along with their promising performance in flexible electronics,[6] these characteristics make them a prime candidate for high-frequency applications such as 5G communication systems.[7]

Another example of significant relevance is the use of GaN-based devices in extreme-radiation environments such as space or high energy experimental facilities. For instance, the European Space Agency's new generation of experimental satellites, such as the Proba-V[8] and the Biomass (2022),[9] are using GaN-based HEMT to communicate with Earth. Another example is the possibility of using GaN-based detectors on CERN's Large Hadron Collider, a change that would significantly reduce maintenance and

downtime.[3,10] The ionic part of the radiation spectrum present in these environments comprises two energy regimes: low to medium, where the ions slow down mainly by colliding with the atoms of the material, and high, where the ions lose their energy mostly by ionization processes. While the former regime is relatively well understood (e.g., Refs. [11–13]), the latter is still under debate.

In space, the strongly ionizing radiation is of particular interest since the Galactic Cosmic Radiation (GCR) and the

M. C. Sequeira
IBC
Helmholtz-Zentrum Dresden-Rossendorf
Dresden, Germany
E-mail: miguel.sequeira@tecnico.ulisboa.pt

The ORCID identification number(s) for the author(s) of this article can be found under https://doi.org/10.1002/smll.202102235.

© 2022 The Authors. Small published by Wiley-VCH GmbH. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

DOI: 10.1002/smll.202102235

M. C. Sequeira, E. Alves, K. Lorenz
IPFN
Instituto Superior Técnico
Campus Tecnológico e Nuclear, Lisbon, Portugal
F. Djurabekova, K. Nordlund
Department of Physics
University of Helsinki
Helsinki, Finland

J.-G. Mattei, I. Monnet, C. Grygiel
CIMAP

Normandie University
CEA, CNRS, UNICAEN, ENSICAEN-BP5133, Caen Cedex 5, France

K. Lorenz
INESC MN
Rua Alves Redol 9, Lisbon 1000-029, Portugal
solar radiation spectrum contain a substantial amount of high energy heavy ions (HZE—High atomic number Z and high Energy). Although with lower fluxes than lighter particles such as protons and α-particles, HZE ions have energies in the order of 10 MeV amu⁻¹, [14-16] which induce extremely high ionization spikes in most materials (e.g., [17-20]). These can cause malfunction of the electronics as well as premature degradation. [21,22] In space, as in any extreme-radiation environment, unreliability and obsolescence should be avoided.

Therefore, there is some urgency in understanding the interaction of strongly ionizing radiation with GaN and establishing a fast and reliable method to predict the radiation effects in related devices. On Earth, HZE radiation effects can be studied using Swift Heavy Ions (SHI) with energies per nucleon in the MeV range and above. In semiconductors, the ionization spikes produced by the SHI often lead to the formation of permanent tracks [17,23,24] and surface modification. [17,25,26] Similar effects have been observed in dielectrics. [27-29] Moreover, it has been reported that SHI radiation can lead to the recovery of pre-damaged material, [17,19,23,30-33] including the recrystallization of amorphous materials. Concerning GaN, its response to SHI is still unsettled. Sall et al. [34], using Transmission Electron Microscopy (TEM), found that SHI can amorphize GaN films. In contrast, Kucheyev et al. [35] used Rutherford Backscattering Spectrometry/Channeling (RBS/C) to reveal that even under high irradiation fluences (above $1 \times 10^{12}$ cm⁻²), GaN films, although damaged, still maintained a crystalline structure. Moreover, the published electronic energy loss thresholds above which permanently damaged tracks begin to form range from just 8 keV nm⁻¹ near the surface [36] to 15-22 keV nm⁻¹ in the bulk. [34,36] Thus, fundamental characteristics of the interaction, such as the type and concentration of defects produced by a given SHI energy, and its dynamics remain questionable.

Recently, Sequeira et al. proposed a Two Temperature Model—Molecular Dynamics (TTM-MD) simulation scheme to describe the interaction of SHI with GaN. [17] Within this model, the ionization spike is depicted as a hot electron gas that quickly relaxes its energy into the target lattice via electron-phonon interactions, originating a sudden increase of the lattice temperature. The TTM is used to study the electronic relaxation (i.e., the first stage of the process) and to obtain the temperature profile in the lattice. [37] Afterward, the corresponding energy profile is deposited in a MD cell that is relaxed until room temperature. Besides successfully describing the effects of SHI in GaN, variations of TTM-MD have been applied to study SHI interactions with other materials, namely with semiconductors, insulators and metals, (e.g., [19,20,32]) with different degrees of agreement with experiments.

However, so far, in GaN, the model has only been tested for SHI inducing mostly congruent melting. [17] Therefore, the applicability of the TTM-MD model to describe the effects of SHI within the entire spectrum of energy loss remains to be confirmed. GaN, alike other group-III nitrides, has a complex phase-transition diagram. [38] In particular, congruent melting (a reversible transition) and molecular dissociation (an irreversible transition) can both occur at high temperatures depending on the pressure in the crystal. [39] SHI with different energies and ionization powers induce distinct temperature-pressure states within their tracks, and both phase-transitions may occur either individually or concurrently.

Here, we aim at understanding the full extent of the effects of SHI irradiation in GaN by combining TTM-MD simulations and Transmission Electron Microscopy (TEM). GaN-on-Al₂O₃ thin films are irradiated with SHI whose electronic energy losses, $ϵ_e$, depict most of the spectrum of experimentally available values (see Table 1; Note S1, Supporting Information). This allows for inspecting the behavior of the material in multi-ion radiation environments (e.g., space) as well as establishing the universality of the TTM-MD model. The same model parameters obtained in Ref. [17], including the temperature-dependent electronic heat capacity and conductivity as well as the electron-phonon coupling, are used. See Experimental Section and Note S2 (Supporting Information) for further details on the simulation scheme and track radii determination. We simulate the interaction at different depths to address the impact of the surface in the defect dynamics. The obtained results reveal an impressive agreement between simulation and experiment. They allow for an understanding of the dynamic processes that lead to the observed defective structures such as amorphous regions, point and extended defects, voids and N₂ bubbles. Furthermore, the simulations clarify the ongoing debate on the energy loss threshold for track formation, including the significant disparities in the values obtained at different depth regions. Based on the present results, we propose thresholds for the formation of discontinuous as well as continuous amorphous tracks, and discuss the existence of a third one above which N₂ bubbles are formed. Despite the complex defect dynamics, we found an empirical relation between $ϵ_e$ (readily obtained using codes such as SRIM, Stopping and Range of Ions in Matter. [40] or CasP, Convolution approximation for swift Particle [41]) and the interaction cross sections for melting and amorphization, which can simplify the planning of future experimental studies.

## 2. Results

### 2.1. Xe Irradiations

According to the TTM-MD simulations, the three Xe irradiations (with $ϵ_e = 14.9 - 19.9$ keV nm⁻¹, Table 1) cause thermal spikes that increase the temperature of the lattice to values above 9000 K. Such temperatures are clearly above the solid-liquid phase transition temperature, and thus a cylindrical molten track is formed along the ion path. However, the phase transition is not immediate since the atomic motion occurs in

<table><caption>**Table 1.** Details of the SHI radiations studied here. The electronic energy losses $ϵ_e$ are obtained using the SRIM code. [40] The samples were irradiated with a fluence of $1 \times 10^{11}$ cm⁻².</caption><thead><tr><th>Ion</th><th>Specific Energy [MeV amu<sup>-1</sup>]</th><th>Energy Loss [keV nm<sup>-1</sup>]</th></tr></thead><tbody><tr><td rowspan="3"><sup>129</sup>Xe<sup>23+</sup></td><td>0.35</td><td>14.9</td></tr><tr><td>0.43</td><td>17.2</td></tr><tr><td>0.54</td><td>19.9</td></tr><tr><td rowspan="2"><sup>208</sup>Pb<sup>28+</sup></td><td>0.6</td><td>28.3</td></tr><tr><td>5.8</td><td>46.2</td></tr><tr><td><sup>238</sup>U<sup>55+</sup></td><td>3.8</td><td>55.2</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Temperature evolution in the cell after being impacted by a 0.54 MeV amu⁻¹ Xe ion. The color code illustrates the temperature distribution. After 0.1 ps following the impact, an area with a radius of 5.5 nm has a temperature above 2500 K, the melting temperature at atmospheric pressure. However, the structural changes in the crystal are still minimal, and the solid–liquid phase transition is just beginning to develop. During the first 6 ps, the melted region increases in a radial direction until it reaches its maximum radius of r_t = 2.8 nm. This area is smaller than the area initially above the melting temperature due to the pressure induced by the thermal spike. When the system returns to room temperature, ≈155 ps after impact, only a small region, magnified in the inset, remains defective, with a radius of r_c = 0.7 nm.</figcaption></figure>

a considerably longer time frame than the fast electronic relaxa-
tion. In **Figure 1**, we show the simulated cross-sectional images
of the pivotal moments occurring after a 0.54 MeV amu⁻¹ Xe
ion impact, the more energetic of the three Xe radiations.
At 0.1 ps after ion impact, only minimal structural changes
occurred, even though the lattice temperature is well above
2500 K (the melting temperature of GaN at atmospheric pres-
sure[³⁸]) in a cylinder with a radius of 5.5 nm. The molten region
increases radially during the first 6 ps of the simulation when
its radius reaches its maximum value, r_t, of 2.8 ±~0.2 nm. Due
to the high pressure within the track, the molten region is sig-
nificantly smaller than the region above melting temperature at
atmospheric pressure. Furthermore, the pressure state ensures
a^ congruent melting, which permits most of the molten mate-
rial to recover its crystalline structure during the subsequent
cooling process.

After 6 ps, recrystallization starts at the edge of the molten track and moves inwards. When room-temperature is reached, a substantial part of the molten region is almost completely recovered. Inspecting the final state of the tracks in the deep regions of the GaN layers, shown in the cross-sectional images of Figure 2a–c (bottom row), reveals that all three Xe irradiations produce discontinuous tracks consisting essentially of small clusters of point defects. The damaged core of the track formed by the 0.54 MeV amu⁻¹ Xe ion (19.9 keV nm⁻¹) has a maximum radius of only r_c = 0.7 ± 0.2 nm. The damaged regions are even smaller in the irradiations with lower energies, being almost inexistent in the case of the 0.35 MeV amu⁻¹ Xe ion (15 keV nm⁻¹) as evident in Figure 2a bottom (see also the plan-view images shown in Figures S2, S4, S6, S8, S10, and S12 of Supporting Information Note S2).

These results are in excellent agreement with TEM results
obtained by Sall et al.[³⁴] Using cross-sectional images of GaN
samples irradiated with 0.45 MeV amu⁻¹ U ions, the authors
observed that at a depth of 2.5 µm, the ion tracks disappeared,
suggesting that the electronic energy loss threshold for track
formation, ε_c^th, is 15 keV nm⁻¹. This value corresponds to
the simulation shown in Figure 2a, where indeed almost no
damage is observed after irradiation. Moreover, in the same
work, it is revealed by plan-view TEM images that the 0.51 MeV

amu⁻¹ Pb ion (e_e = 24 keV/nm) produces discontinuous tracks with a radius of around 0.8 nm, close to the 0.7 nm predicted for the 0.54 MeV amu⁻¹ Xe ion (19.9 keV nm⁻¹). In another study, Karlušić et al., supported by RBS/C measurements of samples irradiated with 0.7 MeV amu⁻¹ Xe ions, proposed that SHI with e_e <22 keV nm⁻¹ do not form tracks.[³⁶] This higher estimation of ε_e^th is likely due to the SHI in this range of e_e creating damaged structures with concentrations too low to be undoubtedly detected in an RBS/C spectrum. Nonetheless, RBS/C measurements of the samples irradiated with Xe ions reported in Ref. [17] reveal that the less energetic 0.35 MeV amu⁻¹ Xe ion (15 keV nm⁻¹) produces no measurable damage while the 0.54 MeV amu⁻¹ Xe ion (19.9 keV nm⁻¹) induces a small increase of the backscattering yield.

In the shallow regions, the surface influences the relaxa-
tion dynamics significantly since the lattice can expand along
the direction normal to the surface. This additional degree
of freedom lowers the pressure near the surface and contrib-
utes to the formation of a pressure gradient along the track,
which ultimately originates the sputtering of Ga and N atoms.
Consequently, voids appear within the track near the surface,
and a nano hill is formed, as shown in Figure 2a–c top. Similar
structures were also observed by Karlušić et al.[³⁶] in Atomic
Force Microscopy (AFM) images of tracks formed by 0.7 MeV
amu⁻¹ Xe ions (22.8 keV nm⁻¹) impinging GaN at a grazing
angle. The 0.54 MeV amu⁻¹ Xe ions create voids until a depth
of 15 nm while in the case of the 0.35 MeV amu⁻¹ Xe ions this
value decreases to 6 nm.

High-Resolution (HR) and Fresnel TEM images of a GaN
sample irradiated with 0.54 MeV amu⁻¹ Xe ions, shown in
Figure 2d,e, supports these results. The HR-TEM image
(Figure 2d) reveals that even after a fluence of 2 × 10¹² cm⁻²,
almost no damage is formed in the film. Only a very obscure
track near the surface could be found (along the arrow) without
any evidence of containing low density structures. However,
inspecting the Fresnel contrast of the same region, shown in
Figure 2e, the presence of low-density structures near the sur-
face of the irradiated sample becomes obvious. Note that the
dimension of structures observed in Fresnel images is magni-
fied by the out-of-focus condition inherent to the technique.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** a–c) Cross-sectional images of the final state of a simulated track morphology in the shallow and deep regions (top and bottom row, respectively) created by a) 0.35 MeV amu⁻¹ ($\epsilon_e$ = 15 keV nm⁻¹), b) 0.43 MeV amu⁻¹ ($\epsilon_e$ = 17 keV nm⁻¹), and c) 0.54 MeV amu⁻¹ Xe ($\epsilon_e$ = 19.9 keV nm⁻¹) ions. d) HRTEM image of a GaN film irradiated with 0.54 MeV amu⁻¹ Xe ions near the surface. The arrow indicates a possible SHI track. The thick dashed line marks the border between surface and bulk regions, as seen by the simulations. e) Under-focused Fresnel image from the sample shown in (d). The sample was tilted by 20° before measurement. Over-focused and focused images are shown in Figure S21 in Supporting Information Note S5. Several low-density structures are visible near the surface. Their size and density decrease with depth, as predicted by the simulations. Note that dimensions in Fresnel analysis are enlarged due to its out-of-focus condition.</figcaption></figure>

Therefore, TTM-MD can describe the damage formed by radiation with $\epsilon_e$ neighboring $\epsilon_e^{th}$, both in the bulk and near the surface. The model explains the morphology of the tracks observed experimentally, in particular, their discontinuity in the deep regions and the formation of voids in the shallow region. I also clarifies the observed energy loss thresholds. A technique with low spatial resolution, such as RBS/C, can hardly detect the low concentration of defects created by SHI within this range of $\epsilon_e$. In contrast, the same structures are detectable by TEM, a technique characterized by its very high spatial resolution. Based on the TTM-MD simulations, we recommend assuming $\epsilon_e^{th}$ = 15 keV nm⁻¹ as the threshold value for permanent damage formation in regions below 25 nm. In the shallow regions, the simulations predict a significant amount of sputtering along with a decline in recrystallization efficiency. Both effects contribute to reducing the value of $\epsilon_e^{th}$, as is particularly evident in the surface simulations of the 0.35 MeV amu⁻¹ Xe ion shown in Figure 2a. Despite barely creating any damage in the deeper regions, the surface is clearly affected by the SHI impact, indicating that $\epsilon_e^{th}$ is lower than 15 keV nm⁻¹. This limit is in good agreement with previous AFM observations of samples irradiated at a grazing incidence, which showed that SHI with $\epsilon_e$ = 8 keV nm⁻¹ would already modify the surface.[36] See Note S2 (Supporting Information) for further details on the recrystallization efficiency in bulk and surface.

## 2.2. 0.6 MeV amu⁻¹ Pb Irradiation

To understand the transition between the discontinuous tracks produced by ions with $\epsilon_e \le 20$ keV nm⁻¹ (Xe irradiations) and the continuous tracks formed by ions with $\epsilon_e \ge 33$ keV nm⁻¹ (0.9 MeV amu⁻¹ Au irradiation[17]), we simulate the 0.6 MeV amu⁻¹ Pb ion ($\epsilon_e$ = 28.3 keV nm⁻¹). The results presented in Figure 3 show the formation of a near-continuous track throughout the entire cell. Alike the Xe and 0.9 MeV amu⁻¹ Au irradiations, a strong recrystallization effect is observed in the molten tracks produced by the 0.6 MeV amu⁻¹ Pb ion, especially in the deeper regions of the cell, as shown in Figure 3a. In this region, the molten track reaches a maximum radius of $r_t$ = 4.1 ± 0.2 nm, 13 ps after ion impact (see Figure S8 of Supporting Information Note S2). After the cooling of the crystal, most of the initial damage is recovered with only a small amorphous core with a radius of $r_c$ = 1.3 ± 0.2 nm remaining, as evidenced in Figure 3b,c. Mansuri et al. used plan-view TEM images of GaN films irradiated with the same Pb ion to obtain the value of $r_c^{exp}$ = 1.5 nm,[42] in excellent agreement with our results.

Although the amorphous core is nearly continuous, the simulations show that there can be some small discontinuities. An example of such discontinuity, with a length of 1.5 nm, can be seen in the bottom part of the magnified image shown in Figure 3c. Their presence suggests that the energy deposited by the 0.6 MeV amu⁻¹ Pb ion, $\epsilon_e$ = 28.3 keV nm⁻¹, should be
0.6 MeV/amu Pb

<figure><img src="image_4.png" /><figcaption>**Figure 3.** a) Cross-section and b) plan-view images of a track produced by a 0.6 MeV amu⁻¹ Pb ion ($e_e$ = 28.3 keV nm⁻¹). The color scheme represents the magnitude of the displacement of the atoms from their original positions. The plan-view image b) and the magnification shown in c) are taken from the region between the dashed lines in (a). A recrystallized shell, with $r_t$ = 4.1 ± 0.2 nm, surrounds a small but almost continuous amorphous core with $r_c$ = 1.3 ± 0.2 nm.</figcaption></figure>

close to a threshold value for continuous track formation. Thus, we propose that SHI with $\epsilon_e$ between 15 and ≈28.3 keV nm⁻¹ induce discontinuous tracks. Below 15 keV nm⁻¹, the crystal remains entirely crystalline with only a few point defects and small point defect clusters, while above 28.3 keV nm⁻¹, the ion induces tracks with continuous amorphous cores. The threshold for continuous track formation is also supported by experimental results. For instance, the 0.45 MeV amu⁻¹ U ion, with $\epsilon_e$ = 24 keV nm⁻¹, produces discontinuous tracks, even in the regions near the surface.[34] In contrast, the 0.9 MeV amu⁻¹ Au ion, depositing energy at a rate of $\epsilon_e$ = 33 keV nm⁻¹, induces continuous and well-defined tracks.[17]

For the same reasons discussed for the case of the Xe irra-
diations, we found voids near the surface, with the deepest
one located at approximately 25 nm, as shown in Figure 3a.
Moreover, the amorphous core has a slightly higher radius
in this region than in the bulk. We obtain $r_c$ = 2 ± 0.3 nm in
the layer extending the first 10 nm from the surface instead of
the 1.3 ± 0.2 nm obtained in the deeper regions. This enlarge-
ment strengthens the suggestion that the surface enhances the

formation of tracks, as discussed in the context of the Xe irra-
diations and in Note S2 (Supporting Information).

2.3. 5.8 MeV amu⁻¹ Pb Irradiation

The 5.8 MeV amu⁻¹ Pb ion, with $\epsilon_e$ = 46.23 keV nm⁻¹, is in the upper limit of experimentally accessible energy losses in GaN (see Figure S1 of Supporting Information Note S1). In Figure 4a,b, we show the TTM-MD outcome along with the TEM images in Figure 4c–f concerning this irradiation. The sputtering and concomitant formation of voids and nano hills is found in the simulations, similar to the results shown for lower energy losses. However, the density of voids, as well as their dimensions, increases significantly and reaches deeper layers. Remarkably, this SHI still induces a substantial amount of recrystallization, as shown in the plan-view image in Figure 4b. Comparing the size of the amorphous cores with those found in the experimental TEM images of samples irradiated with 4.5 MeV amu⁻¹ Pb ions ($\epsilon_e$ = 46 keV nm⁻¹), shown in Figure 4c (adapted from Ref. [34]), reveals a striking agreement. The TTM-MD predicts $r_c$ = 2.2 ± 0.2 nm, while the experimental reported value is, on average, $r_c^{exp}$ = 2.4 nm. Regarding $r_t$, we obtain 4.7 ± 0.1 nm (see Figure S10 of Supporting Information Note S2 and Movie S1, Supporting Information, to visualize the entire relaxation process). Although $r_t$ is not obtainable from the TEM image presented in Figure 4c, one can identify a loss of contrast in a region with a comparable radius surrounding the amorphous cores.

The main features of the tracks produced by the 5.8 MeV amu⁻¹ Pb ions near the surface follow the ones already described above for less energetic SHI. As shown in Figure 4a, TTM-MD predicts the formation of a track with a conical shape near the surface (see Figure S11 of Supporting Information Note S2) with a radius of 2.2 ± 0.1 nm (see Experimental Section for a discussion on the track radii determination). A similar conical structure is identifiable in the near-surface cross-sectional HRTEM image shown in Figure 4d. At a depth of 5 nm, this structure has a diameter of ≈4.1 ± 0.6 nm, matching the predicted value of 4.4 ± 0.2 nm shown in Figure 4a. Such surface structure is also observed in the under-focused Fresnel image, Figure 4e (the over-focused image is shown in Figure S22 in Supporting Information Note S5), which indicates a sharp difference in the density between the structure and the surrounding crystal. We recall that the dimensions and morphology of the defects observed under the Fresnel condition are not directly comparable with the simulations nor TEM images.

However, the most notable difference between the inter-
action of the 5.8 MeV amu⁻¹ Pb ion and SHI with lower ener-
gies is the appearance of N₂ bubbles in the tracks created by the
former ion (see Figure 4a and Figure S17 of Supporting Infor-
mation Note S2). These N₂ bubbles, found in the deep regions,
are a direct product of the dissociation of GaN that occurs
under certain temperature–pressure states.[39] For SHI with low
to medium energies, despite the extremely high temperatures,
the coexistent pressure within the track keeps GaN mostly in
a congruently melted state, apart from the formation of a tiny
amount of N₂ molecules that eventua form split-interstitials.
In contrast, the temperature-pressure induced by high energy
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Simulated a) cross-sectional and b) plan-view images of the track induced by a 5.8 MeV amu⁻¹ Pb (46.23 keV nm⁻¹) ion. The dashed lines in a) indicate the region from which the image in (b) was taken. c) Detail of the N₂ bubbles within the track. Experimental plan-view TEM image of a GaN film irradiated with 4.5 MeV amu⁻¹ Pb ion (46 keV nm⁻¹) adapted from Ref. [34]. Experimental cross-sectional d) HRTEM, e,f) Fresnel analysis of a GaN film irradiated with 5.8 MeV amu⁻¹ Pb ions and a fluence of 1 × 10¹¹ cm⁻². The tracks observed experimentally correspond with great detail to those simulated, including surface pits with different shapes and the presence of N₂ bubbles uniformly distributed with depth. See Movie S1 (Supporting Information) to visualize the entire relaxation process.</figcaption></figure>

irradiations is such that a part of the heated material dissociates into Ga and N₂ molecules. During the cooling of the crystal, these molecules aggregate in the form of bubbles. Due to the strong binding of the N₂ molecule, this process is irreversible, and the bubbles remain stable when room-temperature is reached.

Similar low-density pockets forming tracks along the direction of the ion path are observable experimentally in the cross-sectional Fresnel images shown in Figure 4e,f. Excepting the surface pit mentioned above, the distribution of these pockets is relatively uniform with depth. This uniformity strongly suggests the formation of N₂ bubbles in the deeper regions, as predicted by TTM-MD, instead of voids since the latter would require a substantial upward flow and consequent sputtering to create such a distribution.

Such low-density structures are also observed in TEM images of tracks induced by 4.5 MeV amu⁻¹ Pb ions ($\epsilon_e$ = 46 keV nm⁻¹)[34] and by 7.4 MeV amu⁻¹ Bi ions ($\epsilon_e$ = 45.9 keV nm⁻¹).[43] In the latter work, it is estimated that the bubble-like structures have a radius of 2.9 and 2.2 nm in the shallow and deep regions, respectively. Both values are close to those taken from Figure 4a, 2.2 and 1.8 nm, respectively. Furthermore, a clear sign of the presence of N₂ within the tracks was observed in GaN samples irradiated with fullerene projectiles using Electron Energy Loss Spectroscopy.[44] Despite being

formed by a different mechanism, N₂ bubbles have also been reported previously in ion-implantation studies (i.e., ions with much lower energies that affect the crystal predominantly via nuclear interactions).[45,46] Thus, it is expectable that impacts of SHI with $\epsilon_e$ above a certain threshold produce N₂ bubbles.

### 2.4. 3.8 MeV amu⁻¹ U Irradiation

Focusing on the SHI with the highest energy loss studied here, we analyze the 3.8 MeV amu⁻¹ U irradiation ($\epsilon_e$ = 55.15 keV nm⁻¹). The simulated track morphology, shown in Figure 5a-c, is similar to that produced by the 5.8 MeV amu⁻¹ Pb ion (Figure 4). However, the surface voids are significantly larger in the case of the 3.8 MeV amu⁻¹ U ion, creating a pit with a depth of ≈22 nm and a radius of 2.9 nm, as evidenced in Figure 5a,c. The cross-sectional TEM image of a track induced by the same SHI, presented in Figure 5d, shows a remarkable resemblance with the predicted track. The experimentally observed structure near the surface has a shape and dimension that, when compared with the simulated image in Figure 5a, strongly suggests the existence of a surface pit. This pit has an apparent depth of 25 nm and $r_{exp} = 2.9 \pm 0.4$ nm, both values closely neighboring the predicted ones of 22 and 2.9 nm, respectively. Pits with similar depths are also seen in several tracks reveal the
<figure><img src="image_6.png" /></figure>

**Figure 5.** Simulated a) cross-sectional, b) plan-view images and c) 3D of a track produced by a 3.8 MeV amu⁻¹ U ion ($e_e = 55.15$ keV nm⁻¹). The orange balls in (c) represent N₂ molecules. d) Experimental TEM cross-sectional image of a track produced by the same U ion. The dashed lines in (a) indicate the region from which the image in (b) was taken. The simulated and experimental track morphologies are in excellent agreement, in particular near the surface. The simulated and experimental $r_c$ have an identical value of $2.9 \pm 0.2$ nm. Regarding $r_t$, the predicted value is $5.5 \pm 0.3$ nm, while the experimental one is unknown from the present experimental data. The parallelism between the structure found near the surface in (a), (c), and (d) strongly suggests the existence of a surface pit. Not only its shape is comparable, but also its dimensions are almost identical. The structure seen experimentally reaches a depth of ≈25 nm, very close to the predicted value of 22 nm.

TEM images obtained with less magnification presented in **Figure 6a,a'**.

The recrystallization effect continues to play an essential role in the track morphology. From the simulations, we calculate $r_t = 5.5 \pm 0.3$ nm and $r_c = 2.9 \pm 0.2$ nm in the deeper regions of the track (see Figure S12 of Supporting Information Note S2). That this effect is found in all SHI energy regimes demonstrates the unusual propensity for GaN to recover from irradiation damage. However, we also found a substantial amount of N₂ bubbles within the tracks, as shown in Figure 5a,c (the N₂ molecules are represented by orange balls in Figure 5c). In all experimental images, the track morphology below the pit becomes similar to that produced by the 5.8 MeV amu⁻¹ Pb ions. In particular, the cross-sectional TEM images presented in Figure 6b,b' reveal tracks containing several dark structures that mirror the predicted N₂ bubbles. As expected from the simulations, the supposed bubbles appear to be larger and more pronounced in the tracks induced by the 3.8 MeV amu⁻¹ U ion than by the 5.8 MeV amu⁻¹ Pb ion. Notably evident in Figure 6b', these structures appear elongated along the track axis, similar to the shape of the bubbles in the simulated cross-sectional image presented in Figure 5a,c (see also the magnified images in Figure S17 of Supporting Information Note S2).

<figure><img src="image_7.png" /></figure>

**Figure 6.** Cross-sectional TEM images of tracks induced by 3.8 MeV amu⁻¹ U ions in GaN films obtained at different magnifications. a') and b') show magnifications of the areas delimited by the dashed lines in a) and b), respectively. The solid line in (a) marks the 25 nm depth. The presence of structures resembling the simulated surface pits is evident in all the tracks shown in (a'). Below the pits, the morphology of the track becomes similar to that observed in the tracks shown in **Figure 5a**. The dark structures in (b'), as those within the circles, have shapes and dimensions analogous to the predicted N₂ bubbles. Furthermore, their distribution along the track follows the distribution of the bubbles observed in the simulations.

### 3. Discussion

We summarize the track dimensions and morphologies obtained from the simulations and experiment for the various SHI in **Figure 7**. The TTM-MD results are represented by arrows with their starts and heads indicating $r_t$ and $r_c$, respectively. Besides the SHI discussed here, we also include the dimensions of the tracks formed by the 0.9 MeV amu⁻¹ U ion reported in Ref. [17]. All of the $r_c$ predictions are remarkably close to their corresponding experimental $r_c^{exp}$ values taken from TEM images and also included in the figure. Even for SHI in the upper energy loss range, despite the more complex defect dynamics including, for instance, a significant amount of N₂ nucleation, the final track radii are exceptionally close to the experimental ones.

To emphasize the superior prediction power of TTM-MD with respect to pure TTM simulations, we include the melting radii obtained from the latter model in Figure 7. The TTM parameters are the same in both simulations. In pure TTM studies, the melting radius of the tracks is typically calculated by comparing the lattice temperature profile with the melting temperature of the material at atmospheric pressure (e.g., see Ref. [34] and Note S3 in Supporting Information for further details). Furthermore, in many cases, the melting radius is assumed to correspond to the amorphous radius (i.e., $r_c$), an approach that ultimately excludes the recrystallization effect from the analysis. This can lead to a considerable overestimation of the predicted radius, as evidenced in Figure 7.

Therefore, only the melting radius should be comparable between both models. However, even there, the $r_t$ values
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Compilation of the track radii obtained with TTM-MD (arrows), TTM (filled symbols) and experimentally (hollow symbols). The experimental $r_c$ are obtained from the following sources: the 0.35 and 0.54 MeV amu⁻¹ Xe ion was confirmed using RBS/C in Ref. [17]; the 0.64 and 5.8 MeV amu⁻¹ Pb ions were measured by TEM in Refs. [42] and [34] respectively, (the latter reference used the 4.6 MeV amu⁻¹ Pb with $\epsilon_e = 46$ keV nm⁻¹); the 0.9 MeV amu⁻¹ Au ion was obtained using TEM in Ref. [17]; the 3.8 MeV amu⁻¹ U ion was measured from TEM image (Figure 5d). The arrows illustrate the recrystallization effect, with their starts and heads indicating $r_t$ and $r_c$, respectively. The error bar at the start indicates the uncertainty in $r_t$ and the length of the head equals the uncertainty in $r_c$. The radii of the cores predicted by TTM-MD, across all SHI energies, are exceptionally close to the experimental values.</figcaption></figure>

diverge considerably between the TTM and TTM-MD, as shown in Figure 7. The TTM overestimates $r_t$ by an average of 71%, reaching as high as 161% in the case of the 0.54 MeV amu⁻¹ Xe ion. This overestimation arises from neglecting the influence of the pressure in the melting temperature when computing the TTM $r_t$. As mentioned previously, the pressure increases the melting temperature significantly above 2500 K, resulting in smaller melted regions. Yet, since there is no trivial method to calculate the pressure within the track from the TTM output, the analysis often defaults to assuming a constant melting temperature taken at atmospheric pressure.

To better understand the efficiency of the recrystallization induced by each SHI, we plot in **Figure 8a** the cross-sectional area of the recrystallized volume (i.e., $\pi(r_t^2 - r_c^2)$) as a function of $\epsilon_e$. This area increases linearly with $\epsilon_e$ until 33 keV nm⁻¹. It remains relatively constant between this value and until $\epsilon_e = 46.2$ keV nm⁻¹, when it begins increasing again.

Interestingly, as shown in Figure 8b, the ratio between the recrystallized and melted volumes, i.e. $(r_t^2 - r_c^2)/r_t^2$, decreases linearly with $\epsilon_e$. This linearity is unexpected since SHI with different energies produce different phenomena, including two distinct phase transitions (congruent melting and dissociation) that originate multiple types of defects. Fitting the data we obtain,

$$ \frac{r_t^2 - r_c^2}{r_t^2} = 1.061 - 6.181 \times 10^{-3} \epsilon_e \quad (1) $$

Perhaps the most exciting application of this empirical relation between $r_t$, $r_c$ and $\epsilon_e$ is the estimation of $r_t$ from the

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Study of the recrystallization efficiency as a function of $\epsilon_e$. a) Shows the absolute cross-sectional area of the recrystallized volume $(\pi(r_t^2 - r_c^2))$ while b) shows the ratio between the recrystallized and melted volumes ($(r_t^2 - r_c^2)/r_t^2$). The recrystallized ratio has a linear dependence with $\epsilon_e$. Due to the defects formed, especially the N₂ bubbles, the ratio decreases with energy.</figcaption></figure>

experimentally obtained $r_c$. Knowing both radii allows for an improved understanding of the damage evolution with fluency since it has been shown that impact overlapping can recover damage formed by previous impacts.[17] Furthermore, we can use the relation to estimate $r_c$ from the $r_t$ obtained using TTM, hence introducing the recrystallization effect into the model analysis. However, note that Equation (1) was only validated here for energy losses on the low energy side of the Bragg peak and up to its maximum (see Figure S1 of Supporting Information Note S1).

The same linear relationship is observed in the ratios calculated near the surface (see Figure S16 of Supporting Information Note S2). In the shallow region, we obtain $(r_t^2 - r_c^2)/r_t^2 = 1.020 - 7.135 \times 10^{-3} \epsilon_e$. Comparing with the ratio in the bulk region, we can infer that the decrease in the interception parameter indicates a lower threshold for track formation, while the higher slope points to quicker degradation of the recrystallization efficiency.

The high accuracy of the TTM-MD (with a unique set of SHI-independent parameters) predictions for $r_c$, across a vast range of energy losses, supports the conclusion that $r_t$ is also exact. This means that the extent of the recrystallization observed in a single track is indeed occurring. The recrystallization effect is undoubtedly the reason behind the unusual resistance of GaN to strongly ionizing radiation. However, for high energy irradiations, the dissociation of GaN, and subsequent formation of N₂ bubbles, likely plays a significant role in the electrical properties of devices. Notwithstanding, the TTM-MD can quantify the concentration of the bubbles as well as estimate the formation of point and extended defects by applying the Dislocation Extraction Algorithm (DXA).[17,47] Such insight into the defect dynamics allows for an improved understanding of the response of real devices to strongly ionizing irradiation, such as the study on HEMT reported in Ref. [43]. Moreover,
TTM-MD can describe the effects of the fluence, either in single or multiple ion radiation environment, since it permits the simulation of impact overlapping.[17] We remark that such experiments are often complicated and expensive to perform, making TTM-MD a valuable tool to analyze, plan, and design future experiments and devices.

## 4. Experimental Section

**Materials.**: Commercial c-plane GaN thin films (Lumilog), 3 µm thick, grown by Metal Organic Chemical Vapor Deposition on Al₂O₃ substrates were used.

**Radiations.**:- The samples were irradiated at the Grand Accélérateur National d'Ions Lourds (GANIL). Details on irradiation parameters are summarized in Table 1. In the case of the Xe irradiations, the initial beam energy of 92 MeV was degraded to 70, 55.5, and 45 MeV Xe ($e_e = 19.9, 17.2, \text{and } 14.9 \text{ keV nm}^{-1}$, respectively) using Al foils of different thickness. This process also ensured that the ions had approximately acquired the equilibrium charge state before entering the GaN crystal. The irradiations were performed with a zero angle incidence with respect to the surface normal. Note that $e_e$ is calculated using the SRIM code, which considers the equilibrium charge state of the SHI. Hence, their values for 104 MeV Pb, 1.2 GeV Pb, and 904 MeV U ions may be different from those given in Table 1 in the region very near the surface. To focus on the effects of single ion impacts, the samples were irradiated with a fluency of $1 \times 10^{11} \text{ cm}^{-2}$.

**Transmission Electron Microscopy.**:- High-resolution TEM, scanning TEM and Fresnel experiments were performed on a double-corrected JEOL ARM 200F operating at 200 keV. In scanning TEM mode, this microscope provides a spatial resolution of 0.078 nm. The Fresnel images were obtained in different focusing conditions with the samples tilted by 20°. The samples were prepared for cross-sectional observations with the focused-ion-beam method using 30 keV Ga⁺ ions in a Helios Nanolab 660 system (Thermofisher). The ion energy was decreased to 2 keV in several steps to obtain samples with thicknesses lower than 100 nm. Finally, the samples were cleaned with a plasma cleaner immediately before being measured.

**Two-Temperature Model—Molecular Dynamics Simulations (TTM-MD).**:- The TTM part of the model was solved using its most usual formulation, where the electrons and lattice atoms were depicted as a set of two heat diffusion equations coupled via an exchange term. The energy deposited by the SHI was introduced by adding a source term in the electronic subsystem equation[37,48,49]

$$C_e(T_e) \frac{\partial T_e}{\partial t} = \nabla(K_e(T_e)\nabla T_e) - g(T_e)(T_e - T_l) + A(\vec{r}, t) \quad (2)$$

$$C_l(T_l) \frac{\partial T_l}{\partial t} = \nabla(K_l(T_l)\nabla T_l) + g(T_l)(T_l - T_l) \quad (3)$$

The subscripts e and l designate the electronic and lattice subsystems, $T_{e,l}$ is the temperature, and $C_{e,l}$ and $K_{e,l}$ are the heat capacity and heat conductivity, respectively. $g(T_e)$ in the exchange term is the electron-phonon coupling. $A(\vec{r}, t)$ is the source term, here computed using the Waligoski formula.[50] Note that this formula includes the velocity effect.[51]

The TTM simulations were performed using the temperature-dependent parameters (i.e., the heat capacity and conductivity of the electrons and lattice as well as the electron-phonon coupling, Figures S19 and S20 of Supporting Information Note S4) and methodology described in Ref. [17]. Recall that for GaN, the parameters for the lattice had been measured[52,53] while the electronic parameters are yet unknown. The temperature-dependent parameters previously derived using the effective mass theory for semiconductors to a free electron gas in the GaN Γ-valley were used. As a consequence, the electron relaxation time $\tau_e$ (i.e., the time at which $T_e$ is reduced to 1/e)

remained the only free parameter of the model, estimated to be 85 fs in Ref. [17]. Note that $\tau_e$, like the entire set of parameters, is independent of the SHI under study. The validity of the TTM was assumed until 90% of the initial electronic energy was transferred into the lattice (i.e., the first ≈85 fs). The obtained temperature profiles at this instant are shown in Figure S18 of Supporting Information Note S3.

The resultant energy profile of the lattice was then deposited instantaneously in a MD cell. This cell was relaxed until reaching room temperature (≈155 ps after the ion impact) using the Albe-Nord potential[54] implemented in the PARCAS code.[55] Typically, SHIs have ranges longer than the thickness of our GaN films. Since for the SHIs discussed here, $e_e$ has a relatively slow decrease with depth,[40] its value is considered constant. To account for surface effects in the interaction dynamics, MD simulations were performed for bulk and surface. In the former, periodic boundary conditions in all three directions were imposed, while in the latter, the periodicity normal to the SHI trajectory was removed. In both cases, a Berendsen temperature control[56] was applied at the boundaries parallel to the SHI trajectory. Additionally, to replicate a real size thin film in the surface simulations, yet keep the computation cost reasonable, a fast thermostat was set at the bottom of the cell to prevent surface effects at this boundary.

The simulations of the Xe irradiations were performed using $30 \times 30 \times 50 \text{ nm}^3$ cells. The simulations of the more energetic SHI (i.e., 0.6 MeV amu⁻¹ Pb, 5.8 MeV amu⁻¹ Pb, and 3.8 MeV amu⁻¹ U ions) were performed in cells with $30 \times 30 \times 100 \text{ nm}^3$. Comparing bulk and surface simulations performed on these cells showed that the latter was capable of describing both bulk and surface regions. Hence, only the surface simulations are shown here. The visualization of the simulations was done using the software OVITO.[57]

The predicted radii were determined by visually inspecting slices with a thickness of 25 nm in the bulk region. This thickness was chosen to match that of the TEM samples, which are between-nm. Additionally, to inspect the stochastic behavior inherent to track dynamics, $r_t$ and $r_c$ were computed along the track using slices with a thickness of 5 nm. The radii obtained using thinner slices were smaller than those observed in thicker ones (see Figure S14 of Supporting Information Note S2). This variation occurred due to the melted and final tracks being irregular cylinders, leading to an apparent enlargement of the tracks when the thickness of the sample increases. Since such an effect was likely presented in the TEM images as well, the radii obtained using the 25 nm thick slices were considered to be the appropriate value to compare with the experimental data. Moreover, it was noted that the effect was noticeable until the thickness of the slices reached 15 nm, above which the apparent radii became constant (see Figure S15 of Supporting Information Note S2). The uncertainties associated with the radii determination were estimated to be half of the difference of the values observed in the 25 nm slices, and the average values obtained in the 5 nm slices rounded to one significant digit. See Note S2 in Supporting Information for further details on the subject.

## Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

## Acknowledgements

Financial support by FCT, Portugal and FEDER is acknowledged (PTDC/CTM-CTM/28011/2017, LISBOA-01-0145-FEDER-028011, UID/05367/2020). M.S. thanks FCT Portugal for his PhD grant (SFRH/BD/111733/2015) and the EU H2020 Project No. 824096 "RADIATE." The authors acknowledge the GANIL for SHI beamtime available under the project CIMAP/IPAC2016/LB/P1110-M-S. The authors thank I. S. Roqan (KAUST) for the GaN samples.

Open access funding enabled and organized by Projekt DEAL.
## Conflict of Interest

The authors declare no conflict of interest.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Keywords

defects, GaN, molecular dynamics, radiation, recrystallization

Received: December 1, 2021

Revised: September 14, 2022

Published online: October 30, 2022

[1] S. Nakamura, S. Pearton, G. Fasol, *The Blue Laser Diode*, 2 ed., Springer Berlin Heidelberg, Berlin, Heidelberg **2000**.

[2] R. S. Pengelly, S. M. Wood, J. W. Milligan, S. T. Sheppard, W. L. Pribble, *IEEE Transactions on Microwave Theory and Techniques* **2012**, *60*, 1764.

[3] S. J. Pearton, F. Ren, E. Patrick, M. E. Law, A. Y. Polyakov, *ECS J. Solid State Sci. Technol.* **2016**, *5*, Q35.

[4] U. Mishra, S. Likun, T. Kazior, Y.-F. Wu, *Proc. IEEE 2008*, **96**, 287.

[5] F. Roccaforte, P. Fiorenza, G. Greco, R. Lo Nigro, F. Giannazzo, F. Iucolano, M. Saggio, *Microelectron. Eng.* **2018**, *187-188*, 66.

[6] N. R. Glavin, K. D. Chabak, E. R. Heller, E. A. Moore, T. A. Prusnick, B. Maruyama, D. E. Walker, D. L. Dorsey, Q. Paduano, M. Snure, *Adv. Mater.* **2017**, *29*, 1.

[7] G. Nikandish, R. B. Staszewski, A. Zhu, *IEEE Access* **2019**, *7*, 57138.

[8] Y. Richard, M.-A. Fernandez, J.-L. Issler, H. Guillon, P. Holsters, A. Barnes, K. Mellab, *Proceedings of the 4S (Small Satellites Systems and Services) Symposium* **2014**, 1.

[9] ESA, European Space Agency - Key components of tree-counting Biomass radar cleared for space, **2018**, http://www.esa.int/Our_Activities/Space_Engineering_Technology/Key_components_of_treecounting_Biomass_radar_cleared_for_space.

[10] M. Rahman, A. Al-Ajili, R. Bates, A. Blue, W. Cunningham, F. Doherty, M. Glaser, L. Haddad, M. Horn, J. Melone, M. Mikuz, T. Quinn, P. Roy, V. O'Shea, K. M. Smith, J. Vaitkus, V. Wright, *IEEE Transactions on Nuclear Science* **2004**, *51*, 2256.

[11] S. Kucheyev, J. Williams, S. Pearton, *Materials Science and Engineering: R: Reports* **2001**, *33*, 51.

[12] P. Ruterana, B. Lacroix, K. Lorenz, *J. Appl. Phys.* **2011**, *109*, 013506.

[13] K. Lorenz, E. Wendler, A. Redondo-Cubero, N. Catarino, M.-P. Chauvat, S. Schwaiger, F. Scholz, E. Alves, P. Ruterana, *Acta Mater.* **2017**, *123*, 177.

[14] E. S. Agency, http://www.spenvis.oma.be/.

[15] A. Tylka, J. Adams, P. Boberg, B. Brownstein, W. Dietrich, E. Flueckiger, E. Petersen, M. Shea, D. Smart, E. Smith, *IEEE Transactions on Nuclear Science* **1997**, *44*, 2150.

[16] D. M. Hassler, C. Zeitlin, R. F. Wimmer-Schweingruber, B. Ehresmann, S. Rafkin, J. L. Eigenbrode, D. E. Brinza, G. Weigle, S. Böttcher, E. Böhm, S. Burmeister, J. Guo, J. Köhler, C. Martin, G. Reitz, F. A. Cucinotta, M. H. Kim, D. Grinspoon, M. A. Bullock, A. Posner, J. Gómez-Elvira, A. Vasavada, J. P. Grotzinger, *Science* **2014**, *343*, 1244797.

[17] M. C. Sequeira, J.-G. Mattei, H. Vazquez, F. Djurabekova, K. Nordlund, I. Monnet, P. Mota-Santiago, P. Kluth, C. Grygiel, S. Zhang, E. Alves, K. Lorenz, *Communications Physics* **2021**, *4*, 51.

[18] P. Kluth, C. S. Schnohr, O. H. Pakarinen, F. Djurabekova, D. J. Sprouster, R. Giulian, M. C. Ridgway, A. P. Byrne, C. Trautmann, D. J. Cookson, K. Nordlund, M. Toulemonde, *Phys. Rev. Lett.* **2008**, *101*, 175503.

[19] W. J. Weber, D. M. Duffy, L. Thomé, Y. Zhang, *Current Opinion in Solid State and Materials Science* **2015**, *19*, 1.

[20] H. Vázquez, E. Áhlgren, O. Ochedowski, A. Leino, R. Mirzayev, R. Kozubek, H. Lebius, M. Karlušić, M. Jakšić, A. Krasheninnikov, J. Kotakoski, M. Schleberger, K. Nordlund, F. Djurabekova, *Carbon* **2017**, *114*, 511.

[21] K. F. Galloway, A. F. Witulski, R. D. Schrimpf, A. L. Sternberg, D. R. Ball, A. Javanainen, R. A. Reed, B. D. Sierawski, J. M. Lauenstein, *Aerospace* **2018**, *5*, 67.

[22] A. Nagamatsu, K. Murakami, A. Yokota, J. Yamazaki, M. Yamauchi, K. Kitajo, H. Kumagai, H. Tawara, *Radiat. Meas.* **2011**, *46*, 205.

[23] A. Debelle, M. Backman, L. Thomé, W. J. Weber, M. Toulemonde, S. Mylonas, A. Boulle, O. H. Pakarinen, N. Juslin, F. Djurabekova, K. Nordlund, F. Garrido, D. Chaussende, *Physical Review B - Condensed Matter and Materials Physics* **2012**, *86*, 100102(R).

[24] A. Kamarou, W. Wesch, E. Wendler, A. Undisz, M. Rettenmayr, *Physical Review B - Condensed Matter and Materials Physics* **2006**, *73*, 184107.

[25] O. Ochedowski, O. Osmani, M. Schade, B. K. Bussmann, B. Ban-Detat, H. Lebius, M. Schleberger, *Nat. Commun.* **2014**, *5*, 3913.

[26] M. Karlušić, R. Rymzhanov, J. O'Connell, L. Bröckers, K. T. Luketić, Z. Siketić, S. Fazinić, P. Dubček, M. Jakšić, G. Provatas, N. Medvedev, A. Volkov, M. Schleberger, *Surf. Interfaces* **2021**, *27*, 101508.

[27] R. A. Rymzhanov, N. Medvedev, J. H. O'Connell, A. Janse van Vuuren, V. A. Skuratov, A. E. Volkov, *Sci. Rep.* **2019**, *9*, 1.

[28] P. Zhai, S. Nan, L. Xu, W. Li, Z. Li, P. Hu, J. Zeng, S. Zhang, Y. Sun, J. Liu, *Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms* **2019**, *457*, 72.

[29] R. A. Rymzhanov, J. H. O'Connell, A. Janse Van Vuuren, V. A. Skuratov, N. Medvedev, A. E. Volkov, *J. Appl. Phys.* **2020**, *127*, 1.

[30] Y. Zhang, R. Sachan, O. H. Pakarinen, M. F. Chisholm, P. Liu, H. Xue, W. J. Weber, *Nat. Commun.* **2015**, *6*, 8049.

[31] M. Backman, M. Toulemonde, O. H. Pakarinen, N. Juslin, F. Djurabekova, K. Nordlund, A. Debelle, W. J. Weber, *Comput. Mater. Sci.* **2013**, *67*, 261.

[32] M. Lang, R. Devanathan, M. Toulemonde, C. Trautmann, *Current Opinion in Solid State and Materials Science* **2015**, *19*, 39.

[33] R. A. Rymzhanov, N. Medvedev, A. E. Volkov, J. H. O'Connell, V. A. Skuratov, *Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms* **2017**, *0*.

[34] M. Sall, I. Monnet, F. Moisy, C. Grygiel, S. Jublot-Leclerc, S. Della-Negra, M. Toulemonde, E. Balanzat, *J. Mater. Sci.* **2015**, *50*, 5214.

[35] S. O. Kucheyev, H. Timmers, J. Zou, J. S. Williams, C. Jagadish, G. Li, *J. Appl. Phys.* **2004**, *95*, 5360.

[36] M. Karlušić, R. Kozubek, H. Lebius, B. Ban-d'Etat, R. a. Wilhelm, M. Buljan, Z. Siketić, F. Scholz, T. Meisch, M. Jakšić, S. Bernstorff, M. Schleberger, B. Šantić, *J. Phys. D: Appl. Phys.* **2015**, *48*, 325304.

[37] M. Toulemonde, C. Dufour, A. Meftah, E. Paumier, *Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms* **2000**, *166*, 903.

[38] S. Porowski, B. Sadovyi, S. Gierlotka, S. J. Rzoska, I. Grzegory, I. Petrusha, V. Turkevich, D. Stratiichuk, *J. Phys. Chem. Solids* **2015**, *85*, 138.

[39] W. Utsumi, H. Saitoh, H. Kaneko, T. Watanuki, K. Aoki, O. Shimomura, *Nat. Mater.* **2003**, *2*, 735.

[40] J. F. Ziegler, J. P. Biersack, U. Littmark, *The Stopping and Range of Ion in Solids*, Springer, Boston, MA **1985**.
[41] G. Schiwietz, P. Grande, *Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms* **2012**, *273*, 1.

[42] S. Mansouri, P. Marie, C. Dufour, G. Nouet, I. Monnet, H. Lebius, *Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms* **2008**, *266*, 2814.

[43] P. Hu, J. Liu, S. Zhang, K. Maaz, J. Zeng, P. Zhai, L. Xu, Y. Cao, J. Duan, Z. Li, Y. Sun, X. Ma, *Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms* **2018**, *430*, 59.

[44] J. Mattei, M. Sall, F. Moisy, A. Ribet, E. Balanzat, C. Grygiel, I. Monnet, *Materialia* **2021**, *15*, 100987.

[45] Y. Zhang, M. Ishimaru, J. Jagielski, W. Zhang, Z. Zhu, L. V. Saraf, W. Jiang, L. Thome, W. J. Weber, *J. Phys. D: Appl. Phys.* **2010**, *43*, 085303.

[46] S. O. Kucheyev, J. S. Williams, J. Zou, C. Jagadish, G. Li, *Appl. Phys. Lett.* **2000**, *77*, 3577.

[47] A. Stukowski, V. V. Bulatov, A. Arsenlis, *Modelling and Simulation in Materials Science and Engineering* **2012**, *20*, 085007.

[48] I. Lifshits, M. Kaganov, L. Tanatarov, *Journal of Nuclear Energy. Part A. Reactor Science* **1960**, *12*, 69.

[49] S. I. Anisimov, Kapeliov.BI, T. L. Perelman, *Zhurnal Eksperimentalnoi I Teoreticheskoi Fiziki* **1974**, *66*, 776.

[50] M. Waligórski, R. Hamm, R. Katz, *International Journal of Radiation Applications and Instrumentation. Part D. Nuclear Tracks and Radiat. Meas.* **1986**, *11*, 309.

[51] W. Wesch, A. Kamarou, E. Wendler, *Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms* **2004**, *225*, 111.

[52] H. Shibata, Y. Waseda, H. Ohta, K. Kiyomi, K. Shimoyama, K. Fujito, H. Nagaoka, Y. Kagamitani, R. Simura, T. Fukuda, *Mater. Trans.* **2007**, *48*, 2782.

[53] K. T. Jacob, S. Singh, Y. Waseda, *J. Mater. Res.* **2007**, *22*, 3475.

[54] J. Nord, K. Albe, P. Erhart, K. Nordlund, *J. Phys.: Condens. Matter* **2003**, *15*, 5649.

[55] K. Nordlund, M. Ghaly, R. Averback, M. Caturla, T. D. Rubia, J. Tarus, *Physical Review B - Condensed Matter and Materials Physics* **1998**, *57*, 7556.

[56] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, *J. Chem. Phys.* **1984**, *81*, 3684.

[57] A. Stukowski, *Modelling and Simulation in Materials Science and Engineering* **2010**, *18*, 015012.