Quantitative modeling of point defects in β-Ga₂O₃ combining hybrid functional energetics with semiconductor and processes thermodynamics

Abstract

β-Gallium oxide (β-Ga₂O₃) is of high interest for power electronics because of its unique combination of melt growth, epitaxial growth, n-type dopability, ultrawide bandgap, and high critical field. Optimization of crystal growth processes to promote beneficial defects and suppress harmful ones requires accurate quantitative modelling of both native and impurity defects. Herein we quantitatively model defect concentrations as a function of bulk crystal growth conditions and demonstrate the necessity of including effects such as bandgap temperature dependence, chemical potentials from thermochemistry, and defect vibrational entropy in modelling based on defect formation energies computed by density functional theory (DFT) with hybrid functionals. Without these contributions, grossly-erroneous and misleading predictions arise, e.g. that n-type doping attempts would be fully compensated by Ga vacancies. Including these effects reproduces the experimental facts that melt-grown Sn-doped β-Ga₂O₃ crystals are conductive with small compensation while annealing the same crystals in O₂ at intermediate temperatures renders them insulating. To accomplish this modeling, we developed a comprehensive modelling framework (KROGER) based on calculated defect formation energies and flexible thermodynamic conditions. These capabilities allow KROGER to capture full and partial defect equilibria amongst native defects and impurities occurring during specific semiconductor growth or fabrication processes. We use KROGER to model 873 charge-states of 259 defects involving 19 elements in conditions representing bulk crystal growth by edge-fed growth (EFG) and annealing in oxygen. Our methodology is transferrable to a wide range of materials beyond β-Ga₂O₃. The integration of thermodynamic and first-principles modelling of point defects provides insight into optimization of point defect populations in growth and processing.

---

Introduction

β-Gallium oxide (β-Ga₂O₃) is of intense current interest for power electronics because of its ultra-wide bandgap, high critical field, controllable n-type doping, and the availability of native melt-grown substrates.^1–3 The properties of bulk single crystals and epitaxial layers are intimately tied to point defects and complexes, necessitating predictive models for defect concentrations resulting from varying impurities and processing histories. While density functional theory (DFT) calculations have offered critical qualitative insights into prevalent defects based on calculated formation energies^4–11 translating these findings to actual concentrations of defects present for specific real-world crystal growth and processing requires incorporation of multiple temperature, atmosphere and pressure dependent factors. For example and as we will detail below, a typical calculation holding DFT-calculated formation energies of native and Sn-related defects constant versus temperature predicts that Sn-doped crystals grown at T_melt = 2068 °C by edge-fed growth (EFG) with p_O2 = 0.02 atm should yield crystals fully-compensated by gallium vacancies (V_Ga) and related complexes (Sn_Ga–V_Ga, V_Ga–V_O, etc.). Yet Sn-doped wafers can be purchased commercially and in actuality display (nearly) degenerate n-type doping with <1% compensation.^12,13 The accuracy of DFT-based defect modelling of charge transition levels compared to experiments brought by large supercells with finite-size correction schemes, hybrid functionals and self-interaction corrections usher in a new era in defect modelling in which it is worth the effort to incorporate the aforementioned thermodynamic effects.^14,15 The case study of systematically incorporating a suite of thermodynamic contributions for β-Ga₂O₃ herein sets a benchmark for achieving quantitative agreement with experiments by combining all of these factors.^16,17

Here we utilize a comprehensive, transferrable, quantitative framework dubbed KROGER¹⁸ for modelling full, partial, and constrained defect equilibria to emulate defects resulting from real-world growth and processing.¹⁹ KROGER,¹⁸ named after F. A. Kroger who gave exhaustive treatments of point defect concentrations,^20,21 allows us to take a set of DFT-computed formation energies for the charge states of defects and complexes and compute their numbers given specified thermodynamic conditions representing real-world processing.¹⁴ β-Ga₂O₃ is used as an exemplar and through this modelling we elucidate new insights specific to this material's defect chemistry, but of course the methodology is transferrable to other materials. We account for the temperature dependencies of the bandgap and densities of states, temperature- and pressure-dependent chemical potentials from the Ga–O binary system, degenerate statistics, self-trapped holes (STHs), and a minimal quantum harmonic oscillator model for vibrational entropy change for defects which can exceed 1.5 eV per defect near T_melt. We model EFG-grown unintentionally-doped (UID), Sn-, Fe- and Mg-doped crystals by calculating defect equilibrium based on their reported total dopant and impurity concentrations rather than assuming fixed chemical potentials. We illustrate how assuming kinetic trapping of dopants in combination with the aforementioned temperature dependencies yields agreement with experiments. The ability to model various constrained equilibria gives KROGER added capabilities for modelling semiconductor processing beyond the state of the art.^22–29

We illustrate the scalability of our modelling framework and ability to explicitly handle trace impurities by including a total of 19 elements distributed across 259 defects having 873 charge states, as well as self-trapped holes localized on O_I and O_II atoms;^30,31 limited only by the availability of a self-consistent set of DFT-calculated formation energies for all these defects’ charge states. Such calculations at dozens of temperatures run within minutes on a personal computer, allowing rapid exploration of various effects. We accurately model defects in Sn, Fe, and Mg-doped edge-fed growth (EFG) crystals, Bridgeman crystals grown under different p_O2, and annealing experiments in pure O₂ at 1 atm. We compare our results with the well-established parameters of Sn-doped Ga₂O₃, such as the electron and donor concentration ratio and DLOS data on V_Ga-related defect densities.³² This work provides novel insights into defects and complexes in β-Ga₂O₃ while also illustrating the capabilities of KROGER which was built for generality and transferability to defect systems in other materials. This comprehensive approach yields insights beyond those typically possible from DFT defect formation energetics evaluated in certain rich and poor chemical potential conditions.

Results & discussion

Establishing a baseline model for defects in donor-doped EFG β-Ga₂O₃

In this paper, we primarily focus on the EFG method,³³ as the properties of commercial wafers grown by this method are consistent and well-characterized. We also provide some results in the ESI† for Bridgman³⁴ processes which utilize different noble metal crucibles. Both EFG and Bridgman growth of β-Ga₂O₃ utilize 1 atm total pressure, with EFG using 2% oxygen p_O2 and Bridgman 20% p_O2.^35,36 All of the melt-growth techniques occur at T_m = 2068 K, and despite the different p_O2 values converge on very close predictions because for gasses, chemical potential depends logarithmically on pressure. Additionally, we have carried out annealing experiments on 500–700 µm thick Sn-doped wafers in 1 atm p_O2 at 1300–1400 K for durations up to 2 weeks, which was the time required for elimination of free carrier absorption.

Fig. 1 illustrates three different models of temperature-dependent defect concentrations for EFG-grown Sn-doped β-Ga₂O₃ resulting from different assumptions about chemical potentials and impurity concentrations. To simplify plots, we have grouped similar defects and complexes, as discussed in the ESI.†Fig. 1(a) corresponds to kinetic trapping of [Sn] = 4.5 × 10¹⁸ cm⁻³, [Si] = 2 × 10¹⁷ cm⁻³ and [Fe] = 7 × 10¹⁶ cm⁻³ and O and Ga chemical potentials fixed by eqn (2) and (7) (see in Computational Methods) with p_O2 = 0.02 atm. Subfigure (b) presents results assuming the O and Ga chemical potentials do not vary with temperature for O-rich conditions, in other words for the “O-rich” conditions presented in Varley et al.^8–10 and Frodason et al.^4–6 Part (c) is the same as (a) but with Δμ_Sn set by equilibrium with SnO₂ using measured thermochemical data. The different temperature dependencies of Δμ_Ga, Δμ_O, and Δμ_Sn for these cases are plotted in Fig. 2. In Fig. 1(a) and (c), the total concentrations of other impurities are fixed at representative values from Kuramata et al.³³ In (a) and (c) the temperature dependence of bandgap is taken from Lee et al.³⁷ with a representative assumption of 40% of ΔE_g(T) occurring in the conduction band, as discussed in depth below in the context of Fig. 4. While not shown here, it is rather a simple task to combine models like (a) and (c) to make models accounting for dopant equilibration with 2nd phases at high temperatures but freezing in at some intermediate temperature.

Fig. 1 Defect concentrations calculated for Sn-doped β-Ga₂O₃ in different scenarios. (a) Fixed concentrations of Sn and other impurities and incorporating the full set of temperature-dependent effects for p_O2 = 0.02 atm used in EFG growth.³³ Both equilibrium and quenching scenarios are consistent with Sn-doped wafers being n-type with <1% compensation as long as the native defect system freezes-in by ∼1950 K. Overlap of equilibrium and quenched results predicts insensitivity to cooling rate. (b) “O-rich, Sn-doped” conditions³⁸ without any temperature dependencies predicting that Sn-doped wafers should be insulating with >99% compensation (n in the 10¹⁵ range) for all temperatures, which disagrees with reality. (c) Similar to (a) in terms of other impurities but with μ_Sn(T) set by equilibrium with SnO₂. In this scenario, agreement with real EFG Sn-doped crystals could occur only if the native defect system and [Sn] simultaneously freeze-in at ∼1850 K, which would be more likely for larger crystals. This would be expected to depend sensitively on cooling rate, which is improbable given the widely-reported ease of reliable, low compensation bulk doping with Sn. Note that (c) represents the temperature-dependent solubility limit for Sn in β-Ga₂O₃ from thermochemistry; the fact that [Sn] > 10¹⁸ cm⁻³ is achievable without SnO₂ precipitation is direct evidence for kinetic trapping of Sn. Note, [Si] was omitted in (c) to allow Sn to be the dominant donor for all temperatures and that the equilibrium solubility becomes dependent on the residual [Mg_Ga] + [Ca_Ga] − [Zr_Ga] at low temperatures.

Fig. 2 (a) Temperature dependent chemical potentials for Ga, O and Sn under different thermodynamic scenarios with Δμ_O solid and Δμ_Ga dashed: typical O- or Ga-rich (black), based on thermochemical functions from the Ga–O system for p_O2 = 0.02 and 1 atm (red, blue). For Sn (green), solid denotes the case where the total [Sn] is fixed while dash-dot is derived from the SnO₂ phase boundary for p_O2 = 0.02 atm. (b) Plot of Δμ_Ga and Δμ_O in Ga–O chemical potential space. At 0 or 300 K, the experimental ΔH_F = ΔG° (grey diagonal lines) is within 15% of the DFT-calculated value, but as T_m is approached this error increases to 50% and the locus of conditions for which p_O2 is held constant moves quite far from either O-rich or Ga-rich conditions (blue or red circles). This illustrates why using chemical potentials from thermochemistry gives a more realistic prediction for defects present.

We conclude that the model in (a) is the most likely to represent real-world crystals because it agrees with the physical expectation that the diffusion of Sn will freeze in somewhere near T_m when we consider large crystals having appreciable dimensions (multiple centimeters) and that three important experimental observations are satisfied. These are (1) that real-world EFG-grown Sn-doped crystals are nearly degenerately doped with negligible compensation such that [Sn] = n within experimental limits of a few %; (2) that annealing in p_O2 = 1 atm at 1300–1400 K will change such samples towards insulating as shown in Fig. 3 is widely observed;³⁹ and (3) that the concentration of compensating V_Ga-related defects and complexes responsible for the E_c-2 eV signal measured using deep level optical spectroscopy (DLOS) is in the low 10¹⁶ cm⁻³ range. We note that points (3) and (1) are related since V_Ga and their complexes are the dominant compensating native defects. The suppression of V_Ga and related complexes at intermediate to high temperature is driven primarily by the strong temperature dependence of E_g, as shown in Fig. 4 and the surrounding discussion. As shown further herein and in the ESI,† a limited range of models can yield such agreement with the 3 constraints mentioned. Model (b), which represents the O-rich point in Fig. 2(b), clearly disagrees with as-grown real-world Sn-doped wafers since it predicts that >99% of the Sn_Ga donors in Sn-doped wafers grown by EFG should be compensated with V_Ga and thus Sn-doped wafers should be insulating regardless of the temperature at which defects freeze in. Model (c) represents the extension of typical phase boundary mapping of allowed chemical potential ranges but using measured thermochemical data and for experiments carried out at constant p_O2 vs temperature (as opposed to constant Δμ_O). Note that the electron concentration follows [Sn] closely down to the temperature where [Sn] ≈ [Mg] then it starts to decrease again following Vⁱ_Ga–Sn_Ga, since Ir offers a deep donor transition and Cr offers no transitions in the upper part of the bandgap. We consider it possible but improbable that the scenario of part (c) represents actual samples. This is because (1) the [Sn] in the melt is intentionally kept below the solid solubility limit and (2) no evidence of precipitation of SnO₂ or macroscopic fluctuations in [Sn] and n depending on position or cooling rate have been reported for bulk crystals of β-Ga₂O₃:Sn.

Fig. 3 Calculated concentrations of defects and holding [Sn] = 4.5 × 10¹⁸ cm⁻³ constant for (a) O-rich conditions holding Δμ_O and Δμ_Ga as well as E_g constant vs. T as is typical in much existing literature. Note that the predicted carrier density is not even visible, thus samples would be predicted to be insulating. (b) For EFG crystal growth (p_tot = 1 atm, p_O2 = 0.02 atm, f = 0.40), and (c) O₂ annealing (p_tot = p_O2 = 1 atm, f = 0.40). Solid circles represent equilibrium concentrations, while open circles show results for quenching from each equilibration temperature to 300 K.

Fig. 4 Effects of band edge variation with temperature for conditions modeling EFG of Sn-doped wafers (p_O2 = 0.02 atm, [Sn] = 4.5 × 10¹⁸ cm⁻³) assuming (a) f = 0 such that only the valence band edge varies or (b) f = 1 such that only the conduction band edge varies. If the conduction band and Fermi energy move towards mid gap, the dominant V_Ga native acceptor-like defects are dramatically suppressed. (c) Plots of band edges for f = 0 or 1 along with the highest charge transition levels for V_Ga and its complexes with Sn, and Fermi levels versus temperature. For f = 1 and intermediate temperatures, ε_F is pushed downwards at nearly the same rate as E_C (∼1 eV/1000 K) resulting in 3, 2, and 1× times increases in formation energies for the most negative charge states of V_Ga, V_Ga–Sn_Ga, and V_Ga–2Sn_Ga complexes respectively, as seen in (d). (d) Depicts the most stable charge states of these defects and complexes versus ε_F along with the ε_F values near T_m from part (c).

In Fig. 1(a), the constraint that total [Sn] is constant manifests as the dominant Sn-containing defect or complex changing as temperature is lowered. At high temperatures, isolated Sn_Ga are the dominant species, then below ∼1300 K complexes with V_Ga begin to condense in order of decreasing maximum negative charge state with V_Ga–3Sn_Ga having only a neutral charge state for n-type doping and related complexes being stable at room temperature. We note that complexes containing 2 and 3 Sn would require significant Sn diffusion over distances of 6 nm (the mean Sn–Sn distance at mid-10¹⁸ cm⁻³ concentration) at these lower temperatures which would be kinetically hindered during boule cooling.⁵ In this work we neglect these effects, but we plan to incorporate future improvements into KROGER to account for such defect kinetics. Similar trends for the Si known to provide unintentional doping in most sample types is also seen, but for Fe and other impurities known to be in EFG-grown crystals the stability of single charge states dominate for all temperatures such that they appear as horizontal lines. At the very highest temperatures³⁹ near T_m, V_O and related complexes offering ∼1 eV ionization energies provide equilibrium n > [Sn] (but at room temperature are not ionized). This encourages formation of compensating acceptor defects that, when frozen in, result in lower N_D–N_A at room temperature. To state this plainly, although their donor levels are too deep to be ionized at room temperature, the behavior of V_O at high temperature influences [V_Ga], and thus indirectly does influence n at room temperature depending on where [V_Ga] freezes in which in turn depends upon cooling rate, sample dimensions via diffusion, and the kinetics of complexes. We also investigated the effects of unintentional hydrogen up to densities ∼10¹⁷ cm⁻³ (the upper bound estimated to be unintentionally incorporated into CZ-grown crystals^39,40). We found that this level of hydrogen only minorly changes the results for Sn-doped as-EFG-grown crystals and for O₂ annealing (in which H would all be removed);⁴¹ these results are shown in the ESI.† Finally, we tested the effects of STHs – near the melting point accounting for STHs results in much higher intrinsic carrier density such that carriers rather than defects may dominate charge balance, which again can modify the prediction of dominant defects frozen in during processing. The including or excluding STHs did not change any qualitative results, although they do tend to promote the formation of acceptors like V_Ga and slightly modify the details of predicted carrier density for samples annealed in O₂ at 1000–1100 °C.

We do note that our method of estimating ΔS_vib significantly penalizes vacancy complexes such as V_O–V_Ga divacancies (see Fig. 5). Since changes in mode frequencies may mitigate this effect, we consider our predictions for these complexes at the highest temperatures to carry some increased uncertainty. The computational costs for DFT calculations using hybrid functionals of defective supercells large enough to approach the dilute limit are nearly prohibitive; when high quality calculations become available KROGER can incorporate this effect. Fortuitously, this effect will be largest at high temperatures above which impurities and defects equilibrate for slow growth rates and will decrease in importance linearly with temperature. For these reasons we expect that our predictions for temperatures below which the defect system is frozen-in are on rather sound footing.

Fig. 5 Effects of omitting one temperature dependence at a time to investigate their importances. Conditions kept constant are [Sn] = 4.5 × 10¹⁸ cm⁻³ and f = 0.40. (a) Adopting T-independent O-rich conditions as in Fig. 1(b) predicts O_i are the dominant compensating acceptors at high temperature and at no temperature does n = [Sn]. (b) Effect of E_g(T), N_C(T), and N_V(T) held constant at 300 K values, again there is no temperature where n = [Sn]. (c) Effect of omitting ΔS_vib(T) – here there is only one temperature at which n = [Sn] so coincidence would be required for agreement with experiments. Also, the concentrations of V_Ga and divacancy complexes increase to at% concentrations, which would be possible to measure. (d) Replication of Fig. 3(b) including all effects to facilitate side-by-side comparison.

Summarizing, we find that agreement of defect modelling for β-Ga₂O₃ with real-world observations requires including at least (1) kinetic trapping of dopants like Sn rather than assuming equilibrium with competing phases, (2) temperature dependent chemical potentials taken from thermochemistry, and (3) temperature dependent band edges. The other effects discussed help to tune the details of the results herein, but may be more important for other materials and growth situations.

Importance of “real-world” chemical potentials

We now proceed to analyze and discuss the contributions of these temperature-dependent factors and related assumptions about the thermodynamic environment on predicted defect concentrations. Fig. 2(a) shows the differences in chemical potentials vs temperature for oxygen and gallium under different thermodynamic assumptions. The black solid (dashed) lines show Δμ_O (Δμ_Ga) for β-Ga₂O₃ in equilibrium with “O-rich” conditions discussed in the papers reporting our set of defect formation energies;^4–6,8–10 this corresponds to holding the chemical potentials and formation energy of Ga₂O₃ () constant vs. temperature. Of course, in experiments it is much more common to hold the oxygen partial pressure constant since O₂ is a gas. This results in a roughly linearly-decreasing dependence on T which can be estimated as −1 eV per 700 K. The solid red (solid blue) curves show the experimental temperature dependencies for holding p_O2 = 0.02 atm constant as in EFG growth and p_O2 = 1 atm as during annealing experiments in a tube furnace.^42,43 The Δμ_Ga required to maintain Ga₂O₃ in equilibrium (neither growing nor decomposing) is shown as dashed red or blue lines, taking into account also the temperature dependence of , which is driven by its specific heat capacity (C_p).

The dash-dot green series of Fig. 2(a) shows the temperature dependence of the maximum Δμ_Sn allowable without precipitation of SnO₂ for p_O2 = 0.02 atm calculated analogously to Δμ_Ga. Assuming the concentration of Sn is held constant at [Sn] = 4.5 × 10¹⁸ cm⁻³ for p_O2 = 0.02 atm as for Fig. 1(a) produces the non-monotonic solid green curve. This dependence arises from self-consistently solving for both Δμ_Sn and ε_F at each temperature thus is a function of both phase thermochemistry and the set of defects and complexes included in the calculation. The crossing point of the green solid and dash-dot curves near 1850 K represents the temperature at which the Sn solubility limit falls below 4.5 × 10¹⁸ cm⁻³, as can be seen in Fig. 1(c). In EFG crystal growth of Sn-doped substrates, the concentration of Sn added to the melt is less than the solubility limit in the solid at T_m, thus at high temperatures [Sn] in the solid is limited by supply in the melt but as the crystal cools the Sn is kinetically trapped at concentrations exceeding the solubility limit simply because the diffusion kinetics become too slow to precipitate 2nd phases.

Fig. 2(b) translates Fig. 2(a) into chemical potential space by parametrically plotting Δμ_Ga and Δμ_O at each temperature. O-rich and Ga-rich conditions are single points, shown for the DFT-calculated ΔH_F = −10.22 eV FU⁻¹.⁵ The diagonal line connecting these two points gives all combinations of Δμ_Ga and Δμ_O satisfying 2Δμ_Ga + 3Δμ_O = ΔH_F and moving between the endpoints would correspond to traversing the single-phase β-Ga₂O₃ field on the T-x phase diagram while remaining in equilibrium. Analogous lines based instead on experimental thermochemistry are presented for 300 and 2100 K; it can be seen that the standard DFT result is a minor underestimate at 300 K but that the discrepancy approaches 50% at high temperatures because of the neglect of phonon specific heat capacity. We note that crystal growth requires at least a slight departure up and right from the equilibrium lines in order to provide a driving force; analogously, decomposition/etching occurs for conditions slightly down and left. The blue (red) data series parametrically plot the temperature dependence of Ga and O chemical potentials for β-Ga₂O₃ in equilibrium with p_O2 = 1 atm (0.02 atm), in correspondence with Fig. 2(a). These do not fall on one diagonal line since becomes increasingly negative (by >4 eV FU⁻¹) from 300–2100 K. While the locus (Δμ_Ga, Δμ_O) for 1 atm pure O₂ is <15% away from “O-rich conditions” at 0 or 300 K, at elevated temperatures constant p_O2 experiments have loci many eV away from either O-rich or Ga-rich conditions, translating to enormous differences in predicted defect concentrations (factors >10⁷!). Achieving the “O-rich” conditions presented in DFT papers during bulk crystal growth from the melt would require (1) enormously elevated p_O2, (2) use of alternate O sources such as ozone,^44,45 or (3) dissociation and or excitation of O₂ molecules as in plasma-assisted molecular beam epitaxy (PAMBE). Summarizing, real-world bulk crystal growth or processing conditions are probably quite far from “O-rich” or “Ga-rich” conditions presented in many DFT computation papers, thus large differences in concentrations of defects predicted are to be expected when one or more components of the crystal is not a condensed phase. Growth methods using excited reactants can have higher chemical potentials, but their chemical potential values are difficult to quantify.

Fig. 3 investigates the effect of assumptions about O chemical potential while holding [Sn] = 4.5 × 10¹⁸ cm⁻³ for (a) O-rich conditions, (b) p_O2 = 0.02 atm, p_tot = 1 atm and (b) p_O2 = p_tot = 1 atm. To clarify the effects, all native defects were included in calculations but Sn is the only impurity included. The solid red circles show the free electron density n at the given temperature, while open circles represent predicted results for instantaneous quenching from each temperature to 300 K. It is immediately clear from Fig. 3(a) that the assumption of O-rich conditions at crystal growth and cooling temperatures would yield a prediction that Sn-doped crystals should be insulating rather than the commonly-experienced nearly degenerate doping with nearly unmeasurable compensation. Using experimental thermochemical data for p_O2 = 0.02 or 1 atm as shown in (b) and (c) results in predictions for both equilibrium and quenching compatible with that real-world observation as long as the defect system freezes-in (for large crystals) somewhere between T_m and ∼1950 K. Please note that the transition temperature at which the ratio n/[Sn_Ga] begins to deviate below 1 changes with the assumption of how much of the temperature dependent bandgap occurs in the conduction band, as discussed in detail in the following section. Fig. 3 assumes 40% (f = 0.40), but plots for f between 0 and 1 are shown in the ESI.† As noted before, we suspect that the native defect system probably freezes-in at the lower end of this range given that DLOS finds V_Ga-related defects in the low 10¹⁶ cm⁻³ range³² and annealing above ∼900 °C is required to observe changes in conductivity of doped samples and diffusion of substitutional impurities.^{42,43,46–48} We have performed (unpublished) experiments annealing Sn-doped wafers in 1 atm O₂ at 1000–1100 °C. Samples become transparent to the eye and the free carrier absorption becomes undetectable using Fourier transform infrared (FTIR) after 2 weeks. Comparing Fig. 3(b) and (c), this effect is predicted fairly well since the predicted n(T) for both equilibrium and quenching cases is lower at all temperatures for p_O2 = 1 atm compared to the 0.02 atm used for EFG growth and reaches the 10¹⁶ cm⁻³ range 1300–1400 K. However, our FTIR experiments indicate that the true carrier density is probably lower, thus some further improvements in our modeling such as details of ΔS_vib, fraction f of E_g(T) in the conduction band, hydrogen incorporation, and freeze-out of different defects at different temperatures should be investigated more fully in the future.

Importance of temperature-dependent bandstructure

Having examined the effects of different thermodynamic and kinetic constraints on the elements incorporated into the modeling, we now turn attention to the effects of temperature dependence of the bandgap, which in β-Ga₂O₃ decreases by nearly 2 eV from 300 K to T_m without even accounting for STHs.³⁷ Typical experiments like optical absorption cannot resolve the absolute band edge positions independently; only their difference E_g(T) is accessible. Analogously, DFT using pseudopotentials that exclude deep core levels can result in ambiguity in that different correction schemes for electric potential self-interactions within the supercell may result in conflicting predictions.³⁷ We thus have investigated a range of scenarios in which the temperature dependent change in bandgap ΔE_g(T) ≡ E_c(T) − E_v(T) is apportioned such that fraction f occurs in the conduction band and fraction (1−f) occurs in the valence band. With the definitions E_v(0) ≡ 0 and E_c(0) ≡ E_g(0), the absolute positions of the band edges E_c(T) and E_v(T) are given as:in which E_C (E_V) is the absolute conduction (valence) band energy and f represents the fraction of the total change in bandgap occurring in the conduction band. We note that the lighter-mass oxygen atoms comprising the valence band edge states do appear to dominate the overall band gap reduction at higher temperatures, and that smaller f values (e.g. f = 0.40) are focused on more below.

Fig. 4 shows the dramatic effect exerted by the conduction band's temperature dependence on V_Ga and its complexes, which are the dominant native defects that compensate shallow n-type dopants. In (a) and (b), we show computed defect concentrations representing EFG growth of crystals doped with [Sn] = 4.5 × 10¹⁸ cm⁻³, assuming ΔE_g(T) occurs entirely in the valence band (f = 0) in (a) and entirely in the conduction band (f = 1) in (b). The conditions are the same as in Fig. 1(a). Since ε_F is in the top half of the bandgap for n-type doping, changes in E_V with temperature have no effect on the net doping N_D–N_A and charge balance demands large concentrations of compensating V_Ga and complexes with Sn. However, changes in E_C with temperature dramatically suppress these native defects. This is one of the strongest effects we infer must be occurring in real crystals; without it, Sn-doped wafers grown by EFG should be heavily compensated as in (a) and DLOS experiments should reveal E_C-2 eV trap numbers approaching 10¹⁸ cm⁻³. Large values of f suppress V_Ga compensating defects sufficiently for the appearance of a wide temperature range over which the defect system can freeze-in and result in n = [Sn] at room temperature (ESI†). Larger values of f (>0.50) predict less sensitivity to the freezing-in temperature of compensating defects for Sn-doped EFG crystals and thus the as-grown behavior. However, smaller values of f are required to simultaneously give agreement that 1300–1400 °C O₂ annealed samples should be insulating – thus the annealing experiments provide more stringent constraints on the overall model. The ESI† shows cases for p_O2 = 1 atm for f values of 0, 0.5 and 1. We determined our most-likely estimate for f = 0.40 (±0.05 or so) as the value that best satisfies the two experimental constraints simultaneously (given all the other model assumptions). Lee et al.³⁷ concluded from DFT calculations that most of the temperature dependence of the bandgap occurs in the valence band, although an alternate absolute potential alignment scheme discussed in that work predicted the opposite. Fig. 4(a) shows that setting f = 0 is inconsistent with conducting EFG Sn-doped samples, which is one of the main experimental observations. Thus, our modelling provides circumstantial evidence for the valence band contributing slightly more than half, but direct experimental investigations using e.g. X-ray techniques would be desirable. Our mode-counting of simple quantum oscillators treatment of ΔS_vib suppresses complexes including V_O at high temperature, which will affect the predicted compensation if crystals were quenched rapidly. Fig. 5 delves further into the range of models yielding agreement with experimental observations.

Fig. 4(c) shows the temperature dependent band diagram including Fermi energy corresponding to parts (a) and (b) and, with (d), explains the mechanism behind V_Ga suppression. The black thin (thick) lines show E_C and E_V for f = 0 (f = 1), while red lines show the solution for ε_F in both cases. For f = 1, the fact that we impose [Sn] = constant becomes important, since E_C − ε_F = k_BT ln(n/N_C(T)), in the absence of strong compensation ε_F will follow the temperature dependence of E_C. Examining part (d) makes the mechanism apparent; as E_C and ε_F shift down at high temperatures (by nearly 2 eV from 0 K for f = 1), the formation energies of the dominant V_Ga,ic³⁻, (V_Ga,ic–Sn_Ga)²⁻, and (V_Ga,ic–2Sn_GaII)¹⁻ increase by approximately 6, 4, and 2 eV respectively, resulting in their near-total suppression and prevention of native defect compensation at high temperatures. The case of growth of β-Ga₂O₃ from the melt provides an extreme case of the importance of including temperature dependent band edges in defect modelling because of its high T_m and strong electron–phonon coupling,^3,37 but we recommend that this effect should be more routinely included in defect concentration modelling. If ΔS_vib and E° (see eqn (2) in Computational methods) differ for different charge states, the charge transition level between them will become temperature-dependent.⁴⁹ KROGER has the capability to incorporate such temperature-dependent formation energies, which represent one of the next frontiers in defect energy calculations.

Relative importances of temperature dependent factors

Fig. 5 investigates the criticality of including various temperature-dependent parameters – Δµ(T), E_g(T), and ΔS_vib(T) – in achieving agreement with experimental observations. In each subfigure unless noted, [Sn] = 4.5 × 10¹⁸ cm⁻³, p_O2 = 0.02 atm as for EFG, and f = 0.40 are used. Fig. 3(b) provides the results when all three temperature-dependent effects are included, and is replicated as Fig. 5(d) to facilitate side-by-side comparison. Fig. 5(a) includes E_g(T) with f = 0.40 (and the T^3/2 dependencies of N_c and N_v) and the quantum estimate for ΔS_vib but keeps chemical potentials constant at the DFT-estimated 0 K values (Δμ_Ga = −5.11 eV and Δμ_O = 0 for O-rich conditions). Despite the fact that the V_Ga are suppressed at high temperatures, agreement with reality could not be found as n ≪ [Si] for all temperatures. Notably also, the unrealistically-high Δμ_O compared even to pure oxygen conditions (see Fig. 2) causes O_i to provide compensation at high temperatures rather than V_Ga related defects including V_Ga–V_O divacancies.⁴ O_i is believed to have extremely low migration energies,⁵⁰ thus is expected remain equilibrated even at room temperature, however if other compensators like V_Ga are suppressed at high temperature because of the presence of O_i but frozen in, the presence of O_i can still change the defects present at 300 K. Thus, it would not be proper to build reduced models in which O_i was omitted. This is one example why we favor avoiding a-priori assumptions regarding dominant defects. KROGER is sufficiently computationally efficient that dominant defects and reduced models may be identified after the fact without significant delays.

Fig. 5(b) includes temperature dependent chemical potentials for O-rich conditions (p = 0.02 atm pure O₂ as in Fig. 2), but assumes a temperature-independent bandgap and density of states values. There is no possibility of agreement with real-world Sn-doped wafers since n ≪ [Sn] for all temperatures for equilibrium or full quenching. Finally in Fig. 5(c), temperature dependent band parameters and realistic chemical potentials are incorporated but ΔS_vib is omitted. Here agreement with experiment is possible, but only if the defect system happens to coincidentally freeze in at one specific temperature. Notably, concentrations of V_Ga and divacancy complexes increase to at% concentrations which we believe would have been observed if actually present in real-world samples.

Thus, at least for the case of heavily Sn-doped EFG grown crystals, the real-world temperature dependencies of chemical potentials and bandstructure are absolute requirements for quantitative defect modelling. The role of ΔS_vib is also clearly very important, but we hesitate to draw absolute conclusions because of the oversimplicity of our treatment herein; we eagerly await developments in efficient but high-fidelity computation of defect vibrational spectra. Versions of Fig. 5(a) and (c) for a range of f from 0 to 1 are presented in the ESI.†

Modeling unintentionally-doped (UID), Fe-doped, and Mg-doped crystals

Fig. 6 investigates defects predicted for cases of unintentionally-doped (UID), Fe-doped and Mg-doped β-Ga₂O₃ grown by EFG, with non-dopant-related model parameters fixed at the values determined in Fig. 1(a) unless otherwise noted. The results agree with real-world observations that UID are conducting while Fe and Mg-doped β-Ga₂O₃ are insulating. In Fig. 6(a), the UID case is assumed to have dominant impurities [Si] = 2.3 × 10¹⁷ cm⁻³, [Mg] = 5 × 10¹⁵ cm⁻³ and [Fe] = 3 × 10¹⁶ cm⁻³,³³ and a wide range of temperatures yield equilibrium n = [Si] down to 1700 K. Note that equilibrium n and STH concentrations exceed the low, unintentional [Si] at the highest temperatures; since both electrons and STHs can equilibrate from any processing temperature this is of no consequence at room temperature. In (b), [Fe] = 2.5 × 10¹⁸ cm⁻³ and [Si] = 2.3 × 10¹⁷ cm⁻³ for all temperatures while in (c) [Mg] = 2 × 10¹⁸ cm⁻³ and [Si] = 2.3 × 10¹⁷ cm⁻³.^1,51 Due to the lack of detailed data on other unintentional impurities in such crystals, (b) and (c) can be considered as hypothetical, assuming that the inclusion of Ir and other impurities remain constant. For Fe-doped (b) and Mg-doped (c), electron concentrations fall to immeasurable values by 300 K for both cases whether cooled slowly or quenched, thus both would be insulating for any cooling rates. The main difference compared to the nearly-degenerately Sn-doped cases above is the greater role played by STHs because their numbers either exceed the doping in the case of UID, or because of the acceptor dopants in the Fe and Mg cases. Thus, our modelling reproduces the major known aspects of defect quantities in Sn-, Fe-, Mg- and UID EFG-grown crystals. In the ESI,† we additionally model Bridgeman growth in Pt/Rh crucibles at p_O2 = 1 atm. Follow-on modelling of epitaxially-grown crystals by various methods would help to provide additional constraints on our quantitative defect modelling.

Fig. 6 Calculated defect concentrations for EFG grown (p_O2 = 0.02 atm) (a) UID, (b) Fe-doped Ga₂O₃ and (c) Mg-doped Ga₂O₃. Background impurities have concentrations as in Fig. 1(a), while [Fe] and [Mg] were set to 3 × 10¹⁶ cm⁻³ and 5 × 10¹⁵ cm⁻³. The predicted carrier densities are in good agreement; n = [Si] for UID while immeasurable carrier densities should be present for Fe and Mg doped which is consistent with them being insulating.

Conclusions

We quantitatively modeled point defect concentrations in β-Ga₂O₃ by combining formation energies from DFT using hybrid functionals with advanced thermodynamic modelling using a new program we dub KROGER. We find that, to achieve agreement with real-world experience that Sn-doped wafers are conductive with n = [Sn], many oft-neglected temperature-dependent effects must be included and that constant concentration rather than constant-chemical potential thermodynamic conditions are more appropriate for describing [Sn]. Extensive thermochemical data is available for many semiconductors; using this data to compute temperature- and pressure-dependent chemical potentials allows high-fidelity modelling of specific growth and processing. Including E_g(T) is also critical because it suppresses V_Ga that would otherwise strongly compensate the n-type doping. We included a simple, minimal estimate of defect ΔS_vib based on counting quantum oscillator modes per atom which is especially important for divacancies and other large complexes. Especially for oxides like β-Ga₂O₃ for which crystal growth occurs at very high temperatures, TΔS_vib amounts to 1–2 eV thus while computationally-demanding calculations are required, accurately computing this factor will be critical to advancing defect modelling. We have demonstrated a framework and workflow for high-fidelity modelling of defect concentrations in semiconductors, which we anticipate will be extremely useful in moving defect computations from qualitative insights to quantitative process-dependent predictions.

Calculation methods

Formation free energy of defects

The formation energy of isolated, dilute defects in crystals can be computed under either isobaric (Gibbs energy) or isochoric (Helmholtz energy) constraints. The total free energy of the qth charge state of the jth defect or complex is given as:in which is the finite-size corrected formation energy in the dilute limit for a Fermi level (ε_F) evaluated at the position of the valence band maximum (ε_v)⁵² herein calculated using DFT with Heyd–Scuzeria–Ernzerhof (HSE)^53,54 hybrid functionals, Δε_F represents the Fermi energy (chemical potential for electrons) relative to the electrochemical energy reference (ε_F − ε_v) used for , Δμ_i is the chemical potential of element i referenced to its standard state, m_i,j is the number of i atoms added (+) or removed (−) to form the jth defect or complex, and ΔS^vib_j,q is the associated change in vibrational entropy of the crystal. ΔS^vib_j,q may be sensitive to isochoric versus isobaric conditions for anharmonic bonding and remains computationally demanding to compute accurately for isolated defects at useful levels of theory (e.g. using HSE and very large supercells).

The defect formation energies used in this work were taken from a consistent set of hybrid functional (HSE06) calculations reported previously.^4–6,8–10 These studies follow established practices for first-principles defect modeling, including the use of large supercells, appropriate k-point sampling, full atomic relaxation, and charge correction schemes such as the Freysoldt–Neugebauer–Van de Walle method. While specific computational parameters may vary slightly among different defect types, the dataset provides a uniform and reliable foundation for the thermodynamic modeling presented here.

Defect concentrations

In the dilute limit, the number of each charge state q of defect j (n_j,q) is given by:

We use the number density of appropriately-sized unit cells N_j,q as the basis for each defect or complex since some large complexes may require multiple primitive cells (or equivalently, lattice sites or formula units). For example, we consider complexes containing up to 3Sn_Ga donors bound to a V_Ga, which would require between 2 and 4 formula units depending on whether each of these 4 elementary defects entities occupy Ga₁ or Ga₂ sites. In other words, the numerical prefactor commonly taken as N_sites should vary for large complexes. Additional degeneracy factors for configurational, electronic, and spin degrees of freedom are combined into θ_j,q.⁵⁵ All calculations in this paper included a set of 873 charge states comprising 259 defects and complexes and 19 elements in β-Ga₂O₃. Most (but not all) have been previously published in Varley et al.,^8–10 and Frodason et al.,^4–6 the formation energies are provided for all of these defects in the ESI.† The formation energy of β-Ga₂O₃ was calculated for the same supercells as −10.22 eV per formula unit (FU), which is close to but lower than the experimentally measured value at 300 K.

We also implemented self trapped holes (STH) as new categories of electronic defects in addition to the usual free electrons and holes whose numbers are calculated from effective conduction band and valence band densities of states, respectively. The self-trapping energies for STH on O_I and O_II sites are taken as 0.53 & 0.52 eV.^30,31 An important detail is that, in thermal equilibrium, both band holes and self-trapped holes are equally-accessible microstates for electron–hole pair excitation of the perfect crystal. Since every O_I and O_II atom is capable of localizing a STH and in the fully-localized limit they will form a dispersionless band, the prefactor for each STH type is the primitive unit cell density 1.91 × 10²² cm⁻³ rather than the effective valence band density of states =1.71 × 10²⁰ cm⁻³ we adopt herein for band holes.⁵⁶ While optical excitations measure the “bandgap” of β-Ga₂O₃ to be 4.8–5 eV (strongly modified by Urbach tails³), from the standpoint of calculating the intrinsic carrier density the “bandgap” is arguably 0.53 eV lower because thermal excitations can create STH's while photons cannot. Because the intrinsic carrier density including STH (n_i,STH) approaches ∼10¹⁷ cm⁻³ near T_melt, this slightly affects defect equilibrium at high temperatures but does not change any major findings herein although it may for cases of lower impurity doping concentrations.

Charge neutrality and quenched concentrations

Within the usual grand-canonical formalism, the simplest defect equilibrium problem is full equilibrium at a given temperature with chemical potentials specified for all elements. Each charge state included in modeling adds one unknown concentration but also one equation of the form eqn (3). However, E_j,q for nonzero charge states depends on the additional unknown ε_F, thus the final equation allowing implicit solution of the system of equations comes from charge balance, e.g. solving the charge-neutrality condition:in which provision is made for anonymous ionized donors (N_D) and acceptors (N_A), and the final term is the total charge in all the charge states included in the model. The conduction electron (n), valence band hole (p), and STH densities are computed using Fermi–Dirac statistics.

Eqn (3) can be overridden with fixed values if the concentrations of some charge states should be held constant. If the total concentration n_j,tot of a defect j composed of multiple charge states (given by ) should be held constant, this is accomplished by overriding eqn (3) with values from the Gibbs distribution:

This allows modelling temperature dependent defect concentrations in cases where certain defects stop equilibrating below some temperature but the charge in defect j still changes with temperature. Full quenching a system of defects from a processing temperature to a final temperature (usually 300 K) is accomplished by using eqn (5) for all defects in the model with the n_j,tot determined at the high temperature. The total concentration of impurity i contained in all defects in the model is of course simply the number of i atoms per defect j times the total defect j concentration, summed over all defects in the model:

For the host material elements, this sum must be added to the numbers of atoms in non-defective unit cells.

Specification of chemical potentials or concentrations for elements

In a grand canonical formulation of point defect formation, it is assumed that atoms are exchanged between the crystal and reservoirs setting chemical potentials for each element until equilibrium is reached. It is common (but not universal) in DFT-based work on defects to compute temperature-independent formation energies of all phases also using DFT and then evaluate temperature dependent defect concentrations in rich/poor chemical potential limits for each element. Using temperature-independent chemical potentials assumes that the specific heats of all relevant compounds are negligible. For gaseous reactants, this is corrected for the ideal gas translational kinetic energy contribution , but for multiatomic molecules and high temperatures, the contributions of rotovibrational and electronic excitations also contribute. These can be challenging to compute accurately, so we adopted a pragmatic approach recognizing that for many materials of interest, experimental temperature- and pressure-dependent thermochemical data such as standard Gibbs energies (G°(T,P)) are available that have been validated to reproduce the materials’ phase diagrams. Zinkevich and Adlinger assessed the thermodynamics of the Ga–O binary system (referenced to 0 K) and we adopt their temperature-dependent parameterizations.³⁵ Values for competing phases involving impurities were taken from the Fact Sage database.⁵⁷ The fact that oxygen is a gas for all growth conditions means that temperature and partial pressure significantly change its chemical potentialin which is the standard Gibbs energy as function of temperature, P° is the reference pressure (e.g. 1 atm), p_O2 is the oxygen partial pressure, and in eqn (6a) we have avoided cancelling out P_tot to clarify the roles of total and partial pressures (which can become confusing especially when equilibrium vapors such as Ga₂O + O₂ resulting from Ga₂O₃ decomposition are considered^35,36,58). The temperature dependencies from rotovibrational and electronic degrees of freedom are contained in . If the reference temperature for the thermochemical data is not 0 K then the defect formation energy from the elements should be corrected by . For ternary and more complex materials and for novel materials, this data may not exist thus computation of the formation free energy is required.⁵⁹

Each growth or annealing process we considered occurs under well-defined temperature, total pressure, and oxygen partial pressure (p_O2), allowing us to determine the excess chemical potential for oxygen Δμ_O in eqn (6b) in equilibrium with β-Ga₂O₃. Calculations may be made for values of chemical potentials anywhere within the boundaries of the β-Ga₂O₃ single phase field. At high temperatures, β-Ga₂O₃ may be in equilibrium with Ga_(vap), GaO_(vap) or Ga₂O_(vap) on the Ga-rich side and with O_2(g) on the O-rich side, depending on the experiment. For all growth processes herein, p_O2 is controlled directly so Δμ_Ga is determined from the formation reaction

leading to

Only for p_O2 lower than all cases herein does decomposition into Ga₂O and O₂

Ga₂O + O₂ ⇌ Ga₂O₃ (7c)

become a limiting reaction. To model specific growth processes, we carry out calculations with either the concentration or chemical potential specified for each element, representing whether or not it equilibrates with the growth environment. This is determined by kinetics like surface exchange or diffusion and thus will be dependent on dimensions and cooling rates. In polycrystals or thin films heated or cooled slowly, interstitial diffusing elements are better described via chemical potentials, while substitutional diffusers in rapidly cooled bulk crystals may be better described by fixed concentrations.^60,61 Maximal values of chemical potentials are set by equilibria with 2nd phases, but lower chemical potentials or concentration may be fixed by the particular growth process; for example when a dopant like Sn is supplied at concentrations lower than the solubility limit at T_melt for liquid phase growth. Thus, process modelling of defects requires augmenting the typical grand canonical approach with the ability to simultaneously solve for the Fermi level and for chemical potentials of elements of fixed concentration.

Vibrational entropy

In lieu of directly computing the vibrational energy changes associated with each chargestate, we adopt a simple quantum oscillator model⁶² to calculate the vibrational entropy (ΔS^vib_j,q) using one representative phonon frequency ω₀ for the solid rather than the detailed phonon band structure.⁵⁵ Each atom in the perfect lattice adds 3 degenerate modes at ω₀; therefore forming a vacancy subtracts 3, an interstitial adds 3, a substitutional defect does not change the number of modes, and complexes add or subtract modes according to the net change in number of atoms. Because of our lack of information regarding bond stiffness changes, herein we only count changes in numbers of “average” modes rather than taking into account mode frequency changes. Such information can readily be incorporated, as can integration over phonon bandstructures for perfect and defective supercells.

Within our simple treatment, the vibrational entropy and isochoric specific heat per atom are:

in which ω₀ is the ‘average’ phonon frequency, which we determined by setting the Debye temperature of the associated phonon-only specific heat to 872 K,⁶³k_B is Boltzmann's constant, T is absolute temperature and ℏ is the reduced Planck's constant.

Mass and bonding changes will indubitably result in different magnitudes and even different signs of ΔS_vib^64,65 for different charge states. Improved calculations of ΔS_vib using hybrid functionals or higher levels of theory for dominant defects and complexes is one of the final frontiers for refinement of the modelling herein. We note that crystal growth experiments are nearly always done at controlled by pressure rather than volume; thus computations should reflct consistent assumptions. The difference between the two is related to thermal expansion and constant pressure conditions would in general tend to lower phonon frequencies since interatomic energy vs. distance curves soften at larger distances. Thus our treatment can also be improved by computing ΔS_vib for different lattice constants taking into account (anisotropic) thermal expansion.

Temperature dependence of bandgap

Both bandgap and effective densities of states are strong functions of temperature; many semiconductors exhibit a change in bandgap ΔE_g(T) on the order of −0.1 eV per 300 K, which can be parameterized using e.g. the Varshini or Einstein–Debye equations^37,66 (we used the Varshini values from Lee et al.³⁷). We establish a shared absolute energy scale with 0 fixed at E_v(0) such that both band edges and chargestate formation energies and associated charge transition levels may all have individual temperature dependencies. Because the individual contributions of the valence and conduction bands are not definitively known for β-Ga₂O₃, we introduce the unknown parameter f as the fraction of ΔE_g(T) arising from the conduction band, as shown in eqn (1).^1,37,67

This allows the individual band edges, which determine the effective formation energies⁶⁸ of band holes and electrons, to change with temperature on the same energy scale as the formation energies (and thus densities) of defect chargestates. The extremal cases for f = 0 and f = 1 are depicted in Fig. 4(c). The effective density of states for isotropic, parabolic bands has T^3/2 dependence; users may define other more detailed functions additionally, for example including band non-parabolicity or even direct integration of a computed density of states times occupation function. For generality we utilize Fermi–Dirac statistics for band carriers and also can include degenerate statistics for site occupation by defects⁶⁹ if required. Based on computed effective masess for band electrons and holes, we adopted values of N_C = 3.33 × 10¹⁸ cm⁻³ and N_V = 1.71 × 10²⁰ cm^−3 56 at 300 K and the bandgap at 0 K as E_g = 5.1 eV⁷⁰

Solution methods

For full equilibrium calculations in which all elements are specified by chemical potentials, only eqn (4) must be solved; KROGER's default method is Matlab's fzero function. An initial ε_F guess is generated by computing the net charge on a grid spanning from E_v to E_c, plus 5 k_BT on either side and taking the point where the sign changes. The direct search is done because commonly the net charge barely changes with ε_F over much of the bandgap so traditional downhill optimization methods may stall out unless good initial guesses close to the solution are given. Calculating defect equilibrium with some element concentrations specified is more challenging as it requires simultaneously solving for the ε_F and corresponding chemical potentials (e.g. Δμ_Sn herein) that satisfy charge balance and yield the specified total element concentration across all defects. In KROGER this is achieved by minimizing a composite objective function composed of the absolute value of net charge (right side of eqn (4)) plus the absolute value(s) of the deviation(s) of concentration(s) of the fixed element(s) from their target number(s) (for each element, eqn (6) minus the target value). Because of the exponential dependence of concentrations on ε_F and chemical potentials, the solution is typically located in extremely narrow hypervalleys of width <k_BT surrounded by large hyperareas having ∼0 gradient; thus self-consistent solution requires rather exhaustive searches of the parameter space. KROGER retains a brute-force grid search option but we have found that particle swarm optimization followed by a traditional simplex optimization balances accurate convergence with speed. Calculations versus temperature are done efficienctly by exploiting continuity by proceeding from high to low temperature using the prior solution as the guess for ε_F and Δμ. Solutions are accepted when the composite objective function is of order 10³–10⁶ cm⁻³, which is at least 6–8 orders of magnitude less than practically-measurable concentrations.

Uncertainty estimates

Despite thorough optimization with the model described above, some uncertainties in parameter values remain. Based on sensitivity analysis and the magnitude of free energy terms, we estimate the uncertainties in TΔS_vib at T_m for each charge state to be ±1 eV because of our crude approximate model, the uncertainty in the thermal bandgap energy E_g to be ±0.3 eV (including the issues of STHs mentioned) and the uncertainties related to the HSE-calculated charge state formation energies to be ±0.5 eV.⁴⁹ We have carefully assessed that the impact, even when these are compounded, on the major findings herein is minimal. For example electron density n at 1300–2100 K changes by only a small factor 0.5–2× for Sn-doped samples because of the constraint of charge balance. Addressing these uncertainties represents frontiers of defect computation and will lead to even more accurate predictions.

Data availability

The defect properties including HSE-computed formation energies are available in the ESI.† The KROGER code is available at https://github.com/mikescarpulla/KROGER.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by funding from the Air Force Office of Scientific Research under MURI award no. FA9550-21-0078 (Program Manager: Dr Ali Sayir). The work of J. B. V. was performed under the auspices of the US DOE by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 and partially supported by LLNL LDRD Funding under Project 22-SI-003 and by the Critical Materials Institute, an Energy Innovation Hub funded by the U.S. DOE, Office of Energy Efficiency and Renewable Energy, Advanced Materials and Manufacturing Technologies Office. Financial support is acknowledged from the Research Council of Norway through the GO-POW project (grant no. 314017). Computations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway.

References

1. T. Onuma, S. Fujioka, T. Yamaguchi, M. Higashiwaki, K. Sasaki, T. Masui and T. Honda, Correlation between Blue Luminescence Intensity and Resistivity in β-Ga₂O₃ Single Crystals, Appl. Phys. Lett., 2013, 103(4), 041910,  DOI:10.1063/1.4816759.
2. N. Ueda, H. Hosono, R. Waseda and H. Kawazoe, Synthesis and Control of Conductivity of Ultraviolet Transmitting β-Ga₂O₃ Single Crystals, Appl. Phys. Lett., 1997, 70(26), 3561–3563,  DOI:10.1063/1.119233.
3. A. Islam, N. D. Rock and M. A. Scarpulla, Quantification of the Strong, Phonon-Induced Urbach Tails in β-Ga₂O₃ and Their Implications on Electrical Breakdown, J. Appl. Phys., 2024, 136(3), 035701,  DOI:10.1063/5.0211588.
4. Y. K. Frodason, C. Zimmermann, E. F. Verhoeven, P. M. Weiser, L. Vines and J. B. Varley, Multistability of Isolated and Hydrogenated Ga–O Divacancies in β-Ga₂O₃, Phys. Rev. Mater, 2021, 5(2), 025402,  DOI:10.1103/PhysRevMaterials.5.025402.
5. Y. K. Frodason, P. P. Krzyzaniak, L. Vines, J. B. Varley, C. G. Van De Walle and K. M. H. Johansen, Diffusion of Sn Donors in β-Ga₂O₃, APL Mater., 2023, 11(4), 041121,  DOI:10.1063/5.0142671.
6. Y. K. Frodason, J. B. Varley, K. M. H. Johansen, L. Vines and C. G. Van De Walle, Migration of Ga Vacancies and Interstitials in β-Ga₂O₃, Phys. Rev. B, 2023, 107(2), 024109,  DOI:10.1103/PhysRevB.107.024109.
7. H. Peelaers, J. L. Lyons, J. B. Varley and C. G. Van De Walle, Deep Acceptors and Their Diffusion in Ga₂O₃, APL Mater., 2019, 7(2), 022519,  DOI:10.1063/1.5063807.
8. J. B. Varley, First-Principles Calculations 2, in Gallium Oxide: Materials Properties, Crystal Growth, and Devices, ed. M. Higashiwaki and S. Fujita, Springer International Publishing, Cham, 2020, vol. 293, pp. 329–348 DOI:10.1007/978-3-030-37153-1_18.
9. J. B. Varley, H. Peelaers, A. Janotti and C. G. Van De Walle, Hydrogenated Cation Vacancies in Semiconducting Oxides, J. Phys.: Condens. Matter, 2011, 23(33), 334212,  DOI:10.1088/0953-8984/23/33/334212.
10. J. B. Varley, J. R. Weber, A. Janotti and C. G. Van De Walle, Oxygen Vacancies and Donor Impurities in β-Ga₂O₃, Appl. Phys. Lett., 2010, 97(14), 142106,  DOI:10.1063/1.3499306.
11. X. Zhang, J. Kang and S.-H. Wei, Defect Modeling and Control in Structurally and Compositionally Complex Materials, Nat. Comput. Sci., 2023, 3(3), 210–220,  DOI:10.1038/s43588-023-00403-8.
12. Novel Crystal Technology, Inc. | Creating Gallium oxide will create a new futuristic way of living, https://www.novelcrystal.co.jp/eng/ (accessed 2024-01-25).
13. S. Lany, Defect Phase Diagram for Doping of Ga ₂ O ₃, APL Mater., 2018, 6(4), 046103,  DOI:10.1063/1.5019938.
14. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van De Walle, First-Principles Calculations for Point Defects in Solids, Rev. Mod. Phys., 2014, 86(1), 253–305,  DOI:10.1103/RevModPhys.86.253.
15. D. Wickramaratne and J. L. Lyons, Assessing the SCAN Functional for Deep Defects and Small Polarons in Wide Band Gap Semiconductors and Insulators, Phys. Rev. B, 2024, 109(24), 245201,  DOI:10.1103/PhysRevB.109.245201.
16. V. Stevanović, S. Lany, X. Zhang and A. Zunger, Correcting Density Functional Theory for Accurate Predictions of Compound Enthalpies of Formation: Fitted Elemental-Phase Reference Energies, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85(11), 115104,  DOI:10.1103/PhysRevB.85.115104.
17. Y. Kumagai and F. Oba, Electrostatics-Based Finite-Size Corrections for First-Principles Point Defect Calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89(19), 195205,  DOI:10.1103/PhysRevB.89.195205.
18. GitHub – KROGER: Program for computing full and partial equilibria of point defects in solids, https://github.com/mikescarpulla/KROGER (accessed 2024-09-15).
19. P. C. Bowes, G. H. Ryu, J. N. Baker, E. C. Dickey and D. L. Irving, Fermi Level Pinning in Co-Doped BaTiO₃: Part II. Defect Chemistry Models, J. Am. Ceram. Soc., 2021, 104(11), 5859–5872,  DOI:10.1111/jace.17938.
20. F. A. Kröger and N. H. Nachtrieb, The Chemistry of Imperfect Crystals, Phys. Today, 1964, 17(10), 66–69,  DOI:10.1063/1.3051186.
21. F. A. Kröger, Point Defects in Solids: Physics, Chemistry, and Thermodynamics, Point Defects in Minerals, American Geophysical Union (AGU), 1985, pp. 1–17 DOI:10.1029/GM031p0001.
22. S. R. Kavanagh, A. G. Squires, A. Nicolson, I. Mosquera-Lois, A. M. Ganose, B. Zhu, K. Brlec, A. Walsh and D. O. Scanlon, Doped: Python Toolkit for Robust and Repeatable Charged Defect Supercell Calculations, J. Open Source Softw., 2024, 9, 6433,  DOI:10.21105/joss.06433.
23. J. Buckeridge, Jbuckeridge/Sc-Fermi, 2023, https://github.com/jbuckeridge/sc-fermi (accessed 2024-06-13).
24. T. Ogawa, A. Taguchi and A. Kuwabara, An Extended Computational Approach for Point-Defect Equilibria in Semiconductor Materials, npj Comput. Mater., 2022, 8(1), 79,  DOI:10.1038/s41524-022-00756-0.
25. M. Y. Toriyama, J. Qu, L. C. Gomes and E. Ertekin, VTAnDeM: A Python Toolkit for Simultaneously Visualizing Phase Stability, Defect Energetics, and Carrier Concentrations of Materials, Comput. Phys. Commun., 2023, 287, 108691,  DOI:10.1016/j.cpc.2023.108691.
26. V. Stevanovic, A Python Code for Calculating the Mean-Value (Baldereschi's) Point for Any Crystal Structure, arXiv May 1, 2024, https://arxiv.org/abs/2405.00925 (accessed 2024-06-13).
27. Y. Ma, H. Yu, Y. Zhong, S. Chen, X. Gong and H. Xiang, Transferable Machine Learning Approach for Predicting Electronic Structures of Charged Defects, Appl. Phys. Lett., 2025, 126(4), 044103,  DOI:10.1063/5.0242683.
28. B. Medasani, Mbkumar/Pycdt, 2024, https://github.com/mbkumar/pycdt (accessed 2024-06-13).
29. J. N. Baker, P. C. Bowes, J. S. Harris and D. L. Irving, An Informatics Software Stack for Point Defect-Derived Opto-Electronic Properties: The Asphalt Project, MRS Commun., 2019, 9(3), 839–845,  DOI:10.1557/mrc.2019.106.
30. J. B. Varley, A. Janotti, C. Franchini and C. G. Van De Walle, Role of Self-Trapping in Luminescence and p -Type Conductivity of Wide-Band-Gap Oxides, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85(8), 081109,  DOI:10.1103/PhysRevB.85.081109.
31. T. Gake, Y. Kumagai and F. Oba, First-Principles Study of Self-Trapped Holes and Acceptor Impurities in Ga₂O₃ Polymorphs, Phys. Rev. Mater., 2019, 3(4), 044603,  DOI:10.1103/PhysRevMaterials.3.044603.
32. J. M. Johnson, Z. Chen, J. B. Varley, C. M. Jackson, E. Farzana, Z. Zhang, A. R. Arehart, H.-L. Huang, A. Genc, S. A. Ringel, C. G. Van De Walle, D. A. Muller and J. Hwang, Unusual Formation of Point-Defect Complexes in the Ultrawide-Band-Gap Semiconductor β-Ga₂O₃, Phys. Rev. X, 2019, 9(4), 041027,  DOI:10.1103/PhysRevX.9.041027.
33. A. Kuramata, K. Koshi, S. Watanabe, Y. Yamaoka, T. Masui and S. Yamakoshi, High-Quality β-Ga₂O₃ Single Crystals Grown by Edge-Defined Film-Fed Growth, Jpn. J. Appl. Phys., 2016, 55(12), 1202A2,  DOI:10.7567/JJAP.55.1202A2.
34. K. Hoshikawa, T. Kobayashi, Y. Matsuki, E. Ohba and T. Kobayashi, 2-Inch Diameter (1 0 0) β-Ga₂O₃ Crystal Growth by the Vertical Bridgman Technique in a Resistance Heating Furnace in Ambient Air, J. Cryst. Growth, 2020, 545, 125724,  DOI:10.1016/j.jcrysgro.2020.125724.
35. M. Zinkevich and F. Aldinger, Thermodynamic Assessment of the Gallium–Oxygen System, J. Am. Ceram. Soc., 2004, 87(4), 683–691,  DOI:10.1111/j.1551-2916.2004.00683.x.
36. R. P. Burns, Systematics of the Evaporation Coefficient Al₂O₃, Ga₂O₃, In₂O₃, J. Chem. Phys., 1966, 44(9), 3307–3319,  DOI:10.1063/1.1727229.
37. C. Lee, N. D. Rock, A. Islam, M. A. Scarpulla and E. Ertekin, Electron–Phonon Effects and Temperature-Dependence of the Electronic Structure of Monoclinic β-Ga₂O₃, APL Mater., 2023, 11(1), 011106,  DOI:10.1063/5.0131453.
38. J. Buckeridge, Equilibrium Point Defect and Charge Carrier Concentrations in a Material Determined through Calculation of the Self-Consistent Fermi Energy, Comput. Phys. Commun., 2019, 244, 329–342,  DOI:10.1016/j.cpc.2019.06.017.
39. Z. Galazka, K. Irmscher, R. Uecker, R. Bertram, M. Pietsch, A. Kwasniewski, M. Naumann, T. Schulz, R. Schewski, D. Klimm and M. Bickermann, On the Bulk β-Ga₂O₃ Single Crystals Grown by the Czochralski Method, J. Cryst. Growth, 2014, 404, 184–191,  DOI:10.1016/j.jcrysgro.2014.07.021.
40. Z. Galazka, K. Irmscher, R. Schewski, I. M. Hanke, M. Pietsch, S. Ganschow, D. Klimm, A. Dittmar, A. Fiedler, T. Schroeder and M. Bickermann, Czochralski-Grown Bulk β-Ga₂O₃ Single Crystals Doped with Mono-, Di-, Tri-, and Tetravalent Ions, J. Cryst. Growth, 2020, 529, 125297,  DOI:10.1016/j.jcrysgro.2019.125297.
41. Y. Qin, M. Stavola, W. B. Fowler, P. Weiser and S. J. Pearton, Editors’ Choice—Hydrogen Centers in β-Ga ₂ O ₃: Infrared Spectroscopy and Density Functional Theory, ECS J. Solid State Sci. Technol., 2019, 8(7), Q3103–Q3110,  DOI:10.1149/2.0221907jss.
42. T. Oshima, K. Kaminaga, A. Mukai, K. Sasaki, T. Masui, A. Kuramata, S. Yamakoshi, S. Fujita and A. Ohtomo, Formation of Semi-Insulating Layers on Semiconducting β-Ga₂O₃ Single Crystals by Thermal Oxidation, Jpn. J. Appl. Phys., 2013, 52(5R), 051101,  DOI:10.7567/JJAP.52.051101.
43. J. Jesenovec, M. H. Weber, C. Pansegrau, M. D. McCluskey, K. G. Lynn and J. S. McCloy, Gallium Vacancy Formation in Oxygen Annealed β-Ga₂O₃, J. Appl. Phys., 2021, 129(24), 245701,  DOI:10.1063/5.0053325.
44. H. P. Nair, Y. Liu, J. P. Ruf, N. J. Schreiber, S.-L. Shang, D. J. Baek, B. H. Goodge, L. F. Kourkoutis, Z.-K. Liu, K. M. Shen and D. G. Schlom, Synthesis Science of SrRuO₃ and CaRuO₃ Epitaxial Films with High Residual Resistivity Ratios, APL Mater., 2018, 6(4), 046101,  DOI:10.1063/1.5023477.
45. P. Vogt, F. V. E. Hensling, K. Azizie, J. P. McCandless, J. Park, K. DeLello, D. A. Muller, H. G. Xing, D. Jena and D. G. Schlom, Extending the Kinetic and Thermodynamic Limits of Molecular-Beam Epitaxy Utilizing Suboxide Sources or Metal-Oxide-Catalyzed Epitaxy, Phys. Rev. Appl., 2022, 17(3), 034021,  DOI:10.1103/PhysRevApplied.17.034021.
46. M. Nakanishi, M. H. Wong, T. Yamaguchi, T. Honda, M. Higashiwaki and T. Onuma, Effect of Thermal Annealing on Photoexcited Carriers in Nitrogen-Ion-Implanted β-Ga₂O₃ Crystals Detected by Photocurrent Measurement, AIP Adv., 2021, 11(3), 035237,  DOI:10.1063/5.0031937.
47. R. Sun, Y. K. Ooi, P. Ranga, A. Bhattacharyya, S. Krishnamoorthy and M. A. Scarpulla, Oxygen Annealing Induced Changes in Defects within β-Ga₂O₃ Epitaxial Films Measured Using Photoluminescence, J. Phys. Appl. Phys., 2021, 54(17), 174004,  DOI:10.1088/1361-6463/abdefb.
48. N. D. Rock, H. Yang, B. Eisner, A. Levin, A. Bhattacharyya, S. Krishnamoorthy, P. Ranga, M. A. Walker, L. Wang, M. K. Cheng, W. Zhao and M. A. Scarpulla, Utilizing (Al, Ga)₂O₃/Ga₂O₃ Superlattices to Measure Cation Vacancy Diffusion and Vacancy-Concentration-Dependent Diffusion of Al, Sn, and Fe in β-Ga₂O₃, APL Mater., 2024, 12(8), 081101,  DOI:10.1063/5.0206398.
49. S. Qiao, Y.-N. Wu, X. Yan, B. Monserrat, S.-H. Wei and B. Huang, Temperature Effect on Charge-State Transition Levels of Defects in Semiconductors, Phys. Rev. B, 2022, 105(11), 115201,  DOI:10.1103/PhysRevB.105.115201.
50. M. E. Ingebrigtsen, A. Yu Kuznetsov, B. G. Svensson, G. Alfieri, A. Mihaila, U. Badstübner, A. Perron, L. Vines and J. B. Varley, Impact of Proton Irradiation on Conductivity and Deep Level Defects in β-Ga₂O₃, APL Mater., 2019, 7(2), 022510,  DOI:10.1063/1.5054826.
51. Y. Ueda, T. Igarashi, K. Koshi, S. Yamakoshi, K. Sasaki and A. Kuramata, Two-Inch Fe-Doped β-Ga₂O₃ (010) Substrates Prepared Using Vertical Bridgman Method, Jpn. J. Appl. Phys., 2023, 62(SF), SF1006,  DOI:10.35848/1347-4065/acb55a.
52. C. Freysoldt, J. Neugebauer and C. G. Van De Walle, Fully Ab Initio Finite-Size Corrections for Charged-Defect Supercell Calculations, Phys. Rev. Lett., 2009, 102(1), 016402,  DOI:10.1103/PhysRevLett.102.016402.
53. J. Heyd, G. E. Scuseria and M. Ernzerhof, Hybrid Functionals Based on a Screened Coulomb Potential, J. Chem. Phys., 2003, 118(18), 8207–8215,  DOI:10.1063/1.1564060.
54. J. Heyd, G. E. Scuseria and M. Ernzerhof, Erratum: “Hybrid Functionals Based on a Screened Coulomb Potential” [J. Chem. Phys.118, 8207 (2003)], J. Chem. Phys., 2006, 124(21), 219906,  DOI:10.1063/1.2204597.
55. I. Mosquera-Lois, S. R. Kavanagh, J. Klarbring, K. Tolborg and A. Walsh, Imperfections Are Not 0 K: Free Energy of Point Defects in Crystals, Chem. Soc. Rev., 2023, 52(17), 5812–5826,  10.1039/D3CS00432E.
56. S. Poncé and F. Giustino, Structural, Electronic, Elastic, Power, and Transport Properties of β-Ga₂O₃ from First Principles, Phys. Rev. Res., 2020, 2(3), 033102,  DOI:10.1103/PhysRevResearch.2.033102.
57. C. W. Bale, E. Bélisle, P. Chartrand, S. A. Decterov, G. Eriksson, A. E. Gheribi, K. Hack, I.-H. Jung, Y.-B. Kang, J. Melançon, A. D. Pelton, S. Petersen, C. Robelin, J. Sangster, P. Spencer and M.-A. Van Ende, FactSage Thermochemical Software and Databases, 2010–2016, CALPHAD: Comput. Coupling Phase Diagrams Thermochem., 2016, 54, 35–53,  DOI:10.1016/j.calphad.2016.05.002.
58. C. N. Cochran and L. M. Foster, Vapor Pressureof Gallium, Stability of Gallium Suboxide Vapor, and Equilibria of Some Reactions Producing Gallium Suboxide Vapor, J. Electrochem. Soc., 1962, 109(144) DOI:10.1149/1.2425347.
59. C. Sutton and S. V. Levchenko, First-Principles Atomistic Thermodynamics and Configurational Entropy, Front. Chem., 2020, 8, 757,  DOI:10.3389/fchem.2020.00757.
60. D. Mutter and S. T. Dunham, Calculation of Defect Concentrations in Cu₂ZnSnS₄ from Stoichiometry, In 2014 IEEE 40th Photovoltaic Specialist Conference (PVSC), 2014, pp. 2384–2389 DOI:10.1109/PVSC.2014.6925406.
61. D. E. Sommer, D. Mutter and S. T. Dunham, Defects in Na-, K-, and Cd-Doped CuInSe₂: Canonical Thermodynamics Based on Ab Initio Calculations, IEEE J. Photovolt., 2017, 7(4), 1143–1152,  DOI:10.1109/JPHOTOV.2017.2703958.
62. B. Fultz, Vibrational Thermodynamics of Materials, Prog. Mater. Sci., 2010, 55(4), 247–352,  DOI:10.1016/j.pmatsci.2009.05.002.
63. H. He, M. A. Blanco and R. Pandey, Electronic and Thermodynamic Properties of β-Ga₂O₃, Appl. Phys. Lett., 2006, 88(26), 261904,  DOI:10.1063/1.2218046.
64. X. Wang, T. Liu, Y. Lu, Q. Li, R. Guo, X. Jiao and X. Xu, Thermodynamic of Intrinsic Defects in β-Ga₂O₃, J. Phys. Chem. Solids, 2019, 132, 104–109,  DOI:10.1016/j.jpcs.2019.04.014.
65. K. Reuter and M. Scheffler, Composition, Structure, and Stability of RuO₂ (110) as a Function of Oxygen Pressure, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 65(3), 035406,  DOI:10.1103/PhysRevB.65.035406.
66. Y. P. Varshni, Temperature Dependence of the Energy Gap in Semiconductors, Physica, 1967, 34(1), 149–154,  DOI:10.1016/0031-8914(67)90062-6.
67. A. Mock, J. VanDerslice, R. Korlacki, J. A. Woollam and M. Schubert, Elevated Temperature Dependence of the Anisotropic Visible-to-Ultraviolet Dielectric Function of Monoclinic β-Ga₂O₃, Appl. Phys. Lett., 2018, 112(4), 041905,  DOI:10.1063/1.5010936.
68. J.-H. Yang, W.-J. Yin, J.-S. Park and S.-H. Wei, Self-Regulation of Charged Defect Compensation and Formation Energy Pinning in Semiconductors, Sci. Rep., 2015, 5(1), 16977,  DOI:10.1038/srep16977.
69. M. A. Scarpulla, Statistics of Crystal Defects from Dilute Limit to Alloys, Excluded Volume Exchange Interaction and Effective Energetics of Defect Ensembles, 2024, https://www.mrs.org/meetings-events/spring-meetings-exhibits/2024-mrs-spring-meeting/symposium-sessions/presentations/detail/2024_mrs_spring_meeting/2024_mrs_spring_meeting-4009758 (accessed 2024-07-03).
70. C. Sturm, J. Furthmüller, F. Bechstedt, R. Schmidt-Grund and M. Grundmann, Dielectric Tensor of Monoclinic Ga₂O₃ Single Crystals in the Spectral Range 0.5–8.5 eV, APL Mater., 2015, 3(10), 106106,  DOI:10.1063/1.4934705.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp04817b

---