TCAD Modeling of High-Field Electron Transport in Bulk Wurtzite GaN: The Full-Band SHE-BTE

Gallium Nitride (GaN) High-Electron Mobility Transistors (HEMTs) actually represent one of the best candidates for medium-high power and radio frequency applications. As they operate at large bias and electric fields, a comprehensive analysis of the high-field transport properties is fundamentals, as hot electrons are expected to play a relevant role for the device reliability. In this perspective, Technology Computer-Aided Design (TCAD) simulations can be a very useful tool for the understanding of the phenomena dominating hot-electron degradation mechanisms. The most-accurate modeling approaches are based on the direct solution of the Boltzmann equation, which is not actually available for the GaN material. In this work, the deterministic solution of the Boltzmann transport equation via the spherical-harmonics expansion (SHE-BTE), as incorporated in a commercial TCAD tool, has been extended to the analysis of GaN electrons. To this purpose, the details of the full-band structure has been derived from DFT calculations as in state-of-art literature works, and the electron density of states, $g(E)$ , and group velocity $u_{g}(E)$ , have been calculated for the SHE-BTE for the first time. In addition to this, an accurate calibration of the total scattering rate accounting for nonpolar acoustic and optical carrier-phonon interaction, Coulomb scattering and impact ionization has been carried out against available Monte Carlo data and experiments. The proposed model is also shown to correctly predict the temperature dependence of the electron impact-ionization coefficient and current density up to breakdown.


I. INTRODUCTION
Gallium Nitride HEMTs are nowadays used in medium-high power conversion systems due to their superior electrical performance including high operation frequency. As a wide band-gap semiconductor, it is characterized by high breakdown, which results in superior voltage rating [1], [2]. However, the reliability problems are still the hurdles for a wider application of such technology. The most relevant issue for state-of-the-art GaN-based technologies was the so-called dynamic-Ron, caused by the presence of traps in the GaN channel and at the GaN/AlGaN interface, which has been solved through a careful optimization of the epitaxial stack The associate editor coordinating the review of this manuscript and approving it for publication was Vittorio Camarchia . and of the back end [2]. More recently, stressing a device in the ''semi-ON'' condition led to the generation of defects by hot electrons [3]. This degradation can be caused by a number of mechanisms, often field-driven and strongly correlated with the specific HEMT geometry. Thus, it is difficult to extrapolate insight from the experimental characterizations without an accurate physical modelling. The generation of traps by hot electrons is usually the main cause of permanent degradation, a fundamental limit to ruggedness and longterm reliability. Similar effects have been observed in p-GaN gate HEMTs under repetitive short-circuit stress, where hotelectron-induced defects result in a irreversible shift of V TH and a reduction of the drain current [4]. The past experience with silicon has shown that accurate predictions of hot-carrier-stress degradation can be obtained if the carrier VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ dynamics is described by using a microscopic model, thus solving the Boltzmann transport equation (BTE) for electrons  including the full-band structure and the main scattering rates.  A model was incorporated in the framework of the Synopsys  TCAD tool, based on the deterministic solution of the BTE  based on the spherical-harmonics expansion (SHE-BTE) [5], [6], accounting for the single-electron (SE) and multipleelectron (ME) bond-breaking processes [7], [8]. In general, only electrons with energies greater than a specific activation energy can create or reconfigure defects through a reaction cross-section [8]. The distribution of carriers in energy is thus a key ingredient to this purpose. Although a number of works already addressed the solution of the BTE for GaN by using the Monte Carlo approach [9], [10], the extension of the SHE-BTE solution in the Synopsys TCAD tool for GaN was still lacking and is addressed in this work. The SHE-BTE equation reads [5], [6]: where f (r, E) is the electron density distribution function, g(E) and u g (E) are density of states and group velocity. 1/τ (r, E) is the total scattering rate which accounts for the description of the most relevant scattering mechanisms as usually addressed in Monte-Carlo simulations [9], [10]. s(r, E) is the net in-scattering rate due to inelastic scattering and generation-recombination processes and can be expressed as a function of f and τ . While g(E) has been already reported in different works based on DFT calculations, u g (E) is still not available to the scientific community.
In this work, we calculated u g (E) for electrons in wurtzite GaN (w-GaN) starting from DFT simulations. An accurate calibration of the total scattering rate accounting for non-polar acoustic and optical carrier-phonon interaction, Coulomb scattering and impact ionization has been carried out against available Monte Carlo data and experiments. To this purpose, the SHE-BTE has been solved in two test structures: a GaN resistor and a vertical P-i-N diode.

II. GROUP VELOCITY CALCULATION
In the SHE-BTE equation, the band structure of w-GaN appears through g(E) and u g (E), which can be obtained directly from the full-band system. Available tools for electronic structure provide calculation g(E) but not u g (E). Thus, u g (E) has been determined following [11], [12] from the integration over the angles of the term u 2 g (E)g(E) where 1/(2π) 3 is the density of states in the (r, k) space, i is the band index, u i = (1/h)∇E i . The factor 2 accounts for the spin degeneracy. To this purpose, the electronic structure of w-GaN has been computed using the density-functional-theory (DFT) approach included in the QUANTUM ESPRESSO package [13]. Projector augmented  Density of states of w-GaN calculated according to [11]. Lines corresponding to g 1 , g 2 , g 3 and g 4 show the contribution of the first 4 bands to g(E ), respectively. Black solid line: g(E ) calculated accounting for 8 conduction bands. The vertical dotted line represents the maximum value of energy at which the approximation with NB = 4 is valid..
wave (PAW) pseudo-potential has been used with the energy and the charge cut-off fixed to 80Ry and 640Ry, respectively. Following Monkhorst and Pack [14], a 8×8x8 set of k-points distributed homogeneously in the Brillouin zone, with rows or columns of k-points running parallel to the reciprocal lattice vectors has been used to sample the first Brillouin zone. The exchange correlation potential has been modelled according to the generalized gradient approximation (GGA) as parametrized by Perdew, Burke and Ernzerhof (PBE) [15].
The approach proposed in [11] requires to extract the values of the band energies and their first derivative on a 21 × 21x21 k-points grid in the first octant of the Brillouin zone. In order to reduce the computational cost, the sampling of the Brillouin zone has been carried out by using the Wannier90 tool. It allows, starting from the solution of the Kohn-Sham . Product u 2 g (E )g(E ) of w-GaN calculated according to [11]. Lines corresponding to u 2 g,1 g 1 , u 2 g,2 g 2 , u 2 g,3 g 3 and u 2 g,4 g 4 show the contribution of the first 4 bands to u 2 g (E )g(E ), respectively. Black solid line: u 2 g (E )g(E ) calculated accounting for 8 conduction bands.
equations as that obtained by QUANTUM ESPRESSO on a coarse k-points grid, to obtain the energy eigenvalues in any arbitrary k by using an inverse discrete Fourier transform [16]. In Fig. 1, the validity of the approach has been verified by comparing the band structures calculated with Wannier90 with those obtained by Quantum Espresso. In Fig. 2, g(E) calculated according to [11] is reported along with the separate contributions of the first four bands. As far as the TCAD framework dealing with the BTE solution, the maximum number of conduction bands N B usually adopted for Silicon is 4. The dotted vertical line in Fig. 2 represents the maximum value of the energy at which the approximation with N B = 4 is valid. In w-GaN such value is E = 6.75eV which is well above the barrier height for electrons at the GaN/AlGaN interface meaning that accurate simulations of the most relevant degradation mechanisms can accurately be performed. Similarly, in Fig. 3, the product u 2 g (E)g(E) calculated through Eq. (2) has been reported for the first time. This quantity is directly used in the left-hand side of Eq. (1) along with the total scattering rate.

III. MODELLING OF TOTAL SCATTERING RATE
The total scattering rate 1/τ (E) plays a key role in determining the microscopic carrier mobilities in a semiconductor since it accounts for collisions with charge carriers, phonons and impurities. In the case of the SHE-BTE equation implemented in the Synopsys TCAD, five scattering mechanisms are modelled and implemented as follows [7]: where 1/τ ac (E) is the scattering rate due to acoustic phonon, 1/τ opa (E) and 1/τ ope (E) are related to the absorption and emission of a optical phonon, 1/τ C (E) accounts for the Coulomb scattering mechanism and 1/τ ii (E) accounts for impact ionization phenomena. Differently from the Monte Carlo approach, in the SHE-BTE formalism the scattering rates are directly given as functions of energy by accounting for the integral over the angles of the scattering matrices [6]. The corresponding models of the scattering rates for acoustic, optical emitted and absorbed phonons read [7]: where k B , is the Boltzmann constant, T is the absolute temperature, ρ is the material density, h is the Planck constant, c L is the sound velocity, ε op is the optical phonon energy, N op is the number of optical phonon while and D ac and D op are the acoustic and optical deformation potential, respectively. All the parameters of the electron-phonon collisions have been fixed to values reported in ab-initio calculations of w-GaN [9], [17], with the exception of the D op , which has been determined by fitting the drift velocity curve as a function of the electric field. To this purpose, TCAD simulations have been performed on a two-dimension n-doped GaN resistor. A schematic representation of the device is shown in Fig. 4 (a). The 2-D Poisson and drift-diffusion equations are solved self-consistently to obtain the field profile, the quasi-Fermi level, and the charge distribution in the semiconductor. The electric field and quasi-Fermi level are then used as input for the solution of the SHE-BTE and get the drift velocity calculated by accounting for the full-band structure and the electron distribution function at each bias. As shown in Fig. 5, by tuning the D op parameter, the drift velocity at large electric fields nicely reproduces the data reported in [9], [18], and [19]. It is worth pointing out that such procedure has no influence on the Coulomb and impactionization parameters, which can be determined independently. The resulting values are reported in Table 1. As a further check, in Fig. 6, the obtained total phonon scattering rate is compared with the corresponding ones adopted in [9] and [10]. The Coulomb scattering rate 1/τ C is modelled according to the Brooks-Herring formulation and depends on the total impurity concentration, N i = N D + N A , through the function ζ (N i ) [7]: where N D and N A are the donor and acceptor concentration, respectively. The function ζ (N i ) is a tabulated fitting function introduced to match the experimental low-field mobility curve. More specifically, it accounts for multiple collisions and impurity-clustering effects [6] and needs to be treated as a fitting parameter. Thus, following the previous procedure, VOLUME 11, 2023 the role of the Coulomb scattering term has been addressed by comparing the calculated drift velocity against experiments and other numerical results in Fig. 5. As expected, the role of the Coulomb scattering is limited to the lower electric fields, leading to the reduction of the drift-velocity slope and consequent reduction of the peak. By tuning the ζ (N i ) function, the Monte Carlo data at low fields can be nicely matched, even if they both overestimate the experimental data due to the effect of dislocation scattering as explained in [9]. It is worth noting that the GaN material shows additional scattering mechanisms with polar phonons, which are similar to longitudinal optical phonons, giving rise to an additional scattering contribution a relatively low energies which cannot be accounted for in the SHE-BTE model. As a consequence, the Monte Carlo data show a change in slope at 50 kV/cm which is not reproduced by the proposed approach, leading to a discrepancy between the curves at medium electric fields. For this reason, by using the SHE-BTE approach, the peak of the drift velocity is shifted toward lower electric field. At large fields the drift velocity is anyway nicely reproduced as it depends on the most relevant optical-phonon scattering rate. The use of a full-band numerical structure allows one to accurately account for high-energy effects, which are the most relevant ones for the hot-carrier dynamics. When the absolute value of the electric field exceed 2·10 6 V/cm, impact ionization phenomena start to play a fundamental role, leading to an exponential increase of the electron concentration at the onset of avalanche regime. From a microscopic point of view, impact ionization is modelled as a scattering rate with the following empirical equation [7]: where is the Heaviside function and ε ii,j , v j and s j are fitting parameters. They have been determined in order to reproduce the Monte-Carlo simulations reported in [9] as shown in Fig. 7. All parameters are reported in Table 1. To validate our approach, numerical simulations of a vertical P-i-N diode have been carried out in reverse-biased conditions up to the avalanche onset at different temperatures. The reference device is taken from [20] along with the  Comparison between the total phonon scattering rate (1/τ ph (E ) = 1/τ ac (E ) + 1/τ ope (E ) + 1/τ opa (E )) calculated through Eqs. (4)- (6) and the rigid ion model reported in [9] and [10].   of the non-punch-through diode is reported in Fig. 4 (b). The same approach used for the GaN resistor has been followed: both electron and hole equations are solved, the Poisson solution is then used as input for the SHE-BTE solution of electrons. The Okuto model for the electron and hole impact-ionization coefficients has been used in the drift-diffusion equations with the parameters reported in [20]. The SHE-BTE has been solved at each bias to get the electron current density and to calculate the electron impactionization coefficients. The current densities calculated with the two approaches are reported in Fig. 8 for different values of the ambient temperature showing that the onset of the avalanche regime is perfectly reproduced by the SHE-BTE data. In Fig.9, the impact ionization coefficient α n extracted from the SHE-BTE results are reported along with the Okuto model implemented in the drift-diffusion solver. Both the field and temperature dependence of α n are nicely predicted by the proposed approach. Finally, the electron density as function of energy at the p-n junction is reported in Fig. 10   FIGURE 9. SHE-BTE Impact Ionization coefficient as function of the reverse electric field reported for different temperature and compared with the Okuto model implemented in the Drift-Diffusion solver extracted from [20]. | ⃗ E | = 2.9 · 10 6 V/cm. at different temperatures. A significant number of carriers at high energies is observed with two peaks at 4.9 and 6.8 eV, respectively. The expected reduction of the energy tail with temperature is also observed. The high kinetic energy reached by electrons in such regime might generate defects leading to degradation.

IV. CONCLUSION
In this work, a setup for the solution of the deterministic Boltzmann equation using the spherical-harmonics expansion (SHE-BTE) for electrons in wurtzite GaN has been provided for the first time by using a commercial TCAD framework. The electron group velocity as required by the tool has been calculated using DFT simulations. The scattering mechanisms with phonons and impurities have been calibrated with first principle calculations of the drift velocity reported in literature and validated against corresponding Monte Carlo scattering rates. The parameters for the impact-ionization scattering rate have been fitted against Monte-Carlo data and simulations of a reverse-biased GaN P-i-N diode have been carried out to verify the predictivity of the SHE-BTE solution at very high electric fields at different temperatures. It has been found that the proposed approach nicely reproduces the reference data. Thus, it is now possible to solve the SHE-BTE for electrons in GaN accounting for the high-energy tail of the distribution function. The reported results are a key starting point for the analysis of hot-carrier stress in GaNbased devices.