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.

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) functional6 underestimates the experimental bandgaps of ZnO (3.4 eV7) and β -Ga 2 O 3 (4.9 eV8) 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 β -Ga 2 O 3.

A system with a concentration gradient will, given sufficient time, equalize the concentration in a process called diffusion. In 1855, Adolf Fick published his work on the mixing of salt and water.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
C t = D 2 C x 2 .
(1)
The diffusion coefficient has an Arrhenius relation to the activation energy E A,
D = D 0 exp ( E A / k T ) ,
(2)
where D 0 denotes a diffusion pre-factor, T is the temperature, and k is the Boltzmann constant.

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 sites19 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).

FIG. 1.

Ball-and-stick model illustrating defect-mediated diffusion mechanisms for substitutional and interstitial impurity atoms (red and pink) in a hypothetical diatomic crystal (host atoms are blue and orange): (1) direct interstitial, (2) interstitialcy, (3) vacancy, (4) kick-out, and (5) dissociative mechanism. For the first step of the interstitialcy diffusion (2) and the kick-out (4), the interstitial atom (blue) replaces the substitutional atom (pink) and the pink atom moves to the interstitial site. The pink interstitial atom can then either continue through a direct interstitial diffusion mechanism (1), resulting in a kick-out mechanism (4) or replace a substitutional blue atom and, thus, follow an interstitialcy mechanism (2).

FIG. 1.

Ball-and-stick model illustrating defect-mediated diffusion mechanisms for substitutional and interstitial impurity atoms (red and pink) in a hypothetical diatomic crystal (host atoms are blue and orange): (1) direct interstitial, (2) interstitialcy, (3) vacancy, (4) kick-out, and (5) dissociative mechanism. For the first step of the interstitialcy diffusion (2) and the kick-out (4), the interstitial atom (blue) replaces the substitutional atom (pink) and the pink atom moves to the interstitial site. The pink interstitial atom can then either continue through a direct interstitial diffusion mechanism (1), resulting in a kick-out mechanism (4) or replace a substitutional blue atom and, thus, follow an interstitialcy mechanism (2).

Close modal

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 β -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.

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.

FIG. 2.

Schematic (a) potential energy surface showing the barriers and energy cost of different steps in vacancy-mediated diffusion. This also illustrates the relation between the binding and dissociation energy of a complex. (b) Formation energy diagram for a vacancy ( V) with three stable charge states q = + 1 , 0 , 1 , 2, a shallow donor impurity ( X) and their complex ( V X). (c) Minimum-energy pathway between initial and final defect positions along the migration pathway, as obtained from a nudged elastic band calculation. The contour plot shows the initial guess (from linear interpolation) and the final path.

FIG. 2.

Schematic (a) potential energy surface showing the barriers and energy cost of different steps in vacancy-mediated diffusion. This also illustrates the relation between the binding and dissociation energy of a complex. (b) Formation energy diagram for a vacancy ( V) with three stable charge states q = + 1 , 0 , 1 , 2, a shallow donor impurity ( X) and their complex ( V X). (c) Minimum-energy pathway between initial and final defect positions along the migration pathway, as obtained from a nudged elastic band calculation. The contour plot shows the initial guess (from linear interpolation) and the final path.

Close modal

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).

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.

The defect is simulated in a periodically repeated supercell that is constructed from a repetition of bulk unit cells and, typically, contains a few hundred atoms depending on the computational demand of the studied system and employed exchange-correlation functional. Once geometry optimization has been performed for the defective supercell, the formation energy is calculated as27,28
E f ( X q ) = E tot ( X q ) E tot ( bulk ) i n i μ i + q ( E VBM + E F ) + Δ q ,
(3)
where E tot(bulk) is the total energy of the pristine supercell and E tot ( X q ) is the total energy of the supercell holding the defect in its equilibrium geometry. The defect is formed by exchanging a number of atomic species n i and charges q with reservoirs of chemical potential μ i and E F, respectively. The latter is the Fermi-level position, which is conventionally referenced to the valence band maximum (VBM) E VBM. The chemical potential of atomic species is referenced to the calculated total energy per atom of the element in its standard state, which represents the upper bound (denoted as i-rich conditions). Lower bounds on the chemical potential are imposed by the thermodynamic stability condition of the target material, e.g., Δ μ Zn + Δ μ O = Δ H f ( ZnO ) for ZnO, where Δ H f(ZnO) is the formation enthalpy of ZnO. The chemical potential of impurities is similarly bound by the formation of mixed secondary phases, e.g., μ H is limited under O-rich conditions due to the formation of H 2O.

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 V X. 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 Δ q is a finite-size correction that removes any spurious interaction between the charged defect and its images under periodic boundary conditions.29,30

Defect-assisted impurity diffusion involves the formation of an impurity-defect complex. Depending on the defects, the formation of this complex may be favorable or unfavorable due to elastic, magnetic, and especially electrostatic interactions between the defects. The binding energy E b of the complex A B is the difference in formation energy between the defect complex and the sum of the formation energies of the constituents A and B11,
E b ( A B ) = E f ( A ) + E f ( B ) E f ( A B ) .
(4)
A positive binding energy implies a stable complex. For charged defects, the binding energy can be Fermi-level dependent. However, in Fermi-level regions where the overall charge state is unchanged upon dissociation of the complex, the binding energy will be well defined.

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 V X dissociation energy ( E d ( V X )) is given by the V migration barrier ( E m ( V )) plus the V X binding energy ( E b ( V X )).

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.

Once the formation energies, binding energy, and migration barriers have been computed, the activation energy for vacancy-mediated diffusion can be calculated as
Q = [ E f ( V ) E b ( V X ) ] + E m ( V X ) ,
(5)
that is, the formation energy of a vacancy reduced by the binding energy of the complex plus the limiting migration barrier for the vacancy-mediated diffusion process.19 Note that this activation energy assumes that the X impurities are already incorporated substitutionally.

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 β -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

The diffusivity or diffusion coefficient D of a defect is often modeled according to the Arrhenius equation D = D 0 exp ( E A / k T ), 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 Γ 0 d 2. Here, Γ 0 is the attempt frequency, and d is the jump distance. If we approximate Γ 0 by a typical phonon frequency of 1 × 10 13 s 1 and assume a cubic system with a lattice distance of 3.25 Å, we obtain D 0 1 × 10 2 cm 2 s 1.19,43 For more accurate estimates of the diffusion pre-factor, the attempt frequency can be obtained using transition-state theory20 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.

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.

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.

FIG. 3.

Introduction of defects and diffusion by (a) ion implantation, (b) film deposition, and (c) through a vapor phase in a sealed ampule. (d) Schematic figure of secondary ion mass spectrometry (SIMS), where an ion source sputters the sample surface, releasing ions that are separated by electric and magnetic fields. (e) The result from a typical depth profiling measurement with SIMS. (f) Schematic figure of deep-level transient spectroscopy (DLTS). (g) A depth profiling of the concentration of a specific defect and a typical DLTS spectrum showing electrically active defects.

FIG. 3.

Introduction of defects and diffusion by (a) ion implantation, (b) film deposition, and (c) through a vapor phase in a sealed ampule. (d) Schematic figure of secondary ion mass spectrometry (SIMS), where an ion source sputters the sample surface, releasing ions that are separated by electric and magnetic fields. (e) The result from a typical depth profiling measurement with SIMS. (f) Schematic figure of deep-level transient spectroscopy (DLTS). (g) A depth profiling of the concentration of a specific defect and a typical DLTS spectrum showing electrically active defects.

Close modal

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.

Equation (1) can be solved with the proper boundary conditions. The solution for a constant surface source C S gives the concentration of a diffusing species A as
C A ( x , t ) = C S erfc ( x 2 D A T ) ,
(6)
where D A is the diffusion coefficient for the diffusing species. In this diffusion model, the shape of the resulting diffusion profile is a result of the initial defect gradient in the material and the diffusion barrier of the specific defect in a given solid but nothing else. Thus, this mode of diffusion is often referred to as free diffusion. A diffusion profile shape given by the error function is, therefore, a typical indication of free diffusion in solids, see Fig. 4(a).
FIG. 4.

Schematic of example diffusion modeling results for (a) free, (b) trap-limited, and (c) reaction diffusion.

FIG. 4.

Schematic of example diffusion modeling results for (a) free, (b) trap-limited, and (c) reaction diffusion.

Close modal

1. Self-diffusion in SiC

Studies on self-diffusion of species in crystalline solids are essential for understanding, e.g., defect formation during material growth and thermal response toward melting. An example of such a self-diffusion study in SiC, a promising material for future power electronic devices, was performed in Ref. 47. The diffusion of 13C and 30Si was studied using SIMS in a material stack consisting of a graphitized carbon cap, an isotope purified 4H-SiC epitaxial layer, and a 4H-SiC substrate with natural abundance of the different Si and C isotopes. Figure 5 shows SIMS profiles demonstrating the diffusion of 13C into the isotope purified epitaxial layer from (left side) the carbon cap and (right side) the substrate. Diffusion-profile fitting was performed using the error function according to
C = C 2 + C 1 C 2 2 erfc { x d 2 D T } ,
(7)
with C being the concentration, C 1 and C 2 boundary conditions, T the annealing temperature, and x and d representing depth and interface position, respectively. Arrhenius analysis resulted in an estimated diffusivity of D C = 8.3 × 10 6 × e 10.4 / k B T cm 2 s 1 for self-diffusion of 13C in 4H-SiC.
FIG. 5.

Self-diffusion of 13C in ( 0001) 4H-SiC in a stack consisting of a carbon cap, an isotope purified epitaxial layer, and a 4H-SiC substrate with naturally occurring isotope abundances. Self-diffusion of 13C is shown to occur in the isotope purified epitaxial layer from (left panel) the C-cap and (right panel) the substrate. The inset shows an Arrhenius plot of the diffusivities obtained from the diffusion modeling (solid curves) for 13C self-diffusion from the substrate. The data are taken from Ref. 47.

FIG. 5.

Self-diffusion of 13C in ( 0001) 4H-SiC in a stack consisting of a carbon cap, an isotope purified epitaxial layer, and a 4H-SiC substrate with naturally occurring isotope abundances. Self-diffusion of 13C is shown to occur in the isotope purified epitaxial layer from (left panel) the C-cap and (right panel) the substrate. The inset shows an Arrhenius plot of the diffusivities obtained from the diffusion modeling (solid curves) for 13C self-diffusion from the substrate. The data are taken from Ref. 47.

Close modal

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 trap49 and an important minority carrier lifetime killer50,51 in 4H-SiC. In fact, V C is present in SoA 4H-SiC epitaxial layers to densities around 5 × 10 12 cm 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 generated53 and was shown to migrate at temperatures of 1100–1600  °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 ¯ 0). The vacancies were generated by implantation of 4 MeV carbon ions (to avoid contaminating with impurities) with a projected depth of around 2.8  μ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 ) 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 ¯ 0 ] directions, respectively. The profiles were obtained using DLTS depth profiling at 285 K coinciding with the Z 1 / 2 peak maximum.

FIG. 6.

Diffusion of V C in 4H-SiC along two different crystallographic directions, (a) the c-axis [ 0001 ] and (b) the a-axis [ 11 2 ¯ 0 ], studied using DLTS depth profiling and fitted using diffusing equation modeling. The data are taken from Ref. 56. (c) The atomic structure of 4H-SiC (Si and C shown in gray and orange, respectively) and minimum energy paths for V C migration along [ 0001 ] and [ 11 2 ¯ 0 ], respectively. Reproduced with permission from Bathen et al., Phys. Rev. B 100, 014103 (2019). Copyright 2019 American Physical Society.

FIG. 6.

Diffusion of V C in 4H-SiC along two different crystallographic directions, (a) the c-axis [ 0001 ] and (b) the a-axis [ 11 2 ¯ 0 ], studied using DLTS depth profiling and fitted using diffusing equation modeling. The data are taken from Ref. 56. (c) The atomic structure of 4H-SiC (Si and C shown in gray and orange, respectively) and minimum energy paths for V C migration along [ 0001 ] and [ 11 2 ¯ 0 ], respectively. Reproduced with permission from Bathen et al., Phys. Rev. B 100, 014103 (2019). Copyright 2019 American Physical Society.

Close modal
FIG. 7.

2H concentration-depth profiles (markers) for (a) three different HT ZnO samples implanted with 2H and subjected to a 400  °C anneal for 30 min: one as-received and two pre-annealed for 1 h at 1100 and 1500  °C to reduce the Li content (indicated on the y-axis). Panel (b) shows 2H concentration–depth profiles for a ( 2 ¯ 01)-oriented β -Ga 2 O 3 sample implanted with 2H after 30 min duration heat treatments at selected temperatures between 350 and 600 °C. The solid curves are fits to the data using the trap-limited diffusion model. The data and fits in (b) are taken from Ref. 70. The insets show the lowest barrier for long-range H i migration, i.e., along [ 010 ] in β -Ga 2 O 3 and in the basal plane of ZnO, as obtained from NEB calculations using the strongly constrained and appropriately normed (SCAN)71 functional. (c) Relaxed structures and binding energies from Refs. 70 and 72 of Li ZnH and V Zn3H in ZnO, and V Ga ib2H in β -Ga 2 O 3.

FIG. 7.

2H concentration-depth profiles (markers) for (a) three different HT ZnO samples implanted with 2H and subjected to a 400  °C anneal for 30 min: one as-received and two pre-annealed for 1 h at 1100 and 1500  °C to reduce the Li content (indicated on the y-axis). Panel (b) shows 2H concentration–depth profiles for a ( 2 ¯ 01)-oriented β -Ga 2 O 3 sample implanted with 2H after 30 min duration heat treatments at selected temperatures between 350 and 600 °C. The solid curves are fits to the data using the trap-limited diffusion model. The data and fits in (b) are taken from Ref. 70. The insets show the lowest barrier for long-range H i migration, i.e., along [ 010 ] in β -Ga 2 O 3 and in the basal plane of ZnO, as obtained from NEB calculations using the strongly constrained and appropriately normed (SCAN)71 functional. (c) Relaxed structures and binding energies from Refs. 70 and 72 of Li ZnH and V Zn3H in ZnO, and V Ga ib2H in β -Ga 2 O 3.

Close modal

As shown in Figs. 6(a) and 6(b), high-temperature annealing at 1200–1500  °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 1, E A c = 4.4 eV and D 0 a = 0.02 cm 2 s 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 and k h, 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.

FIG. 8.

(a) The Zn Ga Zn i + donor complex. (b) Formation energies for defects relevant for Zn diffusion in β -Ga 2 O 3 under intermediate conditions. (c) Concentration of Zn as a function of depth in ( 2 ¯ 01) oriented β -Ga 2 O 3.93 The profiles are simulated using a trap-limited diffusion model. The data and fits in (c) are taken from Ref. 93.

FIG. 8.

(a) The Zn Ga Zn i + donor complex. (b) Formation energies for defects relevant for Zn diffusion in β -Ga 2 O 3 under intermediate conditions. (c) Concentration of Zn as a function of depth in ( 2 ¯ 01) oriented β -Ga 2 O 3.93 The profiles are simulated using a trap-limited diffusion model. The data and fits in (c) are taken from Ref. 93.

Close modal
FIG. 9.

(a) Relaxed structure of the V Ga icSn complex (where the purple Sn occupies the i c position), and three-split vacancy V Ga ibc (the intersitial Ga are indicated in blue). (b) Formation energy of defects relevant for Sn diffusion in β -Ga 2 O 3 under intermediate conditions. (c) Sn concentration–depth profiles for isochronally heat treated (30 min from 1050 to 1250  °C), as measured by SIMS. RD simulations are given by the black lines. Sn diffuses from a bulk substrate (highlighted in gray) into an epitaxial layer. The data and fits in (c) are taken from Ref. 16.

FIG. 9.

(a) Relaxed structure of the V Ga icSn complex (where the purple Sn occupies the i c position), and three-split vacancy V Ga ibc (the intersitial Ga are indicated in blue). (b) Formation energy of defects relevant for Sn diffusion in β -Ga 2 O 3 under intermediate conditions. (c) Sn concentration–depth profiles for isochronally heat treated (30 min from 1050 to 1250  °C), as measured by SIMS. RD simulations are given by the black lines. Sn diffuses from a bulk substrate (highlighted in gray) into an epitaxial layer. The data and fits in (c) are taken from Ref. 16.

Close modal
Trap-limited diffusion (TLD) involves trapping and dissociation of the diffusing species, where trapping slows down and dissociation enhances the transport. In the TLD model, the diffusing species X is assumed to be mobile and the trap A is immobile. Upon association between X and A, X goes from being mobile to being immobile and is considered trapped,
A immobile + X mobile X A immobile .
(8)
The model is based on Eq. (1) for diffusion of some species X,
C X t = D X 2 C X x 2 C X A t ,
(9)
where D X is the diffusion coefficient of X, C X is the concentration of mobile X, and C X A is the concentration of X associated with a trap A. The final term describes the association and dissociation of diffusing species X with the trap A to form the immobile complex X A,
C X A t = k X A C X C A k X C X A ,
(10)
where k X A is the capture rate, k X is the dissociation rate, and C A is the trap concentration. The capture rate has an Arrhenius relation to the migration barrier for the diffusing species E m ( X ),
k X A = 4 π R Γ 0 d 2 e E m ( X ) / k T .
(11)
Here, 4 π R describes the capture cross section, Γ 0 is the attempt frequency, and d is the jump distance. The dissociation rate has an Arrhenius relation to the dissociation energy E d ( X A ) needed for X to escape the trap A,
k X = Γ 0 e E d ( X A ) / k T .
(12)
Trap-limited diffusion often results in diffusion profiles with a sharp decrease in concentration, as illustrated in Fig. 4(b). In this model, the dissociation rate depends on temperature and the binding energy of the resulting complex. For a given capture rate, the angle of the diffusion front is strongly dependent on the dissociation rate. As the dissociation rate increases, the diffusion profile will approach the free diffusion profile as described above and lose the characteristic sharp drop in concentration. This will typically happen as the temperature is increased during the diffusion experiment and can be intuitively understood as reduced stability of the trap. Thus, at sufficiently high temperatures, the diffusion model collapses back into the free-diffusion model. For complexes with a very high binding energy, however, this effect is difficult to observe.

1. Hydrogen diffusion in ZnO and β-Ga2O3

Hydrogen is an ubiquitous impurity that strongly influences the properties of solids. This holds true for both ZnO and β -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 ( +/ ) charge-state transition level lies above the CBM in both ZnO and β -Ga 2 O 3, resulting in shallow-donor behavior.60,61 However, H i is generally quite mobile and reactive58—as exemplified by the low migration barriers calculated for H i in ZnO and β -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 β -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  °C for 1 h, and (iii) heat treated at 1500  °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  °C for 1 h the trap concentration is reduced below the detection limit for deuterium of 1 × 10 16 cm 3.

Interestingly, when evaluating the diffusion profiles in the sample heated at 1100  °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  °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 3 × 10 17 cm 3, which is close to the trap concentration. After 1 h at 1100  °C the Li content was lowered in the region of diffusion by approximately one order of magnitude. Finally, after 1 h at 1500  °C, the measured concentration of Li was reduced to below 10 14 cm 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 ZnH complex shows a binding energy of 1.17 eV.80 

TABLE I.

Binding energies for various hydrogenated defects in ZnO and β − Ga2O3 when the Fermi level is at the CBM.70,72 The corresponding dissociation energies can be estimated by adding the calculated Hi migration barriers of 0.72 and 0.24 eV for ZnO and β − Ga2O3, respectively.

Defect reactionBinding energy (eV)
H O + V O 0 + H i + 1.02 
( V Zn H ) V Zn 2 + H i + 3.04 
( V Zn 2H ) 0 ( V Zn H ) + H i + 2.17 
( V Zn 3H ) + ( V Zn 2 H ) 0 + H i + 1.38 
( Li Zn H ) 0 Li Zn + H i + 1.17 
( Fe Ga2 H ) 0 Fe Ga2 + H i + 0.74 
H O1 + V O1 0 + H i + 1.13 
H O2 V O2 0 + H i + + 2 e  0.82 
( V Ga ib H ) 2 ( V Ga ib ) 3 + H i + 2.94 
( V Ga ib 2 H ) ( V Ga ib H ) 2 + H i + 2.62 
( V Ga ib 3 H ) 0 ( V Ga ib 2 H ) + H i + 0.73 
Defect reactionBinding energy (eV)
H O + V O 0 + H i + 1.02 
( V Zn H ) V Zn 2 + H i + 3.04 
( V Zn 2H ) 0 ( V Zn H ) + H i + 2.17 
( V Zn 3H ) + ( V Zn 2 H ) 0 + H i + 1.38 
( Li Zn H ) 0 Li Zn + H i + 1.17 
( Fe Ga2 H ) 0 Fe Ga2 + H i + 0.74 
H O1 + V O1 0 + H i + 1.13 
H O2 V O2 0 + H i + + 2 e  0.82 
( V Ga ib H ) 2 ( V Ga ib ) 3 + H i + 2.94 
( V Ga ib 2 H ) ( V Ga ib H ) 2 + H i + 2.62 
( V Ga ib 3 H ) 0 ( V Ga ib 2 H ) + H i + 0.73 

Considering that the Li concentration level is below that of traps in the sample heat treated at 1100  °C and the observed reduction in dissociation rate, other traps are likely present with only slightly different binding energies to Li ZnH. 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 Zn3H 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 Zn3H. In this scenario, the complex V Zn2H 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 Zn2H 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 eV81), making it immobile during the 400  °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 ZnLi ( Li Zn) and V Zn2H, respectively. Here, it should be pointed out that, in order to suggest V Zn2H, specifically, as a possible trap for H diffusion, DFT has played a substantial role.

β -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 β -Ga 2 O 3. In Fig. 7(b), the concentration vs depth profiles are shown for a ( 2 ¯ 01)-oriented β -Ga 2 O 3 sample subjected to 2H implantation, followed by heat treatments between 350 and 600  °C. Heat treatments above 500  °C result in 2H diffusion, which is evidenced by the increasing penetration of 2H into the sample. Diffusion toward the surface of the sample is also observed for heat treatments in this temperature range. Importantly, the 2H concentration-depth profiles are similar in shape to the concentration-depth profiles observed for 2H diffusion in ZnO and H diffusion in SiC83 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 2H concentration deeper into the sample. This results in an activation energy for 2H 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 ZnO75 and substantially higher than the 0.24 eV barrier predicted for H i migration along [ 010 ] within the large eight-sided channel of β -Ga 2 O 3, as shown in Fig. 7(b). The [ 001 ] direction is, however, perpendicular to the relevant direction in ( 2 ¯ 01) β -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 β -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 β -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 β -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 ib2H complex, as shown in Fig. 7(c). This very configuration has also been assigned to the dominant O–H-related vibrational line observed in β -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 Zn2H and V Zn3H in ZnO, the O–H + bonds rotate to avoid each other but are still quite close together. In V Ga ib2H, 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 β -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 β − Ga2O3

Zn has been suggested as a possible acceptor dopant to achieve semi-insulating layers in β -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 β -Ga 2 O 3, Zn was introduced into β -Ga 2 O 3 through a vapor phase in sealed evacuated quartz ampules at temperatures in the range 900–1100  °C.93 The resulting SIMS profiles are presented in Fig. 8(c) where the Zn concentrations turned out to be high, around 1 × 10 21 cm 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)].

The TLD simulations are blind to the exact details of the diffusion process and can not reveal the identity of the trap. The first-principle calculations predict that Zn migrates via the interstitialcy mechanism through the β-Ga 2O 3 structure, but would get trapped in a Ga site if encountering a Ga vacancy
Zn i + V Ga Zn Ga ,
(13)
where Zn i is mobile and V Ga and Zn Ga are immobile. Once Zn i occupies a Ga site, E b is too high for Zn i to escape. V Ga is, therefore, an obvious candidate as the trap in the TLD model. However, the E d needed for Zn to escape Zn Ga is found with DFT to be 7.1 eV which is considerably higher than the E d extracted from the simulations of the Zn concentration profiles of 3.2 ± 0.6 eV.
DFT also predicted that two Zn atoms can be trapped in the same Ga site in a donor configuration Zn Ga Zn i
Zn i + Zn Ga Zn Ga Zn i ,
(14)
where Zn i is considered mobile and Zn Ga and Zn Ga Zn i immobile, see Fig. 8(a). Such a donor complex could explain why the samples were still conductive after introducing 1 × 10 20 cm 3 Zn impurities. Upon subsequent heat treatment in air, the samples turned insulating, indicating that Zn was dissociating from Zn Ga Zn i, leaving the acceptor Zn Ga behind. The E d for Zn i to escape Zn Ga Zn i was computed with DFT to be 3.0 eV which corresponded much better to the extracted value from the TLD model. Thus the E d from DFT revealed that Zn Ga was the probable identity of the trap in the implemented TLD model.

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 β -Ga 2 O 3 was studied with photoluminescence (PL) and first-principles calculations.94, β -Ga 2 O 3 samples were heat treated in Zn vapor at 1100  °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 β -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.

The reaction-diffusion (RD) model describes a situation where the species X is initially immobile and the diffusion is mediated through a reaction with a defect B through the formation of a mobile defect complex X B. When the complex dissociates the species X will go back to its immobile configuration:
B mobile + X immobile X B mobile .
(15)
The model is based on Eq. (1) for the diffusion of the X B complex:
C X B t = D X B 2 C X B x 2 C X t ,
(16)
where D X B is the diffusion coefficient for the X B complex, C X B is the concentration of the mobile X B complex, and C X is the concentration of immobile X. C X t describes the time-dependent defect reaction, given by the association and dissociation between X and B:
C X t = k X C X B k X B C B C X ,
(17)
where k X B is the capture rate, k X is the dissociation rate and C B is the concentration of B. It is worth noting that the equation is similar to the one used to describe trap-limited diffusion, Eq. (9), however, the roles of the species and defects are reversed and the assumptions, therefore, differ. e.g., in the following we will, in order to simplify the model, assume that the C B equals the thermodynamic equilibrium concentration of defects in the material. C B can then be calculated from the Fermi-level position. In reality, this means that the migration of the mediating defect is not limiting the defect reaction rate, which is satisfied if the transport of B is much higher than X B. Under this assumption C B is given as a function of the carrier concentration n14:
C B = N s ( 2 ( m e k T 2 π 2 ) 3 / 2 ) q exp ( E f k T ) n q ,
(18)
where N s is the number of atomic sites, m e is the electronic effective mass, q is the charge of the defect B and E f is the formation energy of B.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 β − Ga2O3

Sn prefers to substitute for Ga ( Sn Ga) in β -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 β -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 icSn). As shown in Fig. 9(a), the Sn ion takes up the i c 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 icSn, 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 β -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  °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 μ O for equilibrium with 1 atm O 2 gas at 1150  °C (resulting in 14.0 eV). A background donor concentration ( N d +) of 2 × 10 17 cm 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 × 10 1 cm 2 s 1 and 3.0 ±0.4 eV was obtained, respectively. The activation energy, which represents the V GaSn 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 β -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 3 level into two ( 0001)-oriented HT ZnO samples having different N d +: (i) ZnO with a high residual Li content of 1 × 10 17 cm 3, resulting in a low N d + of 3 × 10 13 cm 3,14 and (ii) ZnO with a lower Li content ( < 1 ×10 15 cm 3) doped with Ga donors to a concentration N d + 1 × 10 19 cm 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 β -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.

FIG. 10.

Experimental Al diffusion profiles in two ( 0001) HT ZnO samples with different background-donor concentration: doped ( 1 × 10 19 cm 3) and undoped ( 3 × 10 13 cm 3).96 Al is diffused into the bulk ZnO from an Al-doped thin film deposited by sputtering ( Al 2 × 10 21 cm 3). The samples were isochronally annealed for 30 min starting at 700  °C (doped) and 800  °C (undoped) in temperature steps of 100  °C, and both profiles shown are for the 1100  °C annealing step. The data and fits are taken from Ref. 96.

FIG. 10.

Experimental Al diffusion profiles in two ( 0001) HT ZnO samples with different background-donor concentration: doped ( 1 × 10 19 cm 3) and undoped ( 3 × 10 13 cm 3).96 Al is diffused into the bulk ZnO from an Al-doped thin film deposited by sputtering ( Al 2 × 10 21 cm 3). The samples were isochronally annealed for 30 min starting at 700  °C (doped) and 800  °C (undoped) in temperature steps of 100  °C, and both profiles shown are for the 1100  °C annealing step. The data and fits are taken from Ref. 96.

Close modal

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 °C annealing step shown in Fig. 10, where the Al concentration profile partially exceeds 1 × 10 19 cm 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 ].

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 β -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 β -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 β -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 β -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.

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.

The authors have no conflicts to disclose.

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).

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

1.
J.
Millan
,
P.
Godignon
,
X.
Perpina
,
A.
Perez-Tomas
, and
J.
Rebollo
, “
A survey of wide bandgap power semiconductor devices
,”
IEEE Trans. Power Electron.
29
,
2155
2163
(
2014
).
2.
G.
Iannaccone
,
C.
Sbrana
,
I.
Morelli
, and
S.
Strangio
, “
Power electronics based on wide-bandgap semiconductors: Opportunities and challenges
,”
IEEE Access
9
,
139446
139456
(
2021
).
3.
H.
Morkoç
,
S.
Strite
,
G.
Gao
,
M.
Lin
,
B.
Sverdlov
, and
M.
Burns
, “
Large-band-gap SiC, III-V nitride, and II-VI ZnSe-based semiconductor device technologies
,”
J. Appl. Phys.
76
,
1363
1398
(
1994
).
4.
F.
Nouketcha
,
A.
Lelis
,
R.
Green
,
Y.
Cui
,
C.
Darmody
, and
N.
Goldsman
, “Detailed study of breakdown voltage and critical field in wide bandgap semiconductors,” in 2019 IEEE 7th Workshop on Wide Bandgap Power Devices and Applications (WiPDA) (IEEE, 2019), pp. 200–207.
5.
J. P.
Perdew
and
M.
Levy
, “
Physical content of the exact Kohn-Sham orbital energies: Band gaps and derivative discontinuities
,”
Phys. Rev. Lett.
51
,
1884
1887
(
1983
).
6.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
, “
Rationale for mixing exact exchange with density functional approximations
,”
J. Chem. Phys.
105
,
9982
(
1996
).
7.
D. C.
Reynolds
,
D. C.
Look
,
B.
Jogai
,
C. W.
Litton
,
G.
Cantwell
, and
W. C.
Harsch
, “
Valence-band ordering in ZnO
,”
Phys. Rev. B
60
,
2340
2344
(
1999
).
8.
C.
Janowitz
,
V.
Scherer
,
M.
Mohamed
,
A.
Krapf
,
H.
Dwelk
,
R.
Manzke
,
Z.
Galazka
,
R.
Uecker
,
K.
Irmscher
,
R.
Fornari
,
M.
Michling
,
D.
Schmeißer
,
J. R.
Weber
,
J. B.
Varley
, and
C. G.
Van de Walle
, “
Experimental electronic structure of In 2 O 3 and Ga 2 O 3
,”
New J. Phys.
13
,
085014
(
2011
).
9.
F.
Oba
,
A.
Togo
,
I.
Tanaka
,
J.
Paier
, and
G.
Kresse
, “
Defect energetics in ZnO: A hybrid hartree-fock density functional study
,”
Phys. Rev. B
77
,
245202
(
2008
).
10.
T.
Zacherle
,
P. C.
Schmidt
, and
M.
Martin
, “
Ab initio calculations on the defect structure of β - Ga 2 O 3
,”
Phys. Rev. B
87
,
235206
(
2013
).
11.
C.
Freysoldt
,
B.
Grabowski
,
T.
Hickel
,
J.
Neugebauer
,
G.
Kresse
,
A.
Janotti
, and
C. G.
Van de Walle
, “
First-principles calculations for point defects in solids
,”
Rev. Mod. Phys.
86
,
253
305
(
2014
).
12.
S.
Lany
and
A.
Zunger
, “
Polaronic hole localization and multiple hole binding of acceptors in oxide wide-gap semiconductors
,”
Phys. Rev. B
80
,
085202
(
2009
).
13.
G.
Miceli
,
W.
Chen
,
I.
Reshetnyak
, and
A.
Pasquarello
, “
Nonempirical hybrid functionals for band gaps and polaronic distortions in solids
,”
Phys. Rev. B
97
,
121112
(
2018
).
14.
K. M.
Johansen
,
L.
Vines
,
T. S.
Björheim
,
R.
Schifano
, and
B. G.
Svensson
, “
Aluminum migration and intrinsic defect interaction in single-crystal zinc oxide
,”
Phys. Rev. Appl.
3
,
024003
(
2015
).
15.
W.
Chen
and
A.
Pasquarello
, “
Accuracy of GW for calculating defect energy levels in solids
,”
Phys. Rev. B
96
,
020101
(
2017
).
16.
Y. K.
Frodason
,
P. P.
Krzyzaniak
,
L.
Vines
,
J. B.
Varley
,
C. G.
Van de Walle
, and
K. M. H.
Johansen
, “
Diffusion of Sn donors in β - Ga 2 O 3
,”
APL Mater.
11
,
041121
(
2023
).
17.
A.
Fick
, “
V. On liquid diffusion
,”
Lond. Edinb. Dubl. Philos. Mag.
10
,
30
39
(
1855
).
18.
H. Mehrer
,
Diffusion in Solids: Fundamentals, Methods, Materials, Diffusion-Controlled Processes, Springer Series in Solid-State Sciences, Vol. 155 (Springer Berlin Heidelberg, Berlin, Heidelberg
,
2007
).
19.
D.
Steiauf
,
J. L.
Lyons
,
A.
Janotti
, and
C. G.
Van de Walle
, “
First-principles study of vacancy-assisted impurity diffusion in ZnO
,”
APL Mater.
2
,
096101
(
2014
).
20.
G. H.
Vineyard
, “
Frequency factors and isotope effects in solid state rate processes
,”
J. Phys. Chem. Solids
3
,
121
127
(
1957
).
21.
G.
Gilmer
,
T.
Diaz de la Rubia
,
D.
Stock
, and
M.
Jaraiz
, “
Diffusion and interactions of point defects in silicon: Molecular dynamics simulations
,”
Nucl. Instrum.
102
,
247
255
(
1995
).
22.
I.
Martin-Bragado
,
R.
Borges
,
J. P.
Balbuena
, and
M.
Jaraiz
, “
Kinetic Monte Carlo simulation for semiconductor processing: A review
,”
Prog. Mater. Sci.
92
,
1
32
(
2018
).
23.
R.
Kuate Defo
,
X.
Zhang
,
D.
Bracher
,
G.
Kim
,
E.
Hu
, and
E.
Kaxiras
, “
Energetics and kinetics of vacancy defects in 4H-SiC
,”
Phys. Rev. B
98
,
104103
(
2018
).
24.
E. M. Y.
Lee
,
A.
Yu
,
J. J.
de Pablo
, and
G.
Galli
, “
Stability and molecular pathways to the formation of spin defects in silicon carbide
,”
Nat. Commun.
12
,
6325
(
2021
).
25.
C. G.
Van de Walle
and
J.
Neugebauer
, “
First-principles calculations for defects and impurities: Applications to III-nitrides
,”
J. Appl. Phys.
95
,
3851
3879
(
2004
).
26.
C. E.
Dreyer
,
A.
Alkauskas
,
J. L.
Lyons
,
A.
Janotti
, and
C. G.
Van de Walle
, “
First-principles calculations of point defects for quantum technologies
,”
Annu. Rev. Mater. Res.
48
,
1
26
(
2018
).
27.
S.
Zhang
and
J.
Northrup
, “
Chemical potential dependence of defect formation energies in GaAs: Application to Ga self-diffusion
,”
Phys. Rev. Lett.
67
,
2339
2342
(
1991
).
28.
C. G.
Van de Walle
,
D. B.
Laks
,
G. F.
Neumark
, and
S. T.
Pantelides
, “
First-principles calculations of solubilities and doping limits: Li, Na, and N in ZnSe
,”
Phys. Rev. B
47
,
9425
9434
(
1993
).
29.
C.
Freysoldt
,
J.
Neugebauer
, and
C. G.
Van de Walle
, “
Fully ab initio finite-size corrections for charged-defect supercell calculations
,”
Phys. Rev. Lett.
102
,
016402
(
2009
).
30.
Y.
Kumagai
and
F.
Oba
, “
Electrostatics-based finite-size corrections for first-principles point defect calculations
,”
Phys. Rev. B
89
,
195205
(
2014
).
31.
H.
Jónsson
,
G.
Mills
, and
K. W.
Jacobsen
, “Nudged elastic band method for finding minimum energy paths of transitions,” in Classical and Quantum Dynamics in Condensed Phase Simulations (World Scientific, 1998), pp. 385–404.
32.
G.
Henkelman
and
H.
Jónsson
, “
A climbing image nudged elastic band method for finding saddle points and minimum energy paths
,”
J. Chem. Phys.
113
,
9901
9904
(
2000
).
33.
G.
Henkelman
and
H.
Jónsson
, “
A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives
,”
J. Chem. Phys.
111
,
7010
7022
(
1999
).
34.
S. K.
Estreicher
,
M.
Sanati
,
D.
West
, and
F.
Ruymgaart
, “
Thermodynamics of impurities in semiconductors
,”
Phys. Rev. B
70
,
125209
(
2004
).
35.
D.
Gomes
,
V. P.
Markevich
,
A. R.
Peaker
, and
J.
Coutinho
, “
Dynamics of hydrogen in silicon at finite temperatures from first principles
,”
Phys. Status Solidi B
259
,
2100670
(
2022
).
36.
S. F. J.
Cox
,
E. A.
Davis
,
S. P.
Cottrell
,
P. J. C.
King
,
J. S.
Lord
,
J. M.
Gil
,
H. V.
Alberto
,
R. C.
Vilão
,
J.
Piroto Duarte
,
N.
Ayres de Campos
,
A.
Weidinger
,
R. L.
Lichti
, and
S. J. C.
Irvine
, “
Experimental confirmation of the predicted shallow donor hydrogen state in zinc oxide
,”
Phys. Rev. Lett.
86
,
2601
2604
(
2001
).
37.
R. L.
Lichti
,
Y. G.
Celebi
,
S. P.
Cottrell
,
S. F. J.
Cox
, and
E. A.
Davis
, “
Diffusion and trapping of Mu in the III–V nitrides
,”
J. Phys. Condens. Matter.
16
,
S4721
S4738
(
2004
).
38.
R.
Hauschild
,
H.
Priller
,
M.
Decker
,
J.
Brückner
,
H.
Kalt
, and
C.
Klingshirn
, “
Temperature dependent band gap and homogeneous line broadening of the exciton emission in ZnO
,”
Phys. Status Solidi C
3
,
976
979
(
2006
).
39.
A.
Mock
,
J.
VanDerslice
,
R.
Korlacki
,
J. A.
Woollam
, and
M.
Schubert
, “
Elevated temperature dependence of the anisotropic visible-to-ultraviolet dielectric function of monoclinic β - Ga 2 O 3
,”
Appl. Phys. Lett.
112
,
041905
(
2018
).
40.
T. N.
Sky
,
K. M.
Johansen
,
Y. K.
Frodason
,
T.
Aarholt
,
H. N.
Riise
,
Ø.
Prytz
,
B. G.
Svensson
, and
L.
Vines
, “
Diffusion of indium in single crystal zinc oxide: A comparison between group III donors
,”
Semicond. Sci. Technol.
34
,
025011
(
2019
).
41.
B.
Monserrat
, “
Correlation effects on electron-phonon coupling in semiconductors: Many-body theory along thermal lines
,”
Phys. Rev. B
93
,
100301
(
2016
).
42.
C.
Lee
,
N. D.
Rock
,
A.
Islam
,
M. A.
Scarpulla
, and
E.
Ertekin
, “
Electron–phonon effects and temperature-dependence of the electronic structure of monoclinic β - Ga 2 O 3
,”
APL Mater.
11
,
011106
(
2023
).
43.
J.
Philibert
, Atom Movements: Diffusion and Mass Transport in Solids, Monographies de Physique (Editions de Physique, 1991).
44.
O. N.
Bedoya-Martínez
and
G.
Roma
, “
Activation entropies for diffusion in cubic silicon carbide from first principles
,”
Phys. Rev. B
82
,
134115
(
2010
).
45.
H.
Bracht
, “
Diffusion mechanisms and intrinsic point-defect properties in silicon
,”
MRS Bull.
25
,
22
27
(
2000
).
46.
R. B.
Fair
and
J. C. C.
Tsai
, “
A quantitative model for the diffusion of phosphorus in silicon and the emitter dip effect
,”
J. Electrochem. Soc.
124
,
1107
(
1977
).
47.
M. E.
Bathen
,
M.
Linnarsson
,
M.
Ghezellou
,
J.
Ul Hassan
, and
L.
Vines
, “
Influence of carbon cap on self-diffusion in silicon carbide
,”
Crystals
10
,
1
11
(
2020
).
48.
X.
Yan
,
P.
Li
,
L.
Kang
,
S.-H.
Wei
, and
B.
Huang
, “
First-principles study of electronic and diffusion properties of intrinsic defects in 4H-SiC
,”
J. Appl. Phys.
127
,
085702
(
2020
).
49.
N. T.
Son
,
X. T.
Trinh
,
L. S.
Løvlie
,
B. G.
Svensson
,
K.
Kawahara
,
J.
Suda
,
T.
Kimoto
,
T.
Umeda
,
J.
Isoya
,
T.
Makino
,
T.
Ohshima
, and
E.
Janzén
, “
Negative-U system of carbon vacancy in 4H-SiC
,”
Phys. Rev. Lett.
109
,
187603
(
2012
).
50.
P. B.
Klein
,
B. V.
Shanabrook
,
S. W.
Huh
,
A. Y.
Polyakov
,
M.
Skowronski
,
J. J.
Sumakeris
, and
M. J.
O’Loughlin
, “
Lifetime-limiting defects in n 4H-SiC epilayers
,”
Appl. Phys. Lett.
88
,
052110
(
2006
).
51.
K.
Danno
,
D.
Nakamura
, and
T.
Kimoto
, “
Investigation of carrier lifetime in 4H-SiC epilayers and lifetime control by electron irradiation
,”
Appl. Phys. Lett.
90
,
202109
(
2007
).
52.
B.
Zippelius
,
J.
Suda
, and
T.
Kimoto
, “
High temperature annealing of n-type 4H-SiC: Impact on intrinsic defects and carrier lifetime
,”
J. Appl. Phys.
111
,
033515
(
2012
).
53.
H. M.
Ayedh
,
V.
Bobal
,
R.
Nipoti
,
A.
Hallén
, and
B. G.
Svensson
, “
Formation of carbon vacancy in 4 H silicon carbide during high-temperature processing
,”
J. Appl. Phys.
115
,
012005
(
2014
).
54.
T.
Umeda
,
J.
Isoya
,
N.
Morishita
,
T.
Ohshima
,
T.
Kamiya
,
A.
Gali
,
P.
Deák
,
N. T.
Son
, and
E.
Janzén
, “
EPR and theoretical studies of positively charged carbon vacancy in 4H-SiC
,”
Phys. Rev. B
70
,
235212
(
2004
).
55.
H. M.
Ayedh
,
A.
Hallén
, and
B. G.
Svensson
, “
Elimination of carbon vacancies in 4H-SiC epi-layers by near-surface ion implantation: Influence of the ion species
,”
J. Appl. Phys.
118
,
175701
(
2015
).
56.
M. E.
Bathen
,
J.
Coutinho
,
H. M.
Ayedh
,
J.
Ul Hassan
,
I.
Farkas
,
S.
Öberg
,
Y. K.
Frodason
,
B. G.
Svensson
, and
L.
Vines
, “
Anisotropic and plane-selective migration of the carbon vacancy in SiC: Theory and experiment
,”
Phys. Rev. B
100
,
014103
(
2019
).
57.
C. G.
Hemmingsson
,
N. T.
Son
,
A.
Ellison
,
J.
Zhang
, and
E.
Janzén
, “
Negative-U centers in 4H silicon carbide
,”
Phys. Rev. B
58
,
R10119
R10122
(
1998
).
58.
Hydrogen in Semiconductors, edited by J. I. Pankove and N. M. Johnson (Academic Press, 1991).
59.
C. G.
Van de Walle
and
J.
Neugebauer
, “
Universal alignment of hydrogen levels in semiconductors, insulators and solutions
,”
Nature
423
,
626
628
(
2003
).
60.
C. G.
Van de Walle
, “
Hydrogen as a cause of doping in zinc oxide
,”
Phys. Rev. Lett.
85
,
1012
1015
(
2000
).
61.
J. B.
Varley
,
J. R.
Weber
,
A.
Janotti
, and
C. G.
Van de Walle
, “
Oxygen vacancies and donor impurities in β - Ga 2 O 3
,”
Appl. Phys. Lett.
97
,
142106
(
2010
).
62.
M. G.
Wardle
,
J. P.
Goss
, and
P. R.
Briddon
, “
First-principles study of the diffusion of hydrogen in ZnO
,”
Phys. Rev. Lett.
96
,
205504
(
2006
).
63.
A.
Janotti
and
C. G.
Van de Walle
, “
Hydrogen multicentre bonds
,”
Nat. Mater.
6
,
44
47
(
2006
).
64.
E. V.
Lavrov
,
F.
Herklotz
, and
J.
Weber
, “
Identification of two hydrogen donors in ZnO
,”
Phys. Rev. B
79
,
165210
(
2009
).
65.
E. V.
Lavrov
,
J.
Weber
,
F.
Börrnert
,
C. G.
Van de Walle
, and
R.
Helbig
, “
Hydrogen-related defects in ZnO studied by infrared absorption spectroscopy
,”
Phys. Rev. B
66
,
165205
(
2002
).
66.
E.
Lavrov
,
F.
Börrnert
, and
J.
Weber
, “
Hydrogen motion in ZnO
,”
Phys. B: Condens. Matter
401–402
,
366
369
(
2007
).
67.
D.
Bastin
,
E. V.
Lavrov
, and
J.
Weber
, “
Metastable state of the V Zn H 2 defect in ZnO
,”
Phys. Rev. B
83
,
195210
(
2011
).
68.
F.
Herklotz
,
A.
Hupfer
,
K. M.
Johansen
,
B. G.
Svensson
,
S. G.
Koch
, and
E. V.
Lavrov
, “
Infrared absorption on a complex comprising three equivalent hydrogen atoms in ZnO
,”
Phys. Rev. B
92
,
155203
(
2015
).
69.
P.
Weiser
,
M.
Stavola
,
W. B.
Fowler
,
Y.
Qin
, and
S.
Pearton
, “
Structure and vibrational properties of the dominant O–H center in β - Ga 2 O 3
,”
Appl. Phys. Lett.
112
,
232104
(
2018
).
70.
V. M.
Reinertsen
,
P. M.
Weiser
,
Y. K.
Frodason
,
M. E.
Bathen
,
L.
Vines
, and
K. M.
Johansen
, “
Anisotropic and trap-limited diffusion of hydrogen/deuterium in monoclinic gallium oxide single crystals
,”
Appl. Phys. Lett.
117
,
232106
(
2020
).
71.
J.
Sun
,
R. C.
Remsing
,
Y.
Zhang
,
Z.
Sun
,
A.
Ruzsinszky
,
H.
Peng
,
Z.
Yang
,
A.
Paul
,
U.
Waghmare
,
X.
Wu
,
M. L.
Klein
, and
J. P.
Perdew
, “
Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional
,”
Nat. Chem.
8
,
831
836
(
2016
).
72.
Y. K.
Frodason
,
K. M.
Johansen
,
T. S.
Bjørheim
,
B. G.
Svensson
, and
A.
Alkauskas
, “
Zn vacancy-donor impurity complexes in ZnO
,”
Phys. Rev. B
97
,
104109
(
2018
).
73.
D. G.
Thomas
and
J. J.
Lander
, “
Hydrogen as a donor in zinc oxide
,”
J. Chem. Phys.
25
,
1136
1142
(
1956
).
74.
A.
Janotti
and
C. G.
Van de Walle
, “
Fundamentals of zinc oxide as a semiconductor
,”
Rep. Prog. Phys.
72
,
126501
(
2009
).
75.
K. M.
Johansen
,
J. S.
Christensen
,
E. V.
Monakhov
,
A. Y.
Kuznetsov
, and
B. G.
Svensson
, “
Deuterium diffusion and trapping in hydrothermally grown single crystalline ZnO
,”
Appl. Phys. Lett.
93
,
152109
(
2008
).
76.
L. E.
Halliburton
,
L.
Wang
,
L.
Bai
,
N. Y.
Garces
,
N. C.
Giles
,
M. J.
Callahan
, and
B.
Wang
, “
Infrared absorption from OH ions adjacent to lithium acceptors in hydrothermally grown ZnO
,”
J. Appl. Phys.
96
,
7168
7172
(
2004
).
77.
E. V.
Lavrov
,
F.
Börrnert
, and
J.
Weber
, “
Dominant hydrogen-oxygen complex in hydrothermally grown ZnO
,”
Phys. Rev. B
71
,
035205
(
2005
).
78.
G. A.
Shi
,
M.
Stavola
, and
W. B.
Fowler
, “
Identification of an OH–Li center in ZnO: Infrared absorption spectroscopy and density functional theory
,”
Phys. Rev. B
73
,
081201
(
2006
).
79.
K. M.
Johansen
,
H.
Haug
,
E.
Lund
,
E. V.
Monakhov
, and
B. G.
Svensson
, “
Thermal stability of the OH-Li complex in hydrothermally grown single crystalline ZnO
,”
Appl. Phys. Lett.
97
,
211907
(
2010
).
80.
Y. K.
Frodason
,
K. M.
Johansen
,
A.
Galeckas
, and
L.
Vines
, “
Broad luminescence from donor-complexed Li Zn and Na Zn acceptors in ZnO
,”
Phys. Rev. B
100
,
184102
(
2019
).
81.
A.
Janotti
and
C. G.
Van de Walle
, “
Native point defects in ZnO
,”
Phys. Rev. B
76
,
165202
(
2007
).
82.
J. L.
Lyons
,
J. B.
Varley
,
D.
Steiauf
,
A.
Janotti
, and
C. G.
Van de Walle
, “
First-principles characterization of native-defect-related optical transitions in ZnO
,”
J. Appl. Phys.
122
,
035704
(
2017
).
83.
M. S.
Janson
,
A.
Hallén
,
M. K.
Linnarsson
, and
B. G.
Svensson
, “
Hydrogen diffusion, complex formation, and dissociation in acceptor-doped silicon carbide
,”
Phys. Rev. B
64
,
195202
(
2001
).
84.
A.
Langørgen
,
C.
Zimmermann
,
Y. K.
Frodason
,
E. F.
Verhoeven
,
P.
Weiser
,
R. M.
Karsthof
,
J. B.
Varley
, and
L.
Vines
, “
Influence of heat treatments in H 2 and Ar on the E 1 center in β - Ga 2 O 3
,”
J. Appl. Phys.
131
,
115702
(
2022
).
85.
J. B.
Varley
,
H.
Peelaers
,
A.
Janotti
, and
C. G.
Van de Walle
, “
Hydrogenated cation vacancies in semiconducting oxides
,”
J. Phys.: Condens. Matter
23
,
334212
(
2011
).
86.
A.
Kyrtsos
,
M.
Matsubara
, and
E.
Bellotti
, “
Migration mechanisms and diffusion barriers of vacancies in Ga 2 O 3
,”
Phys. Rev. B
95
,
245202
(
2017
).
87.
J. M.
Johnson
,
Z.
Chen
,
J. B.
Varley
,
C. M.
Jackson
,
E.
Farzana
,
Z.
Zhang
,
A. R.
Arehart
,
H.-L.
Huang
,
A.
Genc
,
S. A.
Ringel
,
C. G.
Van de Walle
,
D. A.
Muller
, and
J.
Hwang
, “
Unusual formation of point-defect complexes in the ultrawide-band-gap semiconductor β - Ga 2 O 3
,”
Phys. Rev. X
9
,
041027
(
2019
).
88.
Y. K.
Frodason
,
J. B.
Varley
,
K. M. H.
Johansen
,
L.
Vines
, and
C. G.
Van de Walle
, “
Migration of Ga vacancies and interstitials in β - Ga 2 O 3
,”
Phys. Rev. B
107
,
024109
(
2023
).
89.
W. B.
Fowler
,
M.
Stavola
,
Y.
Qin
, and
P.
Weiser
, “
Trapping of multiple H atoms at the Ga(1) vacancy in β - Ga 2 O 3
,”
Appl. Phys. Lett.
117
,
142101
(
2020
).
90.
E.
Chikoidze
,
T.
Tchelidze
,
C.
Sartel
,
Z.
Chi
,
R.
Kabouche
,
I.
Madaci
,
C.
Rubio
,
H.
Mohamed
,
V.
Sallet
,
F.
Medjdoub
,
A.
Perez-Tomas
, and
Y.
Dumont
, “
Ultra-high critical electric field of 13.2 MV/cm for Zn-doped p-type β - Ga 2 O 3
,”
Mater. Today Phys.
15
,
100263
(
2020
).
91.
H.
Peelaers
,
J. L.
Lyons
,
J. B.
Varley
, and
C. G.
Van de Walle
, “
Deep acceptors and their diffusion in Ga 2 O 3
,”
APL Mater.
7
,
022519
(
2019
).
92.
J. L.
Lyons
, “
A survey of acceptor dopants for β - Ga 2 O 3
,”
Semicond. Sci. Technol.
33
,
05LT02
(
2018
).
93.
Y. K.
Hommedal
,
Y. K.
Frodason
,
L.
Vines
, and
K. M. H.
Johansen
, “
Trap-limited diffusion of Zn in β - Ga 2 O 3
,”
Phys. Rev. Mater.
7
,
035401
(
2023
).
94.
Y. K.
Hommedal
,
Y. K.
Frodason
,
A.
Galeckas
,
L.
Vines
, and
K. M. H.
Johansen
, “
Broad luminescence from Zn acceptors in Zn doped β - Ga 2 O 3
,”
APL Mater.
12
,
021109
(
2024
).
95.
S.-D.
Lee
,
K.
Kaneko
, and
S.
Fujita
, “
Homoepitaxial growth of beta gallium oxide films by mist chemical vapor deposition
,”
Jpn. J. Appl. Phys.
55
,
1202B8
(
2016
).
96.
T. N.
Sky
,
K. M.
Johansen
,
V.
Venkatachalapathy
,
B. G.
Svensson
,
L.
Vines
, and
F.
Tuomisto
, “
Influence of Fermi level position on vacancy-assisted diffusion of aluminum in zinc oxide
,”
Phys. Rev. B
98
,
245204
(
2018
).