Utilizing $(Al,Ga)_{2}O_{3}/Ga_{2}O_{3}$ superlattices to measure cation vacancy diffusion and vacancy-concentration-dependent diffusion of Al, Sn, and Fe in $\beta\text{-}Ga_{2}O_{3}$

$\beta{\text{-}}Ga_{2}O_{3}$ is an ultra-wide band gap (UWBG) semiconductor with potential for superior performance in power electronics as well as in solar blind UV photodetectors and in transparent contacts for photovoltaics. Pearton et al. ; Green et al. ; Tsao et al. . To further advance this technology, a thorough understanding of the defects present and their transport is needed including during single crystal growth, epitaxial thin film growth, and during other fabrication steps.

Many groups have observed an increase in resistivity of n-type $\beta{\text{-}}Ga_{2}O_{3}$ after $O_{2}$ annealingOshima et al. (a). This formation of an insulating surface layer allows 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 and thus could not contribute directly to n-type conductivityVarley et al. . Oshima, et al. documented the $\sqrt{t}$ kinetics of the insulating layer thickness formed during $O_{2}$ annealing and demonstrated that an $SiO_{2}$ cap could prevent itOshima et al. (b). This indicates that direct exposure of the $Ga_{2}O_{3}$ 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 filmsTuomisto and Makkonen ; Korhonen et al. , and that increases in their numbers and or charge states accompany $O_{2}$ annealing , which makes them a strong candidate for n-type compensationIngebrigtsen et al. (a); Jesenovec et al. . Many mechanisms, involving oxygen, $V_{Ga}$, $V_{O}$, and n-type dopants (Sn and Si) have been proposed, including those which involve defect complexes instead of bare $V_{Ga}$ and $V_{O}$Frodason et al. (a); Karjalainen et al. ; Kyrtsos, Matsubara, and Bellotti ; Sun et al. . Studying the diffusion of $V_{Ga}$ should elucidate what mechanisms were at play, however vacancies are much more difficult to observe directly than e.g. 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.De Souza and Kilner ; Töpfer, Aggarwal, and Dieckmann ; Howard ; Klugkist and Herzig ; Mehrer ; Tsao et al. . Modern computation and experiment have established Ga vacancies and their complexes as the most plentiful native defects in n-type $\beta{\text{-}}Ga_{2}O_{3}$. Substitutional isovalent alloying cations in oxides having ionic radii not significantly less than the host atom – e.g. Al in $Ga_{2}O_{3}$ – should exhibit primarily vacancy-mediated mechanisms under oxygen-rich conditions favoring cation vacancies and suppressing cation interstitials. Interstitial mediated mechanisms may be possible under metal-rich conditions favoring cation interstitials. The local diffusion constant of isotopic or chemical tracers is proportional to the local concentration of mediating native defects. Thus for vacancy-hopping the tracer diffusion “constant” can vary in space and time with the concentration of mediating vacancies and different crystals may exhibit different diffusion constants by virtue of their vacancy content. This in turn is affected by dopant impurities through Fermi level effects and their thermal and chemical history. For dopant ions that control the local Fermi level, local equilibrium of native defects can result in that dopant’s diffusion constant being dependent on its own concentration Johansen et al. ; Frodason et al. (b); Bracht ; Gösele and Tan . The most general microscopic formulation of tracer diffusion problems is coupled reaction-diffusion of defects and complexes wherein the different tracer-native defect complexes have different diffusion coefficients than the isolated speciesHobbs and Ord ; Krasikov and Sankin . The fraction of time, or population fraction, of complexes relates to the binding energy and connects to the correlation coefficient of random walk theory Mehrer . Interestingly, our assumption that the Al impurity diffusion constant in Fick’s laws is proportional to vacancy concentration is not on its face identical to introducing non-zero but constant off-diagonal coefficients within Onsager’s transport formalismOnsager .

Herein, isothermally anneal superlattices of isovalent Al tracer concentration spikes in $\beta$-Ga2O3 to investigate the interdiffusion of Al, diffusion of cation vacancies, and diffusion of Fe or Sn impurities. Commonly it is assumed (and in careful experiments, ensured) that diffusion-mediating native defects equilibriate with local impurity concentrations rapidly compared to the experimental timescale. Rapid equilibration can be assumed in samples containing dense sources and sinks, however in high-perfection crystals like those in this study, equilibration may be much slower and 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 impurity diffusion constants.

Herein, Al diffusion in superlattices (SLs) epitaxially grown on Fe substrates occurs with at least 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-Nicholsen finite differences 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 and the $D_{V_{Ga}}$ 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 $D_{V_{Ga}}$ and Ga self-diffusion constant for single crystalline $\beta\text{-}Ga_{2}O_{3}$. 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 surreptitiously discovered that many other cation species such as Fe, Mn, Ni, Li, and Cu also diffused from the Sn-doped substrates and accumulated near the free surface.

Superlattices (SLs) consisting of 6 or 9 alternating layers of 200 nm unintentionally doped (UID) $\beta{\text{-}}Ga_{2}O_{3}$ and 10-15 nm $Al_{2x}Ga_{2(1-x)}O_{3}$ (AlGO) where $x\cong 0.05$. were epitaxially grown using organometallic vapor phase epitaxy (OMVPE) in an Agnitron Agilis 100 system on (010) Sn or Fe doped $\beta{\text{-}}Ga_{2}O_{3}$ substrates obtained from Novel Crystal Technology (NCT). The vendor reports Fe or Sn concentrations in the mid $10^{18}cm^{-3}$, and etch pits from nanovoids and dislocations are $<10^{5}/cm^{2}$. The lowest practical Al concentration of $\approx 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_{2}O_{3}$ spacers. UID layers grown under the same conditions exhibit n-type doping in the low $10^{16}/cm^{3}$ range. Annealing was carried out for times between 2 and 80 hours with samples face up in a quartz tube furnace in 1 atm flowing $O_{2}$ after acid cleaning the tube and the quartz or sapphire plates on which samples sat. During annealing, temperature was ramped at $10^{\circ}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. 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 Pendellosung fringes allowed us to quantify the uncertainty in depth scales, resulting primarily from growth rate variation, of at most $\Delta x=\pm 6\%$. This translates to relative uncertainty in extracted diffusion constants of order $(\Delta x)^{2}$=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.

Isovalent aluminum and gallium have similar ionic radii, therefore we assume a substitutional, vacancy mediated diffusion mechanism for Al. Somewhat surprisingly, numerical treatment of the case of non-steady-state vacancy-mediated diffusion, which amounts to assuming that the tracer diffusion “constant” varies in space and time, is not widely considered. Thus, a novel numerical approach based on the Crank-Nicholson scheme ($2^{n}d$ order centered finite differences in space, forward and linear in time) was developed. The local Al diffusion constant $D_{Al}$ is assumed to be Arrhenius activated and proportional to the local concentration of cation vacancies (shorthanded as $[V_{Ga}]$ in $\#/cm^{3}$ acknowledging the small Al content) and thus takes the form

in which $D_{o}(cm^{2}/s)$ is the conventional prefactor, $D^{*}(cm^{5}/s)$ is a constant for all temperatures, $E_{a}$ is the Al hopping activation energy given a neighboring vacancy, $k_{B}$ is Boltzmann’s constant, T is absolute temperature, and $D_{oo}(T)=D^{*}exp\left(\frac{-E_{a}}{k_{B}T}\right)$ is a constant at each temperature. Fick’s $2^{n}d$ law for isothermal tracer Al diffusion becomes

Not assuming spatial uniformity of the diffusion constant introduces the cross-derivative term involving gradients of both [Al] and $[V_{Ga}]$. The $2^{nd}$ order space central difference, first order forward time discretization of this equation on a uniform grid of spacing $\Delta x$ is

in which the space and time grids are indexed respectively by the subscript (j) and superscript (N) and $\Delta x$ and $\Delta t$ are the x and t intervals. Dirichlet and von-Neuman boundary conditions were applied for Al at the deepest points in the substrate and surface respectively. To simplify simulation, which we implemented in Matlabnoa , we used analytical solutions for diffusion of $[V_{Ga}]$ which is tantamount to assuming vacancy diffusion does not depend on Al concentration. This assumption was justified by the chemical similarity of Al and Ga and the small maximum Al concentration. We find the best balance of assumptions and distinguishable differences in modeled profiles by assuming depth-constant (but distinct) initial $[V_{Ga}]$ conditions in the OMVPE grown layers and the substrates. Because of the high dynamic range of the SIMS data with low uncertainty (meaning that both the peaks and valleys of the Al profiles are measured with good signal to noise ratios), we chose to use the square difference of $log_{10}([Al])$ between measured and simulated profiles (summed over the J points in each profile) as our goodness of fit (GOF) parameter:

$D_{V_{Ga}}$, the product$[V_{Ga}]_{subs}\cdot Doo$ and the ratio $[V_{Ga}]_{epilayer}/[V_{Ga}]_{subs}$ are the 3 free parameters (degrees of freedom) of each scenario directly-extractable from the simulation/fitting and robustly determined at each temperature. 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_{V{Ga}}$ at each T. Any difference between actual and assumed $[V_{Ga}]_{subs}$ will affect $[V_{Ga}]_{epilayer}$ and $D_{oo}$; if the actual concentration is later found to be 10x higher than assumed in our analysis, then Doo would actually be 10x lower but $[V_{Ga}]_{epilayer}$ also 10x 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 Supplemental Materials.

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. Intuition for a completely ionic crystal suggests $V_{Ga}$ should exist only in the $3^{-}$ charge state and thus have concentration scaling $\propto(\frac{n}{n_{i}})^{3}$ or $\propto exp\left(\frac{-3E_{f}}{k_{B}T}\right)$Jacobs . Modern density functional theory (DFT) with hybrid functionals accurately predicts that defect chargestates may change if $E_{f}$ moves above or below charge transition levels and can account for complexation with other defectsVarley et al. ; Varley . 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. $p_{O_{2}}$) and temperatures during bulk crystal and OMVPE growth will modify this ratio.

A number of experiments havevon Bardeleben et al. ; Jesenovec et al. ; Tuomisto and Makkonen ; Tuomisto ; Karjalainen et al. ; Korhonen et al. 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 spectrscopy (DLOS). DLOS measurements on different but very similar OMVPE layers report a concentration of $\approx 10^{15}/cm^{3}$ for the $E_{c}$-2.0 eV defectsGhadi et al. , while DLOS measurements of $E_{c}-2.0$ eV defects in NCT Sn-doped bulk crystals are reported at $0.8-1.2x10^{16}/cm^{3}$ for wafers matching the nominal concentration of those used herein, $[Sn]=3.5-5x10^{18}/cm^{3}$Johnson et al. . For analysis we assume $[V_{Ga}]=10^{16}/cm^{3}$ vacancies in the Sn-doped substrates.

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^{18}/cm^{3}$ range indicates that $[V_{Ga}]$ can not be much larger. Unfortunately, junction capacitance based defect spectroscopies like DLTS or DLOS can not be performed on insulating samples e.g. Fe-doped wafers or oxygen annealed ones. All substrates were grown at $T_{melt}(k_{B}T\approx 0.18eV)$ by EFG in reduced $p_{O_{2}}$ conditionsGalazka et al. ; Kuramata et al. . In Fe-doped $Ga_{2}O_{3}$, the Fermi level is pinned at the $Fe^{3^{+}/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 techniquesZhang et al. ; Ingebrigtsen et al. (b). Sn-doped substrates are degenerately doped or very close to it, thus $E_{f}\approx E_{c}$. For both of these $E_{f}$ values, isolated $V_{Ga}$ should be in the $q=3^{-}$ charge state; thus $[V_{Ga}^{3-}]$ in the Sn-doped substrates should be higher than in the Fe-doped by a factor of $exp\left(\frac{3\times 0.8eV}{0.18eV}\right)\approx 10^{6}$ if $[V_{Ga}]$ in both substrate types equilibriate at $T_{1}$=$T_{melt}$. If $V_{Ga}$ stay equilibrated to a $T_{1}$ temperature lower than the freezing temperature, this ratio would increase, while if most of the $V_{Ga}$ are in $(Sn_{Ga}-V_{Ga})^{2-}$ complexes or $(V_{Ga}-V_{O})^{-}$ divacancy complexes, this ratio would be approximately $10^{4}$ or $100$ respectively. While the unintentional doping in our OMVPE layers grown at $800^{\circ}C$ is in the low $10^{16}$, 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. Also, in separate experiments we have carried out, we find that annealing Sn-doped wafers in 1 atm $O_{2}$ at $1050^{\circ}C$ requires approximately two weeks for full equilibration; thus $[V_{Ga}]$ will not equilibrate in the full wafer during OMVPE growth. Thus, we consider ratio of $[V_{Ga}]$ in the SL to that in the substrate to be an unknown parameter to be determined in the fitting. Note that the substrate may not equilibrate its $[V_{Ga}]$ during OMVPE growth, so the ratio between substrate and SL may not have this value. 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.

Figure 1a,b shows Al and dopant SIMS profiles before and after annealing for 20 hours in 1 atm $O_{2}$ at $1050^{\circ}C$ for two identical superlattices grown 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 $AlGa_{2}O_{3}$ layers in the superlattices. Meanwhile, for SLs on Sn-doped substrates the extent of Al diffusion is greater for AlGaO layers closer to the substrate than the surface, and overall the Al diffusion is much faster. This is consistent with a larger $[V_{G}a]$ in the Sn-doped substrates compared to the SLs and to that in the Fe-doped substrates which diffuses through the SLs towards 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 at concentrations sufficient to affect Al diffusion from about 5 at.% down to $10^{18}/cm^{3}$. Such an effect could be taking place, but apparently at concentrations significantly lower than the estimated $10^{16}/cm^{3}$ 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.

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 2nd law. We tried many scenarios of different [VGa] in SL and substrate but the best agreement with experiments was found for the case of the same, low concentrations. Evolving the initial to final Al SIMS profiles for the sample annealed at $1050^{\circ}C$ allows us to extract a value of $D_{Al}=D_{oo}(T)\cdot[VGa]=3.5x10^{-18}cm^{2}/s$, assuming $[V_{Ga}]=10^{14}/cm^{3}$, i.e., $\approx 100$ times less than the substrate value used for Sn doped substrates.

The diffusion of Fe itself appears to be complex, probably involving at least two modes such as both interstitial and substitutional. The difference between the integrated Fe concentration in the epi-layer, before and after annealing, is $\approx 8.4x10^{13}/cm^{2}$. Of this, $8.0x10^{13}/cm^{2}$ accumulated near the surface, which we have shown to be neither a SIMS artifact or a result of vapor transport during annealing (See supplemental for details). This effect was reproduced via dynamic and Time of Flight (TOF) SIMS, across multiple samples, and is also reproduced on Sn doped substrates where [Fe] is initially considerably lower in the substrate. The remaining $4x10^{12}/cm^{2}$ Fe that diffused only short distances into the SL region apparently diffuses via a second mode. Assuming this corresponds to diffusion of Fe into a long half-space without complexation yields the analytical solution to Fick’s $2^{nd}$ lawCrank :

in which $C_{subs},C_{SL}$ are the initial [Fe] concentrations in the substrate and SL respectively, $x_{L}$ is the depth of the SL/substrate interface, and $D$ is the diffusion constant for free Fe. The best fit of this equation, excluding the [Fe] near the surface, yields a diffusion constant of $D_{Fe}=6x10^{-16}cm^{2}/s$. The diffusion of Fe in $Ga_{2}O_{3}$ certainly warrants additional experimental and theoretical scrutiny, including Fermi level effects. Fe is present at appreciable concentrations in most if not all $Ga_{2}O_{3}$ bulk crystals as an impurity introduced from raw materials or Ir crucibles and dies McCloy et al. .

Figure 2 a-c) present initial and final SIMS profiles for [Al] for samples annealed at $1000$, $1050$, and $1100^{\circ}C$ for 20 hours. 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 C was grown in a different run; this is important for understanding the apparently lesser degree of Al diffusion despite this sample’s higher annealing temperature, as discussed in more detail below. Figure 2 d and e reproduce the data from a) but also overlay simulated Al concentration profiles at 2 hour 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 the final SIMS profiles. Each tested scenario is defined by Eqs. 1-2, assumed initial $[V_{Ga}]$ in SL and substrate, $[V_{Ga}]$ diffusion behavior explained below, and values for $D_{oo}$ for Al and $D_{V_{Ga}}$. 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 $D_{V_{Ga}}$ 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 Figure 2 d), we first tested the hypotheses that a) $V_{Ga}$ are conserved and b) that they diffuse freely without interactions from the Sn-doped substrate through the SLs according to Eq. 5Crank :

in which $C_{subs},C_{SL}$ are the initial vacancy concentrations in the substrate and SL respectively, $x_{L}$ is the depth of the SL/substrate interface, $D$ is the diffusion constant for free vacancies, and t is the isothermal annealing time. In the supplemental materials we explore the issue of $V_{Ga}$ conservation and implications for assumptions that vacancies equilibrate in other single crystal diffusion experiments. Our simulations reproduced the Al SIMS data fairly well indicating that these two assumptions capture the gross features of $V_{Ga}$ diffusion in our experiments; no generation of $V_{Ga}$ in the SLs or introduction from the free surface are indicated at levels detectable in these experiments.

However, theoretical and experimental evidenceFrodason et al. (b); Varley suggest that $V_{Ga}$ in Sn doped substrates dynamically complex with a small fraction of the available Sn ([$V_{Ga}$]<<[$Sn_{Ga}$]). Complexed $V_{Ga}$ diffuse slower than free $V_{Ga}$, while $Sn_{Ga}$ diffuse only while complexed or correlated with a vacancy Mehrer . 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. [Sn] is roughly 100-1000 times higher than $[V_{Ga}]$ in the substrate, which ensures that the majority of $V_{Ga}$ is on average bound in $Sn_{Ga}-V_{Ga}$ complexes. However, the small equilibrium rate constant ensures that individual $Sn_{Ga}-V_{Ga}$ complexes are constantly dissociating and reforming. Therefore, $V_{Ga}$ in the vicinity of Al has a high likelihood of assisting in Al diffusion even when the $V_{Ga}$ is diffusing as a complex. Because any ionized donor will have Coulomb attraction to cation vacancies, the effective diffusion constant for $V_{Ga}$ and thus for cation self-diffusion in the presence of donors such as Sn should be smaller the value in intrinsic or acceptor-doped $\beta\text{-}Ga_{2}O_{3}$ (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.

Figure 2 e) shows fitting of one data set using a modified assumption for the evolution of $[V_{Ga}](x,t)$, taking into account its complexation with $Sn_{Ga}$. Frodason et al.Frodason et al. (b) studied the diffusion of Sn from (001) Sn-doped substrates into UID HVPE-grown $Ga_{2}O_{3}$ for temperatures $1050$ to $1250^{\circ}C$ and fit them with a diffusion-reaction model. A key assumption of their work was that $[V_{Ga}]$ equilibrated quickly by unspecified generation or transport mechanisms with the local Fermi level set by the local doping – i.e. that $\int[V_{Ga}]dx$ was not conserved.

In Ref. Frodason et al. (b) (and in some similar intermediate Sn d-SIMS profiles obtained in this work but not shown herein) Sn appears on a log scale to diffuse as a nearly-rigid step with maximum slope set by the temperature but not time. In the Supplemental Materials we develop an analytical approximation for the x and t evolution of both [$Sn_{Ga}$] and [$V_{Ga}$] within for this R-D process during isothermal annealing. Remarkably, this is simply a piece-wise-linear function which is initially a step function but in which the initially vertical section pivots vs. time such that its intercepts with the max and min values far away front the interface move as $\sqrt{Dt}$ (see Eq. 6). This piece-wise-linear function accurately captures most of the dynamic range of the $Sn_{Ga}-V{Ga}$ steps from Ref. Frodason et al. (b) and our re-analysis of those data. It only fails to capture some curvature at low concentrations at the diffusion front leading edge and at high concentrations in the substrate. Since Al diffusion scales with $[V_{Ga}]$, the details of the leading edge are not dominant. The difference between the classical $erf()$ solution may be seen by comparing the dashed green lines in Fig. 2 d) and e). Essentially the best-fit model for Al diffusion assuming $erf()$ diffusion for $[V_{Ga}]$ suggests a significant gradient of $[V_{Ga}]$ (and thus $Sn_{Ga}$) should remain at the end of our experiments, while the RD model and our approximation predict the $[V_{Ga}]$ and [Sn] should be closer to uniform at the end of the experiment. The later is in agreement with the measured, nearly-uniform final [Sn] profiles as exemplified by Fig. 1b). Additionally, using the piecewise-linear RD approximation yields an order of magnitude better GOF values and allows us to recover Arrhenius diffusion behavior for both Al and $V_{Ga}$. We take these facts to be strong indicators that the piecewise RD approximation for $V_{Ga}$ diffusion more accurately reflects the actual $V_{Ga}$ diffusion in these SLs. The magnitude of of $D_{V_{Ga}}$ extracted from under the RD model show remarkable agreement with $D_{Sn}$ in Frodason et al. (b) despite $D_{V_{Ga}}$ being extracted from our [Al] profiles and $D_{Sn}$ in theirs from [Sn] profiles, 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 is the rate limiting step(s).

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 $p_{O_{2}}$ during annealing and growth? The creation and transport of vacancies bears further examination in $Ga_{2}O_{3}$ and in other single crystal diffusion experiments. However, for the immediate purpose of analyzing our Al diffusion data, this question is immaterial since the RD models predict the free $[V_{Ga}]$ able to facilitate Al diffusion to have functional shape closely mimicking the Sn profile.

The Supplemental Materials contain a thorough discussion of whether or not $V_{Ga}$ formation in the bulk or transport from the free surface (where formation presumably has lower activation barrier) is detectable in our experiments and those of Ref. Frodason et al. (b). Succinctly, the lower extent of diffusion of Al close to the free surface appears to rule out large concentrations of $[V_{Ga}]$ being injected from that surface (“large” here being relative to  at% concentrations of Al). Also, our models account for the Al diffusion without requiring generation of $V_{Ga}$ in the SLs, only the diffusion of those initially present in the substrate. In Ref. Frodason et al. (b) such mechanisms are implicit requirements of the assumption that $[V_{Ga}]$ equilibrates with the local Fermi level. Since Sn in Frodason et al. (b) concentrations range from $10^{14}$ to $10^{18}/cm^{3}$, these experiments are sensitive to much lower concentrations of $V_{Ga}$. Thus we can only conclude that any bulk generation or injection of $V_{Ga}$ from the free surface must be below the capability of the present experiments with Al to resolve and further experiments with smaller tracer concentrations, longer times, or higher temperatures should be used to elucidate the kinetics of $V_{Ga}$ introduction.

From a visual comparison alone, the sample annealed at $1050^{\circ}C$ (Fig 2b)) shows a greater extent of diffusion than the one annealed at $1100^{\circ}C$ (Fig 2c)), 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^{\circ}C$ (Fig 2c) ) was grown in an earlier growth run using a different substrate, while the samples annealed at $1000$ and $1050^{\circ}C$ (Fig 2a,b) ) were grown at a later date on pieces of a new Sn-doped (010) $\beta\text{-}Ga_{2}O_{3}$ wafer different than the one for the $1100^{\circ}C$ one. $[V_{Ga}]$ is likely to vary spatially within EFG-grown bulk crystals, thus variation from wafer to wafer may occur. 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}]$. Thus in our final analysis, we assumed $[V_{Ga}]=10^{1}6/cm^{3}$ in the substrates for the $1000$ and $1050^{\circ}C$ annealed samples, but explored possibilities that the substrate for the $1000^{\circ}C$ sample had a 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 smaller value of $D_{Al}=D_{oo}\cdot[V_{Ga}]$ for the sample annealed at $1100^{\circ}C$ is best explained by lower initial $[V_{Ga}]$ in that sample’s substrate compared to the samples (from a different OMVPE growth and cut from a different wafer) annealed at $1000$ and $1050^{\circ}C$. After accounting for this difference in $[V_{Ga}]$, and rerunning simulations with initial $[V_{Ga}]$ in the substrate reduced from $10^{16}/cm^{3}$ to $2x10^{15}/cm^{3}$, more sensible behavior for $D_{oo}(T)$ compatible with an Arrhenius law with activation energy $3.9\pm 0.3$ eV is obtained. This emphasizes the extreme sensitivity to the local $[V_{Ga}]$ concentration for self- and impurity-diffusion, whilst also demonstrating the efficacy of our modeling at extracting physically-reasonable behavior even from samples having heterogeneous initial conditions. The 3.9 eV activation energy for $D_{oo}(T)$ is comparable to the value of 4.2 eV obtained for $Al_{Ga}$ diffusion in (-201) $\beta{\text{-}}Ga_{2}O_{3}$ in ref Uhlendorf and Schmidt . Table 1 summarizes all of the parameters from best-fit scenarios under both the $erf()$ and piece-wise-linear RD models for the $V_{Ga}$ diffusion (our final analysis concludes that the later is the correct model).

In Fig. 1a,b) we compare two SLs annealed at $1050^{\circ}C$ on Fe and Sn doped substrates. On the Sn doped substrate, the best-fit scenario within the piece-wise-linear model of $V_{Ga}$ diffusion yields $D_{oo}=2.7x10^{-30}cm^{5}/s$ at an initial value of $[V_{Ga}]=10^{16}/cm^{3}$ 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- or Sn-doped substrates, we estimate the initial $[V_{Ga}]$ in the substrate for samples grown on Fe-doped substrates would be $10^{14}/cm^{3}$, o r $\approx 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.

Returning to the issue of vacancy conservation, our data and analyses point to these experiments lasting on the order of a day or more at $1000-1100^{\circ}C$ $V_{Ga}$ diffusion occurring over a distance of 1.5-2.5 um within a transient, not equilibrated regime. The data show no indication for in-diffusion from the surface, the density of vacancy sources like dislocations is small ($10^{4}-10^{5}/cm^{2}$), spontaneous vacancy pair formation has a large activation energy. These are interesting and impactful results as most analyses of similar diffusion experiments assume that mediating native defects equilibrate quickly compared to impurity or tracer diffusion timescales. Herein, our evidence indicates cation vacancies remaining in a transient diffusion regime over the 1.5-2.5 $\mu m$ SL thickness for most or all of experiments at $1000-1100^{\circ}C$ typically lasting about a full day (at least for concentrations necessary for diffusion of Al at $10^{18}-10^{21}/cm^{3}$ concentrations). This may be a unique experimental case in that the homoepitaxial samples have low densities of vacancy sources/sinks (c.f. grain boundaries or dislocations) and our diffusing species were not introduced as reactive surface layers nor were they ion implanted (both of which can introduce or consume mediating native defects). Higher temperatures and different chemical boundary conditions at free surfaces may change this, but these experiments provide a valuable existence proof of cases where this typically foundational assumption of analysis of native defect-mediated diffusion is not fulfilled.

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). The revealed that many other elements besides Fe or Sn diffused out of a Sn-doped substrate and through the SLs. Figure 3 shows concentration-calibrated SIMS profiles for five illustrative metals present in the SL shown in Fig. 2b), annealed at $1000^{\circ}C$. After collecting the 20 hour anneal SIMS profile shown in Fig. 2b), the sample was immersed in aqua regia for 10 minutes to remove any possibility of metal contamination at the surface, and then annealed in $O_{2}$ for an additional 10 hours (total 30 hours) as part of efforts to eliminate possibilities of Fe contamination depositing on the surface and then diffusing into the SL while annealing inside the tube furnace (Supplemental Materials). Given the known presence of many of these elements in wafers grown from melts held in Ir crucibles- e.g. EFGMcCloy et al. ; Kuramata et al. 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 at concentrations, and to depths, which are unlikely to be due to artifacts. [Cu] exhibits a concentration peak around the middle of the superlattice structure, perhaps indicating interactions between its charge states and the Fermi level within a space charge region. All of these transition metals are likely capable of multiple interstitial and substitutional diffusion modes in multiple charge states. Ti was also detected at $10^{16}/cm^{3}$ 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 electro-chemical potential gradient affecting different elements and charge states differently near the surface. A depletion width and built-in potential $\phi(x)$ induced by surface Fermi-level pinning is the simplest example and defects with different relative charges q would accumulate or deplete as $exp\left(\frac{q\phi\left(x\right)}{k_{B}T}\right)$ self-consistently following the Poisson-Boltzmann equation. However, it has been frequently observed that $\beta\text{-}Ga_{2}O_{3}$ surfaces, especially in the presence of additional impurities, can change to the different $\gamma\text{-}Ga_{2}O_{3}$ 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 diffusion of these impurities, some like Cu and Li being notoriously fast diffusers in semiconductors, calls for further quantification of impurities and their diffusion in multiple phases of $Ga_{2}O_{3}$ and in devices from fabrication to long-term degradation under high fields and temperatures.

Figure 4 compares measured and computed diffusivities in $Ga_{2}O_{3}$ from this work and that from literature. Red pentagons and triangles give the diffusivities inferred in this work for $V_{Ga}$ assuming no interaction with Sn (erf) and coupled reaction-diffusion with Sn (piece-wise linear). $D_{V_{Ga}}$ extracted from our Al diffusion data over the investigated temperature range assuming only the piece-wise linear shape of the $[V_{Ga}](x,t)$ evolution (making no assumptions of diffusion constant) have remarkable agreement with the diffusion of $Sn_{Ga}-V_{Ga}$ complexes reported by FrodasonFrodason et al. (b) and follow a sensible Arrhenius trend with activation energy $1.86\pm 1.13eV$ The values extracted from best-fits at different temperatures assuming free diffusion ($erf(x,t)$) are an order of magnitude lower and deviate from Arrhenius behavior. Blue rhombuses plot the modeled values of free diffusion of uncomplexed $D_{V_{Ga}}$ from Ref Frodason et al. (b). It is clear that interactions of acceptor-like $V_{G}a$ especially with donor-like defects significantly slows their diffusion.

Since Al diffusion is proportional to $[V_{Ga}]$ and this varies with depth and time in the SLs on Sn-doped substrates, we present $D_{Al}$ extracted for each temperature, assuming a uniform concentration local $[V_{Ga}]=10^{16}/cm^{3}$. To emphasize the dependence of $D_{Al}$ on local $[V_{Ga}]$ we also give a second dataset for $D_{Al}$ assuming a uniform local $[VGa]=10^{14}/cm^{3}$. This latter set gives good agreement in magnitude with $D_{Al}$ derived from the interdiffusion of a thin $Al_{2}O_{3}$ layer on $Ga_{2}O_{3}$ during oxygen annealingUhlendorf and Schmidt . The data shown as an inverted brown triangle is the value for the extracted diffusion constants for the slow diffusion channel of Fe out of an Fe-doped substrate. This value is similar to those experienced by Al, so we speculate 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 $\mu m$ thickness during $O_{2}$ annealing. The kinetics of layer thickening were consistent with a diffusion-limited process, thus we plot their equivalent diffusion constant as blue circles. It is very interesting to note that these extracted diffusion constants are much closer to the values for impurities like Al or Fe diffusing to the surface in our samples than they are to any of the cation vacancy diffusion constants or oxygen tracer diffusion constants recently reported. Thus, we speculate that at least part of the phenomenology of insulating layers forming on the surface of n-type $Ga_{2}O_{3}$ crystals may be caused by out-diffusion of unintentional cation impurities towards the O-rich boundary condition, possible surface band bending, and possibly phase changes to g-Ga2O3 with impurity pile-up, rather than directly by in-diffusion of O or Ga native defects as the simplest models might assume. 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_{2}O_{3}$. Our results point out many places where unchallenged assumptions regarding native defect processes in diffusion experiments should be revisited and directly tested.

While many papers have shown a correlation between gallium vacancies and n-type carrier concentration after oxygen annealing, the exact dynamics of this process have not previously been explored. We have shown via Al tracer diffusion, that $V_{Ga}$ are available in great numbers in Sn doped substrates, diffusing outward towards the epi-layers in our SL experiments. However, the rate of diffusion is much faster than that indicated for the formation of the semi-insulating layer measured by Oshima. We developed a modified finite difference scheme in which the diffusion of impurities such as Al is dependent on the concentration of $V_{Ga}$, which is also diffusing. The finite differences model herein developed shows no evidence of Ga vacancies being introduced from the surface, though these may be introduced at lower concentrations. The measured diffusion rates extracted for $V_{Ga}$ are far lower than those predicted for free vacancies based on DFT calculations. Instead they agree well with the measured rates of Sn diffusion, indicating that the formation and diffusion of defect complexes is the primary rate controlling step in the diffusion of the charge-neutral Al defects. The study of unintentionally introduced impurities like Fe, Mn, Cu must be explored further, and may yet play an important role in elucidating the formation of the semi-insulating layer. The accumulation of transition metals at the surface suggests that the effect of the surface potential must be taken into account in future work.

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