Trap-state mapping to model GaN transistors dynamic performance

Trapping phenomena degrade the dynamic performance of wide-bandgap transistors. However, the identification of the related traps is challenging, especially in presence of non-ideal defects. In this paper, we propose a novel methodology (trap-state mapping) to extract trap parameters, based on the mathematical study of stretched exponential recovery kinetics. To demonstrate the effectiveness of the approach, we use it to identify the properties of traps in AlGaN/GaN transistors, submitted to hot-electron stress. After describing the mathematical framework, we demonstrate that the proposed methodology can univocally describe the properties of the distribution of trap states. In addition, to prove the validity and the usefulness of the model, the trap properties extracted mathematically are used as input for TCAD simulations. The results obtained by TCAD closely match the experimental transient curves, thus confirming the accuracy of the trap-state mapping procedure. This methodology can be adopted also on other technologies, thus constituting a universal approach for the analysis of multiexponential trapping kinetics.

Charge-trapping is a critical factor in determining the dynamic performance of electronic devices (transistors and diodes). This is relevant for conventional silicon transistors, but more critical for wide bandgap semiconductor devices, which can operate at much higher frequencies and electric fields. As an example, GaN transistors used in power conversion can reach electric fields above 2 MV/cm 1 , and operate at frequencies above 1-10 MHz and very low duty cycles 2 , with dV/dt in the range of 10-100 V/ns 3 .
Such dynamic performance can be exploited only through a minimization of the trapping phenomena, which originate from defects located in the semiconductor material, at the heterointerfaces, or at the surface of the devices.
For an accurate TCAD (Technology computer-aided design) modeling of device characteristics, a reliable technique for the characterization of trap-states must be developed.
Such characterization is normally carried out through current transient measurements 4 , and trap parameters are typically extracted under the oversimplistic approximation that defects are ideal (thus resulting in singleexponential current transients) or have a Gaussian distribution of time constants 5-10 . However, this is an empirical approach, that does not consider the physical origin of traps: (i) surface or interface states may have a broad distribution of energies and cross sections, thus resulting in turn-on transients that substantially deviate from ideal exponentials, and (ii) a Gaussian distribution relies in the (always false) approximation that the single component kinetic is a step-function instead of an exponential decay.
The question is therefore, is it possible to define a methodology to quantitatively extract trap parameters, in presence of non-ideal (non-exponential) de-trapping processes? In this article, we address this issue by proposing a novel methodology (Trap-State Mapping), that provides a quantitative description of the activation energy and cross section distribution of traps.
Before describing our approach, to give a bit of background, we remind that if carriers are trapped at discrete energy levels in the bandgap, an almost ideal exponential recovery is observed, which is described by the simple exponential decay function 11 : where N t is the density of filled traps after a certain electrical stress, n(t) is the density of trapped electrons, and τ 0 is the time constant of the detrapping process. This equation comes from the solution of the linear first order differential equation: On the other hand, especially for non-ideal interface and surface traps, the pure exponential relaxation rarely forms in nature. In such cases the relaxation kinetics show stretched-exponential trends 12 , that is often explained by considering a distribution of activation energies and cross sections for the individual traps involved [13][14][15] .
Considering that (i) surface and interfaces play a key-role in determining the dynamic performance of wide bandgap HEMT devices and, (ii) de-trapping from surface states is characterized by a distribution of activation energies and capture cross sections, which results in a heavily-spread time-constant spectrum, the mathematical description of how these competing mechanisms interact with each other during the recovery assumes great importance.
A commonly used approach for extracting interface/surface trap parameters is based on capacitance 14,16 or conductance 17 measurements on capacitors, to extract the interface state density (D it ). However, this requires specific test structures, whose properties and processing parameters may deviate from those of final devices.
With this paper, we want to overcome this issue, by proposing an effective and reliable methodology for extracting the interface trap properties (activation energy and capture cross section distribution) starting from the analysis of the stretched exponential transients measured on real transistors.
We start by developing the mathematical framework for the analysis of stretched exponential decay functions, and we show that the trap parameters can be univocally extracted through multi-exponential analysis. Specifically, we consider that the density of occupied traps n(t) during relaxation is described as: where τ 1 , τ 2 , . . . τ k is a finite set of discrete time constants, N T = N 1 + N 2 + · · · + N k is the total density of filled traps and, A i = N i /N T is the normalized amplitude of each exponential decay function associated to the relative time constant τ i . Then, we describe a methodology for the identification of the trap parameters, univocally described by the values of A i and τ i . This procedure is made difficult by the fact that we are dealing with a series of non-linear and non-orthogonal functions.
The developed mathematical framework is then used to interpret the relaxation kinetics of AlGaN/GaN HEMTs submitted to hot electron stress. On top of this, for the first time, we demonstrate how, from the extracted time constant distributions at different temperatures, the trap states properties can be univocally described. This is achieved by adapting the Arrhenius law in case of distributed trap levels for the effective extraction of the activation energy E a,i and capture cross section and σ c,i distributions.
Finally, the extracted trap parameters are used as input for TCAD simulations: we show that the results obtained by TCAD closely reproduce the experimental ones, thus demonstrating the effectiveness of the presented approach in quantitatively mapping the properties of trap states.
The results described here are referred to AlGaN/GaN HEMTs; however, the methodology developed here offers a universal approach to extrapolate the time-constant spectrum profile in heavily stretched exponential decays, that can be used also for investigating different technologies and different physical processes.

Theory of trap-state mapping via stretched-exponential functions
A quantitative trap-state mapping requires the extraction of the distribution of activation energies and cross sections for a set of states, which result in stretched-exponential transient for a transistor. In the past decades, many mathematical methods have been developed in order to study the nature of stretched exponential functions [18][19][20][21][22][23][24] . This is because of its exceptionally wide range of application throughout many fields of science 18 .
These methods rely upon the approach formulated in 1959 by Gardner et al. 21 , which is based on the inversion of the Laplace integral equation by a Fourier transforms method. A main concern in this approach arises from the fact that the experimental data is often noisy and truncated, so an accurate extraction of the parameters may be prevented.
In this section, we propose two approaches to identify the normalized amplitude of the individual exponential components ( A i ) and the related time constants ( τ i ) that constitute a stretched exponential transient. The methodologies can be used whenever the quality of the experimental data permits a fit with the multiexponential decay function f (t).
Asymptotic approach. We start considering that the function f (t) in Eq. (3) is in the form of a Dirichlet series, which may be generalized in a continuum of amplitudes A i = g(ν)dν as: where ν i = τ −1 i and g(ν) is the probability density function such that: Therefore, the extraction of the g(ν) profile is a necessary step to find the time-constant distribution of the experimentally determined data. www.nature.com/scientificreports/ A common methodology, when studying the relaxation kinetics of multiexponential functions, consists of fitting the experimental data with the semi-empirical stretched exponential function 22 , and then posing the equality: where the stretching exponent β is a semi-empirical parameter in the range (0,1]. In this specific case, the g(ν) profile can be computed by expanding the stretched exponential function in Taylor series 22,23 as: where Ŵ is the gamma function 24 . This function does not have a closed solution. However, for the special values β = 1/2 and β = 1 , Eq. (7) can be written as: where δ is the Dirac function centered in τ −1 0 . Recently, Johnston calculated the analytical function for β = 1/3 and β = 2/3 25 . However, for a general value of β the function g(ν) can be obtained in asymptotic form using the saddle point method of integration in Eq. (6) as follows 23 : It can be noted that, for β = 1/2 , Eq. (10) coincides with Eq. (9). Finally, Fig. 1a shows the probability density function g(v) vs ν = τ −1 , illustrated for some values of β with τ 0 = 1 s. If β → 1 , then g(v) is described by the Dirac function centered in τ 0 . Furthermore, by decreasing β , the peak τ pk of the probability density function changes and becomes different from τ 0 . The dependency between the peak τ pk and β can be calculated from the derivative of Eq. (10) as: Figure 1b shows the dependency between τ pk /τ 0 and β . On one hand, for high values of β (ideal exponential decay), τ pk and τ 0 coincides. On the other hand, for low values of β (stretched exponential) the ratio between the two increases significantly leading to counterintuitive results. For example, if β = 0.3 , τ pk is two order of magnitude larger than τ 0 meaning that the extracted τ 0 and the peak of the time-constant spectrum do not coincide.
To conclude, the analytical approach presents some main advantages such as a closed analytical function (easy to implement) and a good accuracy. However, an asymptotic solution is not the most accurate solution (see also next section) and this approach only works with a single stretched exponential centered in τ 0 . Therefore, in the next section, the inverse Laplace transform approach is proposed. Inverse Laplace transform approach. A second approach to extract the probability distribution function, consists in recognizing that the integral of Eq. (4) has Laplacian form. Therefore, the stretched exponential function f (t) is the Laplace transform of g(ν).
Using an analytical expression, such as the one in Eq. (6) to define f (t) , with β and τ 0 taken from the experimental data, it is possible to extract g(ν) by means of inverse Laplace Transform, solving the following equation: . Alternatively, the inversion may be computed numerically. In Fig. 1c the numerical solution (blue) and the asymptotic solution (red) are compared. The result shows that the accuracy obtained using the saddle point method is quite good, validating Eq. (11) for the τ pk extraction. However, as can be seen in Fig. 1d,e, by integrating g(ν) in Eq. (6) the error introduced by the asymptotic approach is not negligible. Finally, the amplitude A i of each component can be extracted as A i = g(ν)dν. Figure 2 shows the (a) contour plot and (b) 3D map of the A i distribution for different values of β with τ 0 = 1 s. A high β (ideal exponential decay) tends to the Dirac function. On the contrary, a low β leads to a very broad distribution. www.nature.com/scientificreports/ Since ν = τ −1 , the A i distribution can be plotted together with the stretched exponential function by simply inverting the x-axis. The superimposition between the A i distribution and its relative stretched exponential is shown in Fig. 2c where, f (t) is defined as the sum of three stretched exponential functions. Each A i is plotted as a delta function centered in its respective time-constant τ i .
To verify the result, the extracted A i distribution is used in Eq. (3) which physically describes the relaxation to the ground state occurring through many competing channels. Figure 2d shows that the summation of each exponential decay multiplied for the extracted A i returns the correct f (t) even for a multiple stretched exponential decay, thus validating the extraction procedure.

Experimental trap-state mapping on GaN-based transistors
The mathematical framework presented in the previous section helps to understand the physical nature of the stretched exponential, which may be useful for a broad range of applications. www.nature.com/scientificreports/ Based on these considerations, we demonstrate the applicability of the developed mathematical framework to the trap-state mapping in in normally-off p-GaN HEMTs after semi-on stress. During this stress condition, the device must sustain a simultaneously high drain-to-gate electric field and carrier density 26 . The carriers accelerated by the high field are described as "hot", and their excess kinetic energy may favor trapping effects that are not reachable under normal (off-state) conditions.
In our previous works we demonstrated that the hot-electrons injection mainly involves surface trapping 27 . To summarize our results, we observed that the trapping kinetics follows a logarithmic trend, thus reconducting it to an inhibition effect originated from the Coulombic repulsion between a filled trap and a free electron injected to the surface 28 .
In GaN HEMTs, the charge trapped at the surface depletes the 2DEG channel and can be quantified by monitoring the dynamic-R on increase. Therefore, the measurements were done in the linear region since, for a given voltage, current is inversely proportional to the dynamic-R on . From previous investigations, we observed that after 10 µs of semi-on stress the very fast hot electron trapping already induced a significant current collapse on the device while the off-state contribution is negligible. Consequently, the test is performed by monitoring, through fast I D V G, the DCT of the device under test (DUT) biased in the linear region at V ON = (V G , V D ) = (4, 0.4) V, after a filling time t fill = 10 µs at V fill = (2, 50) V, to inject the hot electrons to the surface. Figure 3a shows the I D V G taken at different moments during the recovery phase at room temperature. The red curve represents the transfer characteristics taken immediately after the stress. From a qualitative perspective, the 10 µs hot-electron injection in the surface induces a considerable current collapse in the device (from 80 mA/mm to 50 mA/mm at V G = 4 V).
Furthermore, despite the short filling time, the recovery process is very slow. To obtain a quantitative evaluation of the relaxation kinetics, Fig. 3b  The detrapping kinetics are fitted with Eq. (6) and show a strongly stretched-exponential behavior. The set of parameters A T , β T and τ T extracted at each temperature is reported in Table 1. www.nature.com/scientificreports/ Figure 4a,b shows the contour map of the time constant distribution numerically extracted by solving the Inverse Laplace Transform function from the experimental data and mapped at different temperatures. Here, the g(ν) profile extraction is obtained by Euler numerical inversion 29 from the fitted function f (t) . Alternatively, for a direct extraction from the measured I DS (t) , other numerical strategies require the solution of a first kind Fredholm integral equation where the kernel is the decaying exponential function 30,31 . The solution of this illposed inverse problem is extremely sensitive to noise; however, substantial improvements can be obtained by the application of Tikhonov regularization 32 -34 .
The increase in temperature affects the profile by changing both the peak position and the width of the curve. Considering that the physical nature of the stretched behavior can be attributed to the intrinsic defectiveness of the surface region, characterized by a distribution of activation energies ( E a ) and capture cross sections ( σ c ), then the mathematical description developed in the previous section can give an extraordinary insight on the time constant distribution. This is because, in semiconductor physics, time constants are strictly linked to the activation energy and capture cross section of a trap level as 35 : where T is the temperature, h is the Planck constant, k is the Boltzmann constant, m 0 m r is the effective mass of electrons, and g is a degeneration factor (typically 1). Given a temperature dependent time constants distribution, such as the one mapped in Fig. 5b, each E a and σ c can be uniquely extracted as: The results of this analysis, applied to GaN HEMTs, is shown in Fig. 4c,d. The two distributions show the values of activation energy and capture cross section for each of the traps responsible for the stretched exponential process. Remarkably, the extracted E a distribution appears to be consistent with the experimental results reported by Stradiotto et al. 14 where, SiN/AlGaN/GaN structure is used.

TCAD simulation
To demonstrate the correctness of the of trap-state mapping procedure, we used the obtained trap-parameters as input for TCAD simulations and compared the simulated recovery traces with the experimental ones. Specifically, the recovery transients were simulated in Sentaurus TCAD by including the obtained distribution of traps (with the related E a,i and σ c,i ) at the passivation/AlGaN interface. Details on the device structure and the employed TCAD model can be found in 36 . Traps are defined according to the activation energy and cross section distributions obtained via trap-state mapping, and are uniformly distributed at the surface as acceptors (see Fig. 5a). In addition, surface donors (0.8 eV from CB), are equally distributed in the access region 37 . The simulation, solving for standard drift-diffusion equations, consists in forcing all acceptors to be occupied before starting the recovery ( t = 0 ), when their occupation state is then unfrozen. The simulation then considers that in t = 0 all acceptors are filled with electrons, and mimics the recovery transient when 0 V is applied at the contacts, thus accounting exclusively for thermal-induced recovery.
Snapshots of the device are saved during the recovery phase and are fed as input to a separate script where the device is biased in the linear region with V G = 4 V, V D = 0.4 V, in order to emulate the experimental recovery traces.
The simulated recovery traces, reported in Fig. 5b, are normalized to the nominal current and display excellent matching with the experimental ones. This demonstrates that the trap distribution obtained by trap-state mapping can effectively reproduce the dynamic behavior of the devices, even in presence of non-exponential recovery traces.
The interface traps density is used as fitting parameter in order to match the initial current degradation after 10 µs. As observable in Fig. 5c,d, the electron emission from SA during recovery leads to a decrease of the electric field in the AlGaN barrier and the repopulation of the 2DEG density in the channel.

Conclusion
We have proposed a general methodology for mapping the properties (activation energy, cross sections) of a distribution of surface/interface states in GaN-based electronic devices. The properties of the traps are extracted through a mathematical approach, based on inverse Laplace transform. The developed methodology is used to univocally extract the properties of distributions of traps filled during hot-electron stress in GaN-based transistors.
To prove the validity and usefulness of the model, the extracted map distributions are used as input for TCAD simulations. The results obtained by TCAD closely match the experimental transient curves, thus confirming the effectiveness of the developed technique. The proposed methodology is generic and can be extended to other semiconductors for high-speed/high-voltage electronics, including AlN, gallium oxide, and diamond, thus representing a powerful tool for quantitatively assessing and modeling the properties of trap states.  www.nature.com/scientificreports/

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. www.nature.com/scientificreports/ Reprints and permissions information is available at www.nature.com/reprints.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.