For the case of n-type doping of β-Ga2O3 by group 14 dopants (C, Si, Ge, Sn), a defect phase diagram is constructed from defect equilibria calculated over a range of temperatures (T), O partial pressures (pO2), and dopant concentrations. The underlying defect levels and formation energies are determined from first-principles supercell calculations with GW bandgap corrections. Only Si is found to be a truly shallow donor, C is a deep DX-like (lattice relaxed donor) center, and Ge and Sn have defect levels close to the conduction band minimum. The thermodynamic modeling includes the effect of association of dopant-defect pairs and complexes, which causes the net doping to decline when exceeding a certain optimal dopant concentration. The optimal doping levels are surprisingly low, between about 0.01% and 1% of cation substitution, depending on the (T, pO2) conditions. Considering further the stability constraints due to sublimation of molecular Ga2O, specific predictions of optimized pO2 and Si dopant concentrations are given. The incomplete passivation of dopant-defect complexes in β-Ga2O3 suggests a design rule for metastable doping above the solubility limit.

The control of the carrier concentration is essential for virtually all semiconductor applications. Ga2O3 is an emerging ultrawide-gap (4.9 eV1) semiconductor with a range of different potential applications requiring low (solar-blind photodetectors2), intermediate (field effect transistors3), and high (ultraviolet-transparent conducting oxides1) electron densities. Controlling the doping requires incorporating dopants with a sufficiently small ionization energy and simultaneously suppressing the formation of compensating defects.4 First principles calculations of the defect formation energies have elucidated the defect physics of Ga2O3, highlighting the non-conductive nature of the O vacancy5 and the importance of cation vacancies as compensating defects6 and showing that the interstitial defects play a minor role.7 Thermodynamically, the balance between dopants and compensating defects is determined by the defect equilibrium,8 which can be solved using the calculated formation energies.7,9,10 Thus, such predictions can aid the material design for specific applications.

Considering the group 14 dopants C, Si, Ge, and Sn, we construct here a defect phase diagram for Ga2O3, i.e., a quantitative prediction of the net doping (difference NDNA between donor and acceptor type defects), as a function of temperature T, oxygen partial pressure pO2, and concentration c of extrinsic dopant atoms. Further considering the thermodynamic limits to pO2 and T due to the volatility of the molecular Ga2O species,11 specific guidance for the optimization of the net doping is provided. The thermodynamic model for the defect equilibrium includes the effect of dopant-defect pair association within the law of mass action.12 This mechanism is rarely included in first principles defect equilibria but can be essential in determining the doping limits.13 In the simpler model based on isolated dopants and defects, the doping always increases with the dopant density, although it eventually saturates due to compensation.9 In contrast, it is shown here that the dopant-defect pair association, i.e., dopant passivation, causes the net doping to decrease beyond an optimal dopant concentration. Thus, approaching the optimum in the defect phase diagram requires careful adjustment of all three variables (T, pO2, and c).

The first question to address is which of the considered group 14 dopants C, Si, Ge, and Sn create a suitable shallow donor level, rather than a deep state that could pin the Fermi level below the conduction band edge. This problem is particularly acute in Ga2O3 because larger bandgaps tend to favor deep levels. The possibility of a close competition between the deep defect level and the conduction band minimum (CBM) requires a careful choice of the computational approach. With appropriate corrections, density functional theory (DFT) at the level of the standard generalized gradient correction (GGA)14 can quite accurately predict the formation energies of fully ionized defects.15 However, the prediction of deep defect levels close to the band edge requires self-consistent bandgap corrected calculations.16 Here, we use the approximately bandgap corrected hybrid functional of Ref. 17 to calculate the electronic wavefunctions and then use a defect GW approach18,19 to predict the charge transition levels. Further computational details are given below.

In order to characterize the dopant behavior, we are considering the transition between the positively and negatively charged dopants (e.g., SnGa+ and SnGa). Similar as in the case of group 15 dopants in SnO2,20 the existence of a negative charge state follows from the possible multivalency of the group 14 dopants, e.g., Sn(iv) and Sn(ii) for the positive and negative charge states, respectively. The ε(+/−) charge transition level is associated with large lattice relaxations and acts as a pinning level for the Fermi energy. If it lies below the CBM, it will limit the equilibrium free carrier concentration at the level corresponding to this Fermi energy. (However, depending on the energy landscape, a larger non-equilibrium carrier density may be observed due to persistent photoconductivity21). If it lies above the CBM, then the negative charge state is metastable and becomes relevant only in the case of degenerate doping when the Fermi level approaches the ε(+/−) pinning level due to the Burstein Moss shift. In the case of non-degenerate doping, the donor can form a charge-neutral hydrogenic state by binding a conduction band electron with an extended effective mass-like wavefunction.

The monoclinic crystal structure of β-Ga2O3 has two non-equivalent Ga sites, an approximately tetrahedrally (Td) coordinated Ga1 and an approximately octahedrally (Oh) coordinated Ga2 site. Table I summarizes the results for the four group 14 dopants, giving the site preference energies EOh−Td and the ε(+/−) transition level. Among the positively charged dopants, Sn strongly prefers the octahedral Ga2 site, whereas the other dopants occupy the tetrahedral Ga1 site, with a clear trend of increasing site preference energies toward tetrahedral preference for the lighter dopant atoms (see Table I), reflecting the expected trend of coordination number versus ionic radii, similar as previously reported in Ref. 5.

TABLE I.

Calculated site preference energy ΔEOh−Td of the group 14 dopants C, Si, Ge, and Sn in β-Ga2O3 in the (+) and (−) charge states, and the resulting ε(+/−) transition level with respect to the CBM.

ΔEOh−Td (+) (eV)ΔEOh−Td (−) (eV)ε(+/−) − ECBM (eV)
CGa +0.89 +1.57 −0.81 
SiGa +0.68 +0.56 +0.80 
GeGa +0.08 −0.31 −0.29 
SnGa −0.96 −0.50 −0.19 
ΔEOh−Td (+) (eV)ΔEOh−Td (−) (eV)ε(+/−) − ECBM (eV)
CGa +0.89 +1.57 −0.81 
SiGa +0.68 +0.56 +0.80 
GeGa +0.08 −0.31 −0.29 
SnGa −0.96 −0.50 −0.19 

As an impurity or dopant in oxides, carbon has been previously considered to occupy anion, cation, and interstitial sites.22–24 Lyons et al.24 have recently reported hybrid functional calculations showing that under most conditions, the carbon impurity in ZnO, In2O3, and Ga2O3 prefers the substitutional cation site rather than the anion or interstitial sites. Thus, it was suggested to act as a positively charged shallow electron donor. For Ga2O3, it is found here, however, that a negative charge state becomes highly favorable when the Fermi level approaches the CBM. This state is formed by an off-site configuration of CGa1, analogous to the well-known DX center in III-V semiconductors,25 denoting a donor atom changing from the usual positive charge into a negative state via large lattice relaxation. Since the DX center becomes increasingly stable when the bandgap is increased, e.g., due to alloying25 or pressure,26 its prevalence may be expected in wide gap systems like Ga2O3. Figure 1 shows the wavefunction localization of the negative charged dopants. The carbon DX center is formed by relaxation of the C atom in a direction opposite to a C–O bond toward an interstitial site, leading to a strong localization of the defect wavefunction. This relaxation effect can be understood as a bond-breaking mechanism due to the occupation of a C–O antibonding state, which also acts to reduce the strain energy associated with the atomic size mismatch between Ga and C. In the case of the larger dopants Ge and Sn, the off-site relaxed configuration exists only as a metastable local minimum, and the negatively charged defects incorporate substitutionally at the octahedral Ga2 site. Despite maintaining the substitutional configuration, GeGa2 and SnGa2 also form localized defect states inside the bandgap [cf. Fig. 1(c)], and only the Si dopant shows the delocalized character expected for a truly shallow donor.

FIG. 1.

[(a)–(d)] Localization of the negatively charged CDX, SiGa1, GeGa2, and SnGa2 defect states in Ga2O3, illustrated by an isosurface (green) of the corresponding electron density. For the gap states of C, Ge, and Sn, the isosurface density is 5 × 10−2 e3. For the delocalized Si state in (b), which lies above the CBM, this threshold is not surpassed anywhere, and a density of 5 × 10−3 e3 is shown instead. In (a), the broken C–O bond is indicated by a dashed line.

FIG. 1.

[(a)–(d)] Localization of the negatively charged CDX, SiGa1, GeGa2, and SnGa2 defect states in Ga2O3, illustrated by an isosurface (green) of the corresponding electron density. For the gap states of C, Ge, and Sn, the isosurface density is 5 × 10−2 e3. For the delocalized Si state in (b), which lies above the CBM, this threshold is not surpassed anywhere, and a density of 5 × 10−3 e3 is shown instead. In (a), the broken C–O bond is indicated by a dashed line.

Close modal

As expected from the delocalized wavefunction, Si is also the only group 14 dopant with a ε(+/−) level above the CBM (Table I), implying that the low-valent Si(ii) state forms a resonance inside the conduction band continuum but not a gap level. In such a situation, the screened potential of the positively charged dopant creates a hydrogenic donor state. Based on the electron effective mass of 0.31 m0 and the static dielectric constant of 10.8 (3.1 and 7.7 for the electronic and ionic contributions, respectively) calculated for Ga2O3, we obtain 36 meV for the donor ionization energy, small enough to allow for full ionization at room temperature. For the C impurity, the DX center behavior causes a deep ε(+/−) level at 0.8 eV below the CBM. Due to this deep level, the C impurity is electrically inactive and acts as an amphoteric defect that can even reduce the electron concentration due to other dopants. This implication is particularly relevant for chemical vapor deposition techniques that can result in considerable carbon contamination.27 

In the case of Ge and Sn, the ε(+/−) pinning level lies below the CBM (Table I) but possibly close enough to allow for ionization and free electron generation, in particular, when considering that the temperature dependence of the bandgap will further reduce the ionization energy.9 However, a certain reduction of the free electron density due to partial self-compensation might occur depending on the conditions (e.g., temperature and dopant concentration). A recent magnetic resonance study28 reported a negative-U center with an energy of less than 50 meV below the CBM in unintentionally doped Ga2O3. The prediction that Sn and Ge have ε(+/−) transitions in close proximity to the CBM is consistent with these observations. These are negative-U transitions similar to that of the DX center, although the involved lattice relaxations do not involve an off-site displacement. Rather, the relaxation pattern is more symmetric, with an increase of the Sn–O and Ge–O bond length by about 0.25 Å in the (−) state, indicating the formation of low-valent Sn(ii) and Ge(ii) oxidation states analogous to the situation in Bi-doped SnO2.20 For Sn doping, it should be noted that the ε(+/−) level at 0.2 eV below the CBM (Table I) pertains to Sn incorporation at the energetically favorable octahedral Ga2 site, but the corresponding level at the Ga1 site is significantly deeper at 0.4 eV below the CBM. This site dependence of the transition level could considerably reduce the doping efficiency of Sn when “wrong-site” occupancy occurs in non-equilibrium deposition techniques.29 In contrast, Si forms a shallow effective-mass level irrespective of the site occupancy. All dopants with a ε(+/−) level inside the bandgap, i.e., C, Ge, and Sn, could cause persistent photoconductivity, which is a typical fingerprint of such negative-U centers.

We proceed by considering Sn as the most commonly utilized1 and Si as the most favorable dopant identified here. Calculation of the defect equilibrium requires as input the formation energies of the dopant, of the compensating defect, and of the dopant-defect binding energies. Following Ref. 15, we give in Table II the defect formation energy ΔHDref for a reference condition, using Δμα = (μα − μαelem) = 0 and ΔEF = (EFECBM) = 0, i.e., the chemical potentials of the elements α are set to those of the standard phase at 0 K, and the Fermi energy is set to the CBM. The actual formation energies for a given thermodynamic condition and a given Fermi level are then simply obtained by
ΔHD,q(Δμα,ΔEF)=ΔHref+αnαΔμα+qΔEF.
(1)
Here, nα and q are the number of atoms and electrons, respectively, being exchanged with the reservoir when the defect is formed, and q is the charge state of the defect. The required corrections to the DFT energies,15 i.e., the GW calculated band edge shifts and the fitted elemental reference energies, are included in the values shown in Table II.
TABLE II.

Reference formation energies ΔHref of the dopants (D = Sn, Si) and the Ga vacancy defect in Ga2O3, and the binding energies ΔEb of the dopant-defect complexes. All values in eV.

D = SnD = Si
ΔHref DGa1+ +3.71 −0.62 
 DGa2+ +2.75 +0.05 
 VGa3− +1.84 
ΔEb (VGa − D)2− −1.17 −0.94 
 (VGa − 2D) −1.79 −1.60 
 (VGa − 3D)0 −2.05 −1.83 
D = SnD = Si
ΔHref DGa1+ +3.71 −0.62 
 DGa2+ +2.75 +0.05 
 VGa3− +1.84 
ΔEb (VGa − D)2− −1.17 −0.94 
 (VGa − 2D) −1.79 −1.60 
 (VGa − 3D)0 −2.05 −1.83 

Finite temperature effects are taken into account by considering the ideal gas free energy contribution30 for the gaseous species (O2 and Ga2O), as well as the configurational entropies associated with defect and defect pair formation.12 The temperature dependence of the bandgap, which affects the balance between charged defects and free carriers via the ΔEF dependence of ΔHD,q in Eq. (1) (see also Ref. 10), was obtained from molecular dynamics calculations as dECBM/dT = 3 × 10−4 eV/K. Additional temperature dependent vibrational effects can in principle be taken into account by performing phonon calculations31,32 for the defect and host supercells, as well as for all secondary phases, and then adding the respective free energy contributions to ΔHref and Δμα in Eq. (1). Such calculations were, however, beyond the scope of the current work.

In agreement with the previous computational literature,5–7 we find that O vacancies are deep defects that do not contribute to the net doping and that O interstitials are consistently higher in energy than Ga vacancies. Therefore, the O vacancies and interstitials are not explicitly included in the defect phase diagram discussed below. The lower formation energy of Ga vacancies is a notable difference compared to the related sesquioxide In2O3, where negatively charged O interstitials are instead the dominant compensating defect.9,33 This difference is likely a consequence of the smaller cation size in Ga2O3, which reduces the open space at interstitial lattice positions, and is also supported by recent positron annihilation experiments.34 Further, it must be noted that contamination by impurities, especially alkali, earth-alkali, and transition metals, is likely to give rise to additional compensation mechanisms35 not considered in the present model.

A direct comparison of the Sn and Si formation energies given in Table II is not meaningful because the actual energies are subject to boundary conditions for the chemical potentials ΔμGa, ΔμSn, and ΔμSi. Specifically, the existence of the Ga2O3 phase implies 2ΔμGa + 3ΔμO = ΔHf(Ga2O3), where ΔHf(Ga2O3) = −11.26 eV is the calculated formation enthalpy of the monoclinic β phase, in good agreement with the tabulated experimental value of −11.16 eV.36 The Sn and Si solubility limits are determined by similar conditions corresponding to the formation of SnO2, SnO, or SiO2. However, it is well known that the thermodynamic solubility limits can often be exceeded, as the exsolution of dopant phases requires extended annealing at high temperatures, as shown, e.g., for Sn doped In2O3.37 The precipitation of such dopant phases is impeded by the required long-range diffusion of the dopant species and by nucleation barriers for the crystallization of the secondary phases. Therefore, supersaturated, non-equilibrium dopant concentrations are considered in the following, but the regions in the defect phase diagram that lie within the solubility limit are delineated from those that lie beyond.

The atomic structure of the Ga vacancy is not described by simple removal of a Ga atom from a crystallographic lattice site. As shown by Varley et al.,6 a so-called split-vacancy configuration is instead the lowest energy structure, in which one Ga atom assumes an interstitial position between two vacant Ga1 sites. We find here that a split-vacancy structure is also strongly preferable for the (VGa − Sn) and (VGa − Si) defect pairs with the dopant atom at the interstitial site. Considering that there exist other cases, notably Cu2O, where the split vacancy configuration is stabilized by the pairing between the ZnCu donor and the substitutional Cu vacancy,38 this atomic structure feature appears to be a common pattern in the passivation of doping through formation of defect pairs and complexes. The binding energies ΔEb given in Table II are given as the formation energy of the defect complexes with respect to the formation energies of the isolated species, using the lower of the two isolated dopant energies, i.e., SnGa2 and SiGa1. The attachment of a second SnGa2 or SiGa1 dopant to the split-vacancy type dopant-defect pair still yields a significant energy gain of about 0.6 eV, but the third dopant is only weakly bound by about 0.2 eV (Table II). The dopant-defect binding is weaker for Si than for Sn, suggesting that Si is slightly less susceptible to dopant passivation.

The defect equilibria are calculated by the numerical solution of a self-consistency condition for the concentrations of the dopants, defects, and free electrons under the constraint of charge neutrality between charged defects and electrons,9,10 while including the temperature dependence of the bandgap. Further, the enthalpy and (configurational) entropy change due to dopant-defect association is taken into account during the self-consistency loop.12 The resulting defect phase diagram, shown in Fig. 2, describes the net doping NDNA, where the donor concentration ND is taken as the total Sn or Si concentration, and the acceptor concentration NA is three times the total concentration of VGa, accounting for their triple negative charge state. The solid, dashed, and dotted lines delineate the different regions (“defect phases”) relating to the thermodynamic solubility, to the compensation, and to the passivation of doping, respectively. We observe that the net doping increases with temperature, which is a result of the decreasing bandgap and chemical potential ΔμO for a given partial pressure pO2. However, as discussed below in more detail, the volatility of the molecular Ga2O species11 poses additional stability constraints at high T and low pO2.

FIG. 2.

Defect phase diagram for Sn and Si doping of Ga2O3. The contour plot shows the net doping (NDNA) as a function of the dopant concentration and the O2 partial pressure pO2 for three temperatures. (Additional constraints limiting the stability at low pO2 are discussed below.) The solid lines delineate dopant concentrations above and below the thermodynamic solubility limit (th-sol). Higher concentrations are metastable (meta) with respect to precipitation of the indicated dopant phases (SnO2, SnO, Sn, SiO2, or Si). The dashed line gives the 50th percentile for dopant activation (act) vs compensation (comp). Analogously, the dotted line gives the 50th percentile for dopant passivation, i.e., isolated (isol-dop) vs defect paired (def-pair) dopants.

FIG. 2.

Defect phase diagram for Sn and Si doping of Ga2O3. The contour plot shows the net doping (NDNA) as a function of the dopant concentration and the O2 partial pressure pO2 for three temperatures. (Additional constraints limiting the stability at low pO2 are discussed below.) The solid lines delineate dopant concentrations above and below the thermodynamic solubility limit (th-sol). Higher concentrations are metastable (meta) with respect to precipitation of the indicated dopant phases (SnO2, SnO, Sn, SiO2, or Si). The dashed line gives the 50th percentile for dopant activation (act) vs compensation (comp). Analogously, the dotted line gives the 50th percentile for dopant passivation, i.e., isolated (isol-dop) vs defect paired (def-pair) dopants.

Close modal

At low dopant concentrations c and low pO2, the dopants are fully activated, i.e., the net doping approximately equals c. Increasing dopant concentrations raise the Fermi energy and lower the formation energy of the VGa3− defects [cf. Eq. (1)], thereby increasing the compensation and lowering the doping efficiency. Increasing pO2 reduces the chemical potential ΔμGa and also lowers the net doping by increasing VGa formation. Figure 2 further shows that once dopant passivation sets in (formation of close dopant-defect pairs), the net doping decreases with increasing dopant concentration. This is a new effect arising from the inclusion of dopant-defect association in the calculation of the defect equilibrium. If the dopants and defects are only considered as isolated species, NDNA eventually saturates due to compensation but never decreases.9 As a consequence of passivation, there is an optimal dopant concentration copt for any (T, pO2) condition. Depending on the conditions, copt can be remarkably small, as low as 0.01% cat. (see Fig. 2). An interesting observation in the passivation behavior is that even at high dopant concentrations, the charged (VGa − 2Si) complex remains more abundant than the charge-neutral (VGa − 3Si)0 complex, which is due to the weak binding of the third Si atom (Table II). The significance of this finding is that most Si atoms are part of mutually repelling charged complexes, which is a situation that inhibits the exsolution of dopants by nucleation and precipitation of secondary phases (e.g., SiO2). This insight about incomplete dopant-defect passivation should serve as a design rule for metastable doping above the solubility limit, generally favoring dopants with weaker binding, but especially for the last dopant required to passivate the complex into the charge neutral state.

Although Fig. 2 suggests that very high electron densities up to n = 1021 cm−3 could be achievable in Ga2O3 when T is high and pO2 is low enough, it will be difficult to stabilize Ga2O3 under such conditions. A recent study of Ga2O3 thin-film growth by molecular beam epitaxy found that desorption of the sub-oxide Ga2O can limit the thin-film growth.11 In order to incorporate the limits to (T, pO2) resulting from the volatility of the molecular Ga2O species into the defect phase diagram, we calculated the partial pressure pGa2O over Ga2O3 as a function of T and pO2 within the ideal gas approximation, using the tabulated thermodynamic data for Ga2O.36 A similar analysis was also made for atomic Ga, but under most conditions, the decomposition according to
Ga2O3(s)Ga2O(g)+O2(g)
(2)
is the stability limiting process. As seen in Fig. 3, a high pGa2O above 1 atm builds up at higher temperatures and lower O2 partial pressures. In a low pressure thin-film deposition system or during post-deposition annealing, Ga2O3 is expected to be unstable under these conditions.
FIG. 3.

Partial pressure of molecular Ga2O over Ga2O3 as a function of T and pO2, calculated within the ideal gas approximation.

FIG. 3.

Partial pressure of molecular Ga2O over Ga2O3 as a function of T and pO2, calculated within the ideal gas approximation.

Close modal

To guide experimental optimization efforts, we analyzed the defect phase diagram to determine the maximal NDNA for Si doping under the constraint of a partial pGa2O ≤ 10−5 atm. The maximal pGa2O for which Ga2O3 can be stable depends of course on the experimental setup and could even exceed 1 atm in a high-pressure system. The value of 10−5 atm is rather conservative and was chosen such that a Ga2O3 film is likely to survive a post-deposition treatment even in an open system under flowing gas. Table III shows the temperature dependence of the optimal doping conditions. The highest net doping is obtained at lower temperatures, where very low O2 partial pressures are possible without Ga2O desorption. As the temperature increases, also pO2 needs to be increased, so to keep pGa2O within bounds (cf. Fig. 3). The Si concentration which optimizes the net doping is in the order of 1 cation per cent but reduces significantly with increasing temperature. The predicted doping limit is about 1020 cm−3, which agrees with the highest carrier densities reported in the literature for Si doped Ga2O3.39–41 In attempting to compare the present model with process parameters in various deposition techniques, e.g., pulsed laser deposition39,40 or metal organic vapor phase epitaxy,41 it should be noted that the elemental chemical potentials in these processes may not follow from simple ideal gas thermodynamics. Generally speaking, a better understanding of the mapping of process parameters onto effective chemical potentials, e.g., as recently done for non-equilibrium deposition of metastable nitrides by sputtering,42 is highly desirable. On the other hand, post-deposition treatments could more directly relate to thermodynamic models, but special attention needs to be paid toward the kinetics, as seen, for example, in the dramatic increase of equilibration times with decreasing pO2 in ZnO:Ga.13 

TABLE III.

Optimization of the net doping NDNA under the constraint of pGa2O ≤ 10−5 atm. Given are the minimum pO2 allowed under this constraint, the resulting maximal net doping, and the required Si concentration.

T (°C)log(pO2/atm)NDNA (1020 cm−3)Si concentration (1020 cm−3)Si % cat.
600 −29.2 1.4 11.0 2.9 
700 −22.6 1.0 6.9 1.8 
800 −17.4 0.8 4.5 1.2 
900 −13.0 0.6 3.0 0.8 
1000 −9.2 0.4 2.2 0.6 
T (°C)log(pO2/atm)NDNA (1020 cm−3)Si concentration (1020 cm−3)Si % cat.
600 −29.2 1.4 11.0 2.9 
700 −22.6 1.0 6.9 1.8 
800 −17.4 0.8 4.5 1.2 
900 −13.0 0.6 3.0 0.8 
1000 −9.2 0.4 2.2 0.6 

To provide computational details, the first-principles calculations were performed by the projector augmented wave (PAW) method43 as implemented in the VASP (Vienna Ab-initio Simulation Package) code for DFT,44 hybrid functional,45 and GW46 calculations. The GGA functional of Ref. 47 was used for DFT exchange and correlation. The defect levels were calculated in 60 atom defect supercells with a Γ centered 2 × 2 × 2 k-mesh, using the single-shot G0W0 approach20 based on an initial calculation of the atomic relaxation and wavefunctions in the hybrid functional17 with the values α = 0.3 and μ = 0.2 Å−1 for the fractional Fock exchange and the range separation parameters, respectively. The GW step increases the CBM energy relative to the occupied defect states, thereby lowering the ε(+/−) level relative to the CBM by up to 0.8 eV. The defect formation energies were performed in DFT-GGA, using 160 atom supercells, also with a Γ centered 2 × 2 × 2 k-mesh. For bandgap corrections according to Ref. 15, we used the GW result of Eg = 4.96 eV for β-Ga2O3 with the individual band edge shifts of ΔEVBM = −1.43 eV and ΔECBM = 1.10 eV relative to GGA.48 Further details on the GW calculations are given in Ref. 49. For an improved description of the thermodynamic boundary conditions,15 we used the fitted elemental reference energies of Ref. 50. The ideal gas law was used for conversion between chemical potentials and partial pressures at a given temperature,30 using tabulated data for O2 and Ga2O.36 The temperature dependence of the bandgap was determined from molecular dynamics simulations in GGA at 1000 and 1200 K over 1 ps simulation time. The atomic structure files for the dopants and the dopant-defect pairs are available in the supplementary material.

In conclusion, a first-principles-based defect phase diagram was constructed for n-type doping of Ga2O3. Under optimal conditions, which are subject to T, pO2, and dopant concentration c, the model predicts electron densities up to the 1020 cm−3 range. Among the group 14 dopants, only Si is a truly shallow donor, whereas the more commonly used Sn dopant creates a defect state close to the conduction band minimum and is susceptible to non-equilibrium disorder. Carbon exhibits a deep level arising from DX center behavior leading to Fermi-level pinning, compensation, and possibly persistent photoconductivity. Thus, in addition to common electron compensating impurities (e.g., alkali, earth-alkali, and transition metals), also carbon contamination deserves special attention, in particular, since residual C incorporation often occurs in chemical vapor deposition. For two reasons, optimizing the dopant concentration is not a simple matter of “more is better.” First, beyond a certain optimum, a further increase in the concentration causes passivation due to formation of defect complexes with Ga vacancies and a corresponding decline of the net doping. This effect is distinct from compensation by isolated defects with the opposite charge, which causes the net doping to level out but not to decline. Second, the maximal doping is achieved at supersaturated dopant concentrations above the solubility limit, i.e., in a situation when exsolution of the dopant is thermodynamically favored. Therefore, too high doping levels would unnecessarily promote the nucleation of secondary phases (e.g., SiO2). Interestingly, within the β-Ga2O3 structure, the dopant-defect association mechanism is such that most complexes remain at least singly charged, thereby affording for an additional nucleation barrier and hinting toward a new design rule for metastable doping.

See supplementary material for the atomic structure files for the dopants and the dopant-defect pairs.

S.L. thanks Lauren Garten for fruitful discussions and preliminary experimental efforts. This work is supported as part on an Energy Frontier Research Center by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Contract No. DE-AC36-08GO28308 to the National Renewable Energy Laboratory (NREL). This work used computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy, located at NREL.

1.
M.
Orita
,
H.
Ohta
,
M.
Hirano
, and
H.
Hosono
,
Appl. Phys. Lett.
77
,
4166
(
2000
).
2.
F.
Alema
,
B.
Hertog
,
O.
Ledyaev
,
D.
Volovik
,
G.
Thoma
,
R.
Miller
,
A.
Osinsky
,
P.
Mukhopadhyay
,
S.
Bakhshi
,
H.
Ali
, and
W. V.
Schoenfeld
,
Phys. Status Solidi A
172
,
1600688
(
2017
).
3.
K. D.
Chabak
,
N.
Moser
,
A. J.
Green
,
D. E.
Walker
,
S. E.
Tetlak
,
E.
Heller
,
A.
Crespo
,
R.
Fitch
,
J. P.
McCandless
,
K.
Leedy
,
M.
Baldini
,
G.
Wagner
,
Z.
Galazka
,
X.
Li
, and
G.
Jessen
,
Appl. Phys. Lett.
109
,
213501
(
2016
).
4.
A.
Zunger
,
Appl. Phys. Lett.
83
,
57
(
2003
).
5.
J. B.
Varley
,
J. R.
Weber
,
A.
Janotti
, and
C. G.
van de Walle
,
Appl. Phys. Lett.
97
,
142106
(
2010
).
6.
J. B.
Varley
,
H.
Peelaers
,
A.
Janotti
, and
C. G.
van de Walle
,
J. Phys.: Condens. Matter
23
,
334212
(
2011
).
7.
T.
Zacherle
,
P. C.
Schmidt
, and
M.
Martin
,
Phys. Rev. B
87
,
235206
(
2013
).
8.
F. A.
Kroger
,
The Chemistry of Imperfect Crystals
, 2nd ed. (
North-Holland
,
Amsterdam
,
1974
).
9.
S.
Lany
and
A.
Zunger
,
Phys. Rev. Lett.
98
,
045501
(
2007
).
10.
11.
P.
Vogt
and
O.
Bierwagen
,
Appl. Phys. Lett.
108
,
072101
(
2016
).
12.
K.
Biswas
and
S.
Lany
,
Phys. Rev. B
80
,
115206
(
2009
).
13.
A.
Zakutayev
,
N. H.
Perry
,
T. O.
Mason
,
D. S.
Ginley
, and
S.
Lany
,
Appl. Phys. Lett.
103
,
232106
(
2013
).
14.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
15.
H.
Peng
,
D. O.
Scanlon
,
V.
Stevanovic
,
J.
Vidal
,
G. W.
Watson
, and
S.
Lany
,
Phys. Rev. B
88
,
115201
(
2013
).
16.
S.
Lany
,
H.
Raebiger
, and
A.
Zunger
,
Phys. Rev. B
77
,
241201(R)
(
2008
).
17.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
118
,
8207
(
2003
);
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
124
,
219906
(
2006
).
18.
P.
Rinke
,
A.
Janotti
,
M.
Scheffler
, and
C. G.
Van de Walle
,
Phys. Rev. Lett.
102
,
026402
(
2009
).
19.
S.
Lany
and
A.
Zunger
,
Phys. Rev. B
81
,
113201
(
2010
).
20.
H.
Peng
,
J. D.
Perkins
, and
S.
Lany
,
Chem. Mater.
26
,
4876
(
2014
).
21.
S.
Lany
and
A.
Zunger
,
Phys. Rev. B
72
,
035215
(
2005
).
22.
C.
Di Valentin
,
G.
Pacchioni
, and
A.
Selloni
,
Chem. Mater.
17
,
6656
(
2005
).
23.
S.
Sakong
and
P.
Kratzer
,
Semicond. Sci. Technol.
26
,
014038
(
2011
).
24.
J. L.
Lyons
,
D.
Steiauf
,
A.
Janotti
, and
C. G.
van de Walle
,
Phys. Rev. Appl.
2
,
064005
(
2014
).
25.
D. J.
Chadi
and
K. J.
Chang
,
Phys. Rev. Lett.
61
,
873
(
1988
).
26.
C.
Wetzel
,
T.
Suski
,
J. W.
Ager
III
,
E. R.
Weber
,
E. E.
Haller
,
S.
Fischer
,
B. K.
Meyer
,
R. J.
Molnar
, and
P.
Perlin
,
Phys. Rev. Lett.
78
,
3923
(
1997
).
27.
G. A.
Battiston
,
R.
Gerbasi
,
M.
Porchia
,
R.
Bertoncello
, and
F.
Caccavale
,
Thin Solid Films
279
,
115
(
1996
).
28.
N. T.
Son
,
K.
Goto
,
K.
Nomura
,
Q. T.
Thieu
,
R.
Togashi
,
H.
Murakami
,
Y.
Kumagai
,
A.
Kuramata
,
M.
Higashiwaki
,
A.
Koukitu
,
S.
Yamakoshi
,
B.
Monemar
, and
E.
Janzen
,
J. Appl. Phys.
120
,
235703
(
2016
).
29.
P. F.
Ndione
,
Y.
Shi
,
V.
Stevanovic
,
S.
Lany
,
A.
Zakutayev
,
P. A.
Parilla
,
J. D.
Perkins
,
J. J.
Berry
,
D. S.
Ginley
, and
M. F.
Toney
,
Adv. Funct. Mater.
24
,
610
(
2014
).
30.
J.
Osorio-Guillen
,
S.
Lany
,
S. V.
Barabash
, and
A.
Zunger
,
Phys. Rev. Lett.
96
,
107203
(
2006
).
31.
R.
Evarestov
,
E.
Blokhin
,
D.
Gryaznov
,
E. A.
Kotomin
,
R.
Merkle
, and
J.
Maier
,
Phys. Rev. B
85
,
174303
(
2012
).
32.
G.
Miceli
and
A.
Pasquarello
,
Phys. Rev. B
93
,
165207
(
2016
).
33.
G.
Frank
and
H.
Kostlin
,
Appl. Phys. A
27
,
197
(
1982
).
34.
E.
Korhonen
,
F.
Tuomisto
,
D.
Gogova
,
G.
Wagner
,
M.
Baldini
,
Z.
Galazka
,
R.
Schewski
, and
M.
Albrecht
,
Appl. Phys. Lett.
106
,
242103
(
2015
).
35.
J. H. W.
DeWit
,
J. Solid State Chem.
13
,
192
(
1975
).
36.
D. D.
Wagman
,
W. H.
Evans
,
V. B.
Parker
,
R. H.
Schumm
,
I.
Halow
,
S. M.
Bailey
,
K. L.
Churney
, and
R. L.
Nutall
,
J. Phys. Chem. Ref. Data
11
(
Suppl. 2
) (
1982
).
37.
G. B.
Gonzalez
,
T. O.
Mason
,
J. S.
Okasinski
,
T.
Buslaps
, and
V.
Honkimaki
,
J. Am. Ceram. Soc.
95
,
809
(
2012
).
38.
V.
Stevanovic
,
A.
Zakutayev
, and
S.
Lany
,
Phys. Rev. Appl.
2
,
044005
(
2014
).
39.
S.
Müller
,
H.
von Wenckstern
,
D.
Splith
,
F.
Schmidt
, and
M.
Grundmann
,
Phys. Status Solidi A
211
,
34
(
2014
).
40.
F.
Zhang
,
K.
Saito
,
T.
Tanaka
,
M.
Nishio
, and
Q.
Guo
,
J. Mater. Sci.: Mater. Electron.
26
,
9624
(
2015
).
41.
M.
Baldini
,
M.
Albrecht
,
A.
Fiedler
,
K.
Irmscher
,
R.
Schewski
, and
G.
Wagner
,
ECS J. Solid State Sci. Technol.
6
,
Q3040
(
2017
).
42.
C. M.
Caskey
,
R. M.
Richards
,
D. S.
Ginley
, and
A.
Zakutayev
,
Mater. Horiz.
1
,
424
(
2014
).
44.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
45.
J.
Paier
,
R.
Hirschl
,
M.
Marsman
, and
G.
Kresse
,
J. Chem. Phys.
122
,
234102
(
2005
).
46.
M.
Shishkin
and
G.
Kresse
,
Phys. Rev. B
74
,
035101
(
2006
).
47.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
48.
See http://materials.nrel.gov for the GW band gaps and band edge shifts used in this work.
50.
V.
Stevanovic
,
S.
Lany
,
X.
Zhang
, and
A.
Zunger
,
Phys. Rev. B
85
,
115104
(
2012
).

Supplementary Material