Unravelling the secrets of the resistance of GaN to strongly ionising radiation

GaN is the most promising upgrade to the traditional Si-based radiation-hard technologies. However, the underlying mechanisms driving its resistance are unclear, especially for strongly ionising radiation. Here, we use swift heavy ions to show that a strong recrystallisation effect induced by the ions is the key mechanism behind the observed resistance. We use atomistic simulations to examine and predict the damage evolution. These show that the recrystallisation lowers the expected damage levels significantly and has strong implications when studying high fluences for which numerous overlaps occur. Moreover, the simulations reveal structures such as point and extended defects, density gradients and voids with excellent agreement between simulation and experiment. We expect that the developed modelling scheme will contribute to improving the design and test of future radiation-resistant GaN-based devices. Gallium nitride is a wide bandgap semiconductor which is generally expected to replace some silicon-based technologies, despite some of its properties still requiring further investigation. Here, using a two-temperature model coupled to molecular dynamics simulations, the authors investigate and predict the effects of strongly ionising radiation in gallium nitride, revealing the mechanism behind its unusual resistance to radiation.

G roup III nitrides have been the semiconductors of choice for some of the most notable technological achievements of the last few decades. They have found application in devices such as light-emitting diodes, laser diodes and highelectron mobility transistors (HEMTs). More recently, they are raising the standards once again, this time in the field of radiation-hard electronics. For instance, when comparing GaNwith Si-based HEMTs, 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 1,2 . Space and high-energy physics are two examples of areas with ongoing intense research that could benefit enormously from these more efficient and enduring devices. The European Space Agency is currently employing a new generation of space communication systems based on GaN HEMTs in its Porba-V satellite, so far delivering greatly improved data transmission rates and working reliability 3 , as well as in its Biomass satellite due for launch in 2022 4 . The extreme radiation present in these environments introduces significant challenges that should be addressed in order to advance such technologies.
While the effect of low-to medium-energy ions in GaN is relatively well understood [5][6][7] , the high-energy radiation regime is still under debate. In space, the latter is of special interest when studying the influence of the so-called HZE (high atomic number Z and high energy) ions, a crucial component of the galactic cosmic radiation (GCR) and solar radiation spectra. Although less abundant than proton and gamma radiation, HZE ions (such as Fe and Si) have a very high ionisation power and relatively high fluxes [8][9][10] . This likely not only accelerates the degradation but also can cause single-event catastrophic failure of devices 11 .
In this communication, we address fundamental aspects of the interaction of strongly ionising particles with GaN, namely its response to the so-called swift heavy ion (SHI) radiation. With typical energies in the MeV and above range per nucleon, SHIs lose their energy in the target primarily via electronic interactions. In semiconductors, this often induces ionisation spikes along their paths that promote the formation of latent tracks 12 and surface effects 13 . Moreover, it has been reported that SHI radiation can lead to the recovery of pre-damaged material [14][15][16][17][18] , including the recrystallisation of amorphous materials.
Recently, conflicting reports have been published on the matter of strongly ionising radiation effects on GaN. Sall et al. 19 used transmission electron microscopy (TEM) to study GaN films irradiated with SHI and concluded that those are capable of amorphising the crystal. This result opposes the one obtained previously by Kucheyev et al. 20 , who used Rutherford backscattering spectrometry/channelling (RBS/C) to show that even when GaN films are irradiated with high fluences (above 1 × 10 12 cm −2 ), they presented a defective but still crystalline structure. The situation becomes more surprising when we realise that Kucheyev et al. studied SHIs with higher energies than some studied by Sall et al. Moreover, Kucheyev et al. observe that their experimental data follow neither Poisson nor Gibbons models 21 , which already points to the existence of complex damage formation and annihilation mechanisms.
Here, we identify the processes leading to this apparent contradiction by applying a two-temperature model-molecular dynamics (TTM-MD) simulation scheme capable of reproducing the experimental results obtained with both techniques, TEM and RBS/C, and also by electron energy loss spectroscopy (EELS). The TTM is one of the most used approximations to describe the first femtoseconds of the SHI interaction with matter. Within this model, the ionisation spike is depicted as a hot electron gas of temperature T e that quickly relaxes its energy into the target lattice via electron-phonon interactions 22 . However, the TTM does not take into consideration the structural dynamics of the crystal, hence hiding compelling phenomena such as phase transitions or defect dynamics. Thus, the TTM output (i.e., the energy profile of the lattice) is deposited into a MD cell. This simulation scheme, often referred to as TTM-MD, has found success in describing the interaction of SHIs with different materials (e.g., 17,18,23 ). A detailed discussion on the implementation of the TTM for GaN is presented in the Methods section and Supplementary Note 1.

Recrystallisation.
We begin by studying 3-μm thick GaN films on a Al 2 O 3 substrate irradiated with 185 MeV Au ions (with an electronic energy loss of ϵ e = 33 keV/nm as calculated using the SRIM code 24 ). The simulations describing this SHI show that the initial spike induces a solid-liquid phase transition in the crystal, forming a cylindrical track along the SHI path. The radius of this molten track increases radially, reaching its maximum value of r t = 4.4 nm, 13 ps after the SHI impact, as Fig. 1a shows. The high lattice temperature (up to 8000 K), along with the strong confinement of the molten track within the crystalline GaN matrix, gives rise to a high-pressure state. This opens the possibility for the track to recover via recrystallisation similar to what was shown experimentally by X-ray diffraction during the cooling of high-pressure melted GaN 25 . Indeed, at 13 ps, the recrystallisation starts at the edge of the track and proceeds inwards as the system cools down. Due to the high-pressure state, extended defects begin to form as well. When the crystal finally reaches room temperature, Fig. 1b, the morphology of the track consists of a damaged, yet highly crystalline, shell with an outer radius equal to r t and a significantly smaller amorphous core with a radius r c = 1.8 nm (see also Supplementary Movie 1 to visualise the entire recrystallisation process).
Simulated and experimental plan-view TEM in a high-angle annular dark-field setup (TEM-HAADF) images for 185 MeV Au ions are shown in Fig. 2a, b, respectively. We found an exceptional dimensional and morphological agreement. The predicted values for r t and r c , have a straightforward correspondence to the outer (less damaged) and the inner (more damaged) radii of the track seen in the experimental TEM-HAADF image. Moreover, the defects in the outer part of the track induce misoriented crystallographic domains in both simulated and experimental images, as evidenced in the magnifications in Fig. 2. This similarity suggests that the simulations reproduce the nature and influence of the extended defects accurately.
Defect formation. To further study the defects created, we apply the dislocation extraction algorithm (DXA) 26 to the simulation cell. As shown in Fig. 3a, the DXA reveals the existence of dislocation edges (green lines) and stacking faults (orange lines) in the recrystallised volume. Comparing the simulated and experimental cross-sectional TEM-HAADF images, Fig. 3a, b, respectively, reveals a striking resemblance between them.
In addition, the incorporation of atoms into the extended defects leads to a distinct density profile, as the overlay in Fig. 3a, b shows. The ion track morphology consists of a low-density core surrounded by a high-density shell. Using EELS, we found the predicted density profile to be in excellent agreement with the experimental one, Fig. 3c. Moreover, the obtained density of the amorphous core in the simulations is of 5.1 g/cm 3 , comparable with that of amorphous GaN predicted by density functional theory 27 . These results strengthen the validity of the predicted dynamics, namely the recrystallisation and the formation of extended defects.
We suggest that the formation of the extended defects and concomitant density gradient is likely responsible for interrupting recrystallisation. After some time, the density in the melted volume reaches a critical point where recrystallisation becomes inviable.
Surface effects. Near the surface, there is an important additional mechanism of pressure relaxation. While in the deeper regions, the pressure within the track can only relax radially, near the surface, it can relax upwards as well. This difference creates a pressure gradient along the SHI path that ultimately leads to the sputtering of a considerable amount of Ga and N atoms. Consequently, as shown in Fig. 4a, near the surface, we found not only voids within the track but also a hillock on the top of it. Despite the sputtering, a considerable amount of recrystallisation still occurs, particularly in the deeper regions of the cell. Comparing the simulated and experimental crosssectional TEM-HAADF images, Fig. 4b, c, respectively, an excellent agreement is found. Not only the presence of very low-density pockets, presumably voids, inside a highly damaged track is evident in the experiment, but also their shape, density and distribution are simulated accurately. The largest lowdensity pocket observed in the experimental image (Fig. 4c) has a lateral dimension of approximately 4 nm, comparable with the predicted one of about 3.6 nm (Fig. 4a, b). Moreover, the voids are elongated along the track direction in both experimental and simulated images. Regarding their density inside the track, we obtain the experimental value of approximately 0.30-0.35 voids/nm, neighbouring the 0.35 voids/nm obtained from the simulations near the surface.
Note that the differences observed in the surface and bulk simulations are solely due to surface effects. The energy loss is assumed to be equal in both simulations, a fair assumption since the thickness of the GaN films is considerably smaller than the range of the SHIs.
Recrystallisation upon overlapping impacts. Remarkably, efficient recrystallisation also occurs when a SHI impacts predamaged material. This has substantial implications in the damage evolution with fluence. As shown in Fig. 5, the amorphous core of an ion track is almost entirely recrystallised upon the overlapping impact of a second SHI. Interestingly, this phenomenon occurs even when the track overlap is not direct, i.e., when the amorphous track core from the first SHI is slightly outside of the melted volume induced by the second SHI. This implies that the recrystallisation cross-section, with radius r r , is greater than the initial melting cross-section (r t = 4.4 nm). By performing simulations for tracks at different distances, we estimate r r to be approximately 6.5 nm. Despite this recovery being substantial for direct and indirect overlaps, the formation of extended defects is enhanced in the latter case.
Surprisingly, even the voids seen in the surface simulation are partially recovered when overlapped by another SHI. This recovery is more efficient in the deeper regions leading to a density gradient of voids with depth. Damage evolution. The pronounced lattice recovery within a single track and upon multiple impacts makes it unlikely for the crystal to become fully amorphised, even for high fluences. In Fig. 6a, we present the RBS/C measurements of samples irradiated with different fluences of 185 MeV Au ions. As expected, the backscattering yield in the aligned spectra increases with the ion fluence as the defect concentration rises. Nonetheless, even for the highest fluence of 1 × 10 13 cm −2 , the backscattering yield stays well below the random level revealing that the single-crystalline nature of the sample is preserved. For this fluence, Poisson statistics indicate that the whole sample has suffered a high number of impact overlaps. This resistance to radiation follows the results of Kucheyev et al. 20 suggesting that for 200 MeV Au irradiation, no amorphisation of GaN occurs even for high fluences.
To further interpret the RBS/C results, it is paramount to take into consideration the full description of the track morphology, in particular, the presence of extended defects, surface voids and density gradients. None of these structures is considered by prevalent RBS/C models, such as the two-beam approximation 28 or the master equation approach 29 , despite impacting the measurements significantly. Extended defects, in particular, are a well-known issue in RBS/C analysis due to the enhanced dechannelling effect they create 6 . If not taken into account, the models often result in inaccurate damage-depth profiles 30 .
Regarding the influence of density gradients and voids in RBS/ C, to the best of our knowledge, they have not yet been studied. To address these issues, we apply the RBSADEC code 31 , a Monte Carlo implementation of a binary collision algorithm, that replicates the RBS/C experiment in cells with atoms in arbitrary positions. Thus, it is suitable to predict the RBS/C spectra from the final state of our MD simulations, taking into account all the morphological features of the track inherently. A realistic cell for the RBSADEC simulations is created by repeating and joining the final cells of the MD surface and bulk simulations. To account for overlap effects, for each fluence, we determine the appropriate lateral size and the number of impacts inside the MD cell using Poisson statistics. The distance between impacts, an essential parameter due to the recrystallisation dynamics, is chosen based on numerical calculations. Refer to Methods section and Supplementary Note 2 for further details on the implementation of the RBSADEC code. Figure 6 shows an excellent agreement between simulated and experimental RBS/C spectra. Comparison with simulations using only the bulk MD cell reveals that surface voids lower the backscattering yield of the aligned spectrum, as observed in the RBS/C spectra shown in Fig. 6 for the first 300 nm (i.e., energies above 1350 keV). Furthermore, we note that the density of voids is not constant with depth and depends on the fluence. Although all the impacts induce voids until approximately the same depth (approximately 320 nm in the case of the 185 MeV Au ion), the deeper voids formed by previous impacts are recovered in the event of an overlap. Therefore, the areas with more overlaps (e.g., two impacts) have a smaller fraction of deep-voids than those with fewer overlaps (e.g., one impact). As the number of overlaps increases, this effect becomes more pronounced and leads to an inhomogeneous depth distribution of voids, which has a significant influence on the aligned spectra.
Note that there is no fitting involved except for the mentioned surface effects. The RBS/C spectra beyond 320-nm depth are directly simulated from the untreated bulk MD cell. In line with the TEM results, the spectra strongly imply that the simulated track morphology is capable of describing the irradiated samples even for high fluences. Consequently, the dynamics leading to this morphology, namely recrystallisation (with and without overlapping) and sputtering, do occur when the SHIs impact the crystal.
Galactic cosmic radiation. To validate the established modelling scheme (including the thermodynamic parameters) also for lower SHI energies, closer to those found in GCR, we apply the same procedure to study GaN films irradiated with 70 and 45 MeV Xe ions (ϵ e = 19.9 and 14.9 keV/nm 24 , respectively) with a fluence of 2 × 10 12 cm −2 .
The simulations demonstrate a track dynamics similar to that described above for the 185 MeV Au ion but with a considerably smaller track radius. After 6 ps, r t reaches its maximum value of 2.8 and 2.2 nm for the 70 and 45 MeV Xe ions, respectively, instead of the 4.4 nm obtained for the more energetic Au ions. The recrystallisation process in both Xe radiations recovers most of the initial damage and forms considerably fewer extended defects than the Au radiation. When the crystal reaches room temperature, the tracks present discontinuous cores consisting mostly of small amorphous pockets (r c < 1 nm) that get almost completely recrystallised when overlapped by another SHI impact. Although sputtering still occurs, the presence of voids is very superficial, with depths of 18 and 5 nm for 70 and 45 MeV  We present the RBS/C results in Fig. 6b. Remarkably, the results reveal that both Xe SHIs induce nearly no damage in the crystal, even for a relatively high fluence of 2 × 10 12 cm −2 . The yields of the aligned spectra of the samples irradiated with 45 MeV Xe ions and as-grown samples are comparable within the experimental uncertainty. The RBSADEC simulations are in good agreement with the experiment. The slight overestimation of the yields is most likely related to the high defect mobility of point defects in GaN, shown to be mobile at 120 K 32 , leading to defect recovery at time scales not represented by the MD simulations.

Discussion
The excellent qualitative and quantitative agreement between the simulations, TEM, EELS and RBS/C results indicates that the employed model has a significant prediction power. The simulations unveil a strong recrystallisation effect that is responsible for recovering most of the initial melted region induced by the SHI. Interestingly, the effect can recover the damage produced by previous ion impacts, including the amorphous cores of the tracks. This dynamical recovery is of paramount relevance since it lowers the damage evolution with the irradiation fluence considerably, preventing the complete amorphisation of the crystal. Moreover, it is also responsible for the unusual resistance of GaN to strongly ionising radiation.
The model can be used to inspect a broad range of applications and predict the behaviour of GaN under extreme radiation conditions, including the formation of defects such as point and extended defects, amorphous cores, density gradients and voids. This eases the burden of going through complicated and expensive experimental work that characterises this radiation regime. Considering space applications, the multiplicity of ion species and energies natural of this environment can further aggravate these studies. We note that the ϵ e of the Fe component of the HZE radiation peaks at 14 keV/nm 8,9,24 , which is close to the one of the 45 MeV Xe ions studied here. Therefore, we anticipate that important device properties, such as rate of failure, dose thresholds, shielding design, among others, can be better estimated and understood. For instance, our results for the 185 MeV Au radiation are comparable with those reported 33 where GaN-based HEMTs were irradiated with 1540 MeV Bi ions (ϵ e = 45.9 keV/nm). Considering the use of SHIs as a tool for material engineering 34 , our model can give new insights into the influence of the radiation parameters on the planned design.

Experimental details
Materials. Commercial c-GaN thin films (Lumilog), 3-μm thick, grown by Metal Organic Chemical Vapour Deposition on Al 2 O 3 substrates were used.
Radiations. The samples were irradiated with 185 MeV Au 13+ (ϵ e = 33 keV/nm according to simulations done with the Monte Carlo code SRIM 24 ) ions at the Australian National University Heavy Ion Accelerator Facility with fluences of 10 11 , 10 12 and 10 13 cm −2 at room temperature. Moreover, an energy study was performed using a 90 MeV Xe 23+ beam at the Grand Accélérateur National d'Ions Lourds (GANIL). The beam energy was degraded to 70 and 45 MeV Xe (ϵ e = 19.9 and 14.9 keV/nm 24 , respectively) using Al foils of different thickness. The fluence for both energies was 2 × 10 12 cm −2 . All irradiation shown here were performed with a zero-angle incidence with respect to the surface normal. Note that ϵ e is calculated using the SRIM code, which does not consider the charge state of the SHIs. Hence, their values may be smaller than the ones given above in the region very near the surface.
Transmission electron microscopy-high-angle annular dark-field (TEM-HAADF). High-resolution TEM and scanning TEM experiments were performed on a double corrected JEOL ARM 200F operating at 200 kV. In scanning TEM mode, this microscope provides a spatial resolution of 0.078 nm. STEM analyses were carried out using an HAADF detector at the inner and outer angles of 68 and 280 mrad to reduce diffraction contrast efficiently. The samples were prepared for plan-view and 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 the samples with thicknesses lower than 100 nm. Finally, the samples were cleaned with a plasma cleaner immediately before being measured.
Electron energy loss spectroscopy (EELS). EELS experiments were performed with a collection angle of 90 mrad. Elemental mapping was recorded with a dispersion of 1 eV/channel. HR-TEM/STEM micrographs and EELS data were acquired and processed under Gatan Digital Micrograph environment.
Rutherford backscattering spectrometry/channelling (RBS/C). The RBS/C measurements were performed at LATR (Laboratório de Aceleradores e Tecnologias de Radiação) IST using a Van de Graaff accelerator. A 2 MeV He + probe beam was used, and the backscattered particles were detected using a pin-diode detector at a backscattering angle of 165 ∘ . The aligned measurements were performed along the c-axis of the wurtzite crystals with the help of a two-axes goniometer. For the random spectra, to avoid channelling, the samples were tilted by 5 ∘ from the channelling direction and the azimuthal angle rotated during the measurement. Due to the similarities between the random spectra, only the spectrum from the asgrown sample is shown here. The presence of defects increases the backscattering yield in the aligned spectra. In the extreme case of complete amorphisation, the aligned spectrum would reach the random one, since the channelling effect becomes impossible. Therefore, comparing the two spectra allows for a quantification of the defect density in the crystals. The energy-depth scale is obtained considering the energy loss of the probe ions in the material.
Simulation details Two-temperature model (TTM). Being a thermodynamic model, the TTM depends on the heat capacity and conductivity, C and K, respectively, of the lattice and electrons as well as on the electron-phonon coupling, g. For GaN, the parameters for the lattice have been measured 35,36 . In contrast, the electronic parameters are yet unknown. In the context of the TTM, the conduction electrons are usually assumed to have a metal-like behaviour and C e and K e are considered constant (e.g., 18,19,22 ) or have a linear dependence with the electronic temperature, T e 13 . However, it has been shown that even for pure metals, this approximation can be misleading 37 , and even more questionable for semiconductors as discussed for SiC 38 . Here we introduce a simple T e dependence in C e , K e and g by applying the effective mass theory for semiconductors to a free electron gas in the GaN Γ-valley. As a consequence, the electron relaxation time τ e (i.e., the time at which T e is reduced to 1/e) remains as the only free parameter of the model, here estimated to be 85 fs. A detailed discussion on the implementation of the TTM for GaN is presented in Supplementary Note 1.
The use of equilibrium thermodynamic quantities to describe a process far from equilibrium must be taken with some care. Yet, the excellent agreement between experimental and theoretical results presented here and previously for other semiconductors and insulators 12,39,40 suggests that the simplicity and low computational costs offset the small deviations that may arise due to nonequilibrium physics. However, it is important to note that the TTM does not take into consideration the structural dynamics of the crystal hence hiding compelling phenomena such as phase transitions or defect dynamics. Thus, in this work, its validity is assumed only until 90% of the initial electronic energy has been transferred into the lattice (i.e., the first~80 fs). We then transfer the energy profile of the lattice, resultant from the TTM, into a MD cell (see Supplementary Note 1B).
MD simulations (MD). The MD cells are relaxed until reaching room temperature (155 ps after the ion impact) using the Albe-Nord potential 41 implemented in the PARCAS code 42 . Typically, SHIs have ranges longer than the thickness of our GaN films. Since for the SHIs under study, ϵ e has a relatively slow decrease with depth 24 , its value is considered constant. To account for surface effects in the interaction dynamics, we perform MD simulations for bulk and surface separately. In the former, periodic boundary conditions in all three directions are imposed while in the latter, the periodicity normal to the SHI trajectory is removed. In both cases, a Berendsen temperature control 43 is applied at the boundaries parallel to the SHI trajectory. In addition, in order to replicate a real size thin film in the surface simulations, yet keeping the computation cost reasonable, a fast thermostat is set at the bottom of the cell to prevent surface effects at this boundary. All simulations are done using 30 × 30 × 50 nm 3 cells. The visualisation of the simulations is done using the software OVITO 44 .
Transmission electron microscopy-high-angle annular dark-field (TEM-HAADF) simulations. To better compare the TTM-MD simulated track with TEM-HAADF images, we note that this technique is highly sensitive to the atomic number of the atoms of the sample (∝ Z 2 ). Thus, we compute the average of Z 2 at each atom site, considering its neighbouring atoms within a cutoff radius of 1 nm. This value is then represented using a greyscale that is linearly adjusted to match the TEM-HAADF image contrast.
Electron energy loss spectroscopy (EELS) simulations. The density profile is obtained by determining the atomic density within vertical slices (each with a width of 0.5 Å and over the entire 30 nm depth of the MD cell) perpendicular to the plane of view (i.e., the experimental scan direction). By calculating the average over a window of multiple slices and plotting it at the centre of each slice, it is possible to introduce the experimental resolution and therefore compare the simulation with the EELS measurements. In Fig. 3c, we show the simulated densities obtained using average windows with a total width of 1.3 nm (orange line) and 6 nm (green line).

RBSADEC simulations.
To create a realistic cell for RBS/C simulations using the TTM-MD results, we need to consider the effects of the fluence and the depth inhomogeneities. A given fluence can be described by the number of impacts in the cell, N, and the distance between them, D. Note that D is an essential parameter due to the observed distinction between direct and non-direct overlap effects in the final track morphologies. N is calculated according to Poisson statistics. Regarding D, we calculate it by simulating the random position of a large number of impacts (corresponding to the fluence) and computing the average distance between each impact and its nearest neighbour. Afterwards, we run TTM-MD simulations (bulk and surface) for each set of N and D parameters. The depth inhomogeneities are introduced by combining the cells from the bulk and surface simulations. For the deeper regions, the entire bulk cell is repeated until a thickness of t b . For the surface region, we slice the surface cell into thinner subcells and repeat each of these to obtain the desired thickness, t s . The thickness of the final cell, t b + t s , is set to 1200 nm, the typical depth seen in RBS/C. As mentioned in the text, t s assumes a distribution instead of a fixed value. This distribution is introduced by performing RBS/C simulations using n different t i s , with n depending on the fluence. The obtained spectra, Sðt i s Þ, are then averaged to produce the total spectrum S total ¼ 1 n P n i Sðt i s Þ. For most fluences, besides using a minimum and maximum t i s (50 and 320 nm, respectively), we add a simulation using t i s ¼ 205 nm to account for the fact that t s assumes values tendentially closer to its maximum.
Except for the surface part, the same process is applied to study the samples irradiated with Xe ions (with a fluence of 2 × 10 12 cm −2 ). Since the surface morphology is fully contained in the surface MD cell (see DXA in Fig. 6), we can simply append this cell to the bulk cells.
For further details on the implementation of the RBSADEC, please refer to Supplementary Note 2.

Data availability
Supplementary Movie 1 shows the evolution of the melted track formed after the impact of a 185 MeV Au ion from which the plan-view images shown in Fig. 1 are taken. Further material that supports the findings of this study is available from the corresponding author upon reasonable request.