Electronic, Thermal, and Thermoelectric Transport Properties of ε-Ga2O3 from First Principles.

The electronic, thermal, and thermoelectric transport properties of ε-Ga2O3 have been obtained from first-principles calculation. The band structure and electron effective mass tensor of ε-Ga2O3 were investigated by density functional theory. The Born effective charge and dielectric tensor were calculated by density perturbation functional theory. The thermal properties, including the heat capacity, thermal expansion coefficient, bulk modulus, and mode Grüneisen parameters, were obtained using the finite displacement method together with the quasi-harmonic approximation. The results for the relationship between the Seebeck coefficient and the temperature and carrier concentration of ε-Ga2O3 are presented according to the ab initio band energies and maximally localized Wannier function. When the carrier concentration of ε-Ga2O3 increases, the electrical conductivity increases but the Seebeck coefficient decreases. However, the figure of merit of thermoelectric application can still increase with the carrier concentration.

Introduction

Ga2O3 is an emerging wide-band-gap material for semiconductor power devices due to its high breakdown voltage.[1] Ga2O3 also exhibits transparency for the UV region, which makes it a promising material for application in optoelectronics, such as a solar blind deep-UV photodetector.[2,3] Although Ga2O3 has broad application prospects, the fundamental properties are still poorly understood due to its polymorphism.[4] Ga2O3 can form several polymorphs, denoted as α-Ga2O3, β-Ga2O3, γ-Ga2O3, δ-Ga2O3, and ε-Ga2O3,[5,6] which makes the characteristics of Ga2O3 very complex. Recently, ε-Ga2O3 thin films were successfully grown on various substrates.[7−10] It has been proved that an epitaxial ε-Ga2O3 thin film can be thermally stable up to nearly 1000 K.[11]Table shows the phase transition temperatures in a dry atmosphere of these five Ga2O3 phases, which was measured by Roy et al.[5] Theoretical calculations by Yoshioka et al. indicated that the formation energy of these phases is in the order β < ε < α < δ < γ.[12] Although β-Ga2O3 is the most stable at room temperature, ε-Ga2O3 can be stable up to 1000 K, which makes it a promising material for application in electronic devices. To date, a lot of the properties of Ga2O3 come from theory, such as density functional theory (DFT), because it is not easy to prepare high-quality pure Ga2O3 samples for measurement.[13]

Table 1

Transformation Relationship among Five Different Phases of Ga2O3[5]

phase transition	temperature
δ-Ga2O3 → ε-Ga2O3	773 K
α-Ga2O3 → β-Ga2O3	873 K
γ-Ga2O3 → β-Ga2O3	923 K
ε-Ga2O3 → β-Ga2O3	1143 K

An essential factor to consider in high-power electronic applications is the heat dissipation of the devices. Compared with other wide-band-gap semiconductors, such as GaN and GaAs, the thermal conductivity of Ga2O3 is much smaller and thus a clear weak point of Ga2O3 in terms of power device application.[14] The details of electrical transport and energy dissipation could help us develop high-performance electronic devices. Hence, it becomes crucial to understand the details of the thermal and thermoelectric properties of ε-Ga2O3. The latest studies on the transport and thermoelectric properties of β-Ga2O3 show numerical values for the Seebeck coefficient and power factor, which encourages more research on further improvements of Ga2O3. To the best of our knowledge, to date, no specific theoretical research on the thermoelectric properties of ε-Ga2O3 based on ab initio electronic structure calculations combined with Boltzmann transport equations has appeared, and only some experiments related to the electron transport and thermoelectric properties are available. In this paper, we attempt and summarize the results of the theoretical derivation of the transport and thermoelectric properties of ε-Ga2O3 using the semiclassical Boltzmann transport theory.

Results and Discussion

Electronic Properties

The optimized unit cell of ε-Ga2O3 with 40 atoms is shown in Figure , which is rendered by VESTA.[15] The optimized lattice parameters of ε-Ga2O3 are a = 5.12 Å, b = 8.79 Å, and c = 9.41 Å, which are 1.5%, 1.0%, and 1.4% larger than the experimental values by Cora et al.[10] but consistent with previous theoretical values by Yoshioka et al.[12]

Figure 1

Optimized atomic structure of ε-Ga2O3, where the large atoms are gallium and the small red atoms are oxygen.

The calculated electron band structure along with the Brillouin zone, the total density of states (DOS), and the projection density of states (PDOS) are presented in Figure . The band structure is interpolated by the Wannier function. From the band structure and DOS shown in Figure , it can be seen that the energy positions of each group of bands (representing O 2s, Ga 3d, and O 2p in the higher part of the valence band and Ga 4s, Ga 4p, and Ga 4d in the lower part of the conduction band) are rather similar to the band structures of β-Ga2O3. This is clearly reflected in the DOS. The band gap of ε-Ga2O3 is about 4.9 eV as determined by optical spectra, which is comparable to that of β-Ga2O3.[7,16] Our calculated direct band gap at the zone center is 2.11 eV, underestimating the expected experimental band gap of DFT.[17] However, the calculated band gap of ε-Ga2O3 is closest to the value of β-Ga2O3 (see Supporting Information Note II). The effective electron mass tensor is defined by[18]where EC is the dispersion of the lowest conduction band and k and k are the ith and jth elements of k, respectively. If the conduction band minimum is E0 and it locates at point Γ then the first-order approximation of EC can be expressed asIn order to obtain the effective electron mass tensor, the band structure of ε-Ga2O3 was first calculated by Quantum ESPRESSO with a coarse k mesh. Then, the band structure of ε-Ga2O3 was interpolated by Wannier90 using a maximally localized Wannier function. Considering that a good quadratic fit requires very dense sampling at different directions in k space close the center of the Brillouin region, the spacing is the reciprocal of the lattice constant of 0.002. Finally, eq was used to fit the interpolated band structure near the conduction band minimum using the Scikit-Learn python package, yielding the inverse mass tensor. Then, we can invert this tensor to obtain the mass tensor itself as followswhere m0 is the free electron mass. The electron effective mass of ε-Ga2O3 is quite isotropic with an average value of 0.217m0, which is quite close to the effective mass of β-Ga2O3 (see Supporting Information eq S1). Due to the strong Ga–O ionic bond, the energy distribution of the empty and occupied electronic states is little affected by the actual arrangement of the Ga and O atoms.

Figure 2

Dispersion relation (left) and densities of states (right) of electrons in ε-Ga2O3. Gray dashed lines indicate the Fermi levels, which are set at the center of the band gap.

Thermal Properties

Figure shows the phonon dispersion curves along with symmetry lines in the Brillouin zone and the corresponding ε-Ga2O3 phonon DOS. There are 120 phonon modes that span frequencies up to a maximum value of 21.9 THz at the Γ point. The diagonal Born effective charges are Ga1 = (2.92, 2.93, 3.06), Ga2 = (3.51, 3.29, 2.88), Ga3 = (3.35, 3.20, 3.49), Ga4 = (3.35, 3.45, 3.22), O1 = (−2.49, −2.12, −2.31), O2 = (−2.11, −2.30, −2.13), O3 = (−2.13, −2.15, −1.93), O4 = (−1.98, −2.06, −2.53), O5 = (−2.17, −2.14, −1.96), and O6 = (−2.24, −2.10, −1.76) in units of electron charge. The off-diagonal components are below 0.3, and none of them are reported. The high-frequency dielectric tensor of ε-Ga2O3 is somewhat isotropic with , , and . The dielectric tensor calculated by DPFT is consistent with the quasi-particle calculation result with ε = 4.3 solved by the Bethe–Salpeter equation of many-body perturbation theory.[4] The dielectric constant of ε-Ga2O3 is slightly greater than that of β-Ga2O3 (see the Supporting Information Note III), demonstrating the dielectric functionality of ε-Ga2O3.

Figure 3

Dispersion relation (left) and densities of states (right) of phonon in ε-Ga2O3.

The ε-Ga2O3 thermal performance prediction is one of the key factors for rational acceleration of electronic devices, including the specific heat capacity at constant volume (CV) and heat capacity at pressure (CP), Debye temperature (ΘD), speed of sound (vs), mode-resolved Grüneisen (γ(q)), average Grüneisen parameters (γ), thermal expansion coefficient (α), Debye temperature (ΘD), and lattice thermal conductivity (κ).[19] The phonon dispersion is essential to determine other thermal properties of a material such as the heat capacity. The heat capacity CV is shown as a green solid line in Figure a. The Dulong–Petit limit of 3R is plotted as a black dashed line in Figure a, where R is the gas constant. CV reaches the Dulong–Petit limit at high temperatures, analogous to the case of β-Ga2O3.[20] Furthermore, if the ab initio-calculated CV data were fitting to the Debye model, the Debye temperature ΘD of ε-Ga2O3 could be predicted to be 673 K. For comparison, the Debye temperature of β-Ga2O3 was also estimated to be 685 K in the same way (see Supporting Information Note IV and Figure S4), which was close to the experimental β-Ga2O3 Debye temperature of 738 K.[21] Therefore, the Debye temperature of ε-Ga2O3 predicted here is reasonable and is close to the Debye temperature of β-Ga2O3. The thermal expansion coefficient of ε-Ga2O3 is shown in Figure b. At room temperature, the thermal expansion coefficient is about 1.4 × 10–5 K–1, which is higher than that of β-Ga2O3 (see Supporting Information Note IV.) The bulk modulus of ε-Ga2O3 is shown in Figure c. At room temperature, the bulk modulus of ε-Ga2O3 is about 162 GPa, which is also higher than that of β-Ga2O3 (see Supporting Information Note IV).

Figure 4

Thermal properties of ε-Ga2O3. (a) Heat capacity of ε-Ga2O3. Red and blue lines represent the heat capacity at constant pressure (CP) and at constant volume (CV), respectively. (b) Thermal expansion coefficient. (c) Bulk modulus.

From the phonon dispersion near the Γ-point in Figure , it is clear that there are three branch acoustic phonons. The lower two branches (denoted as T1 and T2) belong to the transverse acoustic phonon, while the upper branch (denoted as L) belongs to the longitudinal acoustic phonon. In order to estimate the sound speed of ε-Ga2O3, the slopes of each acoustic branch are calculated near the Brillouin zone center along with three orthogonal directions by eq . The group velocities v, v, and v of acoustic phonon branch T1, T2, and L at each orthogonal direction x, y, and z are shown in Table . The group velocity at each acoustic branch is defined as the average of the squareTwo transverse sound speeds and and one longitudinal sound speed vL are 3.1, 3.7, and 7.0 km/s, respectively. The speed of sound is further estimated by the average of the transverse and longitudinal sound speed usingAfterward, the sound speed of ε-Ga2O3 is vs = 4.9 km/s, which is slightly smaller than the sound speed of β-Ga2O3 (see the Supporting Information Note IV). Besides, the sound speeds at the three orthogonal directions are 5.2, 4.6, and 4.8 km/s, respectively.

Table 2

Group Velocities of ε-Ga2O3 Acoustic Phonons

acoustic phonon	vx (km/s)	vy (km/s)	vz (km/s)	v̅ (km/s)
T1	3.0	3.4	2.8	3.1
T2	4.1	3.6	3.5	3.7
L	7.4	6.3	7.1	7.0
vs	5.2	4.6	4.8	4.9

The frequency-dependent mode Grüneisen parameters of ε-Ga2O3 by the Phonopy software package are shown in Figure . The q mesh is set to 20 × 20 × 20. The colors of the mode Grüneinen parameters are set for band indices with ascending order of phonon frequencies, which is indicated by the color bar in Figure . The Grüneisen parameters of the acoustic modes are mostly positive. A negative Grüneisen parameter indicates an increase in phonon frequency with increasing volume. An average Grüneisen parameter, γ(T), can be obtained usingwhere are the mode contributions to the heat capacity. Using eq , we found that the average Grüneisen parameter of ε-Ga2O3 is 1.4 at 300 K. A simple phenomenological expression for the lattice thermal conductivity due to phonon–phonon scattering has been given by Slack and Morelli[22,23]where κL is the lattice thermal conductivity, M̅ is the average atomic mass, n is the number of atoms per unit cell, δ3 is the volume per atom, T is the temperature, γ is the average Grüneisen parameter for the acoustic branches, and ΘD is the Debye temperature. A is a Grüneisen parameter-dependent quantity equal toif the temperature is in Kelvin, δ is in Angstroms, and the mass is in atomic units, κ is in W/(m·K). This phenomenological expression in eq has been widely used for materials with different crystal structures and with thermal conductivities extending over several orders of magnitude.[23] The lattice thermal conductivity of ε-Ga2O3 calculated by eq is shown in Figure . The thermal conductivity of ε-Ga2O3 is estimated to be 12 W/(m·K) at room temperature by eq , which is about 29% of the β-Ga2O3 thermal conductivity (see Supporting Information Note IV). There are two reasons for the lower thermal conductivity of ε-Ga2O3 than β-Ga2O3. First, the number of atoms in a unit cell of ε-Ga2O3 is two times that of β-Ga2O3. Second, the Grüneisen parameter of ε-Ga2O3 is larger than that of β-Ga2O3. The mode Grüneisen parameter is an indicator of the anharmonic properties of materials.[24] Comparing the mode Grüneisen parameters of ε-Ga2O3 (see Figure ) with those of β-Ga2O3 (see Supporting Information Figure S5), the mode Grüneisen parameters of an acoustic phonon in ε-Ga2O3 are much larger than those of β-Ga2O3. Therefore, the lower thermal conductivity of ε-Ga2O3 can be due to the anharmonic properties of an acoustic phonon of ε-Ga2O3. If we assume that all phonon relaxation times are the same, the thermal conductivity along with three orthogonal directions i can be proportional to the square of the sound speed (21)where i = x, y, z. Using the results in Table , the anisotropy of thermal conductivity can be obtained, κ = 1.12κL, κ = 0.90κL, and κ = 0.98κL. The thermal conductivity anisotropy of ε-Ga2O3 is smaller than that of ε-Ga2O3. It is worth mentioning that this phenomenological model does not take into account the effect of doping and defects. Slomski et al. already indicated that a moderate doping level of β-Ga2O3 does not change the thermal conductivity appreciably.[25] Oxygen vacancy is the most common point defect in Ga2O3.[26] Normally, the temperature-dependent thermal conductivity can be approximated by a simple power law κ(T) ∝ T–. When other phonon scattering mechanisms are considered, the slope m will differ from 1. The point defects in Ga2O3 can provide an additional term for the thermal resistance which is proportional to the temperature and unrelated to the phonon–phonon scattering.[27] Therefore, the point defects can reduce the thermal conductivity of ε-Ga2O3, which can cause the heat dissipation problem in ε-Ga2O3-based electronic devices.

Figure 5

Frequency-dependent mode Grüneisen parameters.

Figure 6

Lattice thermal conductivity as a function of temperature.

Thermoelectric Properties

A detailed understanding of the electrical transport and energy dissipation phenomena is crucial for the development of high-performance electronic materials for application. In the past, the thermal and electrical transport properties of ε-Ga2O3 have been studied in detail. Although there have been many detailed studies on the thermal and electrical transport properties of β-Ga2O3 in the past,[28,29] few systematic works have been reported on the thermoelectric coefficients of ε-Ga2O3 over a wide range of doping concentrations and temperatures. It is clear that the most challenging computational task is the determination of the band velocities v(n, k) from an ab initio perspective. To treat doping of ε-Ga2O3 in the transport calculation, one of the simplest approaches is the rigid band approximation.[30] The band structure will be assumed to remain unchanged as the Fermi level moves up and down to simulate electron doping. To fill this gap, we systematically computed the corresponding Seebeck coefficient S, electrical conductivity σ, and power factor as a function of doping concentrations at various values of the chemical potential φ.

The electron relaxation time is an important parameter to determine the electron transport properties under a constant relaxation time approximation. Similar to GaN and GaAs, a polar optical phonon (POP) plays a vital role in the electron scattering of ε-Ga2O3 because of the strong ionic Ga–O bonding.[32−34] Therefore, both acoustic phonon scattering and polar optical phonon scattering were considered in the electron relaxation time calculation. First, the electron relaxation time will be calculated by eq . The material parameters for calculation of the electron relaxation time of ε-Ga2O3 are shown in Table . The mass density of ε-Ga2O3 is based on the experimental lattice constant.[35] The deformation potential of ε-Ga2O3 is about 12.6 eV as calculated by eq , and it is a little higher than the value of β-Ga2O3. The polar phonon energy in Table is estimated from the average of all of the longitudinal optical phonon energies at all three Cartesian directions. The polar phonon energy of ε-Ga2O3 is about 47.2 meV, which is slightly smaller than the value of β-Ga2O3 (see Supporting Information Note V and Table S3). The electron relaxation times due to both acoustic phonon and polar optical phonon scattering are plotted in Figure . The polar optical phonon scattering is always larger than that of the acoustic phonon when the temperature is larger than about 100 K. In particular, when the temperature is above room temperature, the acoustic phonon scattering in ε-Ga2O3 can be neglected. The electron relaxation time at room temperature is about 15.1 fs, which is smaller than that of β-Ga2O3 (Supporting Information Note V and Figure S6). The calculated Seebeck coefficient of ε-Ga2O3 versus carrier concentration (N) and temperature (T) are shown in Figure a. The absolute value of the Seebeck coefficient is evidently decreasing with increasing carrier concentration. Furthermore, the Seebeck coefficient decreases almost linearly with increasing log(N). The electrical conductivity of ε-Ga2O3 is shown in Figure b, which increases with the doping concentration. It can be seen that for a given temperature the electric conductivity of ε-Ga2O3 has the highest value when the value of φ is close to the bottom of the conduction band minimum. The electronic thermal conductivity of ε-Ga2O3 is shown in Figure c, which increases with both the carrier concentration and temperature. The electronic thermal conductivity is much lower compared to the lattice thermal conductivity shown in Figure . The power factor of a thermoelectric material is given bywhere S is the Seebeck coefficient and σ is the electrical conductivity. The power factor involves all of the important electrical properties of the material. Figure d shows the power factor of ε-Ga2O3, which indicates that the power factor increases with increasing doping concentration.

Table 3

Material Parameters of β-Ga2O3 for Electron Relaxation Time Calculations

parameter	symbol	value
electron effective mass (me)	m*	0.217
deformation potential (eV)	D	12.6
polar phonon energy (meV)	ℏωLO	47.2
mass density (kg/m3)	ρ	5.9 × 103
static frequency dielectric constant[31]	εs	13.2
high-frequency dielectric constant	ε∞	4.3
sound velocity (m/s)	vs	4.9 × 103

Figure 7

Relaxation time of ε-Ga2O3.

Figure 8

Thermoelectric properties of ε-Ga2O3 versus carrier concentration (n) on a logarithmic scale and temperature from 300 to 900 K: (a) Seebeck coefficient, (b) electrical conductivity (on a logarithmic scale), (c) electronic thermal conductibility (on a logarithmic scale), and (d) power factor (on a logarithmic scale).

Although the Seebeck coefficient decreases as the doping concentration increases, the power factor of ε-Ga2O3 still increases with increasing doping concentration as the electrical conductivity increases more with increasing doping concentration. To date, no measured Seebeck coefficient of ε-Ga2O3 can be found in the literature. The Seebeck coefficient of β-Ga2O3 has been reported to be −(300 ± 20) μV/K with a 5.5 × 1017 cm–3 carrier concentration.[28] Our calculated Seebeck coefficient of β-Ga2O3 is 373 μV/K with a 5.5 × 1017 cm–3 carrier concentration (for more details see Supporting Information Note V and Figure S7), while the calculated Seebeck coefficient of ε-Ga2O3 is 572 μV/K with a 5.5 × 1017 cm–3 carrier concentration. Compared with β-Ga2O3, the Seebeck coefficient of ε-Ga2O3 is nearly 1.5 times larger at room temperature, when the carrier concentration of both phases of Ga2O3 is the same. However, because of its lower electrical conductivity, the power factor of ε-Ga2O3 is only 18% of β-Ga2O3 when the carrier concentration is 5.5 × 1017 cm–3 (see Supporting Information Note V). Furthermore, the electrical conductivity and the Seebeck coefficient are quite isotropic. The anisotropy of the electrical conductivity and the Seebeck coefficient at the Cartesian directions is less than 1% and 0.1%, respectively.

The validity of a material for thermoelectric applications depends on a dimensionless parameter, figure of merit, which can be described by[36]where S is the Seebeck coefficient, σ is the electrical conductivity, T is the absolute temperature, and κ is the thermal conductivity. Since κ = κL + κ and κL is much higher than κ, κ ≈ κL. Then, the figure of merit of ε-Ga2O3 can be further obtained from previous calculated results, which is shown in Figure a. When the carrier concentration is 5.5 × 1017 cm–3, the figure of merit of ε-Ga2O3 is about 6.4 × 10–4 at room temperature, which is about 36% lower than the value of β-Ga2O3 (see Supporting Information Note V). The main reasons are that (1) the Seebeck coefficient of ε-Ga2O3 is higher than that of β-Ga2O3 and the lattice thermal conductivity of ε-Ga2O3 is lower than that of β-Ga2O3 when the concentration ranges from 1 × 1016 to 1 × 1019 cm–3 at room temperature and (2) the electrical conductivity of ε-Ga2O3 increases more rapidly as the carrier concentration increases at room temperature. Therefore, the figure of merit of ε-Ga2O3 can increase by increasing the doping concentration. Figure b shows the figure of merit at the Cartesian directions, which is expressed using a function of temperature, when the concentration is about 5.5 × 1017 cm–3. Figure b indicates that the figure of merit at the y direction is always largest among the Cartesian directions. The main reason is that the electrical conductivity and the Seebeck coefficient are nearly the same at all Cartesian directions due to the isotropic band structure near the conduction band maximum. Therefore, the anisotropy of the figure of merit is determined by the anisotropy of the lattice thermal conductivity.

Figure 9

(a) Average figure of merit of ε-Ga2O3 (on logarithmic scale). (b) Figure of merit at the Cartesian directions when N = 5.5 × 1017 cm–3.

Conclusion

In conclusion, we obtained a suitable description using ab initio calculation for the electronic structure of ε-Ga2O3. The phonon dispersion and thermal properties of ε-Ga2O3 were obtained by density perturbation functional theory, the finite displacement method, and the quasi-harmonic approximation method. In the study of the ε-Ga2O3 transport properties, the band structure was first interpolated by a maximally localized Wannier function. Then, the results of the ab initio band energies were integrated with the semiclassical Boltzmann transport theory. This study shows that an appropriate description of the band structure together with phonon dispersion facilitates the study of the transport and thermoelectric properties of ε-Ga2O3. The thermoelectric coefficients of ε-Ga2O3 have been investigated systematically and provide predictive data. The electronic, thermal and thermoelectric properties of ε-Ga2O3 are compared with those of β-Ga2O3. The results reveal that (1) the effective mass, dielectric tensor, heat capacity, average Grüneisen parameter, and thermal conductivity of both ε-Ga2O3 and β-Ga2O3 can be comparable, (2) the Seebeck coefficient of ε-Ga2O3 is larger than that of β-Ga2O3, but the electrical conductivity of ε-Ga2O3 is smaller than that of β-Ga2O3, and (3) the thermoelectric figure of merit of ε-Ga2O3 increases as the carrier concentration increases. Our estimated temperature and concentration dependence of the electrical conductivity, Seebeck coefficient, and figure of merit give guidelines for the thermal management and design of ε-Ga2O3-based electronic devices.

Computational Methods

Density Functional Theory Calculation

The Quantum ESPRESSO package is used for all density generalized function theory (DFT) calculations.[37,38] A plane-wave basis set and the projector augmented-wave (PAW) method are used with the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional.[39,40] All the pseudopotentials are taken from SSSP pseudopotential library[41] and pslibrary.[42] Ga 3d10 4s2 4p1 and O 2s2 2p4 are treated as the valence states. The system energy convergence criterion is set as 1 × 10–9 eV. The plane-wave self-consistence field calculation converges with a plane-wave cutoff of 70 Ry and a 7 × 4 × 4 Brillouin zone grid.[43] Optimization of the structure is truncated after the Hellmann–Feynman force up to 3 × 10–4 eV/Å. Afterward, the details of the ε-Ga2O3 band structure are obtained by maximally localized Wannier functions (MLWFs).[44,45] The SeeK-path tool was used to define the k-point and q-point labels.[46] The dielectric constant and Born effective charge are calculated by the PHonon package using density functional perturbation theory (DFPT) with only the Γ-q-point.[37,38,47,48]

Finite Displacement Method

Under the harmonic approximation, the atoms are assumed to move around their equilibrium positions r, where l is the label of atoms in each unit cell. In the finite displacement method,[49,50] the approximation of the equation for the second-order force constant can be expressed aswhere α and β are the Cartesian indices, l and l′ are the indices of atoms in a unit cell, Δrα, is the finite displacement of atoms l at α Cartesian index, and Fβ,(Δrα,) is the force of atom l′ at β Cartesian index due to Δrα,. The atomic force can be obtained from first principles by Quantum ESPRESSO. After the second-order force constant is calculated, the dynamical matrix D can be further calculated as long as the phonon dispersion ω(q). All of the harmonic phonon properties are obtained using Phonopy packages.[51] After we obtain the full relaxation structure of ε-Ga2O3, we generate a series of 2 × 2 × 1 ε-Ga2O3 supercells with atom displacement, resulting in a total 60 different supercells. The displacement length of each atom from its equilibrium position is 0.01 Å. All of the supercells contain in total 160 atoms. Then, the ground state energies of these 60 supercells are calculated by Quantum ESPRESSO. After the energies of these supercells are obtained, the force constant can be calculated by Phonopy. From the force constant we could further calculate the phonon band structure and the phonon density of states. The phonon dispersion is further interpolated with nonanalytical term correction.[47,52−54] Once the phonon band structure is obtained, the phonon mode contribution to the harmonic phonon energy can be calculated aswhere j is the phonon band index and ω is the phonon frequency. Then, we can express the total phonon energy as followsOn the other hand, the Helmholtz free energy can be expressed asAfterward, the constant volume heat capacity CV can be calculated aswhere is mode contributions to the heat capacity, which is defined byThereby, the heat capacity at constant volume CV can be obtained directly as the second derivative of the Helmholtz free energy with respect to temperature. According to the Debye model, the heat capacity at constant volume CV(T) as function of temperature can be approximated byIf CV(T) is calculated by eq with a temperature step of 10 K, which is used to fit eq by the Python package Lmfit, the Debye temperature ΘD can be obtained. Besides, when the phonon dispersion is acquired, the group velocity vg can be calculated by the finite difference method as followswhere q is the phonon vector and j is the phonon band index. Because the acoustic phonon dispersion at the Γ point is divergent, the group velocity of an acoustic phonon at the Γ point will be approximated by the points slightly away from the center in the different directions, i.e., v, v, and v, where x, y, and z are the Cartesian directions. The sound speed will be further approximated by the average of the square root of the group velocity along the three orthogonal directions.

Quasi-Harmonic Approximation

The thermal properties at constant pressure are further obtained by a quasi-harmonic approximation,[55] in which we need to calculate the Gibbs free energy, which is defined as followswhere V is the volume, p is the pressure, Uelectron(V) is the total energy of the electronic structure with different volumes, and Hphonon(T, V) is the Helmholtz free energy with different temperatures and different volumes. In this work, 9 volume points are used with a temperature step of 10 K, and the lattice constant step is 0.025 Å. After we calculate the Helmholtz free energy with different temperatures and volumes using eq , we could further obtain the Gibbs free energy at the temperature and the respective equilibrium volume V(T). The volume thermal expansion coefficient can be acquired from the equilibrium volume asThe heat capacity at constant pressure Cp is given byFurthermore, the equilibrium volume V(T) is used to fit the Vinet equation of state, which is described by[56]where B0 is the isothermal bulk modulus and B0′ is the derivative of bulk modulus with respect to pressure. If the equilibrium volume at zero pressure is V0, η is cube root of the specific volume, which is defined asIf the equilibrium volume V(T) obtained from eq is used to fit the equation of state in eq , the bulk modulus can be obtained. The mode Grüneisen parameter γ(qj) at the wave vector q and band index j is defined as[51]where V is the volume and ω(q) is the phonon frequency. It is easy to calculate γ(q) from the volume derivative of the dynamical matrix once the dynamical matrix D(q) is acquiredwhere e(q) is the eigenvector. The phonon dispersion is calculated at three different volumes to determine the Grüneisen parameters, one at the equilibrium volume and the other two at slightly distorted volumes (±0.5% lattice constant).

Boltzmann Transport Theory

The calculations of the thermoelectric and electronic transport are performed by Boltzmann transport theory within the constant relaxation time (τ) approximation implemented in the BoltzWann code.[44,45] The acoustic phonon scattering rate as a function of T is given by[34,57,58]where m* is the effective electron mass, kB is the Boltzmann constant, T is the temperature, ρ is the mass density, vs is the sound velocity, D is the deformation potential, and E is the electron energy. The deformation potential is calculated by[59,60]where E0 is the conduction band minimum located at Γ and V0 is the equilibrium unit cell volume. To obtain the absolute deformation potential of the conduction band minimum, the energy change is calculated as the difference between the conduction band minimum and a core level, e.g., the anion 2s core level, using a constant of 0.1% to expand and compress the lattice constant because the absolute position of the energy level is not well established in infinitely periodic crystals.[61] DFT calculation will be performed with several volumes by expanding and compressing the lattice constant by a step of 0.5%. The results of the energy and volume change will be fit to eq by the Lmfit python package in order to obtain a better estimation of the deformation potential. On the other hand, the polar optical phonon scattering rate can be estimated by[32]where ωLO is the longitudinal optical phonon frequency, ε0 is the vacuum permittivity, ε is the high-frequency dielectric constant, εs is the static dielectric constant, and n is Bose–Einstein distribution function for longitudinal optical phonons given byWe acquire the net relaxation rate using Matthiessen’s ruleThe average electron relaxation time is given bywhere f(E, T) is the Fermi–Dirac distributionThe average electron relaxation time will be finally evaluated with numeric integration by SciPy.

On the basis of the Boltzmann transport equation, the following expressions are used to calculate the electrical conductivity σ, electronic thermal conductivity κ, and Seebeck coefficient S as a function of the chemical potential φ and of the temperature THere, i and j are Cartesian indices, denotes the matrix product of the two tensors σ and S, and ∂f/∂E is the derivative of the Fermi–Dirac distribution function with respect to the energy. Moreover, the transport distribution function Σ(E) is defined aswhere the summation is over all bands n and over all of the Brillouin zone, E is the energy for band n at k, τ is the scattering time, and v is the ith component of the band velocity at (n, k) which can be computed byWithin the relaxation time approximation, τ is held constant with respect to the electron on band n at wave vector k; therefore, the Seebeck coefficient is independent of τ. This constant relaxation time approximation is based on the assumption that the variation of energy on the scale of kBT does not cause the electron scattering time to vary with it. This approximation is widely adopted in first-principles calculation for bulk materials.[30,36] BoltzWann code uses a maximally localized Wannier function (MLWF) set to interpolate the band structure obtained from first-principles calculations by Quantum ESPRESSO. First, a 8 × 4 × 4 k-points grid is used for construction of 56 MLWFs around the gap region. Then, a 80 × 40 × 40 k mesh is utilized to calculate the transport properties. The band structure of the Wannier function interpolation matches well with the first-principles calculations of Quantum Espresso. In order to verify the computational methods, we also calculate the electronic, thermal, and thermoelectric properties of β-Ga2O3, which agree well with previous studies (see the Supporting Information).