Universal radiation tolerant semiconductor

Radiation tolerance is determined as the ability of crystalline materials to withstand the accumulation of the radiation induced disorder. Nevertheless, for sufficiently high fluences, in all by far known semiconductors it ends up with either very high disorder levels or amorphization. Here we show that gamma/beta (γ/β) double polymorph Ga2O3 structures exhibit remarkably high radiation tolerance. Specifically, for room temperature experiments, they tolerate a disorder equivalent to hundreds of displacements per atom, without severe degradations of crystallinity; in comparison with, e.g., Si amorphizable already with the lattice atoms displaced just once. We explain this behavior by an interesting combination of the Ga- and O- sublattice properties in γ-Ga2O3. In particular, O-sublattice exhibits a strong recrystallization trend to recover the face-centered-cubic stacking despite the stronger displacement of O atoms compared to Ga during the active periods of cascades. Notably, we also explained the origin of the β-to-γ Ga2O3 transformation, as a function of the increased disorder in β-Ga2O3 and studied the phenomena as a function of the chemical nature of the implanted atoms. As a result, we conclude that γ/β double polymorph Ga2O3 structures, in terms of their radiation tolerance properties, benchmark a class of universal radiation tolerant semiconductors.

Long-range periodicity or translation symmetry is a unique property of solids, even though solids may form amorphous phases too. In this context, accelerated particle beam irradiations are known to induce amorphization in many types of crystals, e.g. in semiconductors 1,2 . In its turn, radiation tolerance in semiconductors is determined as an ability to withstand the accumulation of the radiation disorder, otherwise leading to highly disordered lattice and, upon irradiating with sufficiently high fluences, to amorphization 3,4 . The irradiationinduced disordering mechanisms are generic, even though exhibiting material-specific differences, allowing us to classify semiconductors as low-or high-radiation tolerant [5][6][7][8] . Importantly, very recently, it was discovered that in gallium oxide (Ga 2 O 3 ), which is a promising material for the next generation power electronics [9][10][11][12][13] , the amorphization may be prominently suppressed by the formation of a new metastable crystalline polymorph phase 14 . This process occurs in the irradiation interaction volume and results in a new polymorph film, separated from the initial polymorph by a sharp interface.
In this work, we report that such double polymorph Ga 2 O 3 structures exhibit high radiation tolerance. Specifically, for room temperature experiments, these samples tolerated a disorder equivalent to hundreds of displacements per atom (dpa), without severe degradations of the crystallinity. For comparison, other semiconductors studied in literature for comparative dpa either amorphize or exhibit a high degree of lattice disorder. Notably, to induce such high dpa we use high fluence ion irradiations creating high excess of implanted atoms, maximized at the depth of the ion range, and affecting the process depending on the chemical nature of the implanted atoms. Figure 1 illustrates such high radiation tolerance of the double polymorph Ga 2 O 3 structures, tolerating up to 265 dpa (see Supplementary Note 1 for the dpa calculations) without severe degradation in crystallinity (panels a-c), set in a context of the same characteristics in other semiconductors (panel d). Figure 1a plots the Rutherford backscattering spectrometry in channeling mode (RBS/C) spectra of the double polymorph Ga 2 O 3 structures as a function of the fluence in the range of 1 × 10 16 to 1 × 10 17 Ni/cm 2 . As explained in Methods and in more details in Supplementary Note 2, the characteristic shape observed for the 1 × 10 16 Ni/cm 2 RBS/C spectrum is a fingerprint of the high crystallinity of the double polymorph gamma/beta (γ/β) Ga 2 O 3 structure. Thus, the data in the range of 300-400 channel numbers for 1 × 10 16 Ni/cm 2 implants correspond to the least disordered γ-Ga 2 O 3 in Fig. 1a, which we adapt as a "reference" disorder level to compare with the data for higher fluences. Remarkably, the characteristic shape of the RBS/C spectra is maintained for 3 × 10 16 Ni/cm 2 and for 5 × 10 16 Ni/cm 2 implants, corresponding to 80 dpa and 132 dpa, respectively. Moreover, a minor deviation from the trendrelated to an enhanced RBS/C yieldobserved for the 1 × 10 17 Ni/cm 2 implants (dpa = 265), is related to an increase in the Ni content, with no significant changes in crystallinity of the surrounding matrix. It is evident from the comparison of the corresponding channeling and random spectra in Fig. 1a and from the scanning transmission electron microscopy (STEM) data in Fig. 1b and c, that the sample retains exceptionally high crystallinity even after 1 × 10 17 Ni/cm 2 implants. In fact, the selected area electron diffraction (SAED) indexation of the sample (Fig. 1b) confirms its identification as γ-Ga 2 O 3 /β-Ga 2 O 3 double-layer structure (see Supplementary Note 3). In its turn, we observed 3-6 nm diameter Ni precipitates embedded into the γ-Ga 2 O 3 layer. These are shown in Fig. 1c with a magnified annular dark field (ADF)-STEM image taken along [100] γ-Ga 2 O 3 axis, resolving the precipitates in a brighter contrast (see detailed analysis in Supplementary Note 4). Additionally, the γ/β Ga 2 O 3 interface exhibits stacking with rather low lattice mismatch (Supplementary Note 5). Thus, based on the data in Fig. 1a-c, γ-Ga 2 O 3 tolerates up to 265 dpa without severe degradations in the crystallinity of the semiconductor matrix exposed to Ni implants. This is a remarkable result in particular, when compared with literature data, see the data summarized in Fig. 1d. Indeed, as can be seen from Fig. 1d, such materials as Si, SiC, InP have previously been shown to amorphize already at 0.2-0.5 dpa [15][16][17] 15 , SiC (diamonds) 16 and InP (down triangles) 17 ) and radiation tolerant semiconductors (GaN (squares) 16 , and AlN (up triangles) 18 ) for Au implants at room temperature, as well as for Ga 2 O 3 (this work, circles) -the lines are a guide to the eye. Notably, the depth scale in panel (a) is calculated for Ga atoms, so that the Ni related peak appears deeper in the sample (see Supplementary Note 2 for clarity). Source data are provided as a Source data file. radiation tolerant materials, e.g. GaN or AlN, that are capable to accommodate much higher radiation disorder 16,18 and remaining crystalline. However, none of these materials remains such excellently crystalline as γ-Ga 2 O 3 , see Fig. 1d. Notably, β-Ga 2 O 3 belongs to the low radiation tolerant group of materials. However, we observe that the disorder accumulation in β-Ga 2 O 3 lattice does not result in full amorphization but triggers transformation to a new crystalline polymorph 14,[19][20][21] , as illustrated with arrows showing the trend for converting the irradiated β-Ga 2 O 3 volume into radiation tolerant γ-Ga 2 O 3 /β-Ga 2 O 3 double polymorph structure (see also Supplementary Note 2). Notably, there is a gradual increase in the γ-Ga 2 O 3 thickness as a function of fluence, as highlighted in Fig. 1a by the corresponding dashed lines. This thickness increase is consistent with our hypothesis of the disorder induced β-to-γ-Ga 2 O 3 transformations 14 and the data in Fig. 1a may be used to estimate the corresponding disorder thresholds (see Supplementary Note 6). Further, the fact that the Ni content in Fig. 1b, c was sufficient for the precipitation, implies that one has to account for the chemical nature of the implanted atoms, potentially altering the defect accumulation and eventual amorphization processes in γ-Ga 2 O 3 , as it may occur in other materials too 8 . Thus, for comparison, we investigated these phenomena for several other ions, choosing elements having strongly different chemical capabilities to interact with the matrix atoms. For that matter, Fig. 2 shows examples of the STEM data taken upon the implants resulting in the same dpa range (86-88 dpa) for Au and Ga ions. Importantly, as seen from Notably, Fig. 2g shows a high magnification ADF-STEM image of the interface between γ-Ga 2 O 3 and the amorphous phase. The corresponding fast Fourier transforms (FFT)s in the insets of Fig. 2g confirm that γ-Ga 2 O 3 is oriented along the [100] zone axis while the FFT in the amorphous phase shows the features that are characteristic of amorphous materials. Moreover, the difference in atomic coordination of γ-, βand amorphous Ga 2 O 3 phases is illustrated by the fine structure of the oxygen-K edge in the electron energy-loss spectroscopy (EELS) spectra shown in Fig. 2h. In particular, the oxygen K-edge is characterized by two peaks at 538 eV and 543 eV 22 and the relative intensity of these peaks apparently changes depending on the localization of the measurements 23 Supplementary Fig. 10). Apparently, this feature vanishes and a shape characteristic to the γ-Ga PRDF becomes evident with an increasing number of FPs. The observed change manifests the β-to-γ Ga 2 O 3 phase transformation with an increase of Ga-type defects in β-Ga 2 O 3 , while a similar damage level in γ-Ga 2 O 3 does not result in any significant modification of the γ-Ga PRDF. Additionally, Fig. 3b illustrates the structural differences in β-Ga 2 O 3 (up) and γ-Ga 2 O 3 (down) before and after the introduction of 600 FPs, where the dramatic changes-compared to the initial cellare seen only in β-Ga 2 O 3 (for the more detailed transition process, see Supplementary Note 8, Supplementary Fig. 11). From quantitative comparison of the shapes of PRDFs within the 1st and the 2nd shells separately for both phases before and after introduction of Ga FPs (Supplementary Note 8, Supplementary Fig. 12), we deduce that only the damaged β-Ga PRDF within the 2nd shell underwent the most distinct shape modification. To compare these changes to the γ-Ga PRDF, we map in Fig. 3c the β-Ga PRDF values for the different FP numbers against the corresponding values of the pristine γ-Ga PRDF. This analysis confirms that the shape of the damaged β-Ga PRDF with increase of Ga FPs indeed approaches that of the pristine γ-Ga PRDF. In Fig. 3d we plot the Pearson correlation coefficients (Pr) versus the numbers of FPs, comparing the shapes of the β-Ga at different number of FPs for the pristine βand γ-Ga PRDFs (violet and brown Pr curves in Fig. 3d, respectively). The comparison reveals a high degree of positive correlation (similarity) for the damaged β-Ga PRDF with that of the γ-Ga PRDF after a threshold number of FPs, at~200 FPs per cell (~0.15 dpa) when the β-Ga 2 O 3 phase inevitably transforms into γ-Ga 2 O 3 phase. This is in good agreement with the experimental data (Fig. 1d). Moreover, we see that the FPs have only marginal effect on the γ-Ga PRDFs, as shown in Fig. 3b and Supplementary Note 8, Supplementary  Figs. 11-12, perfectly matching the strikingly high radiation tolerance of the γ-Ga 2 O 3 observed in our experiments.

Results and discussion
To verify the insensitive response of the O-O and Ga-O sublattice to the introduction of FPs observed in our MD simulations of damage accumulation, we performed dynamic single-cascade MD simulations, where the O and Ga atoms were naturally displaced in collision cascades. In these simulations, we see that the O sublattice is highly rigid and strongly prone to recrystallization into the face-centered-cubic (fcc) stacking, despite the stronger displacement of O atoms compared to Ga during the active periods of cascades, see Supplementary Note 8, Supplementary Fig. 13.
Furthermore, we study eventual chemical effects on the accumulation of structural disorder in γ-Ga 2 O 3 using ab initio MD (AIMD). Figure 4a illustrates the PRDFs separately for the O-O and the heavyion sublattices compared to the respective initial PRDFs. The heavy-ion sublattice includes the native Ga and the added Ni, Au, or Ga atoms (~10 at.%). In these simulations, the extra Ni, Ga, or Au atoms were added in random locations between the lattice sites imitating the implantation of ions under high-fluence irradiation. After that the structures were thermally equilibrated using AIMD to obtain the most energetically favorable structures (more data in Supplementary Note 8, Supplementary Figs. 14-15). Since the presence of LRO peaks in an RDF can be used as a crystallinity measure, we analyse both the O-O and heavy ion PRDF fluctuations around the unity (grey dotted line at g(r) = 1 in Fig. 4a) beyond the SRO peaks separated by the vertical grey dotted lines at distances 3.6 Å and 4.0 Å for the O-O and the heavy ions pairs, respectively. In Fig. 4b we quantify the degree of amorphization in the lattices with the implanted ions by integrating the total deviation area of the PRDF curves from the dotted lines (g(r) = 1) beyond the SRO peaks in Fig. 4a by the vertical grey dotted lines. The smaller the deviation area the higher the degree of amorphization the structure exhibits. Naturally, the PRDFs of the stoichiometric γ-Ga 2 O 3 with multiple peaks and valleys along the g(r) = 1 line have the largest deviation area. Remarkably, the strongest disordering effect of the implanted atomsthe smallest deviation areais observed in the cell with the Ga excess, which is in excellent agreement with the experiments (Fig. 2). However, the deviation area for the Ga-Au/Ga PRDF is only marginally larger than that of the Ga-Ga/Ga PRDF. Hence, we further analyse the disorder in the implanted lattices by comparing the bond-angle distributions for the O-O bonds for all three distorted structures with the pristine one in Fig. 4c and the Pr similarity analysis in Fig. 4d. The visual inspection of the plots in Fig. 4c  are close to the unity, which shows a high degree of similarity with the corresponding distribution in the pristine γ-Ga 2 O 3 . In contrast, for the Ga-implanted structure Pr is only~0.5, which essentially indicates the amorphization. This fact may be readily interconnected with the disturbance in the ionic charge distribution (charge transfer from the excess Ga atoms to the closest O ions) affecting the Coulombic interaction that maintains the order in an ionic system. Moreover, in Fig. 4a we see additional peaks in the green (Ni) and the purple (Ga) heavy-ion PRDF curves at 2.5 Å and 2.7 Å, respectively. These peaks can be correlated with metallic Ni precipitates observed in Fig. 1, while the Ga-Ga bonds may contribute to amorphization of the layer with the highest concentration of Ga ions in the Ga-implanted Ga 2 O 3 in Fig. 2.
In conclusion, the full set of experimental and theoretical data in Figs. 1-4 may be seen as solid evidence for a discovery of the remarkable radiation tolerance in the γ-Ga 2 O 3 /β-Ga 2 O 3 doublepolymorph structures, practically independent of dpa. Meanwhile, the chemical effect introduced by high-fluence Ga ions leads to a nonstoichiometric disordered layer. This observation is rationalized by the unique combination of the specific features of both γ-Ga and O sublattices of γ-Ga 2 O 3 . Intrinsically defective, the γ-Ga sublattice is nearly insensitive to new point defects produced in collision cascades during ion irradiation, while the O sublattice is prone to rapid postcascade recrystallization into original fcc stacking. The collaborative effect of both features explains macroscopically negligible structural deformations observed in heavily irradiated γ-Ga 2 O 3 .

Methods
We used commercial (010) monoclinic beta Ga 2 O 3 polymorph (β-Ga 2 O 3 ) single crystals wafers from Tamura Corporation as initial polymorph substrates. To start with the samples were converted to double Ga 2 O 3 polymorph structures with the implantation parameters as reported in Ref. 14. For that matter we used 58 Ni + , 69 Ga + , 197 Au + , and 20 Ne + ion implantation at room temperature, in particular adjusting implantation energies and fluences to obtain double polymorph Ga 2 O 3 structures of comparative thickness while using different ions. Notably, all implants were performed at 7°off the normal direction of the wafer to minimize channeling. Furthermore, to avoid any heating of the samples during the implantations, the beam current was not exceeded 1 μA/cm 2 for Ni/Ne and 0.1 μA/cm 2 for Ga/Au implants. Table 1 summarizes the implantation parameters used in the experiments. Notably, the maximum of the nuclear energy loss profile (R pd ), the projected range (R p ), as well as the dpa values for each ion, were calculated using the SRIM code 25 simulations (see Supplementary Note 1). Table 1 also shows the ion fluences corresponding to 1 dpa in order to facilitate the fluence/dpa conversion for the readers. Importantly, upon each fluence collection step, the samples were measured by the RBS/C, while selected samples were also characterized with the STEM.
The RBS/C measurements were performed using 1.6 MeV He + ions incident along [010] β-Ga 2 O 3 direction and 165°backscattering geometry. Importantly, it is known from the literature that upon the double polymorph Ga 2 O 3 structure formation, the RBS/C yield exhibits a characteristic trend, attributed to the channeling conditions in the newly formed γ-Ga 2 O 3 polymorph film -see Supplementary Note 2 for more details. This trend, if maintained as a function of the further fluence accumulation, is a fingerprint of the maintained crystallinity. Moreover, the horizontal scale in the RBS/C plotsthe channel numbermeasures the thickness of the newly formed polymorph. Notably, Ga-parts of the RBS/C data were used in the analysis because of the significantly higher sensitivity of this method for heavier Gasublattice compared to the O-sublattice.
Further, STEM was used for detailed crystal structure and chemical analysis. For cross-sectional STEM studies, selected samples were thinned by mechanical polishing and by Ar ion milling in a Gatan PIPS II (Model 695), followed by plasma cleaning (Fishione Model 1020) immediately before loading the samples into a microscope. High Resolution Scanning Transmission Microscopy (HRSTEM) imaging, SAED, energy dispersive x-ray spectroscopy (EDS), and EELS measurements were done at 300 kV in a Cs-corrected Thermo Fisher Scientific Titan G2 60-300 kV microscope, equipped with a Gatan GIF Quantum 965 spectrometer and Super-X EDS detectors. The STEM images were recorded using a probe convergence semi-angle of 23 mrad, a nominal camera length of 60 mm using three different detectors: high-angle annular dark field (HAADF) (collection angles 100-200 mrad), annular dark field (ADF) (collection angles 22-100 mrad) and bright field (BF) (collection angles 0-22 mrad). The structural model of both phases was displayed using VESTA software 26 .
EBSD was performed on the Ne irradiated sample in a Zeiss NVision 40 scanning electron microscope (SEM) equipped with a field emission electron cathode and a Bruker EBSD system with an e-Flash HR+ detector. To ensure the removal of a possible carbon contamination layer, the sample was cleaned for 45 s in an air plasma cleaner. The acceleration voltage was set to 30 kV, the beam current to about 10 nA using a 120 µm aperture. In order to record low noise high quality EBSD patterns, the detector resolution was set to 800 × 570 pixels and the exposure time to 8 × 122 ms per frame. EBSD was done as mappings of 20 × 15 steps, with a step size of 1.9 µm on the irradiated and the unirradiated surface sections of the sample The ab initio molecular dynamics (AIMD) simulations were conducted using the Vienna Ab-initio Simulation Package (VASP) 27 , employing the projected augmented-wave method 28 . The Perdew-Burke-Ernzerhof version of the generalized gradient approximation was used as an exchange-correlation functional 29 . The electronic states were expended in plane-wave basis sets with an energy cutoff of 400 eV throughout all AIMD runs. The Brillouin zones were sampled with a single-Γ k-point for a 1 × 2 × 4 160-atom β-Ga 2 O 3 supercell, and a Γ-centered 2 × 2 × 1 k-mesh for a 1 × 1 × 3 160-atom γ-Ga 2 O 3 supercell. In these simulations, the increase of experimental fluence was mimicked by introducing implanted atoms (Ni, Au, or Ga) in interstitial and substitutional (specified by the superscript S) lattice sites. Specifically, we added 8, 12, and 16 atoms of a given species, which corresponded to 5 at.%, 7.5 at.%, and 10 at.% concentrations with respect to the initial number of atoms in the cell. Initially, the obtained structures were relaxed to the local energy minimum with and without constraining the volume of the cell. Then, the relaxed cells were used in AIMD simulations to enable the dynamic evolution of the system to accommodate the added atoms in the best possible configurations. These simulations were performed for 5 ps with the step of 2 fs in isothermalisobaric ensemble 30 at 900 K and 1 bar, employing Langevin thermoand barostats 31 .
The large-scale classical MD simulations were conducted using LAMMPS package 32 . The newly developed machine-learning interatomic potential of Ga 2 O 3 system was employed 24 . The potential is developed to guarantee the high accuracy for β/κ/α/δ/γ polymorphs and universal generality for disordered structures. In these simulations, Ga FPs were generated cumulatively in 1280-atom β-Ga 2 O 3 and a 1440-atom γ-Ga 2 O 3 cells by iteratively displacing a random Ga atom following a randomly directed vector with the norm of 10~15 Å. The two systems then firstly were relaxed to the local minimum to avoid initial atom overlapping and secondly were thermalized with NPT-MD for 5 ps at 300 K and 0 bar. In total, 600 Ga-FP iterations were run for both cells. In addition, we have performed MD simulations of single cascades in β-Ga 2 O 3 at 300 K. The initial momentum direction and the position of a primary knock-on atom (PKA) were selected randomly at the center of the simulation cell. The PKA was assigned the kinetic energy of 1.5 keV. Periodic boundary conditions were applied in all directions. The temperature was controlled using a Nosé-Hoover thermostat 33 only at the borders of the simulation cell to imitate the heat dissipation in bulk materials. To avoid the cascade overlap with temperature-controlled borders, the number of atoms in the simulation cell was increased to 160 000. We applied the adaptive time step 34 for the efficiency of MD simulations in the active cascade phase. Electronic stopping as a friction term was applied to the atoms with kinetic energies above 10 eV. The simulation time of the single cascades was 50 ps. 120 simulations with different PKA were carried out for statistical analysis.
The structural modifications due to accumulated damage in the studied lattices were analyzed using radial distribution functions. The RDF is defined as the ratio of the ensemble-average local number density of particles, ρ r ð Þ , at a distance r from a reference particle to the average number density of particles in the system.
where N at. is the total number of particles, and V is the system cell volume. Essentially, the RDF is a fingerprint descriptor of the structural property of a system of particles down the atomic scale. For a crystal structure, this function is characterized by well pronounced peaks at the radial distances corresponding to the radii of coordination shells. While the SRO peaks are practically always present in a structure, the LRO peaks in amorphous structures are indistinguishable, since the number density of the atoms in the spherical shells at long distances in an amorphous structure is the same as the average number density in the structure. This feature of RDF gives a good measure of crystallinity in the studied structures. The PRDF describes the type-specified sublattice in a multispecies system 35 . We performed a detailed analysis of the Ga-Ga, Ga-O, and O-O PRDFs as a function of stochastically generated FPs. All RDF distributions were obtained by averaging the signals between 2 to 5 ps for the frames recorded at every simulation step. For clarity, we discriminated the Ga-Ga PRDF features within so-called 1st and 2nd shells. The division was based on the significance of changes that are observed in PRDF distributions before and after the damage accumulation. The border was selected at the valley at~4.0 Å. The Pearson correlation coefficient, Pr, between any two given curves is calculated with the formula: where A i and B i are the variable samples of the two curves, respectively, and A and B are the mean values of the variable samples, respectively.

Data availability
The data that support the findings of this study are available within the paper (and its Supplementary Information file) and from the corresponding authors upon request. Source data are provided with this paper. The machine-learning potential parameter files used to run classical MD simulations are openly available at https://doi.org/10. 6084/m9.figshare.21731426.v1. The corresponding raw data of the classical and ab initio MD published in this paper are openly available at https://doi.org/10.6084/m9.figshare.23599950. Source data are provided with this paper.

Code availability
The code and software used in this work are LAMMPS, VASP, OVITO, and SRIM which are openly available online from the corresponding developers and maintainers.