1. Generation of Interface Structures

The electronic structure calculation is performed by using the DFT simulation package, Vienna ab initio simulation package (VASP) ⁽⁴¹⁾. The underestimated bandgap is a major issue for the conventional DFT calculation of semiconductor materials ^(42-47). Hybrid exchange-correlation (XC) functionals or GW approximation methods can generate accurate bandgap values for most semi- conductor materials47-51. But these methods require higher computational cost than conventional local density approximation (LDA) or generalized gradient approximation (GGA). Therefore, they cannot be easily applied to large systems like the interface structure, which typically contains more than hundreds of atoms. For example, the interface structure considered in this work contains 280 atoms. For the electronic band calculation, a dense k-points grid is typically required. When applied to this problem, a hybrid functional results in too heavy computational burden. The TB-mBJ functional may be a good compromise to solve the band gap underestimation problem with reasonable computational resources. So, we use the Tran-Blaha modified Becke-Johnson (TB-mBJ) Meta-GGA functional whose results are well matched with the experimental bandgap results ^{(44, 52-55)}. Especially the non-regular TB-mBJ method has been used for accurate bandgap matching for wide bandgap insulators ^(56,57). Some comparison studies are reported for band structures and density of states calculated with GGA, hybrid functionals, and GW method. The results the TB-mBJ functional recreates electronic properties like other methods ^{(50, 53, 57-60)}. In this study, $Al_{2}O_{3}$ and $Ga_{2}O_{3}$ crystal are tested whether TB-mBJ functional is suitable for this simulation. Electronic band structure and PDOS is evaluated, the offset fluctuation stability is confirmed to the bandgap change. In particular, the previous study for non-regular TB-mBJ functional on $\gamma - Al_{2}O_{3}$ confirm the validity of functional ⁽⁵⁶⁾. Table 1 shows the bandgap values of the regular and non-regular TB-mBJ functionals. The conventional GGA-Perdew-Burke-Ernzerhof (PBE) functional underestimates the bandgap considerably. For the accurate bandgap estimation of both $Ga_{2}O_{3}$ and $Al_{2}O_{3}$, we use the non-regular TB-mBJ functional for the interface structure.

The TB-mBJ potential was proposed as ^(44,52)

where $V^{BR}\left(r\right)$ is the Beck-Roussel (BR) exchange potential, \textit{ $\rho \left(r\right)=\sum _{i}^{N}\left| \psi _{i}\right| ^{2}$ } is the electron density, and $t\left(r\right)$ is the Kohn-Sham (KS) kinetic energy density ⁽⁶¹⁾. The coupling parameter $\textit{c}$ between the BR exchange potential and the electron density-related term is given by ^(44,52),

The coupling parameter c is obtained by the self-consistent process during the DFT calculation. It is noted that $\alpha $ = ${-}$0.012 and $\beta $ = 1.023 are usually used for the conventional TB-mBJ simulation ^{(44, 52, 54)}.

In many cases, the TB-mBJ functional can give bandgaps similar to the experimental values ^(44,52). However, as shown in Table 2, the calculated bandgap of $\gamma - Al_{2}O_{3}$ is still lower than the experimental result. In the non-regular TB-mBJ functional, the coupling parameter $c$ is adjusted to make the bandgap similar to the experimental value. In this study, we adopt the coupling parameter $c$ as follows: $c_{Al}$ = 1.82, $c_{Ga}$ = 1.42, and $c_{O}$ = 1.42. By using these values, the bandgap obtained by the DFT calculation and the experimental bandgap agree very well (see Table 2).

To generate the $am- Al_{2}O_{3}$ structures, we conduct two-step MD (molecular dynamics) sim- ulation. The first step is the reactive force-field (ReaxFF) MD simulation ^(63,64). The canonical ensemble is applied until the temperature decreases from 900 K to 300 K for 25 ps with 5 fs time step. In the second step, the total energy is minimized by the DFT-MD simulation ⁽⁶⁵⁾.

Fig. 2. (a) The crystalline structures of both $\gamma - Al_{2}O_{3}$(110) and $\beta - Ga_{2}O_{3},$ (b) The optimized $\gamma - Al_{2}O_{3}$(110)/ $~ \beta - Ga_{2}O_{3}$(010) interface structures by the MD simulation, (c) the optimized crystalline interface structures from two different projection vectors, [100] and [110], (d) The optimized $\mathrm{am}- Al_{2}O_{3}$/$\beta - Ga_{2}O_{3}$(-201) interface structures by the MD simulation.

The first principles calculation is performed in the microcanonical ensemble, using periodic boundary conditions and constant cell parameters. For the DFT calculation, the GGA-PBE XC functional ⁽⁶⁶⁾ is used for the geometry optimization. The energy equilibrium structure is generated during 9 ps at T = 400 K. For the optimization of the interface structures, a single DFT-MD simulation is done. The DFT-MD simulation conditions are same as the $am- Al_{2}O_{3}$ generation process. The generated interface structures satisfy the structural minimization criteria, 0.01 eV/Å.

In order to verify the generated amorphous model, we have calculated the radial distribution function (RDF). The RDF of the generated am-Al2O3 model is shown in Fig. 1 ⁽⁶⁷⁾. The results indicate that the Al-O bond length is 1.82 Å. This result is similar to the MD simulation result of 1.8 Å ⁽⁶⁷⁾ and an experimental value of 1.8 ${\pm}$ 0.21 Å ⁽⁶²⁾. For the bond length larger than 3 Å, the distribution of the RDF function differs from the experimental result ⁽⁶²⁾. It suggests that the amorphous structure generation process through MD simulation is not perfectly matched with the actual experiment ⁽⁶²⁾. Since the most significant bond length is matched, we have adopted the MD generated structure without further modification.

Fig. 3. Schematic drawings of oxygen atoms on $\gamma - Al_{2}O_{3}$(110) and $\beta - Ga_{2}O_{3}$(010). In the schematic diagram (a) and (b), the red circles indicate the positions of the oxygen atoms when viewed from each interface direction. This corresponds to the position of the oxygen atom indicated by the purple circle in the bottom atomistic structures (c) and (d).

Fig. 4. Electronic band structure and PDOS of the $\gamma - Al_{2}O_{3}$(110)/$\beta - Ga_{2}O_{3}$(010) crystalline interface.

2. Interface Analysis and Discussion

Fig. 2 shows the atomistic structures of crystalline $\gamma - Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ and $\mathrm{am}- Al_{2}O_{3}/$ $\beta - Ga_{2}O_{3}\left(- 201\right)$ interfaces. In Fig. 2(a), the two surfaces in the crystalline structure are not perfectly matched, so that the combined interface can be obtained through the structural optimization. The optimized $\gamma - Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ and $\mathrm{am}- Al_{2}O_{3}/$ $\beta - Ga_{2}O_{3}\left(- 201\right)$ interfaces are shown in Fig. 2(b) and Fig. 2(c), respectively. The crystalline interface model consist of 10-layer $\beta - Ga_{2}O_{3}\left(010\right)$ primitive slab cell and 12-layer $\gamma - Al_{2}O_{3}\left(110\right)$ primitive slab cell. A 15 Å-long vacuum region is added to both sides. The overall slab supercell is the monoclinic system with dimensions of a = 12.25 Å, b = 11.62 Å, c = 50.32 Å, and ${\gamma}$ = 76.37◦.

Fig. 3(a) and (b) show the schematic diagrams of the oxygen atoms on the $\gamma - Al_{2}O_{3}\left(110\right)$ and $\beta - Ga_{2}O_{3}\left(010\right)$ planes, respectively. Fig. 3(c) and (d) show actual atomistic structures. The oxygen atoms in the purple circles correspond to the oxygen positions in the schematic diagram. Double-headed green and blue arrows in Fig. 3(a) and (b) indicate the in-plane distances between two neighboring atomic rows. And $d_{a}^{*}$ and $d_{c}$ are their average lengths along $a^{*}$ and $c$ directions of $\beta - Ga_{2}O_{3}$, respectively. Since the mismatch of in-plane cell parameters is small, a crystalline interface can be constructed. The average mismatch of lattice parameters is 2.2%. As shown in Fig. 3, the model is created by matching the positions of oxygen atoms. The lattice constant obtained from the XPS experimental data of the $\gamma - Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ interface structure is used for the interface model ^(25,26). The interfaces are optimized by DFT-MD simulations.

Fig. 4 depicts the electronic band structure of the $\gamma - Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ interface structure and the projected density-of-states (PDOS) for each atom type. A 7 ${\times}$ 7 ${\times}$ 1 $\Gamma$-centered k-points mesh is used for the electronic structure calculation. The cut-off energy is 600 eV. The relaxation convergence criterion is 0.1 eV/Å. Projector-augmented wave method (PAW) datasets are used in simulation. Ga(s2p1), Al(s2p1), and O(s2p4) valance electrons are applied for the simulation68. The lowest conduction band is similar to that of the bulk $\beta - Ga_{2}O_{3}$. The valance band is almost flat. It is in good agreement with the preliminary experiments and calculation data of the valance band structure of $\beta - Ga_{2}O_{3}$ ^{(5, 69-71)}. The electronic band structure is obtained from the crystalline $\gamma - Al_{2}O_{3}/\beta - Ga_{2}O_{3}$ interface structure. The in- plane band path shown in Fig. 4 is obtained from the monoclinic supercell of the interface structure. The indirect bandgap along $\Gamma$ to C ${-}$ X path is 4.64 eV. The bandgap of $\beta - Ga_{2}O_{3}$ in the interface structure has a slightly smaller bandgap than the bulk case (4.76 eV). In the conduction band, only the $Ga_{2}O_{3}$ side contributes states near the conduction band minimum (CBM).

Fig. 5. Valence band maximum (VBM) and CBM levels of the interface structures (a) Band offset result of $\gamma - Al_{2}O_{3}$(110)/$\beta - Ga_{2}O_{3}$(010) interface and band offset result of am $\gamma - Al_{2}O_{3}$/$\beta - Ga_{2}O_{3}$(-201) interface, (b) calculated offset [31], (c) experimental result of $\gamma - Al_{2}O_{3}/\beta - Ga_{2}O_{3}$ [26], (d) experimental result of $\mathrm{am}- Al_{2}O_{3}$/$\beta - Ga_{2}O_{3}$ [27].

Fig. 5 shows schematic diagrams of band offsets of the $Al_{2}O_{3}/\beta - Ga_{2}O_{3}$ interfaces. The band offset characteristic is evaluated from the interface structure. The CBM and VBM at each side are obtained at a mid-point of each $Al_{2}O_{3}$ or $Ga_{2}O_{3}$ layer along the c-axis, whose direction is perpendicular to the interface. Our results of the DFT calculations in Fig. 5(a) show that the offset is slightly larger than the experimental values ^(26-28). For the $\gamma - Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ interface shown in Fig. 5(a), the band offset of the conduction band side $\Delta E_{C}$ is 2.24 eV and the one of the valence band side $\Delta E_{V}$ is 1.13 eV. The band offset of the $\mathrm{am}- Al_{2}O_{3}/\beta - Ga_{2}O_{3}\left(- 201\right)$ interface (see Fig. 2(c)) is also shown in Fig. 5(a). The band off- set of the conduction band side $\Delta E_{C}$ is 1.83 eV and the one of the valence band side $\Delta E_{V}$ is 1.07 eV. From our numerical experience, it has been found that the coupling parameter strongly changes the band gap. However, the band offsets are less sensitively affected by the coupling parameter used in the TB-mBJ functional.

There are some previous studies on the band offset based on the DFT calculation ^(30,31). In those works, the$\theta \left(\textit{monoclinic}\right)- Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ inter- face has been considered. The offset results in31 are presented in Fig. 5(b). According to the results, $\Delta E_{C}$ is 2.30 eV. It is similar to $\Delta E_{C}$ of our $\gamma - Al_{2}O_{3}\left(110\right)/$ $\beta - Ga_{2}O_{3}\left(010\right),$ but $\Delta E_{V}$ has a negative value of ${-}$0.37 eV. Conversely, the offset results of the ALD deposited interfaces which are shown in Fig. 5(c) and (d) have same tendency with our works.

1. Calculation Method

The defect analysis is performed based on the first principles calculation by using the VASP. Although the non-regular TB-mBJ method is used to calculate the band structure of the inter- face models, the validity of the non-regular TB-mBJ method for the defect analysis is not verified. Therefore, the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional, whose results are in good agreement with the experimental ones ^(72-78), is employed in the defect analysis. A mixing parameter of HSE functional ${\alpha}$ is adjusted as 0.32 to match the experimental gap of $\beta - Ga_{2}O_{3}.$ To reduce computational cost, a sparse k-points grid (2 ${\times}$ 2 ${\times}$ 2 and 2 ${\times}$ 2 ${\times}$ 1) is used for the HSE functional. Accuracy of the evaluated formation energy with the sparse k-points grid is confirmed ^(79-81).

The formation energy ($E^{f}$) and the transition level ($\varepsilon \left(q/q'\right)$) are key quantities for determining the electric characteristics of defects in solid materials. These quantities are defined as follows ⁽⁴⁵⁾. The formation energy $E^{f}\left[X^{q}\right]$ of a defect $X$ at a charge state $q$ is evaluated from inter-supercell calculations. It is defined as ⁽⁴⁵⁾,

where $E_{t}\left[X^{q}\right]$ is the total energy derived from a supercell containing the defect $X$ and $E_{t}\left[bulk\right]~ $is the total energy of the defect-free crystalline supercell. The integer $n_{i}$ is the number of atoms for each defect type, and $\mu _{i}$ is the chemical potential. $E_{F}$ is the Fermi energy level $E_{corr}$ is a correction term for the charge interaction between neighboring periodic supercells. The Freysoldt- Neugebauer-Walle (FNV) finite-size supercell correction scheme is applied for the formation energy evaluation ⁽⁷⁹⁾. It has been reported that the dielectric constant affects the calculated results significantly ^{(34, 35, 37)}. The dielectric constant $\varepsilon _{0}$=11.7 is used for the finite-size supercell correction in the bulk system.

Fig. 6. (a) Primitive unit cell of the monoclinic $\beta - Ga_{2}O_{3}$. The labels represent inequivalent sites for the bulk defect analysis, (b) Supercell structure with gallium (Ga) vacancy defect, (c) supercell structure with oxygen (O) vacancy defect. The gray circles in (b) and (c) are vacancy positions.

The thermodynamic transition level $\varepsilon \left(q_{1}/q_{2}\right)$ is defined as the Fermi-level position for which the formation energies of charge states $q_{1}$ and $q_{2}$ become identical ⁽⁴⁵⁾,

The transition levels within the bandgap are considered as localized charge traps.

2. Bulk Defects

We first calculate the defect energy levels in the $\beta - Ga_{2}O_{3}$ bulk supercell model to verify the consistency of the calculation with preceding works. For the bulk $\beta - Ga_{2}O_{3}$ structure with a defect, a 120 (2 ${\times}$ 3 ${\times}$ 1)-atoms supercell is chosen, and a 2 ${\times}$ 2 ${\times}$ 2 $\Gamma$-centered k-points grid is applied for the calculation. The primitive unit cell of the monoclinic $\beta - Ga_{2}O_{3}$ is shown in Fig. 6(a). Fig. 6(b) and (c) show the structures with Ga- and O-vacancy defects, respectively. All inequivalent sites shown in Fig. 6(a) are considered as the defect positions.

Fig. 7. Formation energies versus the Fermi energy level for each types of vacancies. Breakpoints of each line mean the charged transition level ($\varepsilon \left(q_{1}/q_{2}\right)$) according to the charge state of each supercell with a defect. $q_{1}$ and $q_{2}$ are different charged states of the calculated system (a) Bulk $\beta - Ga_{2}O_{3}$, (b) $\gamma - Al_{2}O_{3}$/$\beta - Ga_{2}O_{3}$ interface calculation results are presented.

Unlike highly symmetric crystal structures such as Si or GaAs, the monoclinic structure (C2/m) of $\beta - Ga_{2}O_{3}$ has low symmetry. There are two inequivalent oxygen positions and three inequivalent gallium positions in a primitive unit cell. For gallium types, there are two different lattice sites.$~ Ga\left(1\right)$ is tetrahedrally coordinated and $Ga\left(2\right)$ is octahedrally coordinated by neighbor oxygen atoms.$~ O\left(1\right)$ and $O\left(2\right)$ sites are three-fold coordinated and $O\left(3\right)$ is four-fold coordinated.

The intermediate pressure condition is applied for the formation energy evaluation ^(34,35). In case of intrinsic defects, the intermediate chemical potential of gallium is $\mathrm{\mu}_{Ga}$ = $\mu _{G{a_{Bulk}}}$+ $\frac{1}{5}\Delta H_{f}^{Ga_{2}O_{3}~ }$and the intermediate chemical potential of oxygen is $\mathrm{\mu}_{O}$ = $\frac{1}{2}\mu _{{O_{2}}}$+ $\frac{1}{5}\Delta H_{f}^{Ga_{2}O_{3}~ }$. For the bulk defect calculation, we use the extended supercell model. A (1 ${\times}$ 3 ${\times}$ 2) expanded supercell with 120 atoms is used for the bulk defect calculation to minimize errors due to the finite-size supercell.

The calculated formation energies for $V_{O}$(Oxygen vacancy) and $V_{Ga}$(Gallium vacancy) in the $Ga_{2}O_{3}$ supercell are shown in Fig. 7(a). The defect formation energies of the bulk supercell are plotted as functions of the Fermi energy. Each curve is obtained by connecting the minimum value of the formation energy line of each charged state obtained from Eq.

The three lowest formation curves in Fig. 7(a) depict the lowest formation energies of $V_{O}$ in $Ga_{2}O_{3}$. The breakpoints of the curves are $\varepsilon \left(2+/0\right)$ transition levels of $V_{O}$(1${-}$3). The number of each curve corresponds to the inequivalent oxygen positions. The transition levels of different oxygen positions are $\varepsilon \left(2+/0\right)$ = 3.12 eV for $O\left(1\right)$, 2.30 eV for$~ O\left(2\right)$, and 3.21 eV for $O\left(3\right)$, respectively. The V + charged states and other positively charged states are not stable for any value of the Fermi level. The energy level of $V_{O}\left(3\right)$ is the closest to the conduction band minimum(CBM), but the energy difference is larger than 1 eV.

The energy difference between the CBM and the transition level of $V_{Ga}\left(1\right)\left(2- /3- \right)$ is 0.51 eV, and it acts as a shallow donor trap in ${\beta}$ ${-}$ Ga2O3. The similar values have been reported in several experimental results ^(82-84). In the case of $V_{Ga}\left(2\right)\left(2- /3- \right)$ = 3.74 eV is the nearest transition level to the CBM. Since the energy difference between the CBM and the transition level is larger than 1 eV, it is difficult to operate as a donor trap. Our calculated formation energies of intrinsic defects in the bulk $\beta - Ga_{2}O_{3}$ are similar to the values of previous studies using HSE functionals ^(33-35).

3. Interface Defects

In case of the interface defect calculation, the MD-generated $\gamma - Al_{2}O_{3}\left(110\right)/\beta - Ga_{2}O_{3}\left(010\right)$ interface model is used. A slab structure with 280 atoms is sampled by 2 ${\times}$ 2 ${\times}$ 1 $\Gamma$-centered k-points grid. The in-plane average dielectric constant $\varepsilon _{0}$ = 10.7 is used for the finite-size supercell correction. The value is obtained from the DFT calculation of the interface structure. The native defects(vacancies) are used for the calculation. The relaxation process for structures is performed with the PBE functional. The convergence criteria in the interface structure generation are adopted without modification. Fig. 8 shows actual interface structures with vacancy defects used in defect calculations. Defect calculations are performed with the most symmetric defect sites on the atomistic $\gamma - Al_{2}O_{3}/\beta - Ga_{2}O_{3}$ interface structure. Especially, in the case of an oxygen vacancy($V_{O}$), three different oxygen sites shown in Fig. 8(b) are evaluated. The Ga vacancy ($V_{Ga}$) and Al vacancy ($V_{Al}$) type defects are formed at the interface site.

Fig. 8. $\gamma - Al_{2}O_{3}$/$\beta - Ga_{2}O_{3}$ interface structures for the defect calculation (a) Defect-free interface, (b) interface with an oxygen vacancy ($V_{O}$), (c) interface with a gallium vacancy ($V_{Ga}$), (d) interface with an aluminum vacancy ($V_{Al}$).

Fig. 9. Transition levels of different charged states and defects in the bulk $\beta - Ga_{2}O_{3}$ and $Al_{2}O_{3}$/$Ga_{2}O_{3}$ interface structure. The left section illustrates transition levels of intrinsic defects from the bulk $\beta - Ga_{2}O_{3}$ supercell and the right hand section shows the levels from the $Al_{2}O_{3}$/$Ga_{2}O_{3}$ interface structure.

Fig. 7(b) shows the formation energies from the interface structures with defects. Regardless of the types of defect, the formation energy of a neutral state is formed between 4 and 6 eV. This result is different from the bulk system in which the formation energy of $V_{O}$ is more stable. In the case of the interface calculation, the negatively charged states of $V_{O}$ are observed. The $\varepsilon \left(0/1- \right)$ charged states are the closest to the CBM for all VO defects in the interface structure. The corresponding transition levels of the different oxygen positions are $\varepsilon \left(0/1- \right)$ = 3.71 eV for $V_{O}$ (interface), 3.73 eV for $V_{O}$ (Al side), and 3.93 eV for $V_{O}$ (Ga side) site, respectively. As in the bulk supercell, $V_{Ga}$ contributes an n-type charge trap. The nearest charged state is $\varepsilon \left(2- /3- \right)$ = 4.24 eV. For the $V_{Al}$ defect, $\varepsilon \left(1- /2- \right)$ = 4.01 eV level is observed. Overall, all vacancy defects formed at the interface structure form a shallow trap level for the CBM ({\textless} 1 eV).

Fig. 9 illustrates the charged transition levels from the bulk and interface systems. All transition levels are sorted based on the Fermi energy level of each system. As mentioned above, the vacancy defects in the interface form shallow traps near the CBM level. In the case of $V_{O}$ defects, only $\varepsilon \left(2+/0\right)$ charged deep acceptor traps are observed in bulk $\beta - Ga_{2}O_{3}$. However, for the oxygen vacancies in the $Al_{2}O_{3}/Ga_{2}O_{3}$ interface structures, donor-like $\varepsilon \left(0/- \right)$ transition levels close to the CBM are generated. In particular, we could observe many defect levels close to the CBM in case of the interface structure with a vacancy. These shallow trap levels are similar to the levels observed from $\beta - Ga_{2}O_{3}$ MOS device measurement ^{(25, 85, 86)}