Since the 1940s, it has been known that diffusion in crystalline solids occurs due to lattice defects. The diffusion of defects can have a great impact on the processing and heat treatment of materials as the microstructural changes caused by diffusion can influence the material qualities and properties. It is, therefore, vital to be able to control the diffusion. This implies that we need a deep understanding of the interactions between impurities, matrix atoms, and intrinsic defects. The role of density functional theory (DFT) calculations in solid-state diffusion studies has become considerable. The main parameters to obtain in defect diffusion studies with DFT are formation energies, binding energies, and migration barriers. In particular, the utilization of the nudged elastic band and the dimer methods has improved the accuracy of these parameters. In systematic diffusion studies, the combination of experimentally obtained results and theoretical predictions can reveal information about the atomic diffusion processes. The combination of the theoretical predictions and the experimental results gives a unique opportunity to compare parameters found from the different methods and gain knowledge about atomic migration. In this Perspective paper, we present case studies on defect diffusion in wide bandgap semiconductors. The case studies cover examples from the three diffusion models: free diffusion, trap-limited diffusion, and reaction diffusion. We focus on the role of DFT in these studies combined with results obtained with the experimental techniques secondary ion mass spectrometry and deep-level transient spectroscopy combined with diffusion simulations.

## I. INTRODUCTION

The diffusion of defects and impurities in semiconductors plays a fundamental role in the processing, operation, and stability of semiconductor devices. To make semiconductor materials suitable, e.g., electronics, sensors, power devices, and quantum technology, it is essential to understand and control the atomic transport of dopants and impurities and their interaction with defects. By identifying the diffusion mechanisms for different defect species, we can in many cases also reveal information about the defect’s behavior in relation to device operation and obtain fundamental knowledge about the semiconductor material’s intrinsic properties.

First-principles calculations based on density functional theory (DFT) enable us to study diffusion on an atomic level and calculate properties like formation energies, binding energies, migration barriers, and probable pathways. These predictions provide a starting point for developing a detailed diffusion model. Through numerical simulation of such models, the results may be compared to experimental data. In this way, one can get more detailed insight into the different diffusion models, as well as opportunities to benchmark the DFT results.

Wide bandgap (WBG) semiconductors are predicted to be the next-generation materials for energy-efficient power devices.^{1,2} The critical field is shown to increase with the size of the bandgap, and an increased breakdown voltage is particularly useful in power devices that need to be efficient and reliable in high voltage, high temperature, and harsh environment applications.^{3,4} Historically, WBG semiconductors have been challenging to describe using DFT, even qualitatively, because the semilocal approximations to the exchange-correlation energy tend to severely underestimate the electronic bandgap.^{5} Indeed, the discrepancy typically increases with the size of the bandgap. For example, the Perdew–Burke–Ernzerhof (PBE) functional^{6} underestimates the experimental bandgaps of ZnO (3.4 eV^{7}) and $\beta -Ga 2 O 3$ (4.9 eV^{8}) by 80% and 60%, respectively.^{9,10} Another related issue with semilocal functionals is their tendency to delocalize charge, which, for instance, manifests itself in the failure to describe polarons. With more precise methods, and in particular, after the introduction of screened hybrid functionals, WBG semiconductors can now be studied quantitatively. Hybrid functionals mix semilocal approximations to the exchange-correlation energy with a portion of nonlocal Hartree–Fock exchange, the fraction of which is often adjusted empirically to reproduce the experimental band gap,^{11} or to satisfy exact physical constraints such as the generalized Koopmans theorem.^{12,13} Importantly, the formation energy and energy levels of defects relevant to diffusion can now be described with a degree of accuracy that allows calculations to serve as a highly useful and effective aid in diffusion studies.^{11,14–16}

In this Perspective, we aim to give an overview of how one can study the diffusion of defects in WBG semiconductors using a combination of DFT calculations and selected experimental methods. We focus on the DFT-related methods and how they compare to experimental results. The article aims to describe how to calculate diffusion related properties such as defect formation energies, binding energies, and migration barriers, and the implementation of the theoretical results in diffusion modeling of experimental results. There are multiple approaches to studying diffusion experimentally, but we have limited our study to direct approaches to measure diffusion, e.g., secondary ion mass spectrometry (SIMS) and deep-level transient spectroscopy (DLTS), making it straightforward to compare results directly with DFT predictions. As case studies, we apply these approaches to defect and dopant diffusion in the WBG semiconductors 4H-SiC, ZnO, and $\beta -Ga 2 O 3$.

## II. MECHANISMS FOR DEFECT DIFFUSION

^{17}He introduced the diffusion coefficient $D$ in a linear relation between the gradient of the concentration $C$ and the flux of particles. By pairing this law with the equation of continuity, he formulated Fick’s second law, also called the diffusion equation

Defect diffusion can, e.g., occur by defects diffusing freely, depend on intrinsic defects, be limited by solid solubility or trapping, or may also involve a more complex rearrangement of the matrix structure. Understanding the diffusion process on an atomic level involves determining a diffusion mechanism for the migrating species. The most technologically relevant mechanisms for diffusion in WBG semiconductors include the direct and indirect (interstitialcy) interstitial, vacancy, kick-out, and dissociative mechanisms,^{18} see Fig. 1 for a schematic of diffusion mechanisms for substitutional and interstitial impurity atoms in a crystal. When a diffusing species jumps between interstitial sites, the diffusion mechanism is called the direct interstitial mechanism [case (1) in Fig. 1]. The interstitialcy mechanism [case (2) in Fig. 1] also involves the diffusing species in an interstitial site, but the jump to the next interstitial site is via a temporary displacement of matrix atoms. In the vacancy mechanism [case (3) in Fig. 1], the diffusing species requires a vacant lattice site to induce transport. The following migration can either occur by a rotation around the vacancy followed by a swapping of sites^{19} or by a swapping of sites followed by dissociation (an effect where transport of vacancies in one direction will effectively move all substitutional atoms, including the diffusing impurity, in the opposite direction of the vacancy transport). The mechanism will depend on the binding energy between the given impurity and the vacancy. The kick-out [case (4) in Fig. 1] and dissociative [case (5) in Fig. 1] mechanisms are both substitutional–exchange mechanisms where migration occurs interstitially, and the exchange between substitutional and interstitial atoms happens through replacement (kick-out) or formation of a vacancy (dissociative).

In general, several diffusion mechanisms may occur in parallel, however, the mechanism that provides the highest total transport of a given impurity will dominate the experimentally observed diffusion profile. Interestingly, the calculated formation energy values (see methodology described below) provided by DFT calculations for the different impurity configurations can be used as a guide to point toward the dominating diffusion mechanisms. To exemplify this, consider the two following cases:

If the formation energy of an interstitial defect configuration for a given impurity is low (i.e., high equilibrium concentration), the interstitial diffusion mechanism will be a likely candidate. However, it is important to understand that even when the formation energy is found to be relatively high (low solubility), substantial transport may still occur if the barrier for diffusion is low. Thus, to reveal the dominating diffusion mechanism, both the formation energy and the diffusion barrier need to be calculated. As will be discussed below, identifying all relevant diffusion pathways can, however, be a substantial challenge. The combination of high formation energy and low diffusion barrier becomes especially interesting when the diffusing species can form immobile complexes with other defects (in other words, become trapped). In such cases, transport may occur below the sensitivity limit of, e.g., SIMS, which precludes direct observation. However, in the presence of defects that trap the diffusing specie, thus forming immobile defect complexes, the diffusing species can be measured due to pileup in defect-rich regions. Interestingly, this provides a way to indirectly observe the defects acting as a trap for the diffusing species and, hence, observe the spatial distribution of defects that can otherwise be difficult to measure experimentally. This is the typical behavior of, e.g., hydrogen, exemplified below in the case studies on both ZnO and $\beta -Ga 2 O 3$.

In the case where both the formation energy for the interstitial configuration and the diffusion barrier are high and the occupation of a lattice site is more favorable, the substitutional impurity will be effectively immobile. In this case, other mediating defects, intrinsic or extrinsic, will be necessary to induce impurity transport. This is exemplified by the vacancy mechanisms described above.

## III. COMPUTING DEFECT DIFFUSION

As discussed, there are several different mechanisms by which impurities can diffuse in semiconductors. The activation energy for each mechanism depends on the energy cost associated with the formation and the migration of the diffusing impurity as well as the diffusion-mediating intrinsic defect. This is illustrated in Fig. 2(a) using vacancy-mediated diffusion as an example, where the impurity can migrate only when a vacancy exists next to it. The formation energy of this vacancy, in turn, depends on the Fermi-level position, defect charge state and chemical potential. Different mechanisms can thus be active, or even occur in parallel, depending on the conditions.

First-principles defect calculations based on the supercell approach provide a means to calculate thermodynamic quantities such as defect formation energies, defect complex binding energies, and energy barriers for defect migration, as described below and can, therefore, be used to estimate the activation energy for different mechanisms. This has proven useful to understand, which mechanisms are likely to contribute to diffusion.

The calculations referred to above provide barriers for migration of defects at zero temperature, from which thermally activated jump rates can be extracted using transition-state theory,^{20} but a variety of other methods are also available and often more appropriate for kinetic studies, including, for example, molecular dynamics (MD) (see, e.g., Ref. 21 for diffusion of intrinsic defects in Si), kinetic Monte Carlo (MC) models (see, e.g., Refs. 22 and 23), and more recently, coupling of MD simulations with neural network based enhanced sampling (see, e.g., Ref. 24 for diffusion of spin qubits in SiC).

### A. Formation energy

An established formalism is available for computing thermodynamic defect quantities from electronic structure calculations based on DFT and the supercell method; detailed reviews on the method are available in the literature (see, e.g., Refs. 11, 25, and 26). Here, we shall briefly explain how the defect formation energy $ E f$, which is the energy cost for creating an isolated defect in a solid, can be calculated using this method.

^{27,28}

The formation energy is, thus, a variable that depends on the Fermi-level position within the bandgap (for charged defects) and the chemical conditions. It is convenient to plot the formation energy as a function of the Fermi-level position in the bandgap, as shown in Fig. 2(b) for a vacancy $V$, an impurity $X$, and their complex $VX$. The slopes of the linear segments correspond to the charge state, and kinks occur at Fermi-level positions where the formation energies of the defect in two different charge states are equal. Such plots are commonly shown for the different limits of the chemical potential values, e.g., Zn-rich and O-rich conditions for ZnO.

The term $ \Delta q$ is a finite-size correction that removes any spurious interaction between the charged defect and its images under periodic boundary conditions.^{29,30}

### B. Binding energy

^{11},

To assess the thermal stability of the complex, one must also calculate the dissociation energy $ E d$, which is the actual activation barrier that must be surmounted in order for the complex to break up.^{25} $ E d$ can be calculated, often to good approximation, as the sum of the binding energy of the complex and the lowest migration barrier of the two isolated constituent defects. For example, in Fig. 2(a), the $VX$ dissociation energy ( $ E d(VX)$) is given by the $V$ migration barrier ( $ E m(V)$) plus the $VX$ binding energy ( $ E b(VX)$).

### C. Migration barrier and path

In addition to the formation energy and the energy differences between distinct stable configurations of isolated defects and complexes, we need to study the kinetics of the atomic hopping processes required for transitions between these stable minima. The transition rates can be explored using transition-state theory, where the initial and final structures are separated by an energy barrier. The problem then amounts to finding the minimum-energy pathway (MEP) for the transition. This can be achieved by performing total energy calculations as a function of the position of the hopping atom, where only the surrounding host atoms are allowed to relax for each position.^{25} The resulting adiabatic total energy as a function of the three spatial coordinates of the hopping atom will then map out a potential energy surface (PES), such as the one illustrated in Fig. 2(c), from which the saddle point and corresponding migration barrier ( $ E m$) can be determined. However, accurate results require a large number of atomic positions to be sampled.

A far more efficient approach is the nudged elastic band (NEB) method,^{31} illustrated in Fig. 2(c), which is an automated search method for the MEP. The typical workflow is as follows: (i) An initial guess for the pathway is provided in the form of a number of structures (or images) between the endpoints. For simple hops over a single barrier, a linear interpolation between the endpoint structures will often suffice as a starting point (care must be taken to ensure that the spacing between any two atoms is not too small); this would be the images along the dashed line on the contour plot in Fig. 2(c). For more complicated processes, e.g., involving bond rotations and metastable minima along the path, it can be sensible to divide the transition into steps. (ii) An initial NEB calculation is performed, which involves a constrained relaxation simultaneously for each image. Specifically, by including artificial spring forces between the images, they are kept equally spaced and prevented from sliding into the initial or final structures. The images are then iteratively driven toward the MEP, shown with a solid line in Fig. 2(c). (iii) Finally, the climbing-image mode of the NEB search is turned on,^{32} where the forces along the reaction path for the highest-energy image is inverted [yellow image in Fig. 2(c)], causing it to climb up to the transition state, and the migration barrier is obtained.

Another approach to automatically search for the saddle point is the dimer method.^{33} Two images or different replicas of the system are needed to form the *dimer* which is then rotated to find the lowest modes and identify the saddle point of the reaction pathway. The two images that form the dimer can either start from the initial state of the transition or be configured closer to the transition state. The dimer method is especially suited for more complex reactions where the final state is not known or for when it is challenging to identify a suitable initial guess of images for NEB calculations. Alternatively, the dimer method can be used to obtain a more accurate estimate of the saddle point starting from an initial NEB calculation with less stringent criteria, e.g., energy cutoff and **k**-mesh.

^{19}Note that this activation energy assumes that the $X$ impurities are already incorporated substitutionally.

### D. Temperature dependency

When the temperature, volume, and pressure dependence are taken into account, the formation energy should be replaced by the Gibbs free energy of formation, which includes vibrational and electronic entropy. These free energy contributions are seldom included, as the vibrational part is computationally demanding, and the electronic part is typically negligible.^{34} For our approach, we are mainly interested in considering the relative energetics of different defects to determine likely diffusion mechanisms. In cases where the calculated formation energy is included in the diffusion modeling, it is merely used as an initial value for a fitting parameter.^{16} As discussed in Sec. IV, quantitative results are improved as more cases are studied within a material system. Moreover, several quantities of interest conserve the chemical species, such as binding energies and migration barriers, making them independent of the chemical potential contribution.^{35} Temperature dependence should always be kept in mind, however, especially in cases where its omission is expected to affect the conclusion. We refer to, e.g., the review by Freysoldt *et al.*^{11} for more discussion on the calculation of free energies of formation. As an example, we also mention the study on H diffusion in Si performed by Gomes *et al.*,^{35} where configurational, vibrational, and rotational (for molecular $ H 2$) degrees of freedom were taken into account to understand the condensation of atomic H into molecules during cooling.

Here, it should also be noted that the zero-point energy contribution to the activation energy can be significant for light elements such as H and D, and especially for diffusion of muonium, which can be regarded as a light H isotope.^{36,37}

In narrower-bandgap semiconductors, the temperature dependence of the Fermi level due to the intrinsic carrier density can be significant, but in WBG semiconductors this effect is often negligible. However, the temperature dependence of the bandgap can usually not be ignored. For example, an experimental bandgap shrinkage of 0.3–0.4 eV is reported for ZnO and $\beta -Ga 2 O 3$ at 800 K.^{38,39} Such data can be included in the diffusion modeling, albeit with an assumption regarding the shifts of the individual valence and conduction band edges.^{40} Alternatively, the temperature dependence of the electronic band structure arising from electron–phonon interactions can be calculated from first-principles (including the shifts of the individual band edges).^{41,42}

### E. Diffusion pre-factor

The diffusivity or diffusion coefficient $D$ of a defect is often modeled according to the Arrhenius equation $D= D 0exp(\u2212 E A/kT)$, where $ E A$ is the activation energy of migration and $ D 0$ is the diffusion pre-factor. In this relation, $ D 0$ is usually assumed to be temperature independent. $ D 0$ contains contributions such as the attempt frequency and the geometry of the system and is an important feature that is quantified when studying defect diffusion experimentally. However, obtaining $ D 0$ from first principles is a rather complex problem, and theoretical diffusion studies often focus on the migration barrier and path alone instead.

In a simple first attempt, the random walk model can be used to obtain an approximate expression for the diffusion pre-factor as $ D 0\u2248 \Gamma 0 d 2$. Here, $ \Gamma 0$ is the attempt frequency, and $d$ is the jump distance. If we approximate $ \Gamma 0$ by a typical phonon frequency of $1\xd7 10 13 s \u2212 1$ and assume a cubic system with a lattice distance of 3.25 Å, we obtain $ D 0\u22481\xd7 10 \u2212 2 cm 2 s \u2212 1$.^{19,43} For more accurate estimates of the diffusion pre-factor, the attempt frequency can be obtained using transition-state theory^{20} and necessitates calculating the vibrational spectrum of the system (see, e.g., Refs. 35 and 44). In the following, we will focus on case studies that do not obtain $ D 0$ using theoretical methods but instead compare theoretically obtained migration barriers and paths to experimental studies and diffusion modeling.

## IV. COMPARISON WITH EXPERIMENTAL MEASUREMENTS THROUGH DIFFUSION MODELING

The physical process of diffusion is described by Fick’s law of diffusion [Eq. (1)]. Traditionally, there are at least two different approaches to modeling diffusion, both using Fick’s law as a starting point (i) to generalize the diffusivity $D$ and allow it to depend on the specific experimental conditions and (ii) to treat $D$ as a material property that depends only on the diffusing species and the host material.

By generalizing the diffusivity, e.g., by treating it as a sum of contributions from different defect configurations, or by effectively making $D$ a function of the Fermi level or the equilibrium concentration of the mediating intrinsic defects, one can recreate the diffusion behavior in a range of complex cases (see, e.g., Ref. 45). A widely used model which is based on this approach is Fair’s vacancy model where $D$ is defined as the sum of the individual diffusivities for each of the specific charge-state configurations taking part in the diffusion model.^{46}

An alternative approach is to treat the diffusivity as a material property dependent only on the chemical nature of the diffusing species and the host material that the species is diffusing through, i.e., where $D$ is independent of the equilibrium concentration of mediating defects, Fermi level, etc. In this approach, the details of diffusion and defect reactions are explicitly described as part of Fick’s law of diffusion [Eq. (1)] by the specific inclusion of defect reactions. The following examples will describe various cases of this approach. The major benefit of such an approach is that the input into the model, and also the resulting output of the diffusion model is restricted to be measurable quantities such as the concentration of a specific defect complex or the vacancy concentration. The equilibrium concentration of the given defect complex or the vacancy concentration can also be calculated using DFT. In this way, one can by diffusion simulations perform benchmark studies by combining state-of-the-art (SoA) DFT calculations with a range of different experimental methodologies to measure independently the different constituents taking part in the diffusion process, as will be described in the case studies below.

A major limitation in all diffusion studies is the fact that given an experimental diffusion profile, e.g., for a specific impurity, one can only extract the total impurity transport. That is, there is typically not enough information in a single diffusion experiment to differentiate between the diffusivity and at which concentration the diffusing species is migrating. For example, one can get a unique number for the total transport or flux of a species, however, not the independent contributions to the total transport. Thus, the same profile may be the result of a fast diffusing species moving at a low concentration or a more slowly diffusing species moving at a higher concentration. However, through the combination of different experimental methods and DFT and by systematically working through the diffusion behavior of a range of different species in a given material system, one can build a self-consistent understanding of that material system. Such an iterative approach resembles the iterative approach utilized to perform the DFT calculations. This is an approach that has proven itself as a highly fruitful method to gain insight into complex systems.

## V. CASE STUDIES

### A. Design of diffusion experiments

#### 1. Introduction of defects and dopants

To study impurity diffusion, one needs in most cases to introduce the impurities into the material. There are numerous approaches to introduce the element of interest. Here, we will consider just a few that are utilized in the case studies. With ion implantation, illustrated in Fig. 3(a), one can control the fluence and the energy of the implanted impurities. The implanted samples are annealed for different time periods and at different temperatures. However, ion implantation will often cause the formation of intrinsic defects in addition to the implanted impurity, which can be both advantageous and challenging. A different approach to introduce impurities is by drive-in diffusion, e.g., by depositing an impurity-doped thin film or by using layered samples with a different impurity concentration, see Fig. 3(b). During annealing, impurities will migrate from the layer with a high concentration into the layer with a lower concentration. In some cases, where the introduced diffusing species has a low evaporation temperature, it is possible to introduce impurities through a gas phase by controlling the ambient. This can be done by using sealed ampules as shown in Fig. 3(c) or by using a gaseous flow during annealing.

#### 2. Techniques to measure diffusion profiles

A multitude of techniques can be used to study the diffusion of defects and impurities. However, in this Perspective, we will focus on direct monitoring of impurity and vacancy diffusion via two techniques, namely, SIMS and DLTS.

A widely used approach to study impurity diffusion is by measuring concentration vs depth of the element of interest using SIMS in samples with different thermal history. Figure 3(d) illustrates the principles behind SIMS, where a primary ion beam is directed toward the sample resulting in the ejection of atoms and ions from the sample surface through a sputter process. The secondary ions are accelerated by an electric field and directed toward a mass spectrometer where ions of different mass over charge ratios are separated by electric and magnetic fields. By scanning the magnetic field, one can construct a mass spectrum to identify the different isotopes originating from the sample. To form a concentration vs depth profile, the elements of interest are monitored as a function of time as the primary beam sputters a crater, see Fig. 3(e). The concentration axis is calibrated by implanting a reference sample with a known fluence while the relation between measurement time and depth is found through the sputter rate. In this way, the SIMS measurement can provide a high precision for the measured concentration vs depth. Furthermore, the SIMS measurement is highly suitable for measuring low impurity concentrations, however, the sensitivity depends on both the matrix and the measured impurity element.

DLTS, on the other hand, detects electrically active defects independent of whether they have an intrinsic or extrinsic origin. This approach is particularly useful for detection of defects that are not measurable with SIMS, such as vacancies, interstitials, and intrinsic defect complexes. A drawback is that one relies on indirect identification of such defects, however, there is a tremendous effort in the community to provide such identifications. As will be shown below, DLTS can be employed to study diffusion of, e.g., intrinsic defects such as vacancies.

DLTS measures the transient capacitance after changing the bias over a Schottky barrier diode or pn-junction. The change in capacitance can be related to the emission of electrons or holes from the defect to the conduction or valence band, respectively. To form a conventional DLTS spectrum [see right panel of Fig. 3(g)], the transient capacitance is measured as a function of temperature, and through an Arrhenius plot, one may extract the activation energy, which is related to the defect’s position in the bandgap, and an apparent capture cross section. However, if an isothermal measurement is conducted at the temperature where the DLTS signal is at its maximum for a specific defect, one may extract the depth distribution of the defect in question [left panel of Fig. 3(g)]. Indeed, this is a powerful approach for studying the diffusion of intrinsic defect levels otherwise difficult to probe, e.g., due to a low defect concentration.

In the following, case studies utilizing both DLTS and SIMS will be discussed and combined with DFT and diffusion modeling to extract information about the diffusivity and diffusion mechanism.

### B. Free diffusion

#### 1. Self-diffusion in SiC

The SIMS measurements enable us to extract the diffusivity but not the mechanism for self-diffusion, i.e., the responsible sublattice and mediating species. For this, DFT calculations are essential and can be used to reveal the activation energy for migration as the sum of the formation energy and migration barriers of the mediating defect species. In the case of C self-diffusion in 4H-SiC, the available defects via which carbon can self-diffuse are the carbon vacancy ( $ V C$), the carbon interstitial (C $ i$), and the carbon antisite ( $ C Si$).^{47} Migration through antisites is unlikely, due to high predicted activation energies, which means that C atoms will migrate strictly on the C sublattice.^{48} Instead, the available pathways are self-diffusion via vacancies injected from the SiC surface during high-temperature annealing, or by Frenkel-pair formation to facilitate self-diffusion via interstitials. Combining theory and experiment, it was shown that the mechanism depends on the specific annealing conditions. Carbon self-diffusion was tentatively assigned to Frenkel-pair formation in the SiC bulk when a carbon rich graphitized photoresist (C-cap) is present during annealing to protect the surface and to $ V C$ injection from the surface in the absence of a C-cap.^{47}

#### 2. Carbon vacancy diffusion in SiC

The carbon vacancy is a prominent defect that is known as an efficient carrier trap^{49} and an important minority carrier lifetime killer^{50,51} in 4H-SiC. In fact, $ V C$ is present in SoA 4H-SiC epitaxial layers to densities around $5\xd7 10 12 cm \u2212 3$,^{52} an amount that is already enough to prevent, e.g., operation of SiC bipolar devices up to 10 kV due to the lifetime reduction. $ V C$ can be thermally generated^{53} and was shown to migrate at temperatures of 1100–1600 $ \xb0$C,^{54,55} precisely in the range used for thermal treatments after implantation doping. However, the migration mechanism was unclear, and the impact of the anisotropic crystal field in 4H-SiC was an open question.

The diffusion of $ V C$ in *n*-type 4H-SiC was studied in Ref. 56 using two different types of epitaxial layers: (i) $c$-cut ( $0001$) and (ii) $a$-cut ( $11 2 \xaf0$). The vacancies were generated by implantation of 4 MeV carbon ions (to avoid contaminating with impurities) with a projected depth of around 2.8 $\mu $m. The presence of $ V C$ can be detected using a Schottky diode structure (or another type of rectifying junction) and by monitoring the $ Z 1 / 2$ peak observed in DLTS spectra and assigned to the negative- $U$ double-electron ( $0$/ $2\u2212$) transition of $ V C$.^{49,57} The curves (dots are experimental data, dashed lines are fits) in Figs. 6(a) and 6(b) show the depth dependence of the $ V C$ concentration before and after C implantation along the $[0001]$ and $[11 2 \xaf0]$ directions, respectively. The profiles were obtained using DLTS depth profiling at 285 K coinciding with the $ Z 1 / 2$ peak maximum.

As shown in Figs. 6(a) and 6(b), high-temperature annealing at 1200–1500 $ \xb0$C for different durations was performed to induce migration of the $ V C$, causing broadening of its depth profile. Evidently, the carbon vacancy migrates faster along the $a$-axis as compared to the $c$-axis in 4H-SiC. Furthermore, diffusion modeling according to Eq. (1) was performed using the fit to the pre-diffusion state [purple curves in Figs. 6(a) and 6(b)] as the initial profile and fitting the calculated profiles to the post-diffusion experimental data [solid pink and yellow curves in Figs. 6(a) and 6(b)] to extract diffusion parameters along the different axes. Arrhenius analysis results in $ D 0 c=0.54 cm 2 s \u2212 1$, $ E A c=4.4 eV$ and $ D 0 a=0.02 cm 2 s \u2212 1$, $ E A a=3.6$ eV. The agreement between the experimental data and Fick’s diffusion law affirms that $ V C$ migrates according to fundamental vacancy migration via vacancy-atom exchange on the carbon sublattice.

To determine the diffusion mechanism for the $ V C$ in 4H-SiC and explain the experimentally observed anisotropy, DFT calculations are needed. Using the NEB method, the importance of the inequivalent lattice sites in 4H-SiC [hexagonal (*h*) and pseudo-cubic (*k*)] for $ V C$ motion could be established. Indeed, as shown by the potential energy surfaces for $ V C$ migration computed using NEB along the $c$- and $a$-axes in Fig. 6(c), axial and basal $ V C$ motion occurs along different pathways. There are two possible jumps along the $c$-axis: $k h \u2032$ and $kh$, see Fig. 6(c), with computed migration barriers of 4.2 and 4.1 eV, respectively. Since both jumps are required for long-range migration of $ V C$, the limiting barrier for $ V C$ motion along the $c$-axis is 4.2 eV. Along the $a$-direction, on the other hand, the two available basal jumps are *hh* and *kk*, making it possible for $ V C$ to travel exclusively within either the hexagonal or pseudo-cubic carbon sublattice plane. The hexagonal plane is found to be preferred due to a lower computed migration barrier of 3.7 eV as compared to 4.0 eV within the *k*-plane.

The computed energy barriers for $ V C$ diffusion in 4H-SiC are thus 4.2 eV for axial migration along the $c$-axis, and 3.7 eV for basal migration along the $a$-axis, in excellent agreement with experimental values of 4.4 and 3.6 eV, respectively. However, DFT calculations are not only useful for discerning the migration barriers; NEB and dimer calculations can also be used to identify the transition-state geometry and explain differences in ease of motion for the relevant defects species. In this particular case, the higher energy barrier (DFT predicts a 4.2 eV activation energy) for axial migration ( $ V C$ jump from *k* to *h’* lattice site) is attributed to the transition state involving a high-energy twofold coordinated C atom. The lower barrier basal migration (DFT predicts a 3.7 eV activation energy for the *hh* $ V C$ jump), on the other hand, it is enabled by a transition state with Si–C bond angles that are the closest to perfect sp $ 3$ bonded structures. Structures showing the transition states can be found in Fig. 9 of Ref. 56.

### C. Trap-limited diffusion

#### 1. Hydrogen diffusion in ZnO and *β*-Ga_{2}O_{3}

Hydrogen is an ubiquitous impurity that strongly influences the properties of solids. This holds true for both ZnO and $\beta -Ga 2 O 3$, two WBG oxides that have received considerable attention for their semiconductor properties. Although interstitial hydrogen ( $ H i$) exhibits amphoteric behavior in most semiconductors,^{58,59} theory predicts that the thermodynamic ( $+$/ $\u2212$) charge-state transition level lies above the CBM in both ZnO and $\beta -Ga 2 O 3$, resulting in shallow-donor behavior.^{60,61} However, $ H i$ is generally quite mobile and reactive^{58}—as exemplified by the low migration barriers calculated for $ H i$ in ZnO and $\beta -Ga 2 O 3$ shown in Fig. 7—and tends to end up in a trapped form.^{61,62} H can, for instance, occupy an O site, resulting in an H $ O$ complex acting as an overall shallow donor.^{61,63,64} As a donor, $ H i$ will also form particularly stable complexes with acceptors such as $ V Ga$ in $\beta -Ga 2 O 3$ and $ V Zn$ in ZnO.^{65–69} Thus, H can act as a marker for other defects, including intrinsic ones that are otherwise challenging to measure chemically, e.g., with SIMS.

Thomas and Lander showed already in 1956 that the solubility of H in ZnO is fairly low.^{73} This is also supported by first-principles calculations,^{9,60,74} which show a relatively high formation energy for H $ i$ in ZnO. Regardless, it has been found that interstitial diffusion is the dominating pathway for H migration in ZnO.^{62,75} In Fig. 7(a), diffusion of implanted deuterium was investigated in three different sets of hydrothermally grown (HT) ZnO samples: (i) as-received, (ii) heat treated at 1100 $ \xb0$C for 1 h, and (iii) heat treated at 1500 $ \xb0$C for 1 h. Deuterium is used instead of hydrogen due to its superior detection limit in SIMS measurements. Interestingly, the three deuterium diffusion profiles are quite different. Interpreted within the framework of trap-limited diffusion, the concentration of traps seems to decrease with the heat treatment. Furthermore, it shows that after 1500 $ \xb0$C for 1 h the trap concentration is reduced below the detection limit for deuterium of $1\xd7 10 16 cm \u2212 3$.

Interestingly, when evaluating the diffusion profiles in the sample heated at 1100 $ \xb0$C with the as-received one, the dissociation rate is much higher in the latter, even though the diffusion is carried out at the same temperature (400 $ \xb0$C). The difference in the dissociation rate can be directly observed as a difference in slope in the diffusion front. This indicates the presence of two different traps for H in the as-received vs heat treated material, however, with similar characteristics.

HT ZnO usually has a high residual Li content, and H is known to form complexes with $ Li Zn$ acceptors in such material.^{76–79} The Li concentration in the as-received sample was measured at $\u223c3\xd7 10 17 cm \u2212 3$, which is close to the trap concentration. After 1 h at 1100 $ \xb0$C the Li content was lowered in the region of diffusion by approximately one order of magnitude. Finally, after 1 h at 1500 $ \xb0$C, the measured concentration of Li was reduced to below $ 10 14 cm \u2212 3$, which is consistent with the trap concentration being somewhere below the SIMS detection limit for deuterium. This makes $ Li Zn$ an excellent candidate for the H trap in the as-received sample. Binding energies obtained from hybrid functional calculations lend further support to this. As listed in Table I, the $ Li Zn$H complex shows a binding energy of 1.17 eV.^{80}

Defect reaction . | Binding energy (eV) . |
---|---|

$ H O +\u2192 V O 0+ H i +$ | 1.02 |

$ ( V Zn H ) \u2212\u2192 V Zn 2 \u2212+ H i +$ | 3.04 |

$ ( V Zn 2H ) 0\u2192 ( V Zn H ) \u2212+ H i +$ | 2.17 |

$ ( V Zn 3H ) +\u2192 ( V Zn 2 H ) 0+ H i +$ | 1.38 |

$ ( Li Zn H ) 0\u2192 Li Zn \u2212+ H i +$ | 1.17 |

$ ( Fe Ga2 H ) 0\u2192 Fe Ga2 \u2212+ H i +$ | 0.74 |

$ H O1 +\u2192 V O1 0+ H i +$ | 1.13 |

$ H O2 \u2212\u2192 V O2 0+ H i ++2 e \u2212$ | 0.82 |

$ ( V Ga ib H ) 2 \u2212\u2192 ( V Ga ib ) 3 \u2212+ H i +$ | 2.94 |

$ ( V Ga ib 2 H ) \u2212\u2192 ( V Ga ib H ) 2 \u2212+ H i +$ | 2.62 |

$ ( V Ga ib 3 H ) 0\u2192 ( V Ga ib 2 H ) \u2212+ H i +$ | 0.73 |

Defect reaction . | Binding energy (eV) . |
---|---|

$ H O +\u2192 V O 0+ H i +$ | 1.02 |

$ ( V Zn H ) \u2212\u2192 V Zn 2 \u2212+ H i +$ | 3.04 |

$ ( V Zn 2H ) 0\u2192 ( V Zn H ) \u2212+ H i +$ | 2.17 |

$ ( V Zn 3H ) +\u2192 ( V Zn 2 H ) 0+ H i +$ | 1.38 |

$ ( Li Zn H ) 0\u2192 Li Zn \u2212+ H i +$ | 1.17 |

$ ( Fe Ga2 H ) 0\u2192 Fe Ga2 \u2212+ H i +$ | 0.74 |

$ H O1 +\u2192 V O1 0+ H i +$ | 1.13 |

$ H O2 \u2212\u2192 V O2 0+ H i ++2 e \u2212$ | 0.82 |

$ ( V Ga ib H ) 2 \u2212\u2192 ( V Ga ib ) 3 \u2212+ H i +$ | 2.94 |

$ ( V Ga ib 2 H ) \u2212\u2192 ( V Ga ib H ) 2 \u2212+ H i +$ | 2.62 |

$ ( V Ga ib 3 H ) 0\u2192 ( V Ga ib 2 H ) \u2212+ H i +$ | 0.73 |

Considering that the Li concentration level is below that of traps in the sample heat treated at 1100 $ \xb0$C and the observed reduction in dissociation rate, other traps are likely present with only slightly different binding energies to $ Li Zn$H. Based on hybrid functional calculations, two candidates of intrinsic nature are proposed. As mentioned, $ V Zn$ is a favorable trap for H, and experiments have shown that up to three H atoms can reside in the vacancy, resulting in a $ V Zn$3H complex acting as a shallow donor.^{68,72} As listed in Table I, the calculated binding energy decreases as H atoms are successively added to the vacancy to a value of 1.38 eV for $ V Zn$3H. In this scenario, the complex $ V Zn$2H would act as the trap for a third H. A prerequisite, for the existence of such trap would be $ V Zn$. However, $ V Zn$ is likely to be present, considering the low formation energy predicted under *n*-type conditions.^{9,72,74,81,82} Moreover, $ V Zn$2H is known to exist in ZnO.^{65,77} One of the key assumptions in the TLD model is that the trap is immobile. Although it would require further studies, one can speculate that since the first two H residing in $ V Zn$ results in a highly stable defect complex, they may effectively pin down $ V Zn$ (which otherwise has a calculated migration barrier of 1.4 eV^{81}), making it immobile during the 400 $ \xb0$C anneal. Another candidate for the trap is $ V O$, which would result in the aforementioned $ H O$, with calculated binding and dissociation energies of 1.02 and 1.74 eV, respectively. However, $ V O$ exhibits a rather high formation energy in *n*-type material and is, therefore, less likely to be present at the observed trap concentration level.^{9,74,81,82} Due to this, it is more likely that the defect complexes acting as traps for hydrogen diffusion are, $ V Zn$Li ( $ Li Zn$) and $ V Zn$2H, respectively. Here, it should be pointed out that, in order to suggest $ V Zn$2H, specifically, as a possible trap for H diffusion, DFT has played a substantial role.

$\beta -Ga 2 O 3$ has received massive attention in the last few years. As in the earlier example of H diffusion in ZnO, one expects H to exhibit trap-limited diffusion in $\beta -Ga 2 O 3$. In Fig. 7(b), the concentration vs depth profiles are shown for a ( $ 2 \xaf01$)-oriented $\beta -Ga 2 O 3$ sample subjected to $ 2$H implantation, followed by heat treatments between 350 and 600 $ \xb0$C. Heat treatments above 500 $ \xb0$C result in $ 2$H diffusion, which is evidenced by the increasing penetration of $ 2$H into the sample. Diffusion toward the surface of the sample is also observed for heat treatments in this temperature range. Importantly, the $ 2$H concentration-depth profiles are similar in shape to the concentration-depth profiles observed for $ 2$H diffusion in ZnO and H diffusion in SiC^{83} and can be modeled using the TLD model. In the model, the total concentration of the trap is assumed to remain unchanged during the diffusion process and found by fitting the height of the diffusion shoulder. The effective capture radius of the trap was set to 5 Å, which is on the order of the lattice constant. The fits from the TLD model are in excellent agreement with the data and capture especially well the steep decrease in the $ 2$H concentration deeper into the sample. This results in an activation energy for $ 2$H diffusion of 1.9 eV and a trap dissociation energy of 2.6 eV.

Interestingly, the 1.9 eV migration barrier is higher than the one obtained for H in ZnO^{75} and substantially higher than the 0.24 eV barrier predicted for H $ i$ migration along $[010]$ within the large eight-sided channel of $\beta -Ga 2 O 3$, as shown in Fig. 7(b). The $[001]$ direction is, however, perpendicular to the relevant direction in ( $ 2 \xaf01$) $\beta -Ga 2 O 3$, and a higher migration barrier can be envisioned for H $ i$ to hop between different channels than to move within a single channel. The 2.6 eV dissociation energy is also higher than the one obtained for ZnO. Again, we can propose potential traps based on calculated binding energies of H-related defects in $\beta -Ga 2 O 3$, as listed in Table I. In contrast to what was found for ZnO, impurities are not likely candidates in this case. Indeed, complexes between H and common acceptor impurities such as $ Fe Ga$ exhibit modest binding energies of 0.74 eV. The binding energy of $ H O$ complexes is similarly too low, irrespective of the prediction that $ H O 2$ in $\beta -Ga 2 O 3$ can be further stabilized under *n*-type conditions by undergoing a large structural rearrangement and capturing two electrons.^{84}

Similar to $ V Zn$ in *n*-type ZnO, $ V Ga$ exhibits the lowest formation energy among the intrinsic defects in $\beta -Ga 2 O 3$, and is a highly favorable trap for H. Interestingly, $ V Ga$ shows a rather unusual behavior where the simple vacancies on the Ga1 and Ga2 sites ( $ V Ga 1$ and $ V Ga 2$) are only metastable.^{85} Indeed, a neighboring Ga ion will rather jump into an interstitial position, such that the vacancy is split between two Ga sites connected by the $ Ga i$. As described in Refs. 85–88, there are three favorable split configurations, labeled as $ V Ga ia$, $ V Ga ib$, and $ V Ga ic$, where the latter is the lowest in energy under *n*-type conditions. The most favorable configuration of $ V Ga$ containing H is the $ V Ga ib$2H complex, as shown in Fig. 7(c). This very configuration has also been assigned to the dominant O–H-related vibrational line observed in $\beta -Ga 2 O 3$ annealed in H $ 2$.^{69} Curiously, $ V Ga ib$ binds the first two H with similar binding energies of 2.94 and 2.62 eV, respectively, which can be understood by considering the Coulomb repulsion between different $ H +$ ions in the complex. For $ V Zn$2H and $ V Zn$3H in ZnO, the O–H $ +$ bonds rotate to avoid each other but are still quite close together. In $ V Ga ib$2H, however, the two $ H +$ ions are situated in different vacancies located on opposite sides of the complex, resulting in minimal Coulomb repulsion. The calculated 2.62 eV binding energy (2.86 eV dissociation energy) of this complex is in excellent agreement with the 2.6 eV dissociation energy extracted using the TLD model. Furthermore, a third H would be captured with a much lower binding energy of 0.73 eV, as it must share the vacancy with another H and form an O–H bond with a threefold—rather than twofold—coordinated O. Other configurations with three H exist with higher binding energies. However, as discussed in detail by Fowler *et al.*,^{89} the rearrangements required to form these configurations after sequential trapping of H are improbable.

These studies on trap-limited diffusion of H in ZnO and $\beta -Ga 2 O 3$ clearly demonstrate the benefits of incorporating first-principles calculations. Calculated migration barriers and binding energies point to specific mechanisms for the diffusion process, while revealing additional insights into the defects involved, which in turn could be used to guide further experimental work.

#### 2. Zinc diffusion in β − *Ga*_{2}O_{3}

Zn has been suggested as a possible acceptor dopant to achieve semi-insulating layers in $\beta -Ga 2 O 3$-based devices.^{90} First-principles calculations predicted that Zn can form a stable acceptor on the substitutional Ga site,^{91,92} see Fig. 8(b). To study Zn diffusion in $\beta -Ga 2 O 3$, Zn was introduced into $\beta -Ga 2 O 3$ through a vapor phase in sealed evacuated quartz ampules at temperatures in the range 900–1100 $ \xb0$C.^{93} The resulting SIMS profiles are presented in Fig. 8(c) where the Zn concentrations turned out to be high, around $\u223c1\xd7 10 21 cm \u2212 1$. Furthermore, the samples were still conductive after introducing a considerable amount of what were presumably acceptors. Due to the familiar shape of the diffusion profiles, the TLD model was implemented to simulate the diffusion profiles. The profiles could indeed be fitted using the trap-limited model [see solid lines in Fig. 8(c)].

SIMS measurements can only provide the total Zn concentration and can not distinguish between different defect configurations. Through the combination of simulations and DFT, the dominating limiting reaction was found to be the trapping of $ Zn i$ by $ Zn Ga$. Without the DFT computed dissociation energies, the identity of the trap might not have been revealed.

To gain further insight into the defect dynamics, the luminescence from Zn-related defects in $\beta -Ga 2 O 3$ was studied with photoluminescence (PL) and first-principles calculations.^{94}, $\beta -Ga 2 O 3$ samples were heat treated in Zn vapor at 1100 $ \xb0$C and later heat-treated in $ O 2$ flow. After the heat treatments in $ O 2$ flow a broad emission peaking at around 2.5 eV emerged. Based on the agreement with the luminescence line calculated based on a one-dimensional configuration coordinate diagram derived from hybrid functional calculations, the 2.5 eV emission was attributed to $ Zn Ga$. It was also predicted that $ Zn Ga Zn i$, as a shallow donor, will not give rise to broad luminescence. Importantly, the luminescence from $ Zn Ga$ increased as the samples were annealed at higher temperatures in $ O 2$ flow. For the highest temperatures, the samples turned insulating. The results from this study supported the diffusion mechanism that was suggested in Ref. 93; after the initial heat treatment, $ Zn Ga Zn i$ dominates the Zn defect configuration, while after the subsequent heat treatment in $ O 2$ flow, the complex dissociates, leaving behind $ Zn Ga$. Zn diffusion in $\beta -Ga 2 O 3$ could have been studied with the experimental techniques only. However, the diffusion mechanism would not have been revealed without first-principles calculations.

### D. Reaction diffusion

^{14}:

^{14}

Figure 4(c) shows examples of reaction-diffusion concentration profiles. A key difference between the shape of profiles described by the reaction-diffusion model and the trap-limited diffusion model is that in the former model, the angle of the diffusion front is not dependent on the temperature as is the case in TLD. The shape will also depend on the charge state $q$ of the mediating defect. e.g., as $q$ increases, it results in a sharper diffusion front due to the much stronger dependency on $n$ in Eq. 18. The model will also collapse into free diffusion when the binding energy for the complex is high or if the background carrier concentration is higher than the change in charge carriers provided by the diffusing complex. In the following, the RD model will be utilized to describe vacancy-mediated diffusion, but it can in principle be applied to cases where the diffusion is mediated by other defects as well.

#### 1. Tin diffusion in *β* − *Ga*_{2}O_{3}

Sn prefers to substitute for Ga ( $ Sn Ga$) in $\beta -Ga 2 O 3$, acting as a shallow single donor,^{61} and is one of the most commonly used *n*-type dopants in this material.^{95} Under *n*-type conditions, the migration of $ Sn Ga$ is mediated by $ V Ga$.^{16} Indeed, interstitial Sn and Ga are both highly unfavorable in *n*-type material, as shown in Fig. 9(b), which means that interstitial-substitutional mechanisms are unlikely to play a role. In contrast, a low formation energy is predicted for $ V Ga$, and there will be a strong attractive interaction between the triply negatively charged vacancy and the Sn donor, making $ V Ga Sn Ga$ complex formation likely. Thus, in contrast to the H and Zn impurities discussed above, which migrate interstitially until trapped by a defect, $ Sn Ga$ is in the trapped state as an isolated defect, and migrates only through the formation of a complex with $ V Ga$.^{16}

As already mentioned, $ V Ga$ prefers to form split configurations in $\beta -Ga 2 O 3$, where $ V Ga ic$ is lowest in energy under *n*-type conditions. This holds true when a complex is formed between $ V Ga$ and Sn ( $ V Ga ic$Sn). As shown in Fig. 9(a), the Sn ion takes up the $ic$ interstitial position in the complex. NEB calculations have shown that split configurations play a crucial role in $ V Ga$ migration.^{86,88} Hybrid functional calculations yield migration barriers that depend strongly on the crystal direction: 2.1 eV for the [100] and [001] directions, and 1.0 eV for the [001] direction. For the latter direction, migration occurs by alternation between $ V Ga ic$ and $ V Ga ib$ configurations. The low barrier is enabled by passing through the three-split vacancy configuration $ V Ga ibc$ shown in Fig. 9(a), rather than the unfavorable $ V Ga 1$, as an intermediate state.^{88} NEB calculations have also revealed the pathways for migration of $ V Ga ic$Sn, which proceeds through a series of exchanges (Sn ion jumps into the $ V Ga$) and rotations (an adjacent Ga ion jumps into the $ V Ga$).^{16} The resulting overall migration barrier is 3.4, 3.2, and 3.4 eV for the [100], [010], and [001] directions, respectively.^{16}

In Fig. 9(c), diffusion of Sn from a Sn-doped bulk substrate with surface orientation (001) into an epilayer of $\beta -Ga 2 O 3$ was investigated. The sample was isochronally annealed under $ O 2$ flow for 30 min in the temperature range from 1050 to 1250 $ \xb0$C. Similar to the trap-limited diffusion profiles of H and Zn discussed above, the Sn diffusion profiles are characterized by a steep drop at the diffusion front. Considering a vacancy-mediated diffusion process, the Sn profiles were modeled using the RD model. As shown in Fig. 9(c), the experimental profiles are well described by the simulations. Notably, since the RD model contains several fitting parameters, these simulations relied heavily on input from the hybrid functional calculations. Specifically, the $ V Ga$ migration barrier and $ V Ga Sn Ga$ binding energy was fixed to the calculated values of 1.0 and 1.6 eV, respectively. The $ V Ga$ formation energy was treated as a fitting parameter (resulting in values of 12.4–12.8 eV), but guided by the value obtained from hybrid functional calculations with a tabulated value of $ \mu O$ for equilibrium with 1 atm $ O 2$ gas at 1150 $ \xb0$C (resulting in 14.0 eV). A background donor concentration ( $ N d +$) of $2\xd7 10 17 cm \u2212 3$ was also included in the charge-neutrality condition (based on the measured concentration of impurities in the epilayer of the as-received sample). Finally, the bandgap shrinkage with increasing temperature was extrapolated from previously reported experimental measurements. From an Arrhenius plot of the $ V Ga Sn Ga$ diffusivities extracted from the RD modeling, a diffusion coefficient and activation energy of $ D 0=2\xd7 10 \u2212 1 cm 2 s \u2212 1$ and 3.0 $\xb1$0.4 eV was obtained, respectively. The activation energy, which represents the $ V Ga$Sn $ Ga$ migration barrier, agrees well with the 3.4 eV obtained from NEB calculations for the [001] direction, and the diffusivity is close to the value expected from a random walk model (see Sec. III E).

The steep drop in concentration is well-captured by the simulations in Fig. 9(c). In the employed RD model, this arises due to the drop in Fermi-level position when going from the Sn-doped substrate to the epilayer, which results in a threefold increase in the formation energy of $ V Ga$ due to its triple negative charge state. The resulting lower concentration of $ V Ga Sn Ga$ complexes slows down the Sn migration at the diffusion front. Thus, in contrast to the TLD model, where the trap concentration is fixed and the steep drop arises due to a high binding energy of the complex, the shape of the Sn diffusion profile is, within the RD model, governed by the Fermi-level dependent access to diffusion-mediating $ V Ga$.

#### 2. Aluminum diffusion in doped ZnO

The reaction-diffusion model used to simulate Sn diffusion in $\beta -Ga 2 O 3$ can also be used to describe diffusion of Al in ZnO,^{14}^{,} *mutatis mutandis*. Al substitutes for Zn ( $ Al Zn$) and behaves as a shallow single donor in ZnO.^{74} In this case, the diffusion is mediated by $ V Zn$ in the double rather than triple negative charge state for $ V Ga$ as described above, resulting in a weaker Fermi-level dependence for the $ V Zn$ concentration. In Ref. 96, this relation was shown, where a nearly quadratic dependence was indeed established between the concentrations of $ V Zn$ and the charge carrier concentration given by the shallow Al-donor as measured using positron annihilation spectroscopy (PAS) and SIMS, respectively. This relation between the *n*-type doping and concentration of diffusion-mediating $ V Zn$ provides a way to experimentally verify the diffusion mechanism for Al in ZnO, and verify the assumption that the defect reaction was not limited by transport of $ V Zn$. If the vacancy mechanism is correct, it should be possible to control the profile by adjusting $ N d +$ through doping. Conversely, the diffusivity of Al within a free interstitial or substitutional-interstitial mechanism would not be enhanced with increasing Fermi level.^{96}

In Fig. 10, the Fermi-level dependence of Al diffusion was studied by comparing Al diffusion from a sputtered thin film containing Al to the $ 10 21 cm \u2212 3$ level into two ( $0001$)-oriented HT ZnO samples having different $ N d +$: (i) ZnO with a high residual Li content of $1\xd7 10 17 cm \u2212 3$, resulting in a low $ N d +$ of $\u223c3\xd7 10 13 cm \u2212 3$,^{14} and (ii) ZnO with a lower Li content ( $<1\xd7$10 $ 15$ cm $ \u2212 3$) doped with Ga donors to a concentration $ N d +\u22481\xd7 10 19 cm \u2212 3$. As shown in Fig. 10, the undoped ZnO shows an Al profile with a sharp drop at the diffusion front (even steeper than the ones obtained for Sn diffusion in $\beta -Ga 2 O 3$). In stark contrast, the Al profile in the predoped sample is flattened significantly and the Al diffusivity is enhanced. Hence, the RD model has almost collapsed into free diffusion due to the abundance of $ V Zn$. This result strongly suggests that Al diffusion is indeed mediated by $ V Zn$.

The inset Arrhenius plot in Fig. 10 shows basically identical activation energies for $ V Zn Al Zn$ diffusion of 2.4 and 2.5 eV for the doped and undoped ZnO, respectively. This is sensible, as the migration barrier of $ V Zn Al Zn$ is a material parameter that should not depend on the Fermi-level. Note that, for the four $ V Zn Al Zn$ diffusivities extracted from the doped sample, the Al concentration is below the background concentration of Ga (isoconcentration regime), i.e., the doping is dictated by Ga. This is not the case for the $1100 \xb0$C annealing step shown in Fig. 10, where the Al concentration profile partially exceeds $1\xd7 10 19 cm \u2212 3$. The extracted $ V Zn Al Zn$ migration barrier of about 2.5 eV is in close agreement with NEB calculations reported by Steiauf *et al.*,^{19} showing a 2.55 eV barrier for $ V Zn Al Zn$ migration along $[0001]$.

## VI. CONCLUDING REMARKS

Defects play a crucial role in the migration of atoms in solids. By studying the atomic interaction between the defects, diffusion processes can be well understood. In this Perspective, we have discussed approaches for combining experimental measurements, diffusion modeling, and first-principles calculations to investigate the interplay between impurities, matrix atoms, and intrinsic defects in wide bandgap semiconductors. The main mechanisms for diffusion in semiconductors are introduced, as well as how to predict relevant properties from theory. We have a particular focus on how properties relevant to diffusion can be obtained by DFT calculations. Important for the calculations are finding the energy barriers and, in particular the minimum energy paths, for migration. For more complex structures, like the monoclinic $\beta -Ga 2 O 3$, this can be a cumbersome exercise. Here, more automated search with the aid of machine learning approaches may be a way forward.

Several case studies are presented herein, and organized according to the complexity of the diffusion model and the integration between experiments and DFT calculations. Examples include models for free diffusion, including the study of self-diffusion and vacancy diffusion in 4H-SiC, where the experimentally obtained activation energies can be compared to similar values extracted from DFT. Indeed, DFT has developed into a framework for quantitative comparison with experiments, where excellent agreement is found for the examples shown.

Next, the trap-limited diffusion model is presented; a model that is successfully applied to wide bandgap oxides like ZnO and $\beta -Ga 2 O 3$, where the solubility of the diffusing species is typically low and the diffusion profile is governed by trapping and dissociation from a complex. Here, DFT is not only used to compare activation energies and confirm pathways, but also to identify the trap and understand the trapping mechanism. Hydrogen diffusion in various oxides, here represented by ZnO and $\beta -Ga 2 O 3$, are prime examples where such a diffusion model accurately describes the process.

Finally, we show two examples from an approach to diffusion modeling where results from DFT are essential in the model and part of the input parameters used, namely, the reaction-diffusion model. In the examples of Sn diffusion in $\beta -Ga 2 O 3$ and Al diffusion in ZnO, the formation of cation vacancies is crucial and given by the local Fermi level. Hence, formation energies can be used as input to the diffusion model and, taking the local charge carrier concentration into account, one can excellently describe otherwise challenging diffusion profiles. This model is particularly useful in wide bandgap oxides where cation vacancies tend to be double or even triple acceptors, and their formation energy is consequently highly Fermi-level dependent.

All case studies have given examples of how the approach of combining first-principle calculations and diffusion experiments has resulted in detailed descriptions of the atomic migration of defects. As more defect diffusion systems are studied in the same material, the extracted values improve, and a more complete picture forms. This approach has shown to be powerful in identifying diffusion mechanisms and can be applied to a wide range of defects in wide bandgap semiconductors.

## ACKNOWLEDGMENTS

Financial support was kindly provided by the Research Council of Norway and the University of Oslo through the frontier research projects GO-POW and QuTe (Grant Nos. 314017 and 325573, FriPro ToppForsk-program) and the Norwegian Micro- and Nano-Fabrication Facility, NorFab, Project No. 295864.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Ylva Knausgård Hommedal:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead). **Marianne Etzelmüller Bathen:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (supporting); Writing – original draft (equal). **Vilde Mari Reinertsen:** Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (supporting). **Klaus Magnus Johansen:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal). **Lasse Vines:** Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – original draft (equal). **Ymir Kalmann Frodason:** Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*2019 IEEE 7th Workshop on Wide Bandgap Power Devices and Applications (WiPDA)*(IEEE, 2019), pp. 200–207.

*Classical and Quantum Dynamics in Condensed Phase Simulations*(World Scientific, 1998), pp. 385–404.

*Atom Movements: Diffusion and Mass Transport in Solids*, Monographies de Physique (Editions de Physique, 1991).

*Hydrogen in Semiconductors*, edited by J. I. Pankove and N. M. Johnson (Academic Press, 1991).