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₃

β-Ga₂O₃ is an ultra-wide bandgap (UWBG) semiconductor with potential for a superior performance in power electronics as well as in solar blind UV photodetectors and in transparent contacts for photovoltaics.^1–3 To further advance this technology, a thorough understanding of native and impurity defect processes during, for example, single crystal growth, epitaxial growth, and other fabrication steps must be developed.

Many groups have observed the formation of insulating surface layers on n-type β-Ga₂O₃ during annealing in air or O₂.^4–6 This formation of an insulating surface layer allows for facile formation of UV photodetectors. This phenomenon was originally attributed to the elimination of V_O, which were presumed to act as shallow donors, before refined DFT calculations indicated that V_O was a deep donor unable to contribute directly to n-type conductivity.⁷ Oshima et al. documented the $t$ kinetics of the insulating layer thickness formed during O₂ annealing and demonstrated that an SiO₂ cap could prevent it.⁵ This indicates that direct exposure of the Ga₂O₃ surface to oxygen does control the compensation of n-type conductivity in some way. Meanwhile, positron annihilation experiments have demonstrated the presence of V_Ga in both bulk crystals and thin films.^8,9 These experiments show increases in V_Ga concentration or charge states accompany O₂ annealing, which makes them a strong candidate for n-type compensation.^10,11 Many mechanisms, involving oxygen, V_Ga, V_O, and n-type dopants (Sn and Si), have been proposed. Some of these involve defect complexes instead of bare V_Ga and V_O.^12–15 Despite these existing studies, the specific mechanisms and rates of V_Ga introduction via processes such as bulk defect pair formation, climb of dislocations, injection from other structural defects, or in-diffusion from surfaces remain unknown. At the root, this is because native defects such as vacancies and interstitials are significantly more difficult to observe than, for example, impurity atoms.

With the exception of the direct interstitial impurity mechanism, self- and tracer-diffusion is mediated by native defects such as vacancies or interstitials.^3,16–22 Low concentrations of isovalent substitutional cations in oxides should exhibit primarily vacancy-mediated mechanisms under oxygen-rich conditions, provided that the radius of the cation is not significantly less than that of the host atom. Such conditions will favor cation vacancies and suppress cation interstitials. Interstitial mediated mechanisms may be possible under metal-rich conditions, which favor cation interstitials. The diffusion constant of isotopic or chemical tracers is proportional to the local concentration of mediating native defects and the correlation factor parameterizing their interactions such as complexation. Thus, for vacancy mediated hopping, the tracer diffusion “constant” can vary in space and time with the concentration of mediating vacancies. Different samples may exhibit different diffusion constants by virtue of their vacancy content, an effect that may explain the typical scatter of experimental diffusion rates. This, in turn, can be affected by dopant impurities through Fermi level effects as well as the samples’ thermal and chemical history. For dopant ions that control the local Fermi level, the local equilibrium of native defects can result in that dopant’s diffusion constant being dependent on its own concentration.^23–29 

Commonly, it is assumed (and in careful experiments, ensured) that diffusion-mediating native defects rapidly equilibrate with local impurity concentrations compared to the experimental timescale. Rapid equilibration can be assumed in samples containing dense sources and sinks such as polycrystals; however, in high-perfection crystals like those in this study, equilibration may be limited by injection at or transport from remote surfaces. Experiments in which native defect processes occur on similar timescales to tracer or impurity diffusion translate to time- and space-varying vacancy concentrations and thus time- and space-varying impurity diffusion constants. A technologically important example is ion implantation in which massive supersaturations of native defects are created in distributions spatially correlated with the implanted impurity profile, which provides an alternative mechanism of concentration-dependent diffusion for short timescales and low temperatures^30,31 applicable to both dopants and isovalent impurities. Etching and layer deposition may also inject native defects from the free surface.³² 

Herein, we grow and isothermally anneal superlattices (SLs) consisting of isovalent Al tracer concentration spikes in β-Ga₂O₃. We use these SLs to investigate the diffusion of cation vacancies (V_Ga), the interdiffusion of Al, and the diffusion of intentional dopants (Fe or Sn). Al diffusion in superlattices epitaxially grown on Fe substrates occurs with spatially uniform vacancy density, while in samples grown on Sn-doped wafers, we show strong evidence for transient diffusion of cation vacancies out of the substrate. Using a novel modification of the Crank–Nicolson finite difference scheme, we infer the diffusion of the vacancies mediating Al diffusion. The final Al profiles in this work record a weighted time average of [V_Ga] at each depth. The $DVGa$ we extract should also be understood as a time-averaged value in the presence of certain Sn (donor) or Fe (acceptor) concentrations and small amounts of Al. Our experiments thus set lower bounds on the true $DVGa$ and Ga self-diffusion constant for single crystalline β-Ga₂O₃. We determine the diffusion constants for Al including its dependence on local V_Ga concentration, for V_Ga (likely in the form of Sn_Ga–V_Ga complexes), and for Fe under conditions of low vacancy concentration and low electron density. Finally, we serendipitously discovered that many other cation species such as Fe, Mn, Ni, Li, and Cu also diffused from the EFG-grown Sn-doped substrates and accumulated near the free surface.

Superlattices (SLs) consisting of 6 or 9 alternating layers of 200 nm unintentionally doped (UID) β-Ga₂O₃ and 10–15 nm (Al_xGa_(1−x))₂O₃ (AlGO), where x≅0.05 were epitaxially grown using organometallic vapor phase epitaxy (OMVPE) in an Agnitron Agilis 100 system on (010) Sn or Fe doped β-Ga₂O₃ substrates obtained from Novel Crystal Technology (NCT) grown by edge fed growth (EFG). The vendor reports Fe or Sn concentrations in the mid-10¹⁸ cm⁻³, and etch pits from nanovoids and dislocations are <10⁵/cm². The lowest practical Al concentration of $≈5$ at. % was used, and layers below the critical thickness were used in order to prevent introduction of misfit dislocations. SLs always began and ended with 200 nm Ga₂O₃ spacers. UID layers grown under the same conditions exhibit n-type doping in the low 10¹⁶/cm³ range. Annealing was carried out for times between 2 and 80 h with the SL side of the sample exposed to O₂ in a quartz tube furnace in 1 atm of flowing O₂. Prior to annealing, the annealing tube and the quartz or sapphire plates on which samples sat were etched with aqua regia to remove any metal contamination. During annealing, temperature was ramped at 10 °C/min and then held constant for the indicated time, after which it was allowed to freely return to room temperature.

Dynamic secondary ion mass spectroscopy (d-SIMS) profiles were collected before and after annealing on each individual sample by EAG Eurofins. The depth was calibrated by measured crater depths and concentrations using ion implanted composition standards. Comparing multiple SIMS profiles from each sample and between samples as well as x-ray diffraction Pendellösung fringes allowed us to quantify the uncertainty in depth scales, resulting primarily from growth rate variation, of at most Δx=±6%. This translates to relative uncertainty in extracted diffusion constants of order (Δx)²=±0.4%, which we expect is much smaller than uncertainties from sample to sample variations in doping and defects. The annealing temperature accuracy at the location of the samples was measured to be within a few degrees of nominal.

The most general microscopic formulation of tracer diffusion problems is coupled reaction–diffusion of defects and complexes wherein the tracer-native defect complexes have different diffusion coefficients than the isolated species.^33,34 The fraction of time, or population fraction, of complexes relates to the binding energy and connects to the correlation coefficient of random walk theory.²⁰ 

$DVGa$⁠, the product $[VGa]subs⋅Doo$⁠, and the ratio $[VGa]epilayer/[VGa]subs$ are the three free parameters (degrees of freedom) of each scenario directly extractable from the simulation/fitting and robustly determined for each annealed sample.

Rather narrow ranges of parameter values reproduce the Al profiles especially for the samples grown on Sn substrates. If we assume a value for [V_Ga] in the substrates (discussed below), minimizing the GOF across many diffusion simulations allows us to extract [V_Ga] in the epilayer, D_oo for the Al at each T_anneal, and D_VGa at each T. Any difference between actual and assumed $[VGa]subs$ will affect $[VGa]epilayer$ and D_oo. Thus, if the actual concentration is later found to be 10× higher than assumed in our analysis, then D_oo would be 10× lower but $[VGa]epilayer$ also 10× higher. Uncertainties in the final extracted parameter values were estimated from the curvature of the GOF with respect to each parameter. Further discussion of numerical methods and error propagation are in the supplementary material.

First, we discuss expectations for the initial conditions of concentrations of V_Ga in the substrates and superlattices and then present and initially assess the SIMS data. An intuition for a completely ionic crystal suggests that V_Ga should exist only in the 3⁻ charge state and thus have concentration scaling $∝(nni)3$ or $∝exp−3EfkBT$⁠.³⁶ Modern density functional theory (DFT) with hybrid functionals accurately predicts that defect charge states may change if E_f moves above or below charge transition levels and can account for complexation with other defects.^7,37 At equilibrium with all other factors held constant, the concentration of V_Ga should be proportional to the n-type doping concentration; however, differences in chemical potentials (e.g., $pO2$⁠) and temperatures during bulk crystal and OMVPE growth will modify this ratio.

A number of experiments^{8,9,11,13,38,39} have detected and studied cation vacancies; however, the most credible quantification of V_Ga or its complexes in Sn-doped wafers comes from the E_c-2.0 eV signal measured by deep level optical spectroscopy (DLOS). DLOS measurements on different but very similar OMVPE layers report a concentration of ≈10¹⁵/cm³ for the E_c-2.0 eV defects,⁴⁰ while DLOS measurements of E_c-2.0 eV defects in NCT Sn-doped bulk crystals are reported at 0.8–1.2 × 10¹⁶/cm³ for wafers matching the nominal Sn concentration of those used herein ([Sn] = 3.5–5 × 10¹⁸/cm³).⁴¹ For analysis, we assume [V_Ga] = 10¹⁶/cm³ vacancies in the Sn-doped substrates unless otherwise mentioned.

While it is possible that DLOS may not detect some vacancies bound in complexes, the near-unity doping efficiency measured for Sn in bulk substrates in the 10¹⁸/cm³ range indicates that [V_Ga] cannot be much larger. Unfortunately, junction capacitance based defect spectroscopies such as DLTS or DLOS cannot be performed on insulating samples, e.g., Fe-doped wafers or oxygen annealed ones; thus, comparable quantification is not available for Fe-doped substrates. All substrates were grown at T_melt (k_BT_melt ≈ 0.18 eV) by EFG under reduced $pO2$ conditions.^42,43 In Fe-doped Ga₂O₃, the Fermi level is pinned at the $Fe3+/2+$ charge transition level corresponding to the 0/- transition of the Fe_Ga acceptor located very near E_c-0.8 eV by DLTS and related techniques.^44,45 Sn-doped substrates are degenerately doped or very close to it, thus E_f ≈ E_c. For both of these E_f values, isolated V_Ga should be in the q = 3⁻ charge state; thus, $[VGa3−]$ in the Sn-doped substrates should be higher than in the Fe-doped substrates by a factor of $exp3×0.8eV0.18eV≈106$ if [V_Ga] in both substrate types equilibrate at T_melt. If V_Ga stay equilibrated to a temperature lower than the freezing temperature, this ratio would increase, while if most of the V_Ga are in $(SnGa−VGa)2−$ complexes or $(VGa−VO)−$ divacancy complexes, this ratio would be ∼10⁴ or 100, respectively. While the unintentional doping in our OMVPE layers grown at 800 °C is in the low 10¹⁶ /cm³, it is very difficult to predict the ratio of [V_Ga] in the SLs to those in the substrates because of the very different growth conditions. In addition, in unpublished experiments we have carried out, we find that annealing Sn-doped wafers in 1 atm O₂ at 1050 °C requires at least two weeks for full equilibration; thus, [V_Ga] will not equilibrate in the full wafer during OMVPE growth. Thus, we consider the ratio of [V_Ga] in the SL to that in the substrate to be an unknown parameter to be determined in the fitting. Despite the nominally identical growth of the superlattices, the initial conditions for the V_Ga or complexes mediating Al diffusion may be different from run to run, for those grown on Fe- or Sn-doped substrates, and from individual substrate to substrate.

Figures 1(a) and 1(b) show the Al and dopant SIMS profiles before and after annealing for 20 h in 1 atm O₂ at 1050 °C for two superlattices grown nominally identically on Sn- or Fe-doped substrates. For the SLs on Fe-doped substrates, Al diffuses without any detectable gradient in the extent of diffusion of the (Al, Ga)₂O₃ layers in the superlattices. Meanwhile, for SLs on Sn-doped substrates, the extent of Al diffusion is greater for (Al, Ga)₂O₃ layers closer to the substrate than the surface, and, overall, the Al diffusion is much faster. This is consistent with a larger [V_Ga] in the Sn-doped substrates compared to the SLs and to that in the Fe-doped substrates, which diffuses through the SLs toward the surface on a non-negligible timescale compared to the annealing time. For both types of samples, we find no evidence for vacancy creation or in-diffusion from the free surface (e.g. greater extent of Al diffusion near the surface) at concentrations sufficient to affect Al diffusion from about 5 at. % down to 10¹⁸/cm³. Such an effect could be taking place, but apparently at concentrations significantly lower than the estimated 10¹⁶/cm³ from the Sn-doped substrates. In the case of Fe-doped substrates, V_Ga diffusion could be much faster than the time required for Al diffusion at the given V_Ga concentration such that it does equilibrate quickly compared to the experimental time. As discussed later, Sn_Ga donors form complexes with acceptor-like V_Ga but acceptor-like Fe_Ga do not; therefore, we should expect isolated V_Ga to diffuse faster through SLs grown on Fe-doped substrates. Figure 1(c) is a 3D rendering of a similar SL (this one having nine repeat units ) constructed from Time of Flight (ToF)-SIMS after being annealed for 20 h at 1100 °C in O₂, showing the same behavior, with greater diffusion of Al near the SL/substrate interface.

For SLs on Fe-doped substrates, uniform [V_Ga] in the SL and substrate eliminates the cross-derivative term in Eq. (2) for the Al diffusion, reducing it to the usual form of Fick’s second law. We tried many scenarios of different [V_Ga] in the SL and substrate, but the best agreement with experiments was found for the case where both began with the same low concentration. Evolving the initial to final Al SIMS profiles for the sample annealed at 1050 °C allows us to extract a value of D_Al = D_oo(T) · [V_Ga] = 3.5 × 10⁻¹⁸ cm²/s, assuming [V_Ga] = 10¹⁴/cm³, i.e., $≈100$ times less V_Ga than for Sn-doped substrates, assuming D_oo(T) is the same for both.

The diffusion of Fe itself appears to be complex, probably involving at least two modes such as interstitial and substitutional (which in turn may each involve different charge states). This conclusion is reached because, in addition to the [Fe] smoothly connected to that in the wafer, we established through many exacting experiments (see the supplementary material) that Fe additionally diffuses through the central portion of the SL and piles up near the surface—i.e., the [Fe] distribution near the surface is not a SIMS artifact nor does it arise from vapor phase transport. This accumulation at the surface was reproduced in both d- and ToF-SIMS, across multiple samples, and is also reproduced for Sn-doped substrates in which [Fe] is initially considerably lower in the substrate as discussed in Sec. VI. This requires a significantly faster diffusion mode such as interstitials and could also involve γ-Ga₂O₃.

Figures 2(a)–2(c) present the initial and final SIMS profiles for [Al] for samples annealed at 1000, 1050, and 1100 °C for 20 h. The extent of Al diffusion within each SL is greater near the substrate than near the surface, indicating that the local concentration of [V_Ga] was non-uniform and greater near the substrate for a significant portion of the annealing time. Samples (a) and (b) were grown in the same OMVPE growth run, while sample (c) was grown in a different run; this is important for understanding the smaller degree of Al diffusion despite this sample’s higher annealing temperature, as discussed in more detail below. Figures 2(d) and 2(e) reproduce the data from (a) but also overlay simulated Al concentration profiles at 2 h intervals (red solid) for the best-fit parameters within two assumptions for V_Ga diffusion (green dashed) described in the following.

Diffusion parameter estimation from each experiment was accomplished by evolving the measured initial Al SIMS profile forward in time under many scenarios and fitting to the final SIMS profiles. Each tested scenario is defined by Eqs. (1) and (2), assuming initial [V_Ga] in the SL and substrate, [V_Ga] diffusion behavior explained below, and values for D_oo for Al and $DVGa$⁠. It is clear for all Sn-doped samples that the initial [V_Ga] in the SLs must be at least 10 times lower than in the substrates, which agrees with intuition given the higher substrate n-type doping. The effects of $DVGa$ and the initial [V_Ga] ratio on the broadening and peak-to-valley [Al] ratio for each Al-containing layer in the SLs are sufficiently distinct that we can determine their values with low uncertainty within each scenario.

As shown in Fig. 2(d), we first tested the hypotheses that V_Ga diffuse freely without interactions from the Sn-doped substrate through the SLs according to the error function solution of Eq. (5).⁴⁶ These simulations were able to reproduce the final Al SIMS data extremely well as judged by inspection and by GoF values ≪1 as seen in Table I. Thus, the gradient in the extent of Al diffusion in the SLs on Sn-doped substrates is consistent with this hypothesized mechanism of out-diffusion of V_Ga from the substrate.

Since the [Al] depth profiles are essentially records of the time-integrated [V_Ga] profile, it may not be possible to infer its exact functional shape from small numbers of experiments. However, the [Al] data very clearly show evidence for larger time integrated [V_Ga] near the substrate, ruling out scenarios in which, for example, [V_Ga] are introduced at the free surface and diffuse inward. In addition, it could be hypothesized that V_Ga are introduced under oxidizing conditions via Frenkel pair generation and out-diffusion of Ga_i. Since this mechanism creates electrically neutral pairs, it should occur at the same rate in the SLs, Fe-doped substrates, and Sn-doped substrates. Thus, any diffusion of Al enhanced by this mechanism would be depth-uniform and have an upper bound set by the rate observed in SLs on Fe-doped substrates—at least 2 orders of magnitude slower than for Sn-doped substrates. While this mechanism cannot be ruled out, it certainly cannot explain the experiments on Sn-doped substrates.

Diffusion of V_Ga according to Eq. (5) would imply that (a) the total number of V_Ga is conserved and (b) the functional shape of their concentration profile is unaffected by non-uniformities in [Sn](x, t) or [Al](x, t) (although the value of D_VGa may be affected). Since Al and Ga are chemically very similar, this is reasonable for Al. Frodason et al.²⁴ studied the diffusion of Sn from (001) Sn-doped substrates into UID HVPE-grown Ga₂O₃ at annealing temperatures 1050–1250 °C and fit them with a numerical reaction–diffusion (R–D) model. In their work (and in some similar intermediate Sn d-SIMS profiles obtained in this work but not shown herein), Sn diffuses with a front shape appearing on a log plot as a rounded step function. Sn is most stable as a substitutional Sn_Ga donor, and, thus, it is logical to expect that Coulomb attraction results in larger correlation for Sn than Al with V_Ga. Theoretical and experimental work has established fairly conclusively that Sn and V_Ga interact strongly while diffusing.^24,38,48

During out-diffusion from Sn-doped substrates, V_Ga dynamically complex with (on average) a small fraction of the available Sn (since [V_Ga] ≪ [Sn_Ga]). Complexed V_Ga diffuse slower than free V_Ga, while substitutional Sn_Ga diffuse only while complexed or correlated with a vacancy.²⁰ Critically for assessing how this will affect Al diffusion, the Sn_Ga–V_Ga complexation reaction’s equilibrium rate constant was determined to be small in our temperature range, meaning that each complex is bound for only a short time before dissociating.²⁴ Therefore, V_Ga passing through the SLs has a high likelihood of assisting in Al diffusion even while co-diffusing with the Sn. The effective diffusion constant for V_Ga and thus for Ga and Al diffusion in the presence of donors such as Sn should be smaller than the value in intrinsic or acceptor-doped β-Ga₂O₃ (e.g., Fe-doped). The degree of retardation of cation vacancy diffusion will depend on the concentrations of donors and vacancies and thus be variable from experiment to experiment and vs depth and time within any diffusion experiment.

At t = 0, this is a step function but the initially vertical section pivots vs time such that its intercepts with the max and min concentration values far away from the interface move as $Dt$⁠. This piecewise linear function accurately captures most of the dynamic range of the Sn_Ga–V_Ga steps from Ref. 24 and our re-analysis of those data. It only fails to capture the shape of the low-concentration leading edge of Sn diffusion (2–3 orders of magnitude below the minimum [Al] herein) and at high concentrations in the substrate (which is irrelevant for the SL diffusion). Since Al diffusion scales with [V_Ga], the details of the leading edge will not dominate the time-integrated Al diffusion. The difference between the classical erf() solution and piecewise linear R–D profiles may be seen by comparing the green dashed curves in Figs. 2(d) and 2(e). Figure 2(e) shows fitting of one example dataset using the piecewise linear model for the evolution of [V_Ga](x, t). By inspection, the Al diffusion is also well-described using this model for V_Ga diffusion and Table I shows that statistically similar GoF values are obtained for both models.

The best-fit model for Al diffusion assuming erf() diffusion for [V_Ga] suggests that a significant gradient of [V_Ga] (and thus Sn_Ga) should remain at the end of our experiments, while the RD model and our piecewise approximation thereof predict that [V_Ga] and [Sn] should be closer to uniform at the end of the experiment. The latter is in agreement with the measured, nearly uniform final [Sn] profiles as exemplified by Fig. 1(b). In addition, using the piecewise linear approximation allows us to recover the Arrhenius behavior for both Al and V_Ga. Moreover, the $DVGa$ values extracted using the piecewise linear model smoothly join with those from Ref. 24 despite $DVGa$ being extracted from our [Al] profiles and from [Sn] profiles in theirs, indicating that in both cases the diffusivity of Sn_Ga–V_Ga complexes or at least the interactions between Sn_Ga and V_Ga determine the rate limiting step(s). We take these observations as confirmations that the piecewise R–D approximation for V_Ga diffusion more accurately represents the microscopics of V_Ga diffusion in these SLs.

A subtle point is that our modeling using either the erf() or piecewise linear models explains the Al diffusion assuming that V_Ga originate in the substrate and that their total number is conserved at least to a tolerance of a few percent. However, in Ref. 24, V_Ga were assumed to appear via unspecified generation or transport processes as needed to maintain local equilibrium with the diffusing Sn—in other words that V_Ga may not be conserved. The supplementary material contains a thorough discussion of this issue as far as it can be addressed in the current work. Succinctly, the lower extent of diffusion of Al close to the free surface rules out large concentrations of [V_Ga] being injected there (“large” relative to at. % concentrations of Al). Since the concentrations of Sn in Ref. 24 range from 10¹⁴ to 10¹⁸/cm³, these experiments are sensitive to much lower concentrations of V_Ga. The high-concentration features of the [Sn](x,t) R–D diffusion front shape can be reproduced assuming that V_Ga is conserved; however, the low-concentration leading tail cannot be reproduced. Thus, we can only conclude that any bulk generation or injection of V_Ga must be below the capability of the present experiments to resolve. Resolution of these fine details should be possible in future experiments using smaller tracer concentrations designed to identify and quantify the mechanisms introducing V_Ga during annealing under oxidizing conditions. The lack of evidence for [V_Ga] in-diffusion from the free surface in our [Al] data herein raises a question for future investigations: what sources, sinks, and transport processes for [V_Ga] operate in order to establish equilibration with local Fermi levels and $pO2$ during annealing and growth? However, for the immediate purpose of analyzing our Al diffusion data, this question is immaterial since the R–D models predict the free [V_Ga] able to facilitate Al diffusion to have functional shape closely mimicking the Sn profile.

The sample annealed at 1050 °C [Fig. 2(b)] exhibits a greater extent of diffusion than the one annealed at 1100 °C [Fig. 2(c)], counter to what would be expected for Arrhenius activation if the two samples had identical initial concentrations of [V_Ga]. This was initially a source of great consternation until we realized that the sample annealed at 1100 °C [Fig. 2(c)] was grown in an earlier growth run using a different substrate, while the samples annealed at 1000 and 1050 °C [Figs. 2(a) and 2(b)] were grown at a later date on pieces of a different Sn-doped (010) β-Ga₂O₃ wafer. [V_Ga] is likely to vary spatially even within the same EFG-grown bulk crystal; thus, some variation from wafer to wafer should not be surprising. Likewise, despite using nominally the same growth recipe for all three samples, some variation is reasonably likely in [V_Ga] within the SL between growth runs.

The finite difference model in Sec. II ultimately depends on D_Al(x, t) as defined in Eq. (1); thus, fitting via simulation does not directly measure [V_Ga] but rather the D_oo(T) × [V_Ga](x, t) product. In our final analysis, we conclude that [V_Ga] = 10¹⁶/cm³ in the substrates for the 1000 and 1050 °C annealed samples, but explored possibilities that the substrate and SL for the 1100 °C sample had different [V_Ga]. In all cases, the value of [V_Ga] in the SL was determined independently for each sample through the fitting process. After testing a wide range of scenarios, the data from the 1100 °C sample are best explained by lower initial [V_Ga] = 2 × 10¹⁵/cm³ in that sample’s substrate. A more sensible behavior for D_oo(T) compatible with an Arrhenius law with activation energy 3.9 ± 0.3 eV is obtained, which is comparable to the value of 4.2 eV obtained for Al_Ga diffusion in (−201) β-Ga₂O₃ in Ref. 49. Table I summarizes the final best-fit parameters for both erf() and piecewise linear models for V_Ga diffusion, and Fig. 4 shows our data in the context of the existing literature.

In Figs. 1(a) and 1(b), we compare two SLs annealed at 1050 °C on Fe- and Sn-doped substrates. On the Sn-doped substrate, the best-fit scenario within the piecewise linear model of V_Ga diffusion yields D_oo = 2.7 × 10⁻³⁰ cm⁵/s at an initial value of [V_Ga] = 10¹⁶/cm³ in the substrate. D_oo for Al in the absence of V_Ga complexing with Sn would be expected to be the same or larger, depending on how significantly the presence of Sn limits the ability of V_Ga to mediate Al diffusion. If we do assume that D_oo is the same for samples grown on both Fe- and Sn-doped substrates, we estimate that the initial [V_Ga] in the substrate for samples grown on Fe-doped substrates would be 10¹⁴/cm³, or $≈100$ times lower that for the Sn-doped substrates. This is on the lower end of the range of estimated ratios based on expected Fermi level positions, but given the lack of direct comparison data and uncertainty regarding the effects of Sn and Fe co-diffusing with V_Ga, considered reasonable.

In addition to the d-SIMS of Al, Fe, and Sn presented so far, we also used ToF-SIMS 3D reconstructions and survey d-SIMS detecting all elements present (at lower spatial and concentration resolutions than Al profiles). These revealed that many other elements besides Fe or Sn also diffuse out of Sn-doped substrates, through the SLs, and accumulated at the free surfaces. Figure 3 shows the concentration-calibrated SIMS profiles for five illustrative metals present in the SL shown in Fig. 2(b), annealed at 1000 °C. After collecting the 20 h annealed SIMS profile shown in Fig. 2(b), the sample was immersed in aqua regia for 10 min to remove any possibility of metal contamination at the surface and then annealed in O₂ for an additional 10 h (total 30 h) as part of efforts to eliminate the possibility of Fe surface contamination diffusing into the SL while annealing inside the tube furnace (supplementary material). Given the known presence of many of these elements in wafers grown from melts held in Ir crucibles—e.g., EFG^43,47 and the high purity of OMVPE-grown layers, we infer that these elements diffused into the SLs from the substrate.

Mn, Fe, and Ni accumulate at the surface to depths that are incompatible with typical SIMS artifacts. [Cu] exhibits a concentration peak around the middle of the SL, perhaps indicating interactions between its charge states and the Fermi level within a space charge region. All these transition metals are likely capable of multiple interstitial and substitutional diffusion modes in multiple charge states. Ti was also detected at 10¹⁶/cm³ uniform in depth through the d-SIMS survey scan, while Si and Na were also detected in ToF-SIMS. Li being extremely small is suspected to be capable of fast interstitial diffusion and is depleted rather than enriched close to the surface. These different behaviors are consistent with an electrochemical potential gradient affecting different elements and charge states differently near the surface. A depletion width and built-in potential ϕ(x) induced by surface Fermi level pinning is the simplest example, and defects with different relative charges q would accumulate or deplete as $expqϕxkBT$ self-consistently following the Poisson–Boltzmann equation. However, it has been frequently observed that β-Ga₂O₃ surfaces, especially in the presence of additional impurities, can change to the different γ-Ga₂O₃ phase, which would of course change the chemical potential for each impurity, perhaps causing them to accumulate. Further analysis of these profiles is beyond the scope of this work, but the observation of these impurities, some like Cu, Na, and Li being notoriously fast diffusers in semiconductors, calls for further quantification of impurities and their diffusion in multiple phases of Ga₂O₃ and in devices from fabrication to long-term degradation under high fields and temperatures.

Figure 4 compares the measured and computed diffusivities in Ga₂O₃ from this work and that from the literature. Our values of $DVGa$ extracted using the piecewise linear model for R–D have remarkable agreement with the diffusion of Sn_Ga–V_Ga complexes reported by Frodason²⁴ and follow a sensible Arrhenius trend with activation energy 1.96 ± 1.19 eV. Our values extracted with the erf(x, t) model are an order of magnitude lower and deviate from the Arrhenius behavior. The blue rhombuses plot the free diffusion of uncomplexed V_Ga from Ref. 24 and are at least three orders of magnitude faster than our results for the vacancies responsible for Al diffusion under either scenario, which is consistent with our conclusion that the co-diffusing Sn also slows the progress of V_Ga in our experiments.

Since D_oo,Al(T) cannot be shown directly on this plot, we multiply it by two different [V_Ga] to show the range of D_Al expected in most situations aside from, for example, ion implantation. [V_Ga] = 10¹⁶/cm³ represents equilibrated conditions for heavy n-type doping as in Sn-doped EFG-grown wafers, while [V_Ga] = 10¹⁴/cm³ may be expected in UID to insulating conditions as found for our SLs on Fe-doped substrates. This latter set is in good agreement with D_Al derived from the interdiffusion of a thin Al₂O₃ layer on Ga₂O₃ during oxygen annealing.⁴⁹ 

The inverted brown triangles show D_Fe for the slow diffusion channel from Fe-doped substrates into the SLs. This value is similar to those experienced by Al, so we speculate that it may be vacancy-mediated but modified by the expected Coulomb repulsion between V_Ga and Fe_Ga acceptors compared to the case of Al. Further detailed studies will be needed to understand the diffusion of Fe and other metals prone to charge state changes and possibly diffusing by multiple mechanisms simultaneously.

Oshima et al. used capacitance methods to measure the expansion of an insulating layer on the surface of n-type samples up to a few μm thickness during O₂ annealing. The kinetics of layer thickening were consistent with a diffusion-limited process; thus, we plot their diffusion constants derived from thickness = $Dt$ as blue circles. These diffusion constants are consistent with the vacancy-mediated diffusion of Sn, Si, and Al under n-type conditions found herein and in Refs. 49 and 24. As discussed in Sec. VI, the conversion of surface layers to insulating may also include surface band bending and possibly phase changes to γ-Ga₂O₃. Clearly, there are many questions to be explored with direct experiment and theory in the diffusion and annealing behavior of native defects and impurities in Ga₂O₃.

In this work, we investigated the diffusion of V_Ga using specially designed superlattices that allow gradients in the extent of diffusion and thus history of [V_Ga] to be inferred. V_Ga are available in great numbers in Sn-doped substrates, transiently diffusing outward toward the epilayers in our experiments. This allowed us to measure D_VGa, which we found was significantly affected by its co-diffusion with Sn donors. We developed a novel finite difference scheme, in which the diffusion of impurities such as Al is dependent on the local concentration of V_Ga. We find no evidence for Ga vacancies introduced from the surface significantly affecting the Al diffusion, although these may be introduced at lower concentrations. The diffusion rates extracted for V_Ga are far lower than those predicted for free vacancies based on DFT calculations. Instead, they are retarded by the strong Coulomb attraction to donors such as Sn and Si highlighting that that doping will have indirect effects on Ga self-diffusion and on isovalent Al diffusion. Diffusion and surface accumulation of unintentionally introduced impurities such as Fe, Mn, Cu, and Li bear further exploration, including the effects of surface potential and phase changes.

See the supplementary material for the derivation of the piecewise linear approximation of the reaction diffusion of SnGa–VGa complexes and quantitative comparison with the numerical solutions described in Ref. 24. The supplementary material also contains a discussion reconciling an apparent contradiction in a priori assumptions on the conservation or non-conservation of VGa in our modeling and that described in Ref. 24, and an analysis of the effects of lateral non-uniformity in d-SIMS profiles. It also describes our efforts to validate the apparent accumulation of Fe at the surface during annealing and eliminate the possibility that this is either a SIMS artifact or a result of surface contamination.

This work was funded by the Air Force Office of Scientific Research—Multidisciplinary Research Program of the University Research Initiative (MURI)—AFOSR under Award Nos. FA9550-18-1-0507 and FA9550-21-0078. We also thank Northrop Grumman-Synoptics and John Blevins for β-Ga₂O₃ boule pieces. We thank K. M. Johansen and Y. K. Frodason for the use of their data and modeling of Sn_Ga and V_Ga reaction–diffusion. We thank D. Feezell and his group for high-resolution x-ray diffraction scans on early samples to help characterize the growth rate variations and F. Guervara-Vasquez for informative discussions on finite difference schemes.

The authors have no conflicts to disclose.

Nathan D. Rock: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Writing – original draft (equal); Writing – review & editing (equal). Haobo Yang: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Brian Eisner: Data curation (supporting); Investigation (supporting). Aviva Levin: Data curation (supporting); Investigation (supporting). Arkka Bhattacharyya: Formal analysis (supporting); Resources (supporting). Sriram Krishnamoorthy: Data curation (supporting); Investigation (supporting); Resources (supporting). Praneeth Ranga: Resources (supporting). Michael A Walker: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Larry Wang: Data curation (supporting). Ming Kit Cheng: Data curation (supporting). Wei Zhao: Data curation (supporting); Investigation (supporting). Michael A. Scarpulla: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.