Ab initio calculation of electron–phonon coupling in monoclinic β-Ga₂O₃ crystal

Recently, β-Ga₂O₃ has emerged as an extremely promising candidate for power electronics due to its large bandgap (4.9 eV) and experimental demonstration of β-Ga₂O₃ transistors with very high breakdown voltage (∼750 V).^1–4 Apart from transistors, Schottky barrier diodes based on this material have also been demonstrated.^5–7 A well-developed bulk and epitaxial crystal growth technology further provides cost advantages.^1–4 It is crucial to have a thorough understanding of the electron transport properties to get an insight into the device operation and also optimize the performance. Although there have been several studies on the electronic structure^8–11 and dielectric properties^12,13 of the material, only a few reports^14,15 are available on carrier transport in the material. Unlike the traditional semiconductors, β-Ga₂O₃ has a highly anisotropic monoclinic base centered crystal (C2/m space group) structure with 10 atoms in the primitive cell.¹⁰ This leads to the large number of phonon modes (30 modes) in the material which makes it complex to analyze the interaction between electrons and phonons that plays a pivotal role in determining mobility. The low symmetry of the crystal poses question on the anisotropy of carrier transport. Though an early experiment showed anisotropic carrier mobility,¹⁶ later on it was found that the low field mobility is actually isotropic.¹⁵ Also, given the presence of both non-polar and polar interactions between the electrons and the phonons, it is critical to evaluate which type would dominate electron transport in β-Ga₂O₃. Conventional non-polar scattering rate equations require deformation potentials (DPs) which are usually found out empirically while polar scattering rates require knowledge of relative contribution of the phonon modes to the dielectric constant. None of these are well established theoretically or experimentally for β-Ga₂O₃. A recent work¹⁵ on low field transport reported the dominant optical phonon mode from empirical fitting of experimental data; however, they did not consider all the phonon modes in the calculation. In this work, we theoretically study the low field electron transport in β-Ga₂O₃ from first principles taking into account all the phonon modes. The dominant phonon mode limiting the low field mobility is identified.

First, we carry out the standard density functional theory (DFT) calculation on a β-Ga₂O₃ primitive cell which contains 4 Ga atoms and 6 O atoms using the ab initio package Quantum ESPRESSO (QE).¹⁷ The real space lattice and the Brillouin zone (BZ) are visualized by XCrySDen¹⁸ and are shown in Figs. 1(a) and 1(b), respectively, with the BZ showing the high symmetry points identified as black dots. The reciprocal vectors are also shown. The details about the DFT calculation are given in the supplementary material. The calculated energy eigen values are interpolated using maximally localized Wannier (MLW) functions.¹⁹ For low field transport, we are interested only with the first conduction band which is dominated by Ga 4s orbital.⁹ The Wannier interpolated first conduction band matches well with previously reported calculations (see Fig. S1).

After obtaining the electronic energy dispersion, we focus on finding out the lattice response under the density functional perturbation theory (DFPT) formalism²⁰ using QE. The details of this calculation are again given in the supplementary material. The dynamical matrices (D_q) and the perturbation in the self-consistent potential are calculated for phonon wave-vectors (we refer this as q points) on a coarse mesh. Next, a Fourier interpolation is done to obtain the eigen values at intermediate q points. The resulting phonon spectrum is shown in Fig. 2(a), which shows the phonon modes dispersion along the three reciprocal vector directions. It shows the 30 phonon modes arising from the 10 atoms in the primitive cell. An important point to note here is that the phonon modes shown in Fig. 2(a) do not take into account the contribution from the macroscopic field due to presence of effective polarization charges on the atoms. The long–range interaction arising from this macroscopic field is calculated separately in the next paragraph. In the following paragraph, we discuss the electron-phonon interaction elements in β-Ga₂O₃.

The interaction elements for the νth phonon mode is given by $Mqν=⟨ψk+q|∂Vscf∂q|ψk ⟩$⁠, where $∂Vscf∂q$ is the perturbation in the self-consistent potential (V_scf) due to phonon wave-vector $q$ and $ψk$ is the electronic wave-function. These interaction elements are calculated on the coarse mesh first. However, transport calculations need very fine resolution of the matrix elements in order to achieve convergence in scattering rates. The coarse mesh matrix elements are interpolated on to a fine mesh using a Wannier–Fourier (WF) interpolation²¹ scheme (Electron–Phonon Wannier (EPW)) utilizing the crystal periodicity. However, as pointed out in the literature,^22,23 the polar interactions between electrons and the polar optical phonon (POP) modes due to a macroscopic electric field should be treated separately due to the long range nature of the polar interaction. Hence, we treat the long-range POP contribution separately. First, we calculate the electron–POP interaction on the coarse mesh by the method described by Verdi.²³ The Born effective charge tensor and the dielectric constant are calculated by DFPT for q = 0. The electron-POP interaction are calculated on the coarse q-mesh using Eq. (4) of Ref. 23. These elements are subtracted from the DFPT calculated matrix elements on the coarse mesh. This leaves us with the short range (cell-periodic part; both non-polar optical and acoustic interactions) part of the interaction on the coarse mesh. Next, we carry out the WF scheme to find out the matrix elements on a fine mesh. The real-space formulation needed for the interpolation is described in the supplementary material. Apart from interpolating the electron–phonon matrix elements, the WF scheme also interpolates the phonon dynamical matrices and hence we could obtain the phonon displacement vectors and eigen values on a fine mesh. These are used to recalculate the electron POP interaction elements on the fine mesh using the same the method described in Ref. 23. It is important to note here is the fact that we use two different fine q-meshes for long range (POP) and short range matrix elements to facilitate our transport calculation. The long range elements (POP) are calculated on a 40 × 40 × 40 uniform mesh covering 0.2 times the length of the reciprocal vectors in each direction surrounding the Γ point while the short range elements are calculated on a similar 40 × 40 × 40 uniform mesh 0.6 times the length of the reciprocal vectors. This choice of meshes ensures accuracy in the long-range elements (since they decay fast near the Γ point) without demanding excessive memory for computation.

Figures 2 (b)–2(d) show the calculated electron-phonon long–range (POP) coupling strength $(gqν=ℏ2ωqυ Mqν)$⁠, in three different Cartesian directions in the BZ, $ωqυ$ is the phonon mode frequency and $ℏ$ is the reduced Planck's constant, $Mqν$ is described before, and the mass factor of the quantized vibrational amplitude is actually taken care of within $Mqν$⁠. It is observed that the coupling strengths are direction dependent, with different modes having strong strengths in different directions. This will play an important role in mobility calculations. There are several significantly strong coupling strengths that could be observed, for example, 53 meV and 71 meV (x and z directions), 39 meV and 56 meV (y direction), 21 meV (in z direction) and 29 meV (in x direction). The polar coupling strength of a mode in a given direction is due to the projection of the net dipole strength along that direction (since polar coupling is the dot product of dipole strength and q vector). A discussion on effective dipole strength in different directions (along with LO–TO splitting) and comparison of corresponding infra-red (IR) strength are given in the supplementary material. The most important information evident from the direction dependent POP coupling elements is the fact that the B_u modes, polarized only in the x-z plane, have different coupling strength in x and z directions which is also evident from the direction dependent LO–TO splitting shown in Table SI. On the other hand, the A_u modes are polarized only in the y direction. Hence, one must be careful while carrying out transport calculations since the conventional way of considering a spherical average in the q-space is not a good approximation in this case. Another important issue to consider is the screening of the vibrational modes. We take care of screening by the fact that higher energy modes can efficiently screen the low energy modes. Note that this screening is taken care of for short range modes automatically since the short rage matrix elements follow the self-consistent potential. However, we need to treat this screening separately for the POP. Here, we include the screening by using a frequency dependent dielectric constant based on Lyddane–Sachs–Teller (LST) theory²⁴ (see the supplementary material for more details).

Next, we calculate the electron scattering rates under Fermi golden rule. As discussed previously, we need to consider direction dependent matrix elements and hence the conventional way of carrying out the q-space integration with spherical symmetry is not a good choice. Here, the q-space integration is performed numerically using a Gaussian smearing of the energy conserving delta function, which is a common practice in numerical BZ integrations.²⁵ The resulting scattering rates are shown in Figs. 3(a) and 3(b). Fig. 3(a) shows the three main different scattering mechanisms that affect the electron transport. The POP scattering rate clearly shows the onset of the emission of the dominant POP mode near 21 meV. The non-polar scattering is a result of combined deformation potential (DP) interaction with acoustic and non-polar optical phonons. It is noted that we calculate the ionized impurity (II) scattering rate using a simple Brooks Herring model²⁶ with a static Debye like screening. An n-type doping of 1.1 × 10¹⁷cm⁻³ and an effective mass of 0.3 times free electron mass is taken for the calculation. Incomplete ionization of the dopants at low temperature is taken into account with a dopant activation energy of 31 meV.⁷ Fig. 3(a) also shows the anti-symmetric part of the electronic distribution function obtained after solving the Boltzmann Transport Equation (BTE) through iteration which is discussed later.

Fig. 3(b) shows the individual contribution of the POP modes to the scattering rate. There are several factors, namely, the coupling strength (⁠$gqν$⁠), energy of the interacting electron (E_k), and the temperature, which control which mode will dominate in determining the low-field mobility. We discuss the role of each of these factors here and show how different modes could become important at different regimes of transport. The coupling strength of a mode has been discussed in detail before. The role of electron energy (E_k) comes in the form of available density of states after the scattering occurs. To clarify further, scattering rate due to a phonon mode with higher coupling strength could be suppressed if the available electronic density of states after emission (or absorption) of that mode is low. So under low electric field conditions, when the anti-symmetric electron distribution is close to the Γ point (see Fig. 3(a), emission of a high energy phonon mode is unlikely even though the coupling strength of that mode is high. Finally, the temperature influences the interaction via the Bose-Einstein occupation statistics which is common to all forms of phonon scattering. A higher temperature will increase the occupation of both low and high energy phonon modes and hence the POP modes at higher energies would also start to affect the electron transport.

The three bar plots in Fig. 3(b) show the electronic scattering rate due to POP interaction for three different electron energies. We see that the most dominant phonon mode for a low energy electron (E_k ≈ 60 meV) is 21 meV. As shown in Fig. 3(a), in low field conditions the electron distribution dies off quickly with energy, hence low energy electrons determine the mobility. On the other hand, for relatively higher energy electrons (E_k ≈ 170 meV) the dominant phonon mode is at 71 meV. This is due to the fact that although the emission of a 71 meV phonon by a 60 meV electron is prohibited, the same by a 170 meV electron is feasible. This trade-off between the phonon mode and the electronic kinetic energy is what makes the electron–polar mode interaction interesting in crystals with multiple LO modes. Hence, as we can see that different modes could affect the electron transport based on the electron energy, it is difficult to specify a mode that will determine the transport across the entire temperature and transport regime. However, for analyzing low-field transport using the Frölich model for POP scattering, a POP mode with energy around 21 meV will give approximate scattering rates based on the calculations discussed above. Hence, majority of the electrons will interact with the 21 meV mode. It is noted that a non-polar phonon of 30 meV was shown to be the dominant mode in Ref. 15, similar low energy phonon mode is also seen here for room temperature mobility. However, Ref. 15 did not consider POP scattering due to the low energy phonons, instead an empirical deformation potential was used to fit the experimental data.

We use Rode's method^27,28 to iteratively solve the BTE which takes into account the inelastic nature of POP scattering. Here, we would like to mention that impurity scattering entered the BTE under the well-known relaxation time approximation (RTA). Fig. 4 shows the calculated net mobility as a function of temperature. We find that at low temperature (30 K–80 K) the mobility shows a brief increasing trend which is due to the II scattering dominated transport. As the temperature is further increased, the POP scattering plays the dominant role and the mobility starts dropping. At room temperature, the estimated mobility is around 115 cm²/V s. We also show the mobility data from recent experimental study by Parisini et al. Good agreement between our calculation and the experimental data could be observed for a wide range of temperature (30 K–650 K).

In conclusion, we have theoretically investigated the electron-phonon interaction in β-Ga₂O₃ from first principles. Low field mobility values are predicted for a varying range of temperature. The interesting physics of electron POP interaction is explored for multi-phonon semiconductors. The polarization of the phonon mode, the electron energy, and the temperature play a crucial role in determining which phonon modes limit the electron transport. At room temperature, POP mode with energy of 21 meV limits the low field mobility.

See supplementary material for descriptions of DFT, DFPT, and electron–phonon calculations, details on dipole oscillation and dielectric tensor, Wannier representations needed for interpolation, and high temperature effects on POP scattering rates.

The authors acknowledge the support from the National Science Foundation (NSF) grant (ECCS 1607833). The authors sincerely acknowledge the excellent high performance computing cluster provided by the Center for Computational Research (CCR) at the University at Buffalo (UB).