Multidimensional defect identification of semiconductors in nonequilibrium

Domain: physics
Published: 2025-11-27
Editorial outcome: minor_revision

Abstract
While the static theory for understanding defect properties in semiconductors in equilibrium has been established for decades, it fails to identify the crucial defects in nonequilibrium, such as under irradiation. In this paper, we develop a robust ab initio-driving multiscale modeling framework to identify deep-level defects in irradiated semiconductors with multidimensional defect properties. It overcomes two challenges unsolved in the past studies, that is, unambiguous nonequilibrium defect identification and exact deep-level transient spectroscopy (DLTS) simulation. Our method, verified by identifying the well-known deep-level defects in neutron-irradiated Si, is successfully applied to identify the controversial deep levels in neutron-irradiated wide-bandgap semiconductor, 4H-SiC, contributing to solving the half-century mystery of their atomic origin. Furthermore, we discover that defect origins of the same DLTS peaks vary significantly with annealing temperature, due to different defect types with distinct dynamic behaviors, breaking the long-lasting belief derived from the static defect theory. Our study not only expands the understanding of nonequilibrium defect physics of semiconductors, but also lays a solid foundation for controlling targeted crucial defects to improve material properties and device performances.

Introduction
The needs of deep-space exploration and nuclear-energy application for the future of mankind continue to explosive growth. Semiconductor devices are essential in space and nuclear equipment like spacecrafts and nuclear-power plants, which will be inevitably exposed to the irradiation of energetic particles like electrons, ions, and neutrons1,2. Energetic particles could knock atoms from the lattice and cause displacement damages to form various defects3,4, which would induce extra energy levels within the bandgap of semiconductors. These defect levels could act as carrier trapping, generation or recombination centers, leading to the performance degradation and even permanent failure of devices3. To eliminate the detrimental effects of defects4,5 and improve device radiation resistance, a long-standing goal is to unambiguously identify the crucial nonequilibrium defects in semiconductors. This is not only an important interdisciplinary study between nuclear physics and material science, but is also of great interest to the development of the semiconductor industry.

Defect-identification in semiconductors is always a great challenging issue, plaguing the researchers for more than half a century. While tracking defect formation and evolution is highly challenging in experiments, the existing static defect theories would fail to capture the complex defect properties like structures, charge states and densities in semiconductors, under nonequilibrium conditions like growth, irradiation, ion implantation, or annealing conditions. Thus, without a breakthrough in computational methods, identifying the crucial nonequilibrium defects in semiconductors is like looking for a needle in a haystack. One typical example is the defect identification in 4H-SiC, a crucial wide-bandgap semiconductor in the post-Si age. Among the detected deep-defect levels including EH1, Z1/2, EH3, EH4, EH5, EH6/7, S1 and S2 in 4H-SiC6–8, only Z1/2 and EH7 have been tacitly recognized as VC(=/0) and VC(+/0), respectively9,10 after decades of intensive efforts, leaving other defect levels largely unclear.

Among all the existing defect detection/characterization techniques11,12, the deep-level transient spectroscopy (DLTS)13 is probably the most widely used technique for non-destructive detection of crucial deep levels that affect device performances. It can detect the important properties of defects like electrical nature (donors or acceptors), level positions and carrier capture cross-sections14,15. Unfortunately, it cannot tell us the defect origins of DLTS peaks. Currently, defect types are commonly inferred using available identification methods that rely on single-dimensional defect property parameters like defect levels16–19. However, according to the principle of DLTS, the unambiguous defect identification requires the three-dimensional defect parameters including defect-level position, carrier capture cross-section, and defect concentration, especially for very close levels related to the superimposed DLTS signals. In the past decades, many theoretical studies20–24 have attempted to simulate DLTS but failed, mainly due to the following two challenges: (1) Accurate calculations of multidimensional defect properties, especially under different temperatures, are difficult when using existing methods25–32; (2) It is very hard to accurately solve the complex intra-cascade electrostatic potential24 induced by the carrier capture/emission of defects that plays a key role in renormalizing the Fermi level and thus regulates defect properties33–35.

In this paper, we have successfully solved the two challenges mentioned above, which allows us to innovatively develop a multiscale modeling framework for the unambiguous identification of deep-level defects in semiconductors at different temperatures. First, the modeling framework is validated by the identification of deep-level defects in the collector of a 2N2222 n-p-n type Si bipolar transistor33 (see sample details in Fig. S2a) under pulsed neutron irradiation. Secondly, we have conducted the DLTS experiment and simulation of a 4H-SiC Schottky Barrier diode (SBD) (see sample details in Fig. S2b) after neutron irradiation and annealing, which are consistent with each other. Two overwhelming advantages of our method are presented over previous ones. Finally, by simulating the DLTS of neutron-irradiated 4H-SiC and comparing them with the experimental ones from 4H-SiC Junction Barrier Schottky (JBS) diodes after neutron irradiation and annealing at different temperatures, we unambiguously identified all typical deep levels in 4H-SiC, a key step to further improve device performance via precise defect passivation. Moreover, this work expands the understanding of defect physics of semiconductors far from equilibrium.

Results
Multiscale modeling framework
The developed modeling framework is shown in Fig. 1. First (Fig. 1a), the parameters of basic defect properties, that is, energy levels, carrier capture cross-sections and migration energies, are mainly obtained by density functional theory (DFT) calculations (see DFT calculations in “Methods” section). Secondly (Fig. 1b), based on these defect parameters, long-term (fs-year) defect generation and evolution are simulated in a large system (up to μm) to obtain time/temperature-dependent concentrations of surviving deep-level defects after irradiation and annealing, by coupling the Monte Carlo (MC) and object kinetic MC (OKMC) methods (see Irradiation and annealing simulations in “Methods” section). Here, defect generation under irradiation is simulated with our self-developed open-source MC code, IM3D25, using the energy spectrum of primary knock-on atoms (PKAs) calculated by Geant436 as inputs. For defect evolution, we consider defects and impurities with different charge states and their corresponding dynamic processes. Thirdly (Fig. 1c), based on the obtained multidimensional defect properties, the DLTS of semiconductor devices are obtained by simulating carrier capture and emission of crucial deep-level defects in the sensitive area at different temperatures, using our improved rate theory model based on the Shockey-Read-Hall theory37,38 (see DLTS simulations in “Methods” section). More importantly, based on the annealed defect spatial distributions, we have constructed an effective polarized region model to obtain the improved intra-cascade electrostatic potential without adjustable parameters to correct carrier densities in the rate equations39. Finally (Fig. 1d), defect identification is achieved by using our proposed peak-level matching (PLM) algorithm (see Peak-level matching algorithm in “Methods” section).

Modeling framework verification: neutron-irradiated Si
Si is a benchmark semiconductor with sufficient theoretical and experimental investigations. The typical deep-level defects under irradiation have been clearly determined in previous research. Among them, divacancies (V2) often act as the most important carrier recombination/generation centers, causing the decrease in gain of bipolar transistors and the increase in leakage current of detectors5,40. Here, we present the general process of defect identification and verify our multiscale modeling framework, through simulating the DLTS of neutron-irradiated Si and then matching it with the experimental base-collector of a 2N2222 n-p-n type Si bipolar transistor33 (see Fig. S2a) after pulsed neutron irradiation with the fluence of 1014 cm–2 at 300 K and thermal annealing at 350 K33.

The overall process of identification of deep-level defects in neutron-irradiated Si is shown in Fig. 2. First, staring from the PKA energy spectrum35, we performed MC and OKMC simulations of neutron irradiation and thermal annealing to acquire the evolution of different types of defects (Fig. 2a). Remarkably, the main surviving defects after annealing are defect complexes rather than single-point defects, including VO, I2, V2, I3, V3, OI and CI, with their numbers decreasing in order. The concentrations of these defect complexes in unirradiated Si are negligible, even below the intrinsic carrier density, indicating that the static defect theory is very challenging to catch the key defects in non-equilibrium. The spatial distributions of different types of surviving defects in the steady state are shown in Fig. 2b. Interestingly, opposite to the common assumption of a uniform distribution of defects in the system, the distribution of surviving defects are strongly inhomogeneous. They form separated defect-aggregation regions. The carrier capture of aggregated defects will cause the formation of the electrostatic potential energy barrier to hinder the transport of free carriers, thus resulting in changes in semiconductor electrical properties.

Secondly, according to the numbers of these surviving defects and their energy levels, we obtained the concerned defect levels (Fig. 2c) that may contribute to DLTS signals including CI(−/0), VO(−/0), V2(−/0), V2(=/−), V3(−/0) and V3(=/−). The electron capture and emission for these concerned levels in Processes Ⅱ and Ⅲ (see DLTS simulations in “Methods” section) are simulated, taking the time-dependent charge-state conversion ratios for V2(=/−) as an example (Fig. 2d). The V2(=/−) level is not completely occupied by electrons, due to the existing electrostatic potential energy barrier (see the dashed line in Fig. 2d). The incomplete occupation of V2(=/−) will cause the decrease in initial defect concentration in Process Ⅲ and thus the decrease in intensity of the DLTS signal of V2(=/−), resulting in the stronger DLTS signal of EC − 0.44 eV than that of EC − 0.24 eV33,41,42. It also strongly indicates the important role of the accurate intra-cascade electrostatic potential in describing DLTS, which is largely underestimated in previous research22–24. Finally, according to the experimental rate window (232 s−1)33, t1 and t2 were taken as 2.08×10–4 and 2.08×10–3 s, respectively, as shown by the black dashed line. We thus calculated the DLTS signals of all defect levels (the dashed line in Fig. 2e) according to Eq. (1) in “Methods”. Summing the signals of all levels, we eventually obtained the simulated DLTS of the base-collector of Si bipolar transistor after neutron irradiation and annealing (solid-line in Fig. 2e). Overall, the simulated DLTS is well consistent with the experimental one33 (circles in Fig. 2e). The five experimental peaks can thus be clearly identified as CI(−/0), VO(−/0), V2(=/−), V3(=/−), and the superposition of V2(−/0) and V3(−/0), respectively. It is consistent with the current understanding43,44, indicating that our method can achieve unambiguous identification of crucial deep-level defects in Si without any biased assumptions.

Modeling framework advantages: neutron-irradiated 4H-SiC
SiC is widely recognized as one of the most important wide-bandgap semiconductors for applications in high-temperature, high-frequency, high-power, and extreme irradiation conditions. However, in 4H-SiC, one of the most superior polytypes of SiC, many crucial defect levels have been under debate for decades, preventing the improvement of device performances via defect engineering. Here, the advantages of our modeling framework are given by identifying the crucial deep-levels in a 4H-SiC SBD (Fig. S2b) after 14 MeV neutron irradiation with the fluence of 7.28×1012 cm-2 and annealing at 300 K.

Similarly, based on the C and Si-PKA energy spectra (Fig. S3a), we obtained the evolution behaviors (Fig. 3a) and spatial distributions (Fig. 3b) of different types of defects at a typical time of months in 4H-SiC during neutron-irradiation and annealing at 300 K. The annealed defects mainly include mono-interstitials, di-interstitials, mono-vacancies, di-vacancies, and anti-site defects. According to DFT calculations, the concerned deep levels are shown in Fig. 3c. Then, we simulated electron capture and emission for these concerned levels in Processes Ⅱ and Ⅲ (see DLTS simulations in “Methods” section) of DLTS experiments, respectively. Due to very small electron capture cross-sections (Table S3), CI(−/0), CI(=/−), CC(−/0) and CC(=/−) levels are almost unoccupied in Process Ⅱ, not contributing to any DLTS signals. The capture and emission of other levels in 4H-SiC behave similarly as the levels except for the partially occupied V2(=/−) in Si, thus not presented here. Finally, according to the experimental rate window of 4.1 s-1, we obtained the DLTS signals (the dashed line in Fig. 3d) of all concerned levels in neutron-irradiated 4H-SiC. Summing the output signals of all these levels, we obtained the simulated DLTS of neutron-irradiated 4H-SiC (the black solid line in Fig. 3d). The positions of three simulated peaks slightly deviate from the experimental ones, which is mainly caused by neglecting the temperature-dependence of energy levels and capture cross-sections. By comparing the simulated DLTS with the experimental one, we achieved the identification of deep levels in 4H-SiC, which is one of the key steps to improve device performance via accurate defect passivation in the future. Meanwhile, our method exhibits two advantages over previous ones:

1). Unambiguous identification of typical deep levels. The first advantage is the power of unambiguous identification of typical deep levels. Through DLTS experiments (see Experiments section), we can detect four DLTS peaks at about 200, 300, 370 and 425 K corresponding to EH1/S1, Z1/2, EH3/S2, and EH4 levels6–8, respectively. Previously, through numerous experimental and theoretical studies, Z1/2 was assigned as VC(=/0), EH1/S1 and EH3/S2 were assigned as the different levels of the same defect like VSi or C45–47, and EH4 was assigned as the divacancy or CSiVC (CAV)48,49. Here, we identified Z1/2 as VC(=/0), which is consistent with previous research9. In addition, we newly identified EH1/S1, EH3/S2, and EH4 levels as VCV C(=/−), VCV C(−/0) and VSiSi(=/−), respectively. Due to much smaller electron capture cross-sections (Table S3) of interstitial-type defects (<10-20 cm2) than those of vacancy-type defects (>10-15 cm2), interstitial-type defects could not be detected by DLTS. Moreover, due to the relatively small threshold displacement energy for the C atom, less than half of that for the Si atom50, more C-relative defects like VC, CI and VCV C than Si-relative defects should form in irradiated 4H-SiC, which proves our identification.

2). Prediction of high-temperature DLTS peaks. The second advantage is the prediction ability of peaks that are not detected by DLTS at high temperatures. Here, four DLTS peaks have been detected below 500 K. Using our modeling framework, we can predict one new DLTS peak (green dashed line in Fig. 3d) at about 570 K that is out of the experimental temperature range. This peak is also detected by other DLTS experiments, which corresponds to the EH6 level6,51 with the activation energy of about 1.55 eV. Our simulations showed that EH6 should be contributed by VSiVC(=/0), which is consistent with the current speculation that the EH6 level is associated with a complex involving VC47,52,53, respectively. The unambiguous identification of this level will be presented in the following part.

Multidimensional identification of deep levels in semiconductors
Again, taking 4H-SiC as an example, we further applied our modeling framework to understand deep-level defects in semiconductors in equilibrium and various nonequilibrium conditions with multidimensional parameters. In equilibrium, the dominating deep level in 4H-SiC is commonly considered as Z1/2. With two-dimensional parameters (ΔE, σ) of the defect level and carrier capture cross-section, we can simulate DLTS peaks of 4H-SiC in equilibrium (Fig. 4a). The dominating deep level of VC(=/0) corresponding to VC with the low formation energy of 1.3 eV is the origin of Z1/2 (see Column Ⅰ in Table 1). Assuming the normalized concentrations of different types of defects, we can simulate a group of possible DLTS peaks with (ΔE, σ) in nonequilibrium states (Fig. 4b). Compared to that in equilibrium, much more peaks would appear in DLTS. In particular, some peaks overlap with each other, causing great difficulties in the identification of deep-level defects (see Column Ⅱ in Table 1).

To unambiguously pick out crucial deep levels in nonequilibrium from the large defect-level groups shown in Fig. 4b, a new dimension of defect concentration (C) must be introduced to determine the intensities of different DLTS peaks. Therefore, we further simulated the DLTS of neutron-irradiated 4H-SiC after annealing at 300 K with (ΔE, σ, C), as shown in Fig. 4c. Compared to the peaks in Fig. 4b, some peak signals corresponding to the levels like CAV(−/0) and SiVC(−/0) vanish in the simulated DLTS, due to their very low concentrations. The left peaks with different intensities are contributed by other crucial deep levels with different concentrations. Overall, the simulated DLTS is consistent with the experimental one. By using the PLM algorithm, the deeper EH6 level mentioned above can be unambiguously identified as VSiVC(=/0) (see Column Ⅲ in Table 1).

However, the S1, S2, and EH5 levels are still not determined due to their low defect concentrations and thus weak intensities of DLTS peaks. It is found that the intensities of DLTS peaks corresponding to these levels rely on irradiation and annealing conditions, especially for annealing temperatures (T). To clarify the influence of annealing temperatures on the intensities of DLTS peaks, we further simulated the DLTS of neutron-irradiated 4H-SiC after annealing at 623 K (Fig. 4d). Interestingly, the pattern of DLTS after annealing at 623 K is similar to that at 300 K (Fig. 4c), but the defect origins contributing to these DLTS peaks vary obviously. S1, S2, and EH5 levels appear due to the appearance of CAV and disappearance of VSi, which are contributed by CAV(=/−), CAV(−/0) and CAV(0/+), respectively, while the EH4 level contributed by VSiSi(=/−) disappears. We thus unambiguously identify all typical deep levels in 4H-SiC, including EH1/S1, Z1/2, EH3/S2, EH4, EH5, and EH6. All these deep levels that appear in the DLTS of 4H-SiC devices are thus determined unambiguously, as listed in Column IV of Table 1. Clearly, our multidimensional DLTS modeling framework for unambiguous defect identification is necessary to determine deep levels in complex semiconductors under different growth, irradiation or annealing conditions.

Expanded understanding of defect physics in nonequilibrium
Based on the above discussion, we can eventually expand understandings of defect properties in nonequilibrium in two aspects, dramatically differing from those in equilibrium. First, for defect properties in nonequilibrium, it is generally thought that the dominating crucial deep-level defects should still be the ones with relatively low formation energies in equilibrium. However, our results show that the crucial deep-level defects in 4H-SiC under irradiation (a typical nonequilibrium state) include not only point defects with relatively low formation energies like VC (<5.0 eV) but also vacancy-type clusters with much high formation energies like VCV C, VSiVC and CAV (up to 10 eV). Again, it indicates great limitations of the static defect theory when applied to nonequilibrium states.

Secondly, it is generally thought that the types of crucial deep-level defects remain almost the same at different growth, irradiation, and annealing conditions18,49,54. However, our simulation results show that the types and concentrations of crucial defects in neutron-irradiated 4H-SiC vary significantly at different annealing temperatures. Especially, CAV is almost absent at 300 K, while it appears at 623 K, mainly caused by different defect evolution stages. At 300 K, there are two stages, athermal and thermal, leading to I-V recombination and cluster formation (see Fig. 4e). Whereas there are three stages at 623 K (see Fig. 4f), with an extra formation process of CAV. The details of defect evolution at 300 and 623 K can be found in Fig. S4. Moreover, VSi exists at 300 K but vanishes after annealing at 623 K (see Fig. 4e, f), due to the more violent I-V recombination at the higher temperature. Therefore, annealing at different temperatures will cause the significant difference in defect properties including type, concentration and dynamic behavior.

Methods
DFT calculations
For the DFT calculations of 4H-SiC, the Vienna Ab initio Simulation Package (VASP)56 was employed. The exchange-correlation energy was treated with the screened HSE06 hybrid density functional57,58. The plane-wave cutoff for the wavefunction expansion was set to 400 eV. A sufficient k-points were selected for integrations over the Brillouin zone for the 4×4×1 supercell (128 atoms). During the structural relaxations, a conjugate-gradient algorithm was used until the force on each atom was lower than 0.01 eV Å−1, and the total energy was converged to 1.0×10–5 eV. Moreover, the spin-polarized calculations were considered for all the defects in 4H-SiC. The Climbing Image Nudged Elastic Band method59 was adopted to calculate the migration energies. The total energy and the force on each atom converged to <1.0×10−5 eV and <0.02 eV Å−1, respectively.

We note that given the substantial computational cost associated with defect pairs in SiC, we initially perform structural relaxations using the GGA-PBE functional with a 2×2×2 k-point grid. Subsequently, static calculations of total energy and eigenvalue are carried out using the HSE hybrid functional at the single Γ-point. This approach effectively reduces errors in formation energies and transition levels caused by the bandgap underestimation while keeping the overall computational cost affordable60.

For the electron/hole capture cross-section calculations, the newly developed theoretical calculation method with multi-phonon emissions30,31 and static approximations59 was used. It has been implanted in PWmat60,61, in which force constant matrices, phonon density of states, and electron–phonon coupling coefficients can be calculated to obtain the carrier capture rate constants.

Irradiation and annealing simulations
DLTS experiments generally characterize the properties of the deep-level defects in the sensitive area of semiconductor junctions (such as p-n or Schottky junctions), not in the whole device13,14. For simplification, in the simulation of defect generation and evolution during neutron irradiation and annealing, we thus focus on the sensitive area instead of a complete semiconductor device. The sensitive area is the collector for the base-collector of the n-p-n type Si bipolar transistor and the semiconductor side for the Schottky diodes, as illustrated in DLTS simulations.

Using the energy spectra of primary knock-on atoms (PKAs) as input, the generation of primary point defects under neutron irradiation is simulated with our self-developed open-source MC code, IM3D25, with high accuracy (close to that of Molecular Dynamics) and efficiency (4 orders of magnitude higher than that of SRIM62). The PKA energy spectra are taken from the work of Myers et al.35 for Si and calculated by Geant436 for 4H-SiC, respectively. The PKA energy spectra for 4H-SiC under 14 and 1 MeV neutron irradiation can be found in Fig. S3a and Fig. S3b, respectively. Primary recoil atoms were randomly sampled from the PKA energy spectra, placed in the center of the IM3D simulation box with random directions, and then traced their detailed cascade trajectories. The 3D spatial distributions of primary point defects can thus be acquired through the ‘Full Cascade’ option with the displacement threshold energy of 20 eV for Si63, and that of 20 eV (C) and 35 eV (Si) for 4H-SiC50,64. The spatial distributions of primary point defects are thus provided as an initial database of collision cascades for subsequent modeling of defect evolution.

For defect evolution, the defects in Si include the intrinsic point defects (self-interstitial I and vacancy V), and their clusters (In, Vn and InVm), the substitutional dopant (phosphorus P) and impurities (carbon C and oxygen O), impurity-point-defect pairs (PI, PV, CI, CV, VO and OI), and the complexes of carbon atoms and interstitials (CnmIm). The typical structures and properties of these defects can be found elsewhere35,65–68. The defects in 4H-SiC include the intrinsic point defects (carbon vacancy VC, carbon interstitial CI, silicon vacancy VSi, silicon interstitial SiI, silicon anti-site SiC and carbon anti-site CSi), their typical double-defect clusters (VC CI, VC SiI, VSi SiI, VSi CI, VSi VC, VSi VSi, VC VC, CC, SiC, SiSi, CSi VC, CSi C, CSi Si, SiC VC, SiC C, SiC Si, and SiC SiC) and the dopant (substitutional and interstitial nitrogen, NS and NI). The structures of all these defects can be found in Figure S1. Defect parameters mainly include energy levels, carrier capture cross-sections, migration energies of mobile defects and binding energies of defect clusters. Defect reaction events mainly include I-V recombination, defect-pair formation and break-up, and cluster growth and dissociation. Defect parameters and reaction events for Si and 4H-SiC are listed in Tables S1–S4. Based on these parameters, the evolution of charged defects and carriers is simulated with the open-source OKMC code, MMonCa69, using the spatial distributions of primary point defects as the database of initial collision cascades. The initial collision cascades were placed at the center of the MMonCa simulation box. The dopant and impurities were randomly placed into the simulation box in the substitutional form according to the measured concentrations (3.0×1017 cm–3 for O, 3.0×1016 cm–3 for C and 9.0×1014 cm–3 for P in Si as well as 1.439×1014 cm–3 for N in 4H-SiC). More specifically, a uniform random number between 0 and 1 in each sampling is generated to obtain the position of the dopant and impurities in the simulation box, according to the size of the box. If the position obtained from one sampling coincides with that of the previous positions, resampling is performed until no overlap occurs. Only one cascade was randomly sampled and modeled at each simulation for simplification, because radiation-induced collision cascades nearly do not overlap under the relatively low neutron fluences here. Thousands of cascades are simulated and statistically averaged to reduce errors.

DLTS simulations
According to the principle of DLTS13,14, charge carriers (electrons or holes) can occupy defect energy levels within the bandgap (carrier capture), and can also be emitted from them to the conduction or valence band (carrier emission). When conducting a DLTS experiment, a constant reverse bias voltage and a constant periodic forward pulse voltage are applied to a semiconductor junction (p-n or Schottky junction) with deep-level defects, at different temperatures. Thus, three processes will happen in the semiconductor junction circularly: (Process Ⅰ) The semiconductor junction is in the reverse bias state, and defect energy levels are almost empty. (Process Ⅱ) The forward voltage is applied to the semiconductor junction, and then carriers rapidly occupy defect energy levels. (Process Ⅲ) The forward voltage is released, and then carriers are gradually emitted from defect energy levels and expelled from the junction due to the existing reverse voltage, causing changes in the junction capacitance. When the forward voltage is applied and released circularly, the carrier capture and emission of deep-level defects will also be cycled in the junction. With increasing temperature, deeper defect energy levels will be activated and detected. DLTS experiments measure temperature-dependent differences in junction capacitances between two moments t1 and t2 (t1 < t2) that are determined by the experimental rate window in Process Ⅲ. For a one-sided abrupt p-n junction or Schottky junction under a reverse bias voltage, carrier capture and emission mainly occur in the sensitive area, that is, the space charge region on the side of the semiconductor junction with a high Fermi level. This space charge region is the n-type side for the p-n junction and the semiconductor side for the Schottky junction. Therefore, the output signal of DLTS corresponding to a defect energy level can be generally given by ref.70,

ΔC±(T) = ±Ndef-o(t1, T) ∓ Ndef-o(t2, T) / 2Ndop (1)

Here ΔC±(T) is the difference in capacitance at t1 and t2 after releasing the forward voltage for the minority (+) and majority (–) carrier traps at temperature T, respectively. C0 is the junction capacitance at the steady state after releasing the forward voltage. Ndop is the dopant concentration in the lightly doped side for a p-n junction or in the semiconductor side for a Schottky junction. Ndef-o(t, T) is the concentration of the concerned defect in the sensitive area whose energy levels are occupied by carriers at time t and temperature T in Process III. The only unknown Ndef-o(t, T) can be determined by simulating the carrier capture and emission of different types of deep-level defects at different temperatures in Processes Ⅱ and Ⅲ. The carrier capture and emission can be generally described by the rate equation based on the SRH theory that37,38,

∂N(Xj±1) / ∂t = vth σXj [nN(Xj) - N(Xj±1) Nc exp(-ΔE(Xj±1) / kBT)] (2)

Here, Xj and Xj±1 represent the defect X with the charges of j and j±1, respectively. N represents the defect concentration. νth = (3kBT/m*)1/2 is the thermal velocity of carriers with the Boltzmann constant kB, the temperature T and the effective mass m*. σXj is the carrier capture cross-section of the defect Xj. n is the carrier density. Nc = Nv = 2(2πm*e/h kBT/h2)3/2 is the effective density of states in the conduction/valence band with the electron/hole effective mass m*e/h and the Planck constant h. ΔE(Xj±1) represent the activation energy of carrier emission of the defect Xj±1.

Here, based on Eq. (2), we developed an improved model to simulate carrier capture and emission in Processes Ⅱ and Ⅲ for DLTS simulations. In Process Ⅱ, the rate of carrier capture is larger than that of carrier emission. As mentioned above, the carrier capture (net capture) of deep-level defects within localized cascades will cause the formation of a complex time-dependent surrounding electrostatic potential. A time-dependent electrostatic potential energy barrier Epeb(t) thus form, which can inhibit the carrier capture and cause incomplete occupation of defect energy levels, eventually affecting the intensities of simulated DLTS signals. In previous studies, there are many shortcomings in solving the complex electrostatic potential33,35. Here, we have acquired an improved intra-cascade electrostatic potential without adjustable parameters, via a proposed effective polarized region model based on the more reasonable annealed defect distributions from OKMC simulations39. We have introduced the improved intra-cascade electrostatic potential to correct the carrier density in Eq. (2), that is, n(t) = n0 exp(-Epeb(t)/kBT), where n0 is the carrier density in unirradiated semiconductors. The improved rate equation describing the charge state transition of the defect X from j to j±1 by hole (+)/electron (–) capture at different temperatures is thus given by ref.39,

∂Ncap(Xj±1) / ∂t = vth(T) σXj n0 exp(-Epeb(t)/kBT) Nc(Xj±1) - Ncap(Xj±1)(T) Nc(Xj±1) exp(-ΔE(Xj±1)/kBT) (3)

Here Ncap(Xj) and Ncap(Xj±1) represent the average concentrations of the defect X with j and j±1 charges within cascades, respectively. In Process Ⅲ, carrier capture can be ignored because carriers can almost be expelled from the semiconductor junction, as mentioned above. Thus, Eq. (2) can be simplified to describe carrier emission as,

∂Nemi(Xj±1) / ∂t = -vth(T) σXj Nc(Xj±1)(T) Nemi(Xj±1) exp(-ΔE(Xj±1)/kBT) (4)

where Nemi(Xj±1) represent the concentrations of the defect X within cascades with j±1 charges, which equals Ndef-o in Eq. (1). Generally, the duration of the forward pulse voltage is long enough to reach the equilibrium state of carrier capture. Thus, the initial value of defect concentration in Eq. (4) can be set as Ncap(t→∞) that can be obtained by solving Eq. (3).

Here, based on the improved model for DLTS simulations mentioned above, we have developed a module of C code, IRadMat-DLTS, to simulate DLTS of semiconductors. The lsoda solver71 of the C version was employed for the rate equations. The charge state transitions of the concerned defect energy levels by carrier capture and emission for Si and 4H-SiC were simulated one by one with the temperature interval of 5 K. The simulated DLTS was thus acquired by summing DLTS intensities of all energy levels for Si or 4H-SiC.

Peak-level matching algorithm
In our modeling framework, a peak-level matching algorithm is used to match defect levels with the experimental DLTS peaks for defect identification. The flow of the algorithm is as follows:

1). Input the property parameters (position, carrier capture cross-section and concentration) of all defect levels.

2). Determine the position, Full Width at Half Maximum (FWHM) and intensity of each experimental DLTS peak.

3). Give the confidence intervals for all property parameters to determine the confidence intervals of position, FWHM and intensity for each simulated DLTS peak.

4). Give a search scheme to select possible defect levels matching with one experimental DLTS peak, and a match criterion to determine whether the total DLTS peak summed by the simulated sub-peaks of possible levels is consistent with the experimental one.

5). Search for possible match defect levels for each experimental peak from low to high temperatures according to the search scheme. If only one level matches with one experimental peak, this level is determined to be the only contribution to this peak. Otherwise, an extra most-close defect level is added to the list of possible levels.

6). Calculate the total DLTS peak according to the property parameters of possible levels. If the match criterion is satisfied for this experimental peak, both of these possible levels contribute to this experimental peak. Otherwise, another most-close level is added to the list of possible levels and the new total DLTS peak is recalculated until the match criterion is satisfied. If the match criterion is still not satisfied when all possible levels have been added, adjust the search scheme and match criteria and repeat steps 5–6.

7). Repeat steps 5–6 until the matching of all experimental DLTS peaks is finished.

Based on the peak-level matching algorithm, we have developed a module of C code, IRadMat-PealMat. In this code, the self-developed IRadMat-DLTS (see DLTS simulations section) is used to obtain the simulated DLTS peak according to the property parameters of defect levels.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors report on a multilevel theoretical simulation framework and on its successful application to irradiated silicon and SiC. Their scheme is based on very high level DFT calculations, coupled to sophisticated MC-simulation on one end, and to innovative simulation of DLTS on the other. These coupled tools have led to excellent results, so their importance for any application area involving irradiated semiconductors just cannot be overrated! The novelty of the work is a comprehensive approach to the DLTS spectra, and consideration of the intra-cascade electric potential but, beyond those, the framework contains a high number of innovative algorithmic solutions. Based on these, I believe, the paper deserves publication in a high-visibility journal.\n\nJudging the scientific validity of the individual parts of this paper, which has been authored by experts from various fields, a panel of referees is needed with experience in the respective fields. I only have active, in-depth experience with electronic structure calculation of defects, so – apart from a few general comments – I will restrict my report to the corresponding part of the paper.\n\nApplying HSE06 for Si and SiC is probably the best choice presently. The chosen parameters of the VASP calculation (PW cut-off, convergence criteria) are adequate. The description of the k-point set is missing for Si and is inadequate for SiC. (“Sufficient” is an insufficient description!)The size of the super-cell for the calculations in silicon is not given at all, and the size (128 atoms) in SiC needs discussion in relation to the larger defect aggregates.\n\nConsidering the method of identification by matching theoretical DLTS parameters calculated on the basis of DFT energetics to the experimental spectrum, the accuracy of the charge transition levels is crucial. The authors should provide numerical values in the SI (possibly with literature comparisons). Most of all, however, it should be clarified, what defect structures are being considered, especially since the notation of the complexes are a bit vague. E.g., what kind of a configuration was considered for OI in silicon? Or does CI mean a (C-Si) dumbbell? What does then C mean? The SI should contain figures with the geometrical structures for all defects (in the neutral charge state). Considering aggregate defects, I have a question. These complexes often have metastable configurations, which may play a role in a non-equilibrium sample. Have these been taken into account?\n\nAnother question: what does it mean, “dopant and impurities were randomly placed into the simulation box”? Arbitrarily, or in an SQS?\n\nGeneral comments.\nI consider this paper as an exceptionally high-quality work, still the authors are trying much too hard to sell it. They make statements which are either not justified or at least hard to prove. Below some examples:\n– 4SiC, a crucial wide-bandgap semiconductor in the post-Si age.\nThat age has yet to come – if ever. SiC might completely replace Si in some applications (high-T, -frequency, -power) in the future (!), but Si will remain the basis of micro-electronics. It is very much an open question, whether SiC (or diamond) will be able to successfully compete with Si in quantum information technology as well.\n– it is extremely difficult to understand the formation process of this complex defect in nonequilibrium, resulting in the generation of NV- center in a very low-efficient and uncontrollable way.11-15\nThe formation process of the NV center is well understood by now (see, e.g., https://doi.org/10.1103/PhysRevB.89.075203). It is true though, that controlled production of this center is difficult and the efficiency is somewhat low. However, the same goes for all qubit defects in SiC as well and, as a matter of fact, for any complex defect in any semiconductor.\n– Currently, the available method for defect-identification is nothing but speculation based on single-dimensional defect property parameter like defect level.21-24\nThe authors are unjustly putting down the efforts and successes in defect identification over the last 40 years. This statement is highly offensive and also very untrue!\n– indicating the static defect theory is unable to capture the key defects in non-equilibrium.\nEven admitting that the authors’ work go to a much higher level, this statement is again simply untrue.\n\nSmaller points:\n– The description of the sample and its preparation in Ref.[38] should be shortly summarized in this paper, since it is very much relevant for evaluating the results in silicon.\n– The paper is very well organized (especially considering its complexity), but its English must be improved considerably!!!\nAll in all, I find this paper very important and the part, I am an expert on, scientifically sound. I suggest the authors just omit the aforementioned overzealous (and not really necessary) statements. Provided they also respond to my questions and comments on the electronic structure calculations properly, I recommend publication.", "reviewer_2": "In this manuscript, the authors present an original modelling framework used to identify deep-level defects in semiconductors formed under nonequilibrium conditions, such as those caused by irradiation and annealing. The proposed approach combines DFT calculations of fundamental defect properties, Monte Carlo simulation of defect generation and their evolution following irradiation and annealing, and simulations of DLTS signals. Defect identification is achieved by comparing and matching the calculated DLTS peaks with experimental data. The paper is very well written and addresses important issues in defect physics, which are likely to be of great interest to the scientific community working on the investigation and accurate identification of defects in semiconductors under both equilibrium and nonequilibrium conditions.\nAfter a thorough reading of the paper, I did not find any serious deficiencies; I would like to offer only a few remarks on the text:\n1.The DLTS spectra shown in Fig. 4c and 4d were measured at temperatures exceeding 600K. However, in the Experiments section (page 27), the stated operating temperature range of the DLTS system is 50-500K. Could the authors clarify this discrepancy ?\n2.The authors use numerous abbreviations and symbols to denote various types of native defects or impurities. While some are commonly known and easily understood by readers (e.g., V_C or V_Si for carbon and silicon vacancies, or Si_CV_Si representing a silicon-carbon antisite and silicon vacancy complex), many others – such as V2, I3, OI, CI, VO, CAV, CICI, etc. – are less standard and may be ambiguous or open to misinterpretation. To improve clarity, the authors should clearly define each of these notations in the text.\n3.With regard to the assignment of the dominant Z1/2 level, which is observed under both equilibrium and nonequilibrium conditions: based on the proposed modelling framework, the Z1/2 peak revealed in DLTS experiments is identified as V_C(=/0) level in equilibrium, but V_C(–/0) in nonequilibrium (as shown in Fig. 4 and Table 1). However, on page 11 (line 253), the authors state that this level actually corresponds to V_C(=/–), which was previously assigned to a different peak arising from negative-U behaviour of the V_C (see page 9, line 222). Could the authors clarify this apparent inconsistency ?\n4.Furthermore, previous papers (e.g. Capan et al, JAP 123, 2018) have generally attributed the Z1/2 level to a superposition of two emission peaks originating from the doubly negatively charged V_C defect, exhibiting negative-U properties, as demonstrated by by high resolution Laplace DLTS measurements. Could Authors also comment on this point ?\n5.There are also a few minor issues – mostly typographical errors or inconsistencies – that should be corrected:\n-page 21, line 513, in the phrase “capture section calculations” please add “cross” to clarify it as “capture cross section calculations.”,\n-page 21, line 524: the term “Schottky dioxide” should be corrected to “Schottky diode”,\n-page 26, line 652, the phrase “only one level could math with it” should be corrected to “only one level could match with it.”,\n-page 27, line 666: please change “Schottky Barrier dioxide” to “Schottky Barrier Diode (SBD)”"}, "author_response": "We thank the referees for their careful review on the manuscript and we appreciate very much that both referees hold positive recommendation on our manuscript. We have carefully considered all the comments and made responses as follows. Meanwhile, the changes in the manuscript are marked as blue color.\nResponse for the Referee #1:\nComment 1: The authors report on a multilevel theoretical simulation framework and on its successful application to irradiated silicon and SiC. Their scheme is based on very high level DFT calculations, coupled to sophisticated MC-simulation on one end, and to innovative simulation of DLTS on the other. These coupled tools have led to excellent results, so their importance for any application area involving irradiated semiconductors just cannot be overrated! The novelty of the work is a comprehensive approach to the DLTS spectra, and consideration of the intra-cascade electric potential but, beyond those, the framework contains a high number of innovative algorithmic solutions. Based on these, I believe, the paper deserves publication in a high-visibility journal.\nJudging the scientific validity of the individual parts of this paper, which has been authored by experts from various fields, a panel of referees is needed with experience in the respective fields. I only have active, in-depth experience with electronic structure calculation of defects, so – apart from a few general comments – I will restrict my report to the corresponding part of the paper.\nReply 1: We are grateful for your positive assessment and insightful suggestions, which have helped us improve the manuscript.\nComment 2: Applying HSE06 for Si and SiC is probably the best choice presently. The chosen parameters of the VASP calculation (PW cut-off, convergence criteria) are adequate. The description of the k-point set is missing for Si and is inadequate for SiC. (“Sufficient” is an insufficient description!) The size of the super-cell for the calculations in silicon is not given at all, and the size (128 atoms) in SiC needs discussion in relation to the larger defect aggregates.\nReply 2: We thank you for this important comment.\nThe defect-property parameters in our study:\nFor the case of Si, the defect properties including Fermi level-dependent formation energies, migration energies and carrier capture cross-sections have been extensively studied in the literature. These well-established parameters are not only able to reproduce the experimental results [Liu, J. et al. IEEE Trans. Nucl. Sci. 70, 113-122 (2023); Hehr, B. D. et al. IEEE Trans. Nucl. Sci. 61, 2847-2854 (2014); Myers, S. M., et al. J. Appl. Phys. 104, 044507 (2008)], but also have been verified and adopted in practical simulations [Liu, J. et al. IEEE Trans. Nucl. Sci. 70, 113-122 (2023); Deak P. et al. Nat Commun 14, 361 (2023).]. Therefore, we adopted these well-established parameters for Si. To further clarify this point: (1) In Irradiation and annealing simulations of Methods section, we added “The typical structures and properties of these defects can be found elsewhere35, 68-71.” (2) In Table S1 and Table S2 of Supplementary Note 1, we further added the complete defect-property parameters and defect reaction events of Si, respectively. (3) In Supplementary Note 1 of Supplementary information, we added the sentences of “Defect property parameters mainly include defect energy levels, carrier capture cross-sections, migration barriers of mobile defects and binding energies of defect clusters. These defect parameters are taken from first-principles calculations, molecular dynamics simulations, and experimental measurements, which have been used and verified widely2, 5, 6…. All the property parameters of typical defects and reaction events for Si are listed in Tables S1 and S2, respectively.”.\nDiffering from the case of Si, a few studies have been performed for the calculation of complete defect-property parameters in 4H-SiC, especially for defect clusters. In addition, the deviation in calculated parameters among different studies may cause an inconsistent understanding of defect properties in 4H-SiC. Therefore, we systematically performed DFT calculations for defect-property parameters of 4H-SiC. As described in the DFT calculations of Methods section, we note that “For the energy calculations of 4H-SiC, …The total energy and the force on each atom converged to <1.0×10−5 eV and < 0.02 eV Å−1, respectively.” In particular, we note “We note that given the substantial computational cost associated with defect pairs, we initially perform structural relaxations using the GGA-PBE functional with a 2×2×2 k-point grid. Subsequently, static calculations of total energy and eigenvalue are carried out using the HSE hybrid functional at the single Γ-point. This approach effectively reduces errors in formation energies and transition levels caused by band gap underestimation while keeping the overall computational cost affordable.”.\nThe larger defect aggregates in SiC:\nThe relationship between the supercell size and the larger defect aggregates for 4H-SiC can be explained from the following two aspects. First, the clusters larger than double defects are not considered in our simulations, as described in Irradiation and annealing simulations of Methods section “The defects in 4H-SiC include the intrinsic point defects (carbon vacancy V , carbon interstitial C, silicon vacancy V , silicon interstitial Si, silicon anti-site Si and carbon anti-site C ), their typical double-defect clusters (V C, V Si, V Si, V C, V V , V V , V V , CC, SiC, SiSi, C V , C C, C Si, Si V , Si C, Si Si, and Si C ) and the dopant (substitutional and interstitial nitrogen, N and N).”. In fact, we also tested larger defect clusters like V V V , VV V V and V V V in the simulations, which shows a negligible influence on defect evolution and DLTS signals. For the sake of simplicity, it is reasonable to ignore larger clusters binding more than two intrinsic point defects. Secondly, we compared defect property parameters using the super-cell size of 128 atoms with those of larger super-cell size (512 atoms) from the reference [A. Csóré, et, al, Phys. Rev. B, 96, 085204 (2017)]. In general, the differences within 5% are found in defect property parameters like formation energy and charge state transition energy level of typical defects. Such small differences would have little impact on defect evolution and DLTS simulations, and thus on defect identification, as long as the error caused by the supercell size effect does not cause an intrinsic change in the nature of defects. Considering the calculation cost increasing dramatically with increasing supercell size and numbers of complex defects with diverse configurations and to balance accuracy and efficiency, we selected defect property parameters with sufficient precision using a 128-atom supercell, which is appropriate for the stable double-defect clusters studied in this work.\nBased on this comment, we have carefully revised the Methods section (Pages 19 and 21 in the revised manuscript) by including all the computational details. Meanwhile, we have added Supplementary Note 1 as well as Tables S1 and S2 in the Supplementary Information.\nComment 3: Considering the method of identification by matching theoretical DLTS parameters calculated on the basis of DFT energetics to the experimental spectrum, the accuracy of the charge transition levels is crucial. The authors should provide numerical values in the SI (possibly with literature comparisons). Most of all, however, it should be clarified, what defect structures are being considered, especially since the notation of the complexes are a bit vague. E.g., what kind of a configuration was considered for OI in silicon? Or does CI mean a (C-Si) dumbbell? What does then C mean? The SI should contain figures with the geometrical structures for all defects (in the neutral charge state). Considering aggregate defects, I have a question. These complexes often have metastable configurations, which may play a role in a non-equilibrium sample. Have these been taken into account?\nReply 3: We thank you for this important comment. We apologize for the confusion caused by the notations in our manuscript. We have added defect property parameters including transition energy levels for both Si and 4H-SiC, as in Tables S1 and S3 in the revised Supplementary Information, respectively. We also added the explanations of notations and the structures of all defects in 4H-SiC in Fig. S1 in Supplementary Note 2, and described in Irradiation and annealing simulations in Methods section “The defects in 4H-SiC include the intrinsic point defects (carbon vacancy V , carbon interstitial C, silicon vacancy V , silicon interstitial Si, silicon anti-site Si and carbon anti-site C ), their typical double-defect clusters (V C, V Si, V Si, V C, V V , V V , V V , CC, SiC, SiSi, C V , C C, C Si, Si V , Si C, Si Si, and Si C ) and the dopant (substitutional and interstitial nitrogen, N and N). The structures of all these defects can be found in Fig. S1.”. For Si, we also added defect notations in Irradiation and annealing simulations in Methods and in Supplementary Note 1, as illustrated in Methods section “For defect evolution, the defects in Si include the intrinsic point defects (self-interstitial I and vacancy V), and their clusters (I , V and I V ), the substitutional dopant (phosphorus P) and impurities (carbon C and oxygen O), impurity-point-defect pairs (PI, PV, CI, CV, VO and OI), and the complexes of carbon atoms and interstitials (C I ). The typical structures and properties of these defects can be found elsewhere35, 68-71.” and in Supplementary information “In Si, typical intrinsic and extrinsic defects are considered in the model. Intrinsic defects include self-interstitials (Is), vacancies (Vs) and their clusters (I , V and I V ). For the collector of the 2N2222 n-p-n Si bipolar transistor, we also considered the extrinsic defects including the substitutional phosphorus (P) dopant, carbon (C) and oxygen (O) impurities, defect pairs of PI, PV, CI, CV, VO and OI, and the complexes of C I . These widely proved defects as well as their structures and properties have been extensively studied1-5.”.\nWe agree with you that incorporating different metastable configurations of aggregate/cluster defects is indeed essential for more accurately modeling their evolution behavior. Unfortunately, the evolution of defects in irradiated semiconductors is really a complex multiscale and highly coupled process of multi-micromechanics involving charged mobile defect diffusion and drift, defect–defect reactions, and defect-carrier interactions. Till now, the kinetic Monte Carlo (KMC) method used here has been one of most effective approaches for simulating long-term defect dynamics evolution in large-scale systems after irradiation and annealing. However, to best of our knowledge, there is still no available KMC models that can involve different configurations for aggregate/cluster defects. Even so, we still have considered all their configurations in DFT calculations, and tried to take into account all the key defects and clusters as completely as possible by introducing their most stable configurations into KMC modeling of defect evolution. This approximation should be reasonable enough, as explained in the following three aspects: 1) The formation and evolution of defects most likely to occur in the most stable configuration, which catches the general dynamic rules of defect evolution; 2) According to our DFT calculations, only a small deviation (about 4.2% on average) in the formation energies of cluster defects is observed among different configurations, as shown by the typical examples in Table R1. Such small deviation would have little impact on defect evolution and thus the defect identification; 3) Our simulation results show a good consistency with several experimental measurements, which is also sufficient to indicate the reasonability of this approximation and the reliability of our simulations. We have added these descriptions in the Discussion and outlook Section in the revised manuscript as “Unfortunately, our multiscale modeling cannot involve different configurations for cluster defects, although the KMC method used here has been one of most effective approaches for simulating long-term defect dynamics evolution in large-scale systems after irradiation and annealing. But we still considered all configurations of cluster defects in DFT calculations, and tried to take into account all the key defects and clusters as completely as possible by introducing their most stable configurations into KMC modeling of defect evolution.”\nTable R1. Formation energies of different configurations for three typical cluster defects in SiC.\nConfiguration notations hh-1 hk-1 kk-1\nFormation energies for V Si (eV) 3.94 3.95 3.95\nC I\nFormation energies for C V (eV) 7.13 7.12 7.36\nSi C\nFormation energies for V V (eV) 7.31 7.47 7.12\nSi c\nBased on this comment, we have carefully revised the Method section (Page 21 in the revised manuscript). We have added some of the above discussion (method limitation for metastable configurations) in the Discussion and Outlook section in the revised manuscript. We also added Table S1 and Fig. S1, and revised Table S3 in the revised Supplementary Information.\nComment 4: Another question: what does it mean, “dopant and impurities were randomly placed into the simulation box”? Arbitrarily, or in an SQS?\nReply 4: In our simulations, we randomly placed dopants and impurities into the simulation box to achieve their uniform distribution, neither arbitrarily nor in an SQS. In experiments, impurities are expected to be uniformly distributed, but achieving strict uniformity of dopants at the micrometer scale is indeed challenging. At the nanoscale, however, it is reasonable to approximate the dopant distribution as spatially uniform and random. In this study, we simply focus on the sensitive area of Si and 4H-SiC-based devices with approximately uniform dopants and impurities in nanoscale to simulate their dynamic behaviors. More specifically, a uniform random number between 0 and 1 in each sampling is generated to obtain the position of dopant and impurities in the simulation box, according to the size of the box. If the position obtained from one sampling coincides with that of the previous positions, resampling is performed until no overlap occurs.\nBased on this comment, we have added this explanation in the revised Methods section.\nComment 5: General comments.\nI consider this paper as an exceptionally high-quality work, still the authors are trying much too hard to sell it. They make statements which are either not justified or at least hard to prove. Below some examples:\n– 4H-SiC, a crucial wide-bandgap semiconductor in the post-Si age.\nThat age has yet to come – if ever. SiC might completely replace Si in some applications (high-T, -frequency, -power) in the future (!), but Si will remain the basis of micro-electronics. It is very much an open question, whether SiC (or diamond) will be able to successfully compete with Si in quantum information technology as well.\nReply 5: We thank you for this comment. We have changed the description as “SiC is recognized as one of the most important wide-bandgap semiconductors for the applications in high-temperature, high-frequency, high-power and extreme irradiation conditions.” in the revised manuscript.\nComment 6: – it is extremely difficult to understand the formation process of this complex defect in nonequilibrium, resulting in the generation of NV- center in a very low-efficient and uncontrollable way.11-15\nThe formation process of the NV center is well understood by now (see, e.g., https://doi.org/10.1103/PhysRevB.89.075203). It is true though, that controlled production of this center is difficult and the efficiency is somewhat low. However, the same goes for all qubit defects in SiC as well and, as a matter of fact, for any complex defect in any semiconductor.\nReply 6: We thank you for this comment. We have deleted this overzealous sentence in the revised manuscript. In addition, we have cited this nice PRB paper in the revised manuscript as Ref. 56.\nComment 7: – Currently, the available method for defect-identification is nothing but speculation based on single-dimensional defect property parameter like defect level.21-24\nThe authors are unjustly putting down the efforts and successes in defect identification over the last 40 years. This statement is highly offensive and also very untrue!\nReply 7: We thank you for this comment. We have replaced this unsuitable expression as “Currently, defect types are commonly inferred using available identification methods that rely on single-dimensional defect property parameters like defect levels16-19.”\nComment 8: – indicating the static defect theory is unable to capture the key defects in non-equilibrium.\nEven admitting that the authors’ work go to a much higher level, this statement is again simply untrue.\nReply 8: We thank you for this comment. We revised this unsuitable expression as “indicating the static defect theory is challenging to catch the key defects in non-equilibrium.”\nComment 9: Smaller points:\n– The description of the sample and its preparation in Ref.[38] should be shortly summarized in this paper, since it is very much relevant for evaluating the results in silicon.\nReply 9: We thank you for this helpful suggestion. We added the description of the sample and its preparation in Ref. [38] as “The vertical 2N2222 n-p-n bipolar transistor was made by growing epitaxial silicon on a highly As-doped substrate7. Base and emitter wells were made by diffusing Boron (B) or phosphorus (P) for p-type or n-type layers, respectively. The collector was doped with P to a level of 9×1014 cm-3 with an oxygen impurity level of 3×1017 cm-3 and a carbon level of 3×1016 cm−3.” in the legend of Fig. S2 in revised Supplementary information. We also added “(see sample details in Fig. S2a)” in the sentence “First, the modelling framework is validated by the identification of deep-level defects in the collector of a 2N2222 n-p-n type Si bipolar transistor33 (see sample details in Fig. S2a) under pulse neutron irradiation, a benchmark material system.” in on Page 3 in the revised manuscript.\nBased on this suggestion, we have added the requested description in the legend of Fig. S2 and also mentioned in the main manuscript.\nComment 10: – The paper is very well organized (especially considering its complexity), but its English must be improved considerably!!!\nReply 10: We thank you for this suggestion. We have done our best to improve our English expression.\nComment 11: All in all, I find this paper very important and the part, I am an expert on, scientifically sound. I suggest the authors just omit the aforementioned overzealous (and not really necessary) statements. Provided they also respond to my questions and comments on the electronic structure calculations properly, I recommend publication.\nReply 11: Again, we appreciate very much for your very high recommendation of our manuscript. We have done our best to revise the manuscript based on your comments and suggestion.\nResponse for the Referee #2:\nComment 1: In this manuscript, the authors present an original modelling framework used to identify deep-level defects in semiconductors formed under nonequilibrium conditions, such as those caused by irradiation and annealing. The proposed approach combines DFT calculations of fundamental defect properties, Monte Carlo simulation of defect generation and their evolution following irradiation and annealing, and simulations of DLTS signals. Defect identification is achieved by comparing and matching the calculated DLTS peaks with experimental data. The paper is very well written and addresses important issues in defect physics, which are likely to be of great interest to the scientific community working on the investigation and accurate identification of defects in semiconductors under both equilibrium and nonequilibrium conditions. After a thorough reading of the paper, I did not find any serious deficiencies; I would like to offer only a few remarks on the text:\nReply 1: We are grateful for your highly positive assessment and insightful suggestions, which have helped us improve the manuscript. We have done our best to revise the manuscript according to your comments and suggestions, as follows.\nComment 2: 1. The DLTS spectra shown in Fig. 4c and 4d were measured at temperatures exceeding 600K. However, in the Experiments section (page 27), the stated operating temperature range of the DLTS system is 50-500K. Could the authors clarify this discrepancy?\nReply 2: We apologize for this confusion. We cannot measure the DLTS over 500 K in our experiments. Since we want to further check our calculations at higher temperature, we adopted the DLTS at 623 K using other experimental results, that is, the DLTS spectra shown in Fig. 4 are taken from the literature.\nBased on this comment, we have added description in the legend of Fig. 4 that, “Here, the experimental DLTS of the 4H-SiC JBS diode is taken from Ref. 7” on page 13 in the revised manuscript.\nComment 3: 2. The authors use numerous abbreviations and symbols to denote various types of native defects or impurities. While some are commonly known and easily understood by readers (e.g., V_C or V_Si for carbon and silicon vacancies, or Si_CV_Si representing a silicon-carbon antisite and silicon vacancy complex), many others – such as V2, I3, OI, CI, VO, CAV, CICI, etc. – are less standard and may be ambiguous or open to misinterpretation. To improve clarity, the authors should clearly define each of these notations in the text.\nReply 3: We thank you for this helpful suggestion. We have added the definitions/explanations for notations of all defects included in our simulations for both Si and 4H-SiC in Irradiation and annealing simulations in Methods and in Supplementary Note 1 and Note 2 in the revised manuscript. The detailed revisions can also be found in our responses to Comments 2 and 3 from the Referee #1.\nComment 4: 3. With regard to the assignment of the dominant Z1/2 level, which is observed under both equilibrium and nonequilibrium conditions: based on the proposed modelling framework, the Z1/2 peak revealed in DLTS experiments is identified as V_C(=/0) level in equilibrium, but V_C(–/0) in nonequilibrium (as shown in Fig. 4 and Table 1). However, on page 11 (line 253), the authors state that this level actually corresponds to V_C(=/–), which was previously assigned to a different peak arising from negative-U behaviour of the V_C (see page 9, line 222). Could the authors clarify this apparent inconsistency?\nReply 4: Thank you for your very careful reading of our manuscript and we apologize for the mistakes in the original version.\nFirst of all, V (−/0) was incorrectly written as V (=/0) in Table 1, Column Ⅰ, C C\nand as V (=/-) on page 11 in the original manuscript. We carefully checked it again C\nand found that our previous identification of Z should be incorrect. In fact, the Z 1/2 1/2\npeak is indeed contributed by the V (=/0) level rather than V (-/0), as expected in the C C\nexperiments. According to the suggestion of Referee #1, we have recalculated the Fermi-dependent formation energies of V using a larger supercell of 576 atoms. We C\nfound that only the V (=/0) level appears when using such a very large supercell, in C\ncontrast to the co-existence of V (-/0) and V (=/-) levels obtained with a small C C\nsupercell of 128 atoms in our original manuscript. This indicates that, in the case of few special defects like V with a negative-U property, the formation energy error C\ninduced by the supercell size effect could change the nature of defects. We have revised the identification result for the Z peak as the V (=/0) level rather than V (- 1/2 C C\n/0), as given in Figs. 3c and 3d, Figs. 4a-d and Table 1 in the revised manuscript (also copied here). Now the simulated DLTS of 4H-SiC devices are more consistent with different experimental results. We also added the corresponding description of “Here, we identified Z as V (=/0), which is consistent with previous research” on 1/2 C\npage 8 and revised the sentence “In addition, we newly identified EH1/S1, EH3/S2, and EH4 levels as V V (=/-), V V (-/0) and V (=/-), respectively.” on pages 8-9 in C C C C Si\nthe revised manuscript. Meanwhile, we have removed the second, unrealistic advantage over previous methods that. We also revised the sentence of “The dominating deep level of V (=/0) corresponding to V with the low formation energy C C\nof 1.3 eV is the origin of Z (see Column Ⅰ in Table 1).” on page 11 in the revised 1/2\nmanuscript.\nRevised Figs. 3c and 3d in the manuscript │DLTS simulations of the neutron-irradiated 4H-SiC SBD. c, Energy levels and numbers per cascade of the concerned deep-level defects.\nd, Simulated and experimental DLTS of the 4H-SiC SBD after neutron irradiation and annealing at 300 K, acquired by our modelling framework and experimental measurements, respectively.\nRevised Figs. 4a-4d in the manuscript│ New understanding of deep-level defects in 4H-SiC through defect evolution under different conditions and DLTS simulations with different dimensions of defect parameters. a,b, DLTS of 4H-SiC in equilibrium (a) and nonequilibrium (b) with two-dimensional parameters of (E,). c,d, DLTS of a 4H-SiC JBS diode after 1 MeV neutron irradiation with the fluence of 6.6×1013 cm-2, and annealing at 300 K (c) with three-dimensional parameters of (E,,C) at 623 K (d) with the extra dimensional parameter of T. Here, the experimental DLTS of the 4H-SiC JBS diode is taken from Ref. 7.\nRevised Table 1 in the manuscript│Identification of typical deep levels in 4H-SiC in equilibrium and nonequilibrium with different dimensional parameters. I corresponds to 4H-SiC in equilibrium with (E,). Ⅱ corresponds to 4H-SiC in nonequilibrium with (E,). Ⅲ corresponds to irradiated 4H-SiC with (E,,C) at the room temperature. Ⅳ corresponds to irradiated 4H-SiC with (E,,C) at the high temperature.\nⅠ Ⅱ Ⅲ Ⅳ\nV V (=/−)/CAV(=/-)/\nEH1/S1 Si C V C (=/-)/Si V (−/0) V C V C (=/−) V C V C (=/−)/CAV(=/-)\nC Si C Si\nZ V (=/0) V (=/0) V (=/0) V (=/0)\n1/2 C C C C\nEH3/S2 V C V C (−/0)/V Si (=/-)/CAV(-/0) V C V C (−/0) V C V C (−/0)/CAV(-/0)\nV V (−/0)/V (=/-)/\nC C Si\nEH4 V (=/−)\nCAV(-/0)/CAV(0/+) Si\nEH5 V Si (=/-)/CAV(-/0)/CAV(0/+) CAV(0/+)\nEH6/7 V V (=/0) V V (=/0) V V (=/0)\nSi C Si C Si C\nBased on this comment, we have corrected the mistakes and updated the corresponding Figures 3, 4 and Table I in the revised manuscript.\nComment 5: 4. Furthermore, previous papers (e.g. Capan et al, JAP 123, 2018) have generally attributed the Z1/2 level to a superposition of two emission peaks originating from the doubly negatively charged V_C defect, exhibiting negative-U properties, as demonstrated by high resolution Laplace DLTS measurements. Could Authors also comment on this point?\nReply 5: Thank you for this comment. As mentioned in the above Reply 4, a mistake caused by the size effect of the supercell has been made in the original manuscript. The reasonable interpretation is that only the V (=/0) level appears when using the C\nlarger supercell of 576 atoms in DFT calculations of defect formation energies, rather than the co-existence of V (-/0) and V (=/-) levels in a 128-atom supercell. Again, C C\nwe apologize for making this mistake. We have revised all corresponding expressions and results in the revised manuscript.\nComment 6: 5. There are also a few minor issues – mostly typographical errors or inconsistencies – that should be corrected:\n-page 21, line 513, in the phrase “capture section calculations” please add “cross” to clarify it as “capture cross section calculations.”,\nReply 6: We thank the referee for carefully pointing out these mistakes. We have added “cross-” to in the phrase “capture cross-section calculations.” in the revised manuscript.\nComment 7: -page 21, line 524: the term “Schottky dioxide” should be corrected to “Schottky diode”,\nReply 7: We have corrected the term “Schottky dioxide” to “Schottky diode” in the revised manuscript.\nComment 8: -page 26, line 652, the phrase “only one level could math with it” should be corrected to “only one level could match with it.”,\nReply 8: We have corrected the sentence of “If only one level matches with one experimental peak” in the revised manuscript.\nComment 9: -page 27, line 666: please change “Schottky Barrier dioxide” to “Schottky Barrier Diode (SBD)”\nReply 9: We have corrected the term “Schottky Barrier dioxide” to “Schottky Barrier diode (SBD)” in the revised manuscript.\nList of major changes\nFollowing the referee’s suggestions, improvements in the revised manuscript are listed below:\n1. We delete the sentence “Another famous example is the generation of NV- center in diamond. The NV- center, created by ion implantation, holds great promise for advanced quantum qubit applications 11-13. However, it is extremely difficult to understand the formation process of this complex defect in nonequilibrium, resulting in the generation of NV- center in a very low-efficient and uncontrollable way 14, 15.” in lines 56-60 on page 2 in the original manuscript and corresponding references.\n2. We deleted the paragraph “2) Level resolution of metastable states. The second advantage is overcoming the low-time resolution in DLTS measurements to present more refined defect levels. We predicted a new peak in our simulated DLTS contributed by V C (=/-) without any experimental correspondence, caused by the negative-U properties of V 9. V can capture two electrons, and the C C\nbinding energy of the second electron is higher than that of the first one. Hence, during the reverse process (electron emission), the emission of the first electron from V ( C V 2C − − e − = V −C ) will cause the instantaneous emission of the second one ( V C −−e− =V C 0 )55. Present DLTS techniques with limited time resolution cannot detect the ultrafast second emission process. Thus, two levels of V (-/0) and C V C (=/-) were considered as V C (=/0) in DLTS experiments.” in lines 219-227 on page 9 in the original manuscript.\n3. We revised the identification results of Z1/2 as “V (=/0)”, and accordingly C\nrevised Figs. 3c and 3d, Figs. 4a-d and Table 1. We also revised corresponding wrong expression about the identification results of Z1/2.\n4. We added the sentences “Given the substantial computational cost associated with defect pairs, we initially perform structural relaxations using the GGA-PBE functional with a 2×2×2 k-point grid. Subsequently, static calculations of total energy and eigenvalue are carried out using the HSE hybrid functional at the single Γ-point. This approach effectively reduces errors in formation energies and transition levels caused by band gap underestimation while keeping the overall computational cost relatively low.” in lines 479-484 on page 19.\n5. We added the description about property parameters, structures and notations of defects in Si and 4H-SiC on page 21 in Section Irradiation and annealing simulations in Methods\n6. We added the description of defect-property parameters and reaction events in Si in Supplementary Note 1, and also added physical parameters of typical defects in Si in Table S1 and reaction events among defects in Si in Table S2 in Supplementary Note 1. We also revised the descriptions of defects property parameters and reaction events in 4H-SiC in Supplementary Note 2, added the notations and geometrical structures of defects in 4H-SiC in our simulations in Fig. S1, and revised defect property parameters in 4H-SiC in Table S3.\n7. We added the device sample details in the legend of Fig. S2 in lines 91-102 on page 9.\n8. We added 7 new references as Refs. 56, 61, 68-71 and 74 in the revised manuscript, and eight references as Refs. 1-8 in the revised Supplementary information.\n9. We revised some not necessary/proper expressions mentioned by Referee 1.\n10. We carefully corrected the English grammar and improved the overall English expression.\nThese major changes and other updates are marked in the files."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have properly responded to all of my comments and I find the revised MS publishable as it is.", "reviewer_2": "After a careful reading of the revised manuscript and the authors' honest responses to all my comments and remarks, I consider the article ready for further editorial work and recommend it for publication."}, "author_response": ""}]

Decision letter