Structural, electronic, elastic, power and transport properties of $\beta$-Ga$_2$O$_3$ from first-principles

We investigate the structural, electronic, vibrational, power and transport properties of the $\beta$ allotrope of Ga$_2$O$_3$ from first-principles. We find phonon frequencies and elastic constants that reproduce the correct band ordering, in agreement with experiment. We use the Boltzmann transport equation to compute the intrinsic electron and hole drift mobility and obtain a room temperature values of 258 cm$^2$/Vs and 1.2 cm$^2$/Vs, respectively as well as 6300 cm$^2$/Vs and 13 cm$^2$/Vs at 100 K. Through a spectral decomposition of the scattering contribution to the inverse mobility, we find that multiple longitudinal optical modes of B$_u$ symmetry are responsible for the electron mobility of $\beta$- Ga$_2$O$_3$ but that many acoustic modes also contributes, making it essential to include all scattering processes in the calculations. Using the von Hippel low energy criterion we computed the breakdown field to be 5.8 MV/cm at room temperature yielding a Baliga's figure of merit of 1250 with respect to silicon, ideal for high-power electronics. This work presents a general framework to predictively investigate novel high-power electronic materials.


I. INTRODUCTION
The β allotrope of Ga 2 O 3 has attracted some attention as an ultra wide bandgap transparent semiconducting oxide [1].As a consequence of its large bandgap, β-Ga 2 O 3 possess a very high breakdown electric field of 8 MV/cm [2] and a large Baliga's figure of merit (BFOM) [3], which makes it a promising alternative to GaN and SiC for high-power electronics [4,5].In addition it can be synthesized by melt-growth method which allows for low cost and large-scale production [6,7].Its electronic and optical properties also make it a good candidate for UV transparent conducting oxide (TCO) [8,9].
One property of β-Ga 2 O 3 that makes it so attractive is its high carrier mobility for a material with such a wide bandgap.The electron mobility of β-Ga 2 O 3 has been studied more extensively than the hole mobility due to experimental interest and the fact the hole mobility is two order of magnitude smaller.Given the promise offered by β-Ga 2 O 3 , it is surprising that many basic properties have not been investigated in detail.From a theoretical perspective, this might be due to the fact that β-Ga 2 O 3 has a 10-atom primitive cell, which makes firstprinciples calculations in this material more challenging than for standard tetrahedral semiconductors.In particular, the shape of the conduction band was not well understood until recently.Indeed, Ueda et al. [10] measured a strong anisotropy of the conduction-band effective mass.However, since then many experiments and theoretical studies indicated that the conduction band is nearly isotropic [6,[11][12][13][14][15][16][17].Another question relates to the relative importance of nonpolar optical-phonon, polar optical-phonon, and ionized-impurity scattering at room temperature.Initially it was thought that the dominant scattering mechanism in β-Ga 2 O 3 was due to nonpolar optical phonons with a large deformation potential of 4×10 9 eV/cm [18].However, Ghosh and Singisetti [19] identified a longitudinal-optical phonon mode with energy around 21 meV as the dominant mechanism in the mobility of β-Ga 2 O 3 , and this finding was later confirmed by multiples authors [20][21][22].Finally, there is some debate about the ordering of the zone-centred phonons, namely the Raman-active A g mode and the infraredactive B u TO mode [22][23][24][25].
One crucial material property for high-power electronics is the breakdown field, i.e. the magnitude of the external electric field that a material can sustain before incurring permanent damage.The breakdown field can be computed from first-principles using the von Hippel low energy criterion [26][27][28], and was recently computed ab-initio by Mengle and Kioupakis [22] to be 5.4 MV/cm in β-Ga 2 O 3 considering only the dominant LO phonon mode.They further estimated that considering all modes would increasing the theoretical intrinsic breakdown field by 20% to 6.8 MV/cm.Such calculation assumes total impact ionization for all electrons with energies above the bandgap, and should therefore be seen as a lower bound; it can also be improved by computing impact ionization coefficient from first principles [29].
The BFOM [3] describes the current handling capability of a material and is often given relative to silicon.In addition to the breakdown field, the second material's parameter entering into the BFOM is the intrinsic carrier mobility.The electron room-temperature mobility of β-Ga 2 O 3 was computed to be 115 cm 2 /Vs at a carrier concentration of 10 17 cm −3 , with a temperature dependence in good agreement with experiment [19], by using Wannier interpolation of the electron-phonon matrix elements [30] and Rode's method [31].The mobility was also estimated to be below 200 cm 2 /Vs using k • p perturbation theory [20].
In this context, a careful and detailed analysis of the crystal structure, electronic, optical, vibrational, elastic, and transport properties of β-Ga 2 O 3 is warranted.The manuscript is organized as follows.In Section II we discuss the relaxed crystal structure of the monoclinic β-Ga 2 O 3 and the importance of spin-orbit coupling.Section III is dedicated to the study of the electronic properties including bandgaps, electronic bandstructure and effective masses.In Section IV we analyze the phonon dispersion, infrared and Raman spectra, dielectric constant and Born charges, and elastic properties.The Section V present the computed electron and hole carrier mobility with temperature as well as a mode-resolved analysis of the scattering contribution to the mobility.Finally in Section VI, we discuss and compute Baliga's figure of merit of β-Ga 2 O 3 and compare it with Silicon.

II. CRYSTAL STRUCTURE
The crystal structure of β-Ga 2 O 3 was originally investigated by Geller to be monoclinic with the C 2h (2/m) point group [32] and later refined by Åhman [33] using single crystal diffraction.The measured lattice parameters of the conventional unit cell are a=12.214Å, b=3.037Å, c=5.798Å and β=103.83• [33].The conventional cell vectors are (a, 0, 0), (0, b, 0) and (c cos β, 0, c sin β) while the primitive cell vectors are ( a 2 , − b 2 , 0), ( a 2 , b 2 , 0) and (c cos β, 0, c sin β).Any atomic coordinate expressed in the conventional cell (c x , c y , c z ) can be expressed in the primitive cell by using the transformation (c x − c y , c x + c y , c z ).The primitive and conventional cell are made of 10 and 20 atoms, respectively.
The gallium atom sits in two inequivalent position with octahedral and tetrahedral coordination, respectively.There are three inequivalent oxygen atoms occupying a distorted cubic lattice, with two oxygen atoms being three-fold coordinated and one oxygen atom fourfold coordinated.All the five inequivalent atoms have 4i Wyckoff position which corresponds to symmetry (x,0,z) and (−x,0,−z).The system has four crystal symmetries: the identity, a π rotation around the Cartesian y axis and their inversions.
To determine the atom positions, we relaxed the lattice parameters and atomic coordinates, starting from the experimental data.We used the Quantum ESPRESSO software suite [34] with relativistic LDA pseudopotentials from PseudoDojo [35], including the 3s 2 3p 6 3d 10 4s 2 4p 1 semicore states for gallium and the 2s 2 2p 4 electrons for oxygen.The wavefunctions were expanded in a plane-wave basis set with energy cutoff of 120 Ry (160 Ry for the elastic response) and an homogeneous Γ-centered FIG.1: Relaxed crystal structure of the primtive cell of β-Ga2O3 where the large atoms are gallium and the small red atoms are oxygen.Rendered using VESTA [36].
Brillouin-zone sampling of 8×8×8 points.We converged the structure such that the maximum force was smaller than 2•10 −7 Ry/ Å and the maximum stress component was lower than 0.07 Ry/ Å3 .
The relaxation yielded the lattice parameters a=12.128Å, b=3.016Å, c=5.752Å and β=103.75• , which slightly underestimates the experimental one as expected from LDA.The relaxed primitive cell crystal structure is shown in Fig. 1 and is formed by two distorted octahedra and two distorted tetrahedra.The gallium and oxygen atoms occupy two and three inequivalent sites at the 4i Wyckoff position, respectively, whose coordinates are provided in Table I and are in close agreement with the experimental assignment [33].The inequivalent gallium-oxygen bond lengths are also reported in Table I with the thetrahedra having smaller bond-lengths than the octahedra.Interestingly, due to their distorted nature, there are two inequivalent Ga II -O III bond lengths in the octahedral configuration despite having only one inequivalent oxygen position.
We also report in Table I the volume, density, atomic coordinates, and bond lengths, and compare them with experimental data.The calculations were made without spin-orbit coupling (SOC) but we tested that including this effect modifies the crystal data shown in Table I by less than 0.005%.Hence, this effect is neglected for the rest of this work.We finally note that the primitive cell vectors can equivalently be rotated such that a=b=11.809Å, c=10.869Å, α = β = 103.335• and γ = 27.933• .

III. ELECTRONIC PROPERTIES
The room-temperature optical bandgap of β-Ga 2 O 3 obtained through absorption measurements is estimated to be between 4.54 eV and 4.90 eV [37][38][39].
Our calculated direct bandgap at the zone center is 2.55 eV, strongly underestimating experiments as expected from density functional theory (DFT).In agreement with prior work [15], we find that the valence band maximum (VBM) is located on the I − L high-symmetry lines in the Brillouin zone, and yields a slightly smaller There has been some confusion in the literature about the shape of the Brillouin zone of β-Ga 2 O 3 [11,[40][41][42][43].The first band structure using the correct monoclinic variation was reported in 2015 [15].It is therefore important to pay close attention when constructing the Brillouin Zone of β-Ga 2 O 3 .
We note that, as the definition of two of the primitive cell vectors in the Quantum Espresso software are inverted with respect to prior studies, we had to adapt the definition of the high-symmetry points of the Brillouin zone.We give the conversion for clarity in Appendix Table V as well as the value of the four parameters that define the high-symmetry points.To avoid further confusion, the primitive vectors for the base centered monoclinic Bravais lattice have been modified in Quantum Espresso version 6.5 to use the same definition as in the literature [15].The electronic bandstructure along highsymmetry lines is given in Fig. 2(a), where the highest valence band and lowest conduction band are highlighted in orange color.
We computed the electron effective mass using finite differences, and found 0.267, 0.254 and 0.244 along the Γ-X, Γ-Y and Γ-Z direction, respectively.The electron effective mass is quite isotropic with an average value of 0.255 as reported in Table II, which compares well with prior theoretical work, and is also close to the experimental value of 0.28 [40,42].This level of agreement gives us confidence that our calculations of electronic transport properties will be reliable.
In contrast, the hole effective mass at the zone-center is highly anisotropic, with very heavy masses along the Γ-X and Γ-Y direction, and a small hole mass of 0.35 m e along the Γ-Z direction.As a result, this should be an ideal hole transport direction.However, the VBM is not located at the zone centred but 26 meV higher in energy on the I-L line.The transverse and perpendicular hole effective mass at that point is 3.0 m e and 3.6 m e , respectively in agreement with previous work [11].As transport properties scale inversely with the effective mass, we expect at least an order of magnitude lower hole mobility than the electron mobility.III: Phonon frequencies (meV) of β-Ga2O3 at the zone center using a 16×16×12 k-grid.The three zero-frequency acoustic modes are not reported.The infrared experimental values from Dohy et al. [44] are measured by transmission and the frequency of the maxima lies between the LO and TO frequencies.
The point groups along high-symmetry lines is either C s (m) or C 2 .The C s (m) point group contains two symmetry operations: the identity operation E and a mirror plane σ.This point group possess two irreducible representations: the phonon branches belonging to the A irreducible representation are symmetric with respect to both the identity operation and reflection through the mirror plane while the branches belonging to the A representation are symmetric with respect to the identity but antisymmetric with respect to reflection (coloured in Fig. 2(b) in gray and blue, respectively).The other point group is the C 2 point group which contains the identity (gray) and a π rotation around the [0,1,0] Cartesian axis (displayed with red lines in Fig. 2(b)).Note that some directions in the Brillouin-Zone are less symmetric and only possess the identity (gray).In addition, specific high symmetry points have higher symmetries: (i) the point group at the N and M points is C i (−1) with A g and A u symmetry operation; (ii) the point group at the X, I points is C 2 (2) with identity E and C 2 with π rotation around the [0,1,0] Cartesian axis; (iii) the point group at the

B. Infrared and Raman spectra
The infrared spectrum as well as polarization and temperature-dependent Raman spectra of bulk β-Ga 2 O 3 were first measured by Dohy et al. [44] in 1982.The measured normal modes frequencies are reported in Table III along with more recent measurements and previous ab-initio values, and are compared to the calculated frequencies from this work.Our calculated phonon bandstructure slightly underestimates experiments but are in better agreement than previous calculations.Overall, our calculations agree with previous theoretical work [22,24,47] with a notable difference: in agreement with experiments [23,25] we find that the A g (3) Raman active mode has a lower frequency than the B u (TO1) mode.The highest phonon frequency at the zone center is a LO mode in the Z direction, with a frequency of 97 meV, very close to the experimental value of 100 meV [25].However we note that the highest phonon frequency occurs at the Z point with a value of 99.12 meV (not shown in Table III).Our predicted Raman-active phonon frequencies are 2.5% within the experimental data [23] with the largest difference being attributed to the A g (6) mode.Our predicted infraredactive LO modes are even closer, with deviation of 1.4 % from experimental data [25], while the agreement with LO modes is not as good with a deviation of 5.4%.

C. Dielectric constant and Born charges
The high-frequency dielectric tensor is fairly isotropic, with ε xx = 3.98, ε yy = 4.09 and ε zz = 4.08, slightly overestimating the experimental value of 3.53-3.6[48][49][50] obtained as an isotropic average in thin films.The slight overestimation of the theoretical dielectric tensor is a direct consequence of the underestimation of the bandgap by DFT as the electronic part of the dielectric function is inversely proportional to the bandgap [51].We note one experimental work which obtained a direction-dependent dielectric tensor ε xx = 3.7, ε yy = 3.2 and ε zz = 3.7 [25] using generalized spectroscopic ellipsometry within the infrared and far-infrared spectral region.This anisotropy was not observed in another recent experiment reporting ε xx = 3.6, ε yy = 3.58 and ε zz = 3.54 [52] also using generalized spectroscopic ellipsometry.Our calculations appear to support an isotropic dielectric tensor.β-Ga 2 O 3 also possesses one non-zero off-diagonal component of the dielectric tensor, but the computed value was lower than 10 −4 and therefore is not reported. The

D. Elastic properties
The stiffness C ij and compliance S ij = C −1 ij tensors link the stress tensor to the strain tensor following the generalized Hooke's law: where Einstein's notation is implied.The Young modulus E is the linear response of a material to a uniaxial stress where the response is measured in the direction of the applied stress and the Bulk modulus B is the response to an isotropic stress.The Young and bulk moduli can therefore be expressed as a function of a single unit vector in Cartesian space expressed in spherical coordinates 0 ≤ θ ≤ π and 0 ≤ φ ≤ 2π as u = (sin θ cos φ, sin θ sin φ, cos θ) [60]: where in the case of the Young modulus we have transformed the head of the compliance tensor from the Cartesian basis to a new basis whose first unit vector is u following the transformation: where a ij are the direction cosine specifying the angle between the i th axis of the new basis and the j th axis of the initial basis.The Bulk modulus is simpler because it is obtained by applying an isotropic stress (pressure p) such that ε ij = −pS ijkk .
Other elastic properties such as the shear modulus G or Poisson's ratio ν depend on the direction in which the stress is applied u but also the orthogonal direction in which the response is measured v and can be parametrized with three angles θ, φ and 0 ≤ ξ ≤ 2π: The shear modulus and Poisson's ratio can therefore be obtained as: We note that for the elastic properties studied here we only need up to two vectors (or three angles) in the new basis because the direction of applied stress and measured response are orthogonal but a general elastic property where this was not the case would require three vectors (or four angles) in the transformed basis.These elastic properties can be averaged by direct integration on the unit sphere to give the standard Young modulus, bulk modulus, shear modulus and Poisson's ratio.However, very popular averaging approximations have been developed including the Voigt approximation where the average bulk and shear moduli are given by [61]: In the Reuss approximation, the bulk and shear modulus are defined as [61]: The Voigt approximation provides an upper bound for the bulk and shear moduli, while the Reuss approximation gives a lower bound.We can therefore define the arithmetic mean, referred to as the Void-Reuss-Hill approximation [61], as We then express the effective Young E modulus and Poisson ratio ν as: where the relations apply to the Voigt, Reuss and Hill approximation of the Young, bulk and shear moduli and the Poisson's ratio.We can also define the universal elastic anisotropy as [55]: Finally, we can obtain the bulk sound velocity v B , the compressional velocity v P , shear velocity v G and the av-  erage sound velocity v av as: where ρ is the average mass density.Using the average sound velocity, the Debye temperature can be estimated within the Debye model as: where h, k B and N at are the Planck constant, Boltzmann constant and the number of atoms in the primitive cell, respectively.We studied the elastic properties of β-Ga 2 O 3 using the thermo pw code [62].The stiffness tensor of Laue class C 2h for base centered monoclinic crystals has 13 independent elastic constants written in Voigt notation as follow: and, C 66 .The stiffness matrix C ij was obtained by third-order polynomial fitting using 12 deformations with strain intervals of 0.001 to remain in the linear regime.The strains were applied along the crystal lattice vector of the β-Ga 2 O 3 primitive cell presented in Fig. 1 such that the resulting stiffness matrix is expressed in that basis.For each strain, the ions were relaxed to their equilibrium positions with a very tight convergence threshold of 4•10 −6 Ry/ Å on forces.We used a 160 Ry energy cutoff on planewaves and a 12×12×9 k-point grid.All the elastic coefficients and elastic properties are reported in Table IV.Our calculations compare well with prior theoretical work and with resonant ultrasound spectroscopy coupled with laser-Doppler interferometry [53].We computed all coefficients independently such that we can estimate off-diagonal accuracy when symmetry constrain are not precisely fulfilled.The most sensitive coefficient is the C 12 , with an accuracy of ±3.4 GPa.
Using Eqs. ( 9)-( 14), we obtained a bulk modulus of 184 GPa, a Young modulus of 207 GPa, a shear modulus of 79 GPa and a Poisson's ratio of 0.313.Those numbers agree well with recent experimental elastic constants of B = 183 GPa, E = 210 GPa, G = 80 GPa and ν = 0.31 [53].Finally, using Eqs.( 15), ( 17)-( 19) and (20) we compute the universal elastic anisotropy A U to be 0.84, the average sound velocity to be 4.01 km/s and the estimated Debye temperature Θ D to be 551 K.
Using the ELATE software [59,60], we show in Fig. 3(a) the parametrized Young modulus of Eq. (3) as a parametrized three dimensional surface and in Fig. 3(b,c) the parametrized shear modulus and Poisson's ratio of Eqs. ( 7) and ( 8) where the maximum and minimum value of the third angle is shown in blue and green, respectively.Compared to simple semiconductors where the bulk modulus is spherical, β-Ga 2 O 3 is strongly anisotropic.For example the Young modulus has a minimum value of 134 GPa in the xz plane with a unit vector (0.94, 0, 0.34) while the maximum value of the Young modulus is 293 GPa in the (0.34, 0.93, 0.13) direction.In the case of the shear modulus and the Poisson's ratio presented in Fig. 3, they are also highly anisotropic with values ranging from 50 GPa to 133 GPa for the shear modulus and from 0 to 0.67 for the Poisson's ratio which displays a flower-like shape along the diagonal axes.

V. CARRIER MOBILITY
We now analyze the intrinsic carrier transport properties of β-Ga 2 O 3 .We compute the ab-initio drift carrier mobility through the linear response ∂ E β f nk of the electronic occupation function f nk to the electric field E and where V uc is the unit cell volume, Ω BZ the first Brillouin zone volume and, n c = (1/V uc ) n (d 3 k/Ω BZ )f nk is the carrier concentration.We solve the linearized Boltzmann transport equation (BTE) [63][64][65]: with τ nk being the total scattering lifetime: Here v nk is the electronic velocity of the eigenstates ε nk , f nk is the Fermi-Dirac occupation, and n qν is the Bose-Einstein distribution function.The electron-phonon matrix elements g mnν (k, q) are the probability amplitude for scattering from an initial state nk to a final state mk + q via the emission or absorption of a phonon of frequency ω qν .A common approximation known as the self-energy relaxation time approximation (SERTA) consists in neglecting the second term on the right-hand side of Eq. ( 22).The mobility then takes the simpler form: We used the EPW software [30,66] to interpolate the electron-phonon matrix element g mnν (k, q) from a coarse 8 × 8 × 6 k-points and 4 × 4 × 3 q-points grids to a dense 160 × 160 × 120 k and q grids, as required to converge the electron mobility.The interpolation uses maximally localized Wannier function [67] and the Wannier90 software [68].We used 22 Wannier functions of initial s character centered on the gallium atoms and of p character centered on the oxygen atoms.The Dirac delta function in Eqs. ( 22) and ( 23) were computed using the adaptive smearing method of Refs.[64,69,70].
To reduce computational cost, we computed separately the electron and hole mobility by explicitly interpolating only the matrix elements for which their electronic eigenvalues at k and k + q where within 0.3 eV of the band edges.We also relied on crystal symmetries to decrease the number of k-points.In the case of the electron mobility we explicitly interpolated 13,516 k and 101,346 q-points, instead of the 3,072,000 points that would have been required by computing all the points from the 160 × 160 × 120 grid.In the case of the hole mobility, owing to very flat bands the majority of grid points contribute to the hole mobility, as can be seen in Fig. 2(a).Thus the computational cost is much higher and our densest interpolated grid is 56×56×42 points, which corresponds to 55,892 k-points and 131,712 qpoints explicitly computed.
We obtained the following room-temperature electron and hole drift mobility tensor (cm 2 /Vs) in the SERTA: (26) Interestingly, although the electron effective mass is isotropic (see Table II), we observe about 15% anisotropy for the electron mobility resulting from anisotropic electron-phonon scattering.This result is in line with the recently observed 10-15% anisotropy in the electron mobility of β-Ga 2 O 3 [16].Based on our convergence study with increasing fine grid size, we estimate an accuracy of ±3 cm 2 /Vs for the electron mobility and ±0.5 cm 2 /Vs on the hole mobility.The anisotropy of the hole mobility is within the uncertainty of the calculations.
The temperature dependence of the BTE electron and hole mobility as a function of temperature is presented in Fig. 4, slightly overestimating experimental data.The isotropic average of the electron and hole mobility are 258 cm 2 /Vs and 1.2 cm 2 /Vs, respectively.To our knowledge, this is the first time that the hole mobility of β-Ga 2 O 3 is computed from first-principles.
Our room-temperature value of the electron mobility of β-Ga 2 O 3 is slightly higher than prior theoretical studies: Ref. [19] gives 115 cm 2 /Vs at a carrier concentration of 10 17 cm −3 using Rode's method [31], and 200 cm 2 /Vs using k • p perturbation theory [20].Ref. [21] obtained an electron mobility of 155 cm 2 /Vs using the SERTA, in close agreement with our SERTA value of 167 cm 2 /Vs.
The overestimation with respect to experimental electron mobility can be traced back to the fact that our calculated electron effective mass is 7% smaller than in experiments and that electron-phonon matrix elements are dominated by Fröhlich polar scattering, which in turn scales with the dielectric constant.Our calculated dielectric constant is approximately 11% higher than in experiments.Taken together, these estimates indicate that our calculation underestimate the Fröhlich coupling by approximately 13%.In Ref. [73] we have shown that the mobility is inversely proportional to the Fröhlich coupling and effective mass, therefore we expect that the use of DFT leads to an overestimation of the mobility by approximately 24%.Experimental Hall electron mobilities of 125 cm 2 /Vs [12] and 152 cm 2 /Vs [71] were reported and are consistent with our findings.
Since lattice scattering becomes negligible at low temperature, the mobility computed using Eq. ( 21) diverges when T tends to zero.At low temperature other scattering mechanisms dominate carrier transport includ-ing defect [74] and impurity scattering [63].The impurity scattering may be included using the semi-empirical model developed by Brooks and Herring [75,76].The ionized-impurity limited mobility µ i can be evaluated analytically assuming spherical energy surfaces, negligible electron-electron interactions, and complete ionization of the impurities: where Here m * d = 0.26 m 0 and 3.39 m 0 is the density-of-state effective mass for the electron and hole, respectively, n and n i are the electron or hole densities and the density of ionized impurities, respectively, s = 4.05 0 is the average dielectric constant, 0 is the permittivity of vacuum, and h is Planck's constant.In the above expressions, the concentrations are expressed in cm −3 , and the temperature T is in K.The mobility including phonon (µ) and impurity (µ i ) scattering can be computed using the mixed-scattering formula [76] } , where X 2 = 6µ/µ i and ci(X) and si(X) are the cosine and sine integrals.The resulting combined mobility for a concentration of 10 15 cm −3 of ionized impurity is shown with a dashed line in Fig. 4, improving the agreement with experiment in the low temperature regime.
Finally, to shed light on the microscopical mechanisms driving the electron mobility in β-Ga 2 O 3 we computed the isotropic average of the momentum and mode resolved contribution to the SERTA mobility as where the mode resolved inverse mobility T qν is where w q is the weight of the q-point.
We show in Fig. 5 the mode contribution to the inverse mobility as well as the density of state inverse mobility along with the cumulative integral (dashed red line).The mode contribution spans a region close to the zone center, since as discussed above, larger momenta have negligible contribution to the mobility.The spectral decomposition is separated into three defined energy regions: low energy ( ω < 50 meV), middle energy (50 meV ≤ ω < 71 meV), and high energy ( ω ≥ 71 meV) regions.The high energy phonons alone account for 62% of the inverse mobility at room temperature followed by the low energy phonons (22%) and middle energy phonons (16%).We mention the following 10 modes, in relation with Table III, that contributes significantly to reducing the mobility: the B u (LO z 1-3,8) and B u (LO y 2-3, 5-8) modes.Interestingly, all the dominant modes have B u symmetry and are longitudinal optical modes.
As can be seen on the left-side of Fig. 5, the spectral decomposition of the mode contribution to the inverse mobility is complex, with many modes contributing to the mobility.Such complexity in the phonon spectrum of β-Ga 2 O 3 with 30 crossing and intertwined phonon branches translates into many ways for the electrons to interact with the bosonic continuum yielding increased scattering and reduced mobility.It is worth comparing such behavior of the electron scattering with a related material, wurtzite GaN that possess similar electron effective mass ≈ 0.2-0.3m e .In the nitride compound, the phonon bandstructure is composed of 12 modes clearly separated by a 20 meV gap [77].This translates into a reduced scattering with two dominant scattering at around 2 meV and 92 meV [78] and explains why the electron mobility in wurzite GaN is four times larger than in β-Ga 2 O 3 despites similar effective masses.
Figures of merit have been introduced as a way to quantify the influence of materials parameters on the performance of semiconductor devices.The most common figures of merit include the Johnson figure of merit (JFOM) which assess the quality of a semiconductor for high frequency power transistor application [79], the Keyes' figure of merit (KFOM) which quantifies the thermal limitation of transistors switching frequency [80] and the Baliga's figure of merit (BFOM) [3].In this work we focus solely on the BFOM which is used to identify material's parameters so as to minimize losses in power field effect transistors [1].The BFOM relies upon the assumption that power losses are solely due to power dissipation in the on-state by current flow through the on-resistance of the device.As a result, the BFOM is used for device operating at low frequency where the conduction losses are dominant.
The BFOM is given by where E b is the computed breakdown field, µ the computed mobility from Eq. ( 21) and 0 is the temperaturedependent experimental static dielectric function with the field perpendicular to the (100), (010) and (001) direction, respectively [81], which we reproduce in Fig. 6(c).Importantly, we stress that all the quantities entering in Eq. ( 30) are temperature-dependent.The temperature and direction-dependent mobility has already been obtained in Section V. Therefore we only need to compute the breakdown field to obtain the BFOM.Refs.[84,85] proposed the following model:   The BFOM of Si was obtained using dielectric constants from Ref. [82] and a breakdown field of 0.3 MV/cm, as well as the experimental electron mobility from Norton [83].
where E g is the bandgap of the materials in eV, ω max the phonon cutoff frequency in THz and E b the breakdown field in MV/m.Although successful, the main limitation of this model is that it is independent of temperature.For this reason, we aim at computing the BFOM from first principles while retaining the temperature-dependence.
To do so, in addition to the intrinsic carrier mobility, we need to compute the intrinsic breakdown field.
The most common theory for a material breakdown rely on electron avalanche [86] which occurs when the electron energy reaches the threshold for impact ionization.This is the energy at which an electron generates a second conduction electron by excitation across the electronic energy gap, causing electron multiplication (avalanche) and leading to a breakdown of the material [27].As a result, the threshold for impact ionization is usually taken as the electronic bandgap.The idea behind the theory relies on accelerating the conduction electron with a laser field and taking into account the electron scattering with the lattice during pumping.Indeed the phonon collision reduce the acceleration of the electron by modifying their momentum.
The von Hippel low energy criterion is more stringent and states that breakdown will occur when the rate of energy gain A(E, ε, T ) by an electron of energy ε due to the external field E at temperature T is larger than the energy-loss rate B(ε, T ) to the lattice due to electronphonon interaction [26][27][28]: for energies ε going from the conduction band minimum to the threshold for impact ionization, i.e. the bandgap of the materials.The steady-state solution for the average energy-gain rate from the electric field is [27]: where e 2 is the electron charge, m * = 0.3 [87] the electron effective mass.The energy and temperature-dependent electron-phonon lifetime is given by: where D(ε) is the density of state and τ −1 nk is given by Eq. ( 23).
The field-independent net rate of energy loss B(ε, T ) to the lattice is obtained by subtracting the rate of phonon absorption from phonon emission [27,28]: where n qν are the Bose-Einstein occupation factors in the absence of an electric field.We computed the energy gain and energy loss rates using the EPW software by interpolation on a dense 80×80×60 k-point grid and a 40×40×30 q-point grid with a constant smearing of 20 meV.In Fig. 6(a) we present the change of energy loss rate as a function of energy, starting from the conduction band minimum (CBM).On the same figure, we compare the loss rate with the average energy gain rate for increasing external electric field E at room temperature.We define the intrinsic breakdown field E b as the smallest external electric field such that the energy gain curve is larger than the energy loss curve for all energies between the CBM and the CBM plus the energy of the bandgap (4.5 eV).This value provides an estimate of the electric field range for which the material will not undergo dielectric breakdown.We compute that at room temperature the breakdown field is 5.8 MV/cm including all electron-phonon scattering processes.Using the same approach for different temperatures, we can obtain the change of breakdown field with temperature shown in Fig. 6(b).We find a breakdown field of 6.64 MV/cm at 500 K.
Such calculation was performed by Mengle and Kioupakis [22] for the intrinsic electron breakdown field at 300 K.They obtained 5.4 MV/cm by considering only the dominant LO phonon mode, and estimated that the contribution of other modes would lead to 6.8 MV/cm.We note that the experimental breakdown field in β-Ga 2 O 3 is typically reported around 8 MV/cm [1].This is in line with our calculations, since the von Hippel low energy criterion should be seen as a lower bound for the breakdown field.
Using this information and the experimental dielectric function, we can compute the temperature and directiondependent BFOM.The BFOM is typically given with respect to the BFOM of Silicon.
In this case we computed the reference BFOM of Silicon by using the temperature-dependence dielectric constant of Refs.[82,88] and breakdown field of 0.3 MV/cm as well as the experimental temperature-dependent electron mobility from Norton [83].The resulting change of BFOM is given in Fig. 6(d).The direction-averaged minimum and maximum values are 1130 and 2035, respectively.We see that even though the computed breakdown field underestimates the experiment, this effect is compensated by an overestimation of the mobility.As a result, our calculated BFOM is close to experimental estimates of 2000-3000 [1].
This cancellation suggests that the current level of theory could be sufficient to predict the BFOM of new ma-terials.

VII. CONCLUSION
In this work, we performed an in-depth study of the structural, vibrational, elastic, electrical, and transport properties of β-Ga 2 O 3 using state-of-the art firstprinciples simulation tools.We carefully analyzed the structural properties of the monoclinic variation of β-Ga 2 O 3 and analyzed the effect of spin-orbit coupling on those properties.We studied the electronic structure and carrier effective masses.We made a careful analysis of the vibrational properties including a symmetry analysis of β-Ga 2 O 3 using first-order response function theory including dielectric and Born effective charges study.We calculated many elastic properties by computing the elastic constants tensor including bulk, shear and Young modulus tensor using parametric three dimensional visualization but also Poisson's ratio, universal elastic anisotropy, sound velocities and Debye temperature and found a strong directional anisotropy.We use the Boltzmann transport equation to compute the intrinsic electron and hole drift mobility and obtain a room temperature values of 258 cm 2 /Vs and 1.2 cm 2 /Vs, respectively.We found that the mobility in β-Ga 2 O 3 was limited by a series of longitudinal optic phonons with symmetry character B u at the zone center.Finally we used the von Hippel low energy criterion to compute fully from first-principles the breakdown field which allowed us to compute the direction and temperature-dependent Baliga's figure of merit for high power device.We saw that the predicted figure of merit was in good agreement with experiment and attributed this to an overestimation of the computed mobility compensating an underestimation in the computed breakdown field.
The present analysis may serve as the basis for a general, consistent, and predictive framework to study materials for power electronics from first principles.

VIII. ACKNOWLEDGEMENT
The authors thank Paolo Giannozzi for his help with the monoclinic structure in Quantum Espresso and Andrea Dal Corso for his help with the monoclinic structure in the thermo pw code.Computer time was provided by the PRACE-15 and PRACE-17 resources MareNostrum at BSC-CNS, and the Texas Advanced Computing Center (TACC) at the University of Texas at Austin.S.P. acknowledge support from the European Unions Horizon 2020 Research and Innovation Programme, under the Marie Sk lodowska-Curie Grant Agreement SELPH2D No. 839217.F.G.'s contribution to this work was supported as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award de-sc0020129.

FIG. 2 :
FIG.2: Electronic (a) and phonon (b) bandstructure of β-Ga2O3 along high-symmetry lines of the monoclinic Brillouin Zone.For the phonons, the gray lines denotes bands with the identity point group while the red lines have a π rotation around the [0,1,0] Cartesian axis (C2 point group) and the blue line are lines with inversion symmetry (Cs(m) point group).

FIG. 3 :
FIG. 3: Spatial dependence of the (a) Young modulus, (b) shear modulus and (c) Poisson's ratio of β-Ga2O3 using the ELATE software [59] for visualization of second-order elastic constants.In (b,c) the blue and green surfaces represents the maximum and minimum of the third angle parametrization, see text.The directions x, y, and z represents the increments along the a, b and c directions of the primitive cell shown in Fig. 1.

FIG. 4 :
FIG. 4: (a) Electron and (b) hole drift mobility of β-Ga2O3 using the Boltzmann transport equation along the three principal directions where the dashed line indicate the directionaveraged drift mobility including the effect of 10 15 cm −3 ionized impurity scattering.The experimental data of Oishi et al. [71], Irmscher et al. [12] and Chikoidze et al. [72] are Hall measurement of the mobility.

FIG. 5 :
FIG.5: Direction-averaged momentum and mode-resolved contribution to the inverse electron mobility Tqν at room temperature (left) as well as spectral decomposition of the inverse mobility (right).The dashed red line represent the cumulative integral of each spectrum ∂µ −1 /∂ ω and adds up to the inverse electron mobility in the self-energy relaxation time approximation.

FIG. 6 :
FIG.6:(a) Average energy gain rate A(ε) from an applied external field, and average energy loss to the lattice B(ε).Both quantities are for 300 K.The intrinsic breakdown occurs when the applied electric field is such that the gain rate is larger than the loss rate for all energies between the conduction band minimum (CBM) and the CBM plus the energy of the bandgap (4.5 eV).(b) Variation of computed breakdown field with temperature.(c) Experimental variation of dielectric function with direction and temperature of β-Ga2O3 from Ref.[81].(d) Baliga's figure of merit (BFOM) with respect to the BFOM of Silicon.The BFOM of Si was obtained using dielectric constants from Ref.[82] and a breakdown field of 0.3 MV/cm, as well as the experimental electron mobility from Norton[83].

TABLE IV :
Comparison between our calculated elastic constants Cij, bulk modulus BH, Young modulus EH, shear modulus GH, Poisson's ratio νH, universal elastic anisotropy A U , bulk sound velocity vB, compressional velocity vP, shear velocity vG, average velocity vav, and Debye temperature ΘD, and prior theoretical and experimental work.The subscript H denotes the Void-Reuss-Hill averaging approximation.

TABLE V :
Reciprocal space coordinates of the highsymmetry point in the Brillouin zone of β-Ga2O3.Ψ =