Impact ionization in β-Ga₂O₃

Wide-bandgap semiconductors are attractive for high-power electronics and UV optoelectronic applications. A recently emerged material β-Ga₂O₃ has gained a lot of attention¹ due to its immense potential in both electronics and photonics. High break-down voltage MOSFETs,^2–6 Schottky diodes,^7,8 and deep UV photodetectors^9,10 are experimentally demonstrated. Well-developed bulk and thin film growth techniques^11,12 make this material a strong candidate for future applications. Accurate n-type doping and the difficulty in p-type doping make electrons the dominant carriers in this material. The electronic structure,^13,14 optical absorption,^15,16 and lattice dynamical calculations^17–19 in this material have been reported in the last few years. Theoretical investigation of electron transport in this material is crucial to augment the experimental advancements. There have been a few low-field transport calculation reports^20,21 in this material, revealing that the long-range polar optical phonon (POP) electron-phonon interactions (EPIs) limit the electron mobility. Recently, the authors reported a high-field transport calculation including full-band EPI to predict velocity-field curves in this material²² for an electric field up to 0.4 MV/cm. However, as the electric field is further increased, interband transitions and electron ionization become important. Indeed, in power devices, the electric field reaches up to several MV/cm and the resulting ionization of electrons could lead to breakdown of devices. An empirical estimate of the critical breakdown electric field in β-Ga₂O₃ based on the bandgap was reported.¹ The authors have previously reported²³ on the impact ionization co-efficient in β-Ga₂O₃ using a simple classical calculation. But, given the exponential sensitivity of the ionization co-efficients on the electric field, it is crucial to carry out a much rigorous calculation from first-principles. Proper understanding of the electron-electron interaction is the key to probe impact ionization.

Recent advancement in the electronic structure and lattice dynamics calculations motivates accurate model development for non-equilibrium carrier dynamics. There have been several reports on the first-principles EPI calculation and subsequent carrier relaxation time estimation with high accuracy.^24–26 There are reports on the high-field transport based on empirical pseudo-potential models under rigid-ion approximations²⁷ that include impact ionization. Conventionally, deformation potential based theories for EPI and a Keldysh empirical ionization model are common practice^28,29 in Monte Carlo simulations. Here, the authors explore impact ionization starting from density functional theory (DFT) under local density approximation (LDA). The uniqueness of this work is that the authors utilize maximally localized Wannier functions (MLWFs) to calculate electron-electron interactions (EEIs). While MLWFs have been used in EPI formulation with high accuracy,³⁰ they have not been used to calculate EEI to the best of the authors' knowledge. Using MLWFs in computing EEI helps integrating EPI and EEI under the same theoretical manifold, thus providing a single framework for far-from-equilibrium transport calculations. Using these EEI calculations, ionization rates are obtained using the Fermi-Golden rule. Full-band Monte Carlo (FBMC) simulation is carried out to extract the IICs. The IICs are fitted to a Chynoweth model to help calibrate device simulators. First, the authors discuss the theory and methods for the calculations followed by the results obtained for β-Ga₂O₃.

The primary interaction mechanism that limits transport of hot electrons is EPI. The effect of EPI is shown in a schematic on the left pane of Fig. 1. The type of EPI that controls the transport depends on the regime of the electronic transport. Under a low-electric field, the dominating EPI mechanism in ionic semiconductors is due to polar optical phonons (POPs). The POP scattering rates are calculated using the Fermi-Golden rule as follows:

Here $k, q, m$⁠, and $υ$ are the electronic wave-vector, phonon wave-vector, electronic band index, and phonon mode index, respectively. $wq$ is the weight of the $q$ point arising from the sampling of the phonon Brillouin zone. $Emk$ is the Kohn-Sham energy eigen value of an electron in band $m$ and wave-vector $k$⁠, and $ωυPOPq$ is the phonon energy eigen value of a phonon in mode $υ$⁠. $δ$ is a small smearing in energy. $gPOPυ q$ is the long-range coupling elements calculated using the Vogl model.³¹ The Vogl model requires the Born effective charge tensor and the phonon displacement patterns both of which come from DFPT calculations. Also note that in the current formulation, $gPOPυ q$ is independent of the electronic wave-vector. This is because under the long-wavelength limit (⁠$q→0$⁠), the overlap of the initial and final electronic states from the same band results in unity due to the orthogonality of the Kohn-Sham states. This is also the reason behind not considering any interband scattering mediated by POP. The electronic energies used are the Wannier interpolated energies from the actual Kohn-Sham calculation. Note that this formulation takes into account the anisotropies in the POP matrix elements arising from the low-symmetry of the crystal. For screening of the POP elements, a standard Lyddane-Sachs-Teller model is used where it is assumed that a given mode could only be screened by the modes that are higher in energy compared to that mode. No screening contribution from free carriers is taken into account.

As the electric field is increased, short-range non-polar phonon (NP) scatterings become important and those scatterings provide higher momentum relaxation due to the involvement of long phonon wave-vectors. The scattering rate is calculated in a similar way as in Eq. (1) but with the exception that the sum in Eq. (1) runs over final electronic band indices besides $υ,q$ since interband scatterings are possible by short-range coupling. The short-range coupling elements are computed as

Here, $ψkm$ s are the electronic wave-functions and $∂Vscfυ∂q$ represent the perturbation in the periodic crystal potential for a phonon mode $υ$⁠. The computation of $gNPυ, mn k,q$ on fine electronic and phonon BZ meshes is done using an efficient Wannier function scheme.^30,32 Note that the overlap of the initial and final electronic states in this case cannot be approximated unlike the long-range case and needs to be explicitly taken care of. Computational details of the short-range elements in β-Ga₂O₃ are documented in a previous work by the current authors.²² 

The electron-electron interaction is a two particle process and the matrix element for the interaction term can be written as

Here, $φk1n$ denotes the electronic wavefunction with wavevector $k1$ and on band $n$ and $V(r1−r2)$ is the screened Coulomb interaction which (in atomic units) is given by

where $ε$ is the dynamic dielectric constant and $qD$ is the Debye length. Writing the Coulomb potential as the sum of Fourier components and obeying the momentum conservation, one can rewrite Eq. (3) as

Now using maximally localized Wannier (MLW) functions,³² the Bloch electronic wave-functions can be written as $|nk⟩=∑mReUk†eik.Re|mRe⟩$⁠, where $|mRe⟩$ is an MLW function centered at $Re$ and $Uk$ is the gauge-transforming unitary matrix that rotates the gauge of the Bloch functions in a way such that the subsequent Wannier functions achieve minimum spread in real space. Utilizing the orthonormality relation of the Wannier functions and under a small q limit, one can write

A similar expression is used in Ref. 31 to calculate polar electron-phonon interactions in a Wannier gauge. In the current context, the small q limit is justified due to the long-range nature of the Coulomb interactions.

The advantage of the proposed approach of using MLW functions for EEI is the ability to obtain a very fine sampling of the Brillouin zone in calculating the interaction terms. The $Uk$ elements are the eigen functions of the interpolated Kohn-Sham (KS) Hamiltonian. First, the $Uk$ matrices are calculated on the coarse mesh followed by a Fourier interpolation of the KS Hamiltonian in the Wannier gauge, which will give the $Uk$ matrices on a finely sampled Brillouin zone. Impact ionization is a two electron process which involves relaxation of a hot electron in the conduction band and excitation of a valence electron into the conduction band. This is schematically shown along with a Feynman diagram on the right side pane of Fig. 1. Using the prescription given in Eq. (6), the impact ionization interaction term can be written as

Here, $|mk⟩$ and $|nk+q⟩$ are the initial states of the hot electron before and after ionization, respectively, and $|n′k′⟩$ and $|m′k′−q⟩$ are the states of the electron being ionized before and after ionization, respectively. $Vq=4πΩεq2+qD2$ is the Fourier transformed Coulomb interaction. Screening of the Coulombic interaction is taken into account by considering a dynamic polarizability (frequency dependent)³³ under a long-wavelength limit. The ionization rate is calculated in a similar way to electron-phonon scattering rates using the Fermi-Golden rule enforcing the energy conservation $δEmk+ En′k′−Em′k′−q− Enk+q$⁠. As seen on the electron-electron interaction diagram shown in Fig. 1, calculating the ionization rate for the electron at $|mk⟩$ involves summing over the internal degrees of freedom which implies integrations over $k′$ and q while summations over $m′$⁠, $n′,$ and $n$⁠. So, someone interested in just evaluating the ionization rates need not store the matrix elements in Eq. (5); rather, the rates can be computed on the fly by summing over the internal degrees of freedom. However, for final state calculation in the FBMC simulation, the entire matrix elements $gnk1→n′k1+qmk2→m′k2−q ee$ are required and hence they need to be stored. It is important to note that with $k′$ being confined to the valence band manifold only, a much coarser grid is used for the $k′$ integration since valence bands are flat. To correct the mean-field (LDA) estimated bandgap, the conduction band energies are shifted to match the experimental bandgap. It is noted that the mean-field estimated wave-functions are close to the actual quasiparticle wave-functions,³³ and hence, the formulation of $Uk$ starting from LDA estimated wave-functions is justified. Also, it is to be noted that in the current work, only the direct term of the electron-electron interaction is considered. As pointed out in Ref. 34, the contribution from the exchange term could be maximum in the case of point interactions, and under such a situation, the ionization rates could be at most underestimated by a factor of 0.75 due to neglecting the exchange term. However, the Coulomb interaction is far from behaving as a point interaction, and hence, the factor of underestimation is expected to be much lower. For electron-electron scattering rate calculations, an importance-sampling (with Cauchy distribution) of the Brillouin zone is used with peak density of the q points at the zone center and gradually the density decays out towards the zone edges. This helps in capturing the long-range behavior of the Coulomb interaction within a reasonable computational cost. This type of sampling has been used previously for polar optical phonon scattering calculations in GaAs.²⁴ 

Electrons get energized by the applied electric field and relax momentum and energy through EPI and EEI. Under a low-field, this could be treated by relaxation time approximation or iterative techniques starting from equilibrium distribution. However, the high-field drives the distribution far away from equilibrium and the low-field techniques fail. Monte Carlo simulation stochastically computes the trajectories of the ensemble of electrons. The details of the simulation can be found in Ref. 22. EEI is added as an additional scattering mechanism in a way similar to EPI. The EEI is restricted to be only between valence electrons and conduction electrons. In other words, EEI within the conduction band manifold is not considered and that does not affect the computed transport properties since the ensemble averages remain unchanged by mere exchange of momentum and energy within the conduction band. The valence electrons getting excited to the conduction band are referred to as secondary electrons and the high energy conduction band electron losing energy as the primary electron. As described in the previous work²² by the authors, in the case of EPI mediated scatterings, during the final state selection in the FBMC algorithm, the scattering mechanisms were classified based on the final band index of the scattered electron, phonon mode index, polar/non-polar nature of the scattering, and absorption/emission. In the case of EEI, the corresponding mechanisms are classified as the conduction band index of the primary and secondary electrons after ionization and the k-point index of the secondary electron before ionization (which is defined on a much coarser grid than the actual fine k grid, see Sec. III). The generation rate of the secondary electrons (G) is computed to be as the difference between the total number of electrons at the end of the simulation and that at the beginning divided by the simulation time. The IIC (α) is defined as the reciprocal of the mean free path traversed by an electron before creating an ionization. The IIC is extracted from the generation rate using the relation $GF= α(F)vd(F)$⁠, where $vd(F)$ is the drift velocity for an applied electric field F which is calculated in the FBMC simulation.

First, density functional theory (DFT) calculations are carried out on a β-Ga₂O₃ unit cell¹³ under LDA using norm-conserving pseudopotentials³⁵ in Quantum ESPRESSO.³⁶ The Brillouin zone is sampled with a Monkhorst-Pack³⁷ grid of 8 × 8 × 4 with an energy cut-off of 80 Ry to truncate the reciprocal vectors. The Kohn-Sham eigen values are interpolated on a fine electronic grid, and the rotation matrices are computed on the same grid. The β-Ga₂O₃ conventional unit cell is shown in Fig. 2(a) (visualized by Vesta³⁸), and the interpolated KS eigen values are shown in Fig. 1(b) for two reciprocal crystal directions. The inset in Fig. 2(b) shows the BZ (visualized by XCrySDen³⁹) with the corresponding reciprocal directions. To obtain the screening element, $εω$⁠, the full-frequency epsilon calculation is used as implemented in BerkeleyGW.³³ While the details of EPI could be found in previous reports^20,22 by the current authors, here the methodology is described in a few sentences. Using density functional perturbation theory,⁴⁰ the phonon eigen values, displacement patterns, and EPI elements are calculated on the coarse mesh. Next, the Wannier-Fourier interpolation^30,32 is carried out to calculate the EPI elements and the phonon dynamical matrices. Long-range POP scattering is calculated separately following Ref. 31.

The ionization rates are computed using the theory described in Sec. III. There are a total of 6 conductions bands taken in the transport calculation, and the electron-electron self-energies for conduction bands 5 and 6 are shown in Fig. 2(c) along the two reciprocal vector directions. In the Monte Carlo scheme, ionization rates are included only from bands 5 and 6 since the lower bands (1–4) do not have energy states high enough to cause ionization. It is noted that band 4 has some states away from the zone center which have energy higher than the bandgap, but their computed ionization rates (not shown) are much smaller than the EPI scattering rates. It could be seen that the ionization rates are much higher near the zone center compared to the zone edges. This peculiarity arises from the isolation of the conduction band minimum at the Γ point by a significantly high energy difference. To be more specific, the electrons near the zone edge in bands 5 and 6 which have energies higher than the bandgap cannot ionize as there are no available final states. The long-range nature of the interaction prohibits a zone edge electron to release enough energy (>bandgap) and end up near the zone center. Hence, the ionization rates mediated by the zone edge hot electrons are significantly lower than those mediated by zone center electrons. No phonon-assisted ionization mechanism is considered here, however. Figure 3(a) shows the contributions of the individual electronic bands in impact ionization. It could be seen that qualitatively the ionization rate follows a power law²⁹ near the ionization threshold which is slightly above 5 eV. The abrupt drop of the ionization rates around 6.5 eV does not happen due to truncation of the conduction bands after the sixth band. There are significant available states on the sixth band much beyond 6.5 eV as well. The abrupt drop rather arises from the absence of the final density of states after an ionization event. Since the bandgap is around 4.9 eV, the final energy of the ionizing electron has to be at least 4.9 eV less than the energy it had prior to ionization. However, due to the absence of remote satellite valleys within 0–2.4 eV of the CB minima, the hot electrons can only end up on the Γ valley. On the sixth band, the states which possess more than 6.5 eV occur much away from the Γ point and hence lack final density of states to cause an ionizing event.

The best practice to determine the number of conduction bands to be taken in the Monte Carlo simulation should be “a posteriori.” However, there are a few important aspects to be noted here. The hot electron distribution function dies off quickly after its peak and only a minuscule of the population survives beyond 6 eV as is shown in Fig. 3(b). This is the scenario even for the highest magnitude of the electric field considered in this work which is 8 MV/cm. So, the authors speculate that considering beyond 6 bands might not make any difference in the distribution function and hence the ionization co-efficient will not get affected. The top 4 valence bands were used in the ionization rate calculation considering a 1 eV window from the VB maxima. This translates to around a minimum of 6 eV energy for the ionizing electron. To further address the need for an “a posteriori” convergence study with respect to the number of conduction bands, the authors would like to point out the computational limitation caused by the huge memory requirement of storing the EPI elements. The 30 phonon modes along with fine Brillouin zone sampling require an enormous amount of memory. Some programing strategies²² are taken to circumvent the memory overflow challenges and that is why it was possible to carry out the computation with 6 bands. However, considering beyond 6 conduction bands could not be achieved. Please note that the required memory for storing the matrix elements will scale quadratically with the number of bands (since the number of allowed transitions will scale quadratically).

The FBMC scheme initializes the ensemble of electrons thermodynamically after which the electric field is turned on. Trajectories of the electrons are formed in reciprocal space stochastically by using the electron-phonon scattering rates and electron ionization rates. Six conduction bands are taken into account in the FBMC simulation. The FBMC simulation is run for electric fields ranging from 1 MV/cm to 8 MV/cm. Figure 4(a) shows the transient electron dynamics under an applied electric field of 2 MV/cm. The oscillations that are observed initially result when the electrons cross the Bragg planes and they are subsequently suppressed out due to EPI. Under a high enough electric field, the transit time to reach the Bragg plane (let's call it τ_B) can become comparable to mean free time (τ_s) between successive scattering events (τ_B ≈ τ_s) making it possible for the electrons to reach the plane. This oscillation, often known as Bloch oscillation (BO), is experimentally observed in superlattices⁴¹ at room temperature but very rare in bulk semiconductors. In β-Ga₂O₃, the satellite valleys occur at an energy comparable with the zone edge maxima of the first conduction band. Hence, the onset of intervalley scattering occurs only near the zone edge where the reflection (reversal of electronic group velocity) is also likely to onset. The time-period of oscillation (TB) in Fig. 4(a) is comparable with the analytically calculated BO time-period, $TB=ℏeFd$⁠, where $d$ is the distance of the Bragg plane from the zone-center, further confirming the origin of the oscillations. In bulk semiconductors, there are primarily two prohibitive actions that prevent the observation of Bloch oscillations. First is the very high scattering rates mediated by non-polar optical phonons at high electronic wave-vectors demanding for an unrealistically high electric field to satisfy the criterion τ_B ≈ τ_s. Second, under such a high electric field in most semiconductors, band-to-band tunneling comes into action. This significantly increases the net number of electrons having a positive group velocity compared to the ones having a negative velocity (arising due to reflection on the Bragg plane). thereby prohibiting the observation of Bloch oscillation. In the current work, it is shown that it is possible to satisfy the first criterion in β-Ga₂O₃. But no band-to-band tunneling phenomenon is considered, and hence, it cannot be conclusively said whether Bloch oscillation can actually be observed.

Figure 4(b) shows the occupation of the bands as the electric field is increased in two different directions. In this calculation, the interband transitions occur only via short-range EPI and long-range EEI. It could be seen that the population on the first band drops to about half of the total as the field reaches 8 MV/cm. Although the population on bands 5 and 6 is really low even at very high fields, given the fact that impact ionization is a cumulative process, even a small population can eventually trigger avalanche breakdown. Figure 4(c) shows the calculated IICs along three different directions. The anisotropy in the IIC is attributed to the anisotropy in EEI which in turn originates from the anisotropy of the higher conduction bands (bands 5 and 6 in this case) even near the Γ point as seen in Fig. 2(c). There are some uncertainties with the bandgap of β-Ga₂O₃, and a range of values 4.5 eV–4.9 eV could be found in the literature based on experiments. Hence, the calculations are carried out for two different gaps 4.5 eV and 4.9 eV. The uncertainty in IIC due to the bandgap uncertainty is within an order of magnitude. Given the classical exponential dependence⁴² of IIC on the bandgap, this much uncertainty in the IIC is reasonable. Also, as expected, the IIC is lower for a higher bandgap. For comparison, the calculated IIC is less than that computed in GaN²⁷ for a similar range of electric fields which indicates a higher avalanche breakdown field for β-Ga₂O₃.

Finally, to help facilitate device simulation, the authors perform a Chynoweth fitting of the calculated IIC and attempt to identify favorable transport directions for high power applications. The Chynoweth model⁴³ formulates the IIC as $αF=ae−bF$⁠, where $a$ and $b$ are the fitting parameters. In order to account for the anisotropy in IIC, three sets of the parameters are provided for three different directions in Table I. In electronic devices under high electric fields, impact ionization induced avalanche breakdown occurs if the generation of secondary carriers becomes self-sustainable which in turn requires the ionization integral⁴⁴ to be greater than 1. To have a qualitative understanding of the avalanche breakdown field in β-Ga₂O₃, a hypothetical triangular electric field profile (which is approximately the case in p+-n junctions) of the peak electric field of 8 MV/cm is considered and the base of the triangle to be 1 μm wide. Considering the estimated Chynoweth parameters, the ionization integral values are also tabulated in Table I. The computation of ionization integral (I_α) requires knowledge of the ionization co-efficient of both electrons (α_n) and holes (α_p). However, the hole ionization coefficients for β-Ga₂O₃ have not been estimated. So, the ionization integral is shown for the two limits α_n ≈ α_p, and α_p ≪ α_n. As could be seen, for the considered electric field profile, avalanche breakdown is very likely in the y direction while very unlikely in the x direction. So, for high voltage/power applications, electron transport along the x direction is more favorable than that along the other two directions. Based on the hypothetical electric field profile, the critical electric fields are evaluated under the limit α_n ≈ α_p and are also listed in Table I. As revealed in a previous work⁴⁵ by the current authors and also from experimental observations,⁴⁶ the electron mobility is slightly lower in the x direction compared to that in the y direction. But considering the joint effect of on-resistance and breakdown field, transport along the x direction is supposed to provide a higher Baliga's figure of merit.⁴⁷ Hence, one can conclude that the x direction is the most suitable transport direction for high power applications using β-Ga₂O₃.

Maximally localized Wannier functions are utilized to compute the ionization rates in β-Ga₂O₃ from first-principles. FBMC simulation is done to compute IICs for a wide-range of high electric fields along two different directions. IIC is fitted to an empirical Chynoweth model to calibrate device simulators. A hypothetical estimate using the computed Chynoweth parameters predicts the avalanche breakdown field to be higher than the empirically predicted value of 8 MV/cm in the x direction.

The authors acknowledge the support from the National Science Foundation (NSF) Grant No. (ECCS 1607833) monitored by Dr. Dimitris Pavilidis. The authors also acknowledge the excellent high performance computing cluster provided by the Center for Computational Research (CCR) at the University at Buffalo.