Theory of the Thermal Stability of Silicon Vacancies and Interstitials in 4H–SiC

: This paper presents a theoretical study of the electronic and dynamic properties of silicon vacancies and self-interstitials in 4H–SiC using hybrid density functional methods. Several pending issues, mostly related to the thermal stability of this defect, are addressed. The silicon site vacancy and the carbon-related antisite-vacancy (CAV) pair are interpreted as a unique and bistable defect. It possesses a metastable negative- U neutral state, which “disproportionates” into V + Si or V − Si , depending on the location of the Fermi level. The vacancy introduces a ( − / +) transition, calculated at E c − 1.25 eV, which determines a temperature threshold for the annealing of V Si into CAV in n-type material due to a Fermi level crossing effect. Analysis of a conﬁguration coordinate diagram allows us to conclude that V Si anneals out in two stages—at low temperatures ( T (cid:46) 600 ◦ C) via capture of a mobile species (e.g., self-interstitials) and at higher temperatures ( T (cid:38) 1200 ◦ C) via dissociation into V C and C Si defects. The Si interstitial (Si i ) is also a negative- U defect, with metastable q = + 1 and q = + 3 states. These are the only paramagnetic states of the defect, and maybe that explains why it escaped detection, even in p-type material where the migration barriers are at least 2.7 eV high.


Introduction
Silicon carbide (SiC) is currently a mature semiconductor that supports a variety of technologies, most notably in power electronics [1,2]. Among the merits of this material, particularly of its 4H polytype, we find high crystalline quality, a 3.2 eV wide bandgap, a large breakdown field, the possibility to grow both n-type and p-type layers, and exceptional thermal and mechanical stability.
SiC can be grown with a relatively low residual defect concentration (∼ 10 11 -10 13 cm −3 ). Intrinsic defects in SiC have recently received much attention from the community for their negative effects on device performance and for their promising role as building blocks in quantum technologies (e.g., quantum-bit holders, single-photon emitters, or quantumsensors) [3,4]. Most of these centers involve vacancies; they are usually introduced via electron or ion irradiation, and a precise understanding of their electronic and dynamic properties is of utmost importance, for instance, in order to control their positioning upon dynamic and thermal annealing treatments.
Due to the extreme radiation hardness and the low leakage current of SiC junctions, this material has been proposed for the fabrication of radiation detectors [5,6], including neutrons [7], allowing for operation under harsh conditions, for instance, at high temperatures and under intense radiation fields [8].
During exposure to high-energy radiation or high temperatures, the formation of intrinsic defects can play a decisive role regarding the resulting material specifications. A critical example is the impact of deep carrier traps on the minority carrier lifetime (due to non-radiative recombination at point and extended defects), and therefore on the functioning of bipolar devices.
Defects that are capable of trapping carriers can occur in at least two distinct charge states. These defects are termed acceptors or donors, whether they become negatively or V1/V1' photoluminescence lines (known to arise from transitions of the spin-3/2 state of V − Si ) as a function of the proton implantation fluence [32]. Again, this picture is not free of controversy. An obvious question we may ask is how consistent the annealing of S 1/2 is with the calculated barriers involving V Si migration and reactions. Hence, during annealing, we either have (i) an unknown defect X that becomes mobile and reacts with V Si , (ii) defect X is immobile and traps a fast-diffusing V Si , or (iii) V Si transforms into a different configuration (e.g., the CAV structure). The fastest among these mechanisms must have a barrier of E a = 1.8 eV. Additionally, in order to comply with the observed first-order kinetics, the concentration of defect X must be at least one order of magnitude higher than that of V Si . For that reason, (i) is hardly consistent with the assignment of X to a fast-diffusing self-interstitial, whereas (ii) and (iii) are not compatible with the large barriers calculated for V Si migration (E a > 2.6 eV) and conversion into CAV (E a > 3.2 eV) in 4H-SiC [33]. EH 1 (E c − 0.45 eV) and EH 3 (E c − 0.72 eV) traps have been detected in 4H-SiC irradiated with low-and high-energy electrons [26,27,[34][35][36], protons [30,37], and heavier ions and neutrons [38,39]. In order to clarify whether S 1/2 and EH 1/3 correspond to the same or distinct defects, Alfieri and Mihaila [36] investigated 116 keV electron irradiated material that had been oxidized in order to suppress the as-grown carbon vacancies. With these low-energy electrons, only carbon-displacements (vacancies, interstitials, and Frenkel pairs) could be introduced. EH 1 and EH 3 annealed out simultaneously at about 350 • C with an activation energy E a = 1.13 eV and 1.17 eV. Their concentration depth profiles corresponded to a one-to-one correlation, thus indicating that they belong to the same defect. The activation energies for annealing are about 0.7 eV lower than those of S 1 and S 2 , thus indicating that EH 1/3 and S 1/2 are not the same. Based on their observations, Alfieri and Mihaila [36] suggested that EH 1/3 is a carbon interstitial related defect.
Due to the injection of interstitials during the oxidation treatment, the concentration of Z 1/2 in the pre-irradiated samples of Ref. [36] dropped below the detection limit of DLTS. However, upon irradiation (with a dose of 5 × 10 16 cm −2 ) the introduction of the EH 1/3 the peak was about two times more effective than that of Z 1/2 . Considering that the latter corresponds to a double emission, the concentration of both EH 1/3 and Z 1/2 was essentially identical. Although this points to the assignment of EH 1/3 to a carbon interstitial, further work is necessary to clarify the introduction efficiency of EH 1/3 and Z 1/2 , in particular its dose dependency.
The M-center is a bistable defect observed in two configurations (A and B). It becomes evident in the DLTS spectra of electron-and proton-irradiated samples after gentle annealing at ∼200 • C for 30 mins [40]. Such treatments clean up the crowded spectra of the as-irradiated samples, and it is well possible that the defect is already formed at this stage. Configuration A is formed when the sample temperature is T > 20 • C and reverse bias is applied. This process is thermally activated (E a = 0.90 eV) and corresponds to a B → A transformation. On the other hand, A → B takes place at T > 140 • C with no bias, and the activation energy is now E a = 1.40 eV [41]. Configuration A is responsible for two DLTS peaks, M 1 (E c − 0.42 eV) and M 3 (E c − 0.83 eV), while configuration B has one peak referred to as M 2 (E c − 0.63 eV) [40,41]. The M-center anneals out in the range 310-370 • C according to first-order kinetics with an activation energy of E a = 2.0 eV [41]. The introduction of the M-center is also reported in a paper of this special issue [39], where samples were both ion-implanted and neutron-irradiated, adding further support for its connection to an intrinsic defect.
The bulk of previous theoretical work on migration and reaction mechanisms involving intrinsic defects in SiC was carried out for the simpler 3C allotrope (see for instance [23]). Although elucidative, we now know that these results do not necessarily apply to 4H-SiC [42,43]. In addition, the narrower gap of 3C or the use of too aggressive approximations, such as the local or semi-local treatments of the many-body electronic interactions, implies that they should be revised using more accurate methods. The same applies to results obtained for 4H-and 6H-SiC [33,[42][43][44][45]. This issue was recently investigated by Yan et al [46] using hybrid density functional theory (DFT), where it was found a much stronger charge state dependence of migration barriers of intrinsic defects in 4H-SiC than previously thought [23]. For instance, the migration barrier of V C increased from about 2.5 eV in the double-negative charge state to 5.5 eV in the double-positive charge state. Another interesting aspect is the anisotropic effect found by Bathen et al. [47], showing both theoretically and experimentally that the annealing of V C proceeds preferentially via migration of the defect along the basal plane of the 4H-SiC crystal.
The present paper attempts to consolidate some of the most recent findings and to revisit some of the open issues highlighted above, in particular by analyzing the charge-state dependence of transformation, dissociation, and migration of Si-related intrinsic defects in 4H-SiC. The data were obtained from highly accurate hybrid density functional calculations. After a short description of the methodology (Section 2), results are reported for the Si vacancy and self-interstitial (Section 3). In both cases, charge-state-driven bistability and migration mechanisms are discussed (Section 4).

Theoretical Methods
The electronic structure calculations reported below were carried out within rangeseparated hybrid density functional theory using the Vienna Ab-initio Simulation Package (VASP) [48][49][50][51]. The many-body exchange-correlation interactions were accounted for using the functional proposed by Heyd, Scuseria, and Ernzerhof (commonly named HSE06 functional) [52]. The projector-augmented wave method was used to treat the core electrons [53], whereas valence states were expanded in plane-waves with a cut-off energy of 420 eV. The latter was solved self-consistently until the total energy between two consecutive steps differed by less than 10 −7 eV.
Defects were inserted in 400-atom supercells of 4H-SiC with a hexagonal shape, constructed by replication of 5 × 5 × 2 primitive cells. The calculated basal and axial lattice constants were a = 3.0712 A and c = 10.051 A, respectively. Forces acting on atomic nuclei were obtained within the generalized gradient approximation (GGA) to the exchangecorrelation potential [54]. They were minimized with respect to the atomic coordinates until the maximum residual force was lower than 0.01 eV/Å. For the calculation of forces and structural relaxation (GGA-level), the Brillouin zone (BZ) was sampled on a Γ-centered 2 × 2 × 2 k-point grid (Γ−2 3 ). On the other hand, for total energy calculations (HSE06-level), the sampling was done at the Γ-point (see [47,55] for details).
The formation energy of a Si intrinsic defect in a state described by its net charge and atomistic structure (q, R) is obtained as a function of the Fermi energy (E F ) using the usual expression, where E def is the calculated energy of the defective supercell, E bulk is the energy of a bulk (defect-free) supercell, µ Si is the chemical potential of a Si atom removed from (plus sign) or inserted into (negative sign) the bulk supercell in order to obtain a vacancy or an interstitial, and E v is the valence band maximum energy, here taken as the highest occupied Kohn-Sham state of a bulk calculation. The chemical potential of silicon was calculated for a carbon-rich (silicon-poor) crystal, µ Si = µ 0 Si + ∆H 0 SiC , where µ 0 Si is the energy per silicon atom in pure silicon and ∆H 0 SiC = −0.61 eV the calculated heat of formation of 4H-SiC (to be compared to −0.72 eV from calorimetry measurements [56]). We note that E def includes a charge correction, calculated here according to the recipe of Freysoldt, Neugebauer, and Van de Walle [57].
Migration and transformation mechanisms were investigated using the climbing image nudged elastic band method (CI-NEB) [58]. While already reported elsewhere [47], a summary of the main workflow steps are as follows: the first step consisted of postulating a mechanism, which started and finished at specific minimum-energy configurations (previously obtained). Essentially, it formed a sequence of seven intermediate structures between the end-configurations. On the second step, an exploratory NEB-constrained calculation was performed using the Γ-point for sampling the BZ. This allowed for a relatively fast screening of the configurational space, close to the minimum energy path between the end-structures. In the third step, the relaxed image sequence from the previous step was further refined by turning on the climbing-image mode of the NEB algorithm and improving the forces by increasing the BZ sampling to Γ−2 3 . Finally, a self-consistent calculation within HSE06-level was carried out using the transition state geometry obtained from the previous step.
With the above specifications, quantities derived from the energy of (local) minimum structures (e.g., formation energies, transition levels) have an estimated numerical accuracy of about 10 meV. On the other hand, transition state calculations (activation energies for migration and structural transformation) have a numerical accuracy of 0.1 eV [47,55].

Results
Before moving on to the results, let us introduce a few useful words about the notation employed throughout. A defect state is referred to as S D q (R), where D stands for a stoichiometric label (V Si or Si i for a silicon atom subtracted or added to a perfect crystal, respectively), q = {· · · , −1, 0, +1, · · ·} is a charge state, S = {0, 1/2, 1, · · ·} the total spin, and R is a structural label to indicate a specific atomistic geometry. For instance, 3/2 V −1 Si (k) represents the lack of a pseudo-cubic (k) Si atom in 4H-SiC, with a trapped electron (net charge -1) and total spin S = 3/2. Should a generalization be made, or when the context makes the notation unambiguous, we may reduce the above to S D q , D q , D or simply R. Some defect structures can be described as a combination of elementary defects (e.g., antisites and vacancies). The complex pair formed by first-neighboring carbon-antisite and carbon-vacancy (CAV) can adopt up to four different alignments, depending on the lattice location of its elements. One of its alignments could be represented as C Si (k)V C (h). However, because its stoichiometry is identical to that of a silicon vacancy, we simply write V Si (CAV, kh), where first and second lattice indexes refer to the antisite and vacancy sites (not interchangeable).

The Silicon Vacancy
We investigated several configurations for V Si in 4H-SiC, of which the most stable are represented in Figure 1. The left-most structure represents a pristine piece of crystal for the sake of reference. White and dark grey spheres represent Si and C atoms, respectively. The first configuration is the site vacancy structure, V Si (k), which is simply the lack of a Si atom, in this case, removed from the pseudo-cubic site (represented by a translucent white sphere). The second structure is V Si (CAV), a carbon antisite-vacancy first neighboring pair (in this case, with an off-axis kh alignment). This structure "communicates" with V Si (k) via jumping of the C Si atom into the V C site. Another stable structure investigated, and also shown in Figure 1, is V Si (CAV2). It is also a carbon antisite-vacancy pair, but V C is on the second neighboring C-shell with respect to C Si (both lying on pseudo-cubic sub-lattice sites).
The relative stability of different V Si configurations are reported in Table 1. For each net charge of the defect, q, zero energy indicates the most stable state considering both atomistic and spin configurations. In agreement with previous results [22], the simple site configurations ({k, h}) are more stable in the negative charge states, while the CAV structure is more favorable when positively charged. Both configurations have four unsaturated radicals, and their electronic structure share an identical shape, i.e., a ↑↓ 1 a ↑ 1 e ↑ for V 0 Si within the C 3v point group (the singlet state within brackets is resonant with the valence band, while the other states are located in the bandgap). However, while site vacancies have four nearly equivalent and highly electronegative carbon radicals with states in the lower part of the gap, the CAV defect has a doublet localized on three nearly equivalent Si atoms with states high in the gap. Hence, the deep C-radicals of {k, h} site structures can trap more electrons than CAV, whereas the latter can trap more holes. It is noted that some defect structures mentioned above, notably 1 Si (CAV, {kk, hh}) and 1/2 V −1 Si (CAV, {kk, hh}) do not retain the full C 3v symmetry and display small monoclinic Jahn-Teller distortions of the order of tens of meV. In this case, the doublet splits into close a and a" singlets. The relative stability of different Si configurations are reported in Table 1. For each net charge of the defect, , zero energy indicates the most stable state considering both atomistic and spin configurations. In agreement with previous results [22], the simple site configurations ( { , ℎ} ) are more stable in the negative charge states, while the CAV structure is more favorable when positively charged. Both configurations have four unsaturated radicals, and their electronic structure share an identical shape, i.e., ↑↓ ↑ ↑ for Si within the point group (the singlet state within brackets is resonant with the valence band, while the other states are located in the bandgap). However, while site vacancies have four nearly equivalent and highly electronegative carbon radicals with states in the lower part of the gap, the CAV defect has a doublet localized on three nearly equivalent Si atoms with states high in the gap. Hence, the deep C-radicals of { , ℎ} site structures can trap more electrons than CAV, whereas the latter can trap more holes. It is noted that some defect structures mentioned above, notably Si / (CAV, { , ℎℎ}) do not retain the full symmetry and display small monoclinic Jahn-Teller distortions of the order of tens of meV. In this case, the doublet splits into close ' and '' singlets.
The site vacancy in 4H-SiC was found to be stable in charge states from 0 to −3. The strong localization of the C-radicals favors high spin states [59], namely ↑ ↑ , ↑ ↑↑ , ↑↓ ↑↑ and ↑↓ ↑↑↓ , for = 0, −1, −2 and −3, respectively, where we dropped the resonant ↑↓ level for the sake of simplicity. Conversely, , ↑ , ↑ ↑ , and ↑↓ ↑ were found for charge states = +2, +1, 0, and −1 of CAV, respectively, where the spin = 1 of the neutral state originates from a very small overlap of the ↑ level (strongly localized on the unique C atom) and the ↑ level (localized on the Si atoms and nodal on the C atom).
The CAV2 configuration, comprising a more remote C Si -C defect pair (see righthand side of Figure 1) can be thought of a carbon vacancy perturbed by the tensile strain field of a nearby inert C Si defect. We actually found that like C [14,16], CAV2 is stable in charge states from +2 to −2 with a negative-ordering of the acceptor levels. However, CAV2 is always metastable with respect to both CAV and { , ℎ} configurations by about 1 eV and it is not expected to be found under equilibrium conditions.
The energy of a neutral uncorrelated pair of C Si and C defects were also calculated as 1.1-1.2 eV higher than that of Si (CAV, ℎℎ), depending on the sub-lattice occupation of the defect components. This figure is representative of the binding energy of CAV. We  Table 1. Energies (in eV) of different silicon-vacancy structures in 4H-SiC, within the same charge state, with respect to the corresponding ground state (zero energy). Site (k, h) and axial CAV configurations (kk, hh) were considered. The total spin of each state is given between parentheses. The site vacancy in 4H-SiC was found to be stable in charge states from 0 to −3. The strong localization of the C-radicals favors high spin states [59], namely a ↑ 1 e ↑ , a ↑ 1 e ↑↑ , a ↑↓ 1 e ↑↑ and a ↑↓ 1 e ↑↑↓ , for q = 0, −1, −2 and −3, respectively, where we dropped the resonant a ↑↓ 1 level for the sake of simplicity. Conversely, a 0 1 e 0 , a ↑ 1 e 0 , a ↑ 1 e ↑ , and a ↑↓ 1 e ↑ were found for charge states q = +2, +1, 0, and −1 of CAV, respectively, where the spin S = 1 of the neutral state originates from a very small overlap of the a ↑ 1 level (strongly localized on the unique C atom) and the e ↑ level (localized on the Si atoms and nodal on the C atom).
The CAV2 configuration, comprising a more remote C Si -V C defect pair (see right-hand side of Figure 1) can be thought of a carbon vacancy perturbed by the tensile strain field of a nearby inert C Si defect. We actually found that like V C [14,16], CAV2 is stable in charge states from +2 to −2 with a negative-U ordering of the acceptor levels. However, CAV2 is always metastable with respect to both CAV and {k, h} configurations by about 1 eV and it is not expected to be found under equilibrium conditions.
The energy of a neutral uncorrelated pair of C Si and V C defects were also calculated as 1.1-1.2 eV higher than that of V 0 Si (CAV, hh), depending on the sub-lattice occupation of the defect components. This figure is representative of the binding energy of CAV. We finally note that V Si structures involving Si antisites, e.g., C Si Si C V Si obtained by displacing a Si atom into V C in CAV (see Figure 1) were not stable and relaxed back to CAV.

Bistability of the Silicon Vacancy
We now look at the relative stability of the several vacancy states in equilibrium at a specific temperature T. Figure 2a depicts the formation energy for the relevant states as obtained from Equation (1). Formation energies for site and CAV structures are shown in red and blue, respectively. Solid and dashed lines represent the energies of defects located on hexagonal and cubic sites, respectively. There is a total of four possible alignments for CAV. However, we only considered the axial defects, V Si (CAV, hh) and V Si (CAV, kk). The lower horizontal axis in Figure 2a shows the Fermi energy with respect to the valence band top, The upper horizontal axis marks the location of the Fermi level in the material (assumed to be n-type), for the temperatures indicated. For that, a donor concentration N D = 10 14 cm 3 was assumed. Accordingly, for an extrinsic semiconductor [60], where E i is the intrinsic Fermi level, n is the free electron density, and n i is the intrinsic carrier density,  The formation energy of the Si (CAV) state is rather low in p-type material (about 3-4 eV). Of course, this is still high enough to prevent the thermal generation of these defects at room temperature. However, the formation of CAV defects should be preferred in detriment of the site configuration, for instance in irradiated or ion-implanted p-type SiC. This was actually predicted experimentally by Konstantinov [61], who postulated the  In the above, E c,v are the energies of the conduction band bottom and valence band top denoted by c and v subscripts, respectively. These correspond to a density of states , where g c,v are band edge degeneracy factors (g c = 3 and g v = 1 in 4H-SiC), k B is the Boltzmann constant and h the Plank constant. In the density of states and Equation (3), m * e,h are the free electron and hole effective mass. The results of Figure 2a are in line with previous reports using hybrid density functional theory [22]. While the silicon vacancy adopts the CAV structure in p-type material, it displays the "simpler" site configuration in n-type material. The defect is amphoteric and strongly compensating. It can capture up to three electrons in n-type 4H-SiC, whereas in p-type material it captures up to two holes. The site configuration has three acceptor transi- These are average values obtained from the calculated levels of cubic and hexagonal configurations (c.f. Table 2). On the other hand, the CAV configuration has donor levels at E(+/ + +) = E c − 1.98 ± 0.03 eV and E(0/+) = E c − 1.03 ± 0.01 eV, and one acceptor transition at E(−/0) = E c − 0.44 ± 0.03 eV.  The formation energy of the V +2 Si (CAV) state is rather low in p-type material (about 3-4 eV). Of course, this is still high enough to prevent the thermal generation of these defects at room temperature. However, the formation of CAV defects should be preferred in detriment of the site configuration, for instance in irradiated or ion-implanted p-type SiC. This was actually predicted experimentally by Konstantinov [61], who postulated the formation of antisite-type defects in the context of the formation of shallow and deep boron acceptors in Si and C sites, respectively.
It is noteworthy from Figure 2a that the neutral vacancy is not stable (under thermodynamic conditions), irrespective of is configuration and location of the Fermi level. The reason is that the following disproportionation reaction is exothermic.
Although the energy of the intermediate state with respect to reactants is positive, U = 0.63 eV (which is the Hubbard correlation energy of the neutral CAV state, in this case for the kk alignment), this figure is lower in magnitude than the energy release of the final step due to the V − Si (CAV) → V − Si reconfiguration. In the end the effective U = −0.46 eV. This feature is what distinguishes a negative-U defect [9], thus implying the existence of a (−/+) transition, here calculated at E c − 1.25 ± 0.02 eV.
From the upper horizontal axis of Figure 2a we arrive at the conclusion that in ntype material and above ∼600 • C the CAV structure becomes the most favorable state (positively charged). This suggests that an annealing/quenching process could be used in order to create a non-equilibrium population of CAV defects in n-type 4H-SiC. We note that although being exothermic, the kinetics of Reaction (6) strongly depends on the emission/capture rates of carriers, as well as on transformation barriers along the process. This will be discussed further below. In p-type 4H-SiC the site structure of V Si is invariably unstable against CAV.
The left-hand side of Figure 3 depicts a configuration coordinate diagram of the silicon vacancy in n-type 4H-SiC. It shows the potential along the transformation coordinate V Si (k) ↔ V Si (CAV, kk) only, although that will suffice for our needs. The barriers involving transformations between V Si (h) and V Si (CAV, hh) systematically exceeded those of Figure 3 by ∼0.5 eV. Barriers for all charge states are indicated on the trail of the arrows along the potential energy surfaces. Vertical arrows indicate the energy difference between the respective energy minima. Note that V −1 Si and V −1 Si (CAV) have a different spin in their ground states (S = 3/2 and 1/2, respectively). This means that a spin-flipping event will take place somewhere along the conversion 1/2 V −1 Si (CAV) ↔ 3/2 V −1 Si . Our calculations assume that single atomic jumps occur on a much shorter time scale than that involving the change of the spin state [62]. Hence, transformations are expected to conserve the spin of the initial state.  On the other hand, if a Si (CAV) defect is found in n-type 4H-SiC without an external bias, it will be "eager" for electrons. First of all, it will quickly capture two electrons and become negatively charged, Finally, it is noted that previous semi-local calculations gave a qualitatively similar picture for the silicon vacancy, including the acceptor character and the preference for the CAV configuration in p-type material [59]. However, and despite its merits, the semi-local approach used to describe the exchange-correlation functional lacks accuracy, making it difficult to draw quantitative conclusions. For instance, the calculated electronic transition levels reported in [59] suffer from the same insufficiency that underestimates the bandgap width by about 40%. By employing the experimental bandgap in the construction of the formation energy diagram, the calculated acceptor levels came out too deep (with respect to c ) and actually, a (−4/−3) transition was artificially predicted. For the same reason, the calculated Si ↔ Si (CAV) conversion barriers were about 30-50% smaller than those The V Si → V Si (CAV) transformation essentially involves a jump of a carbon atom from the edge of the Si vacancy into the empty site. This may occur either in the neutral or negative charge states. Only the former leads to an exothermic reaction. Other charge states, namely V −3 Si and V −2 Si must emit electrons into the conduction band or capture holes from the valence band before making the jump. The barriers for jumping into the CAV structure are rather high, namely 4.3 eV and 3.6 eV for charge states q = −1 and q = 0, respectively.
Another distinct feature of a negative-U defect is the double emission (or double capture) of carriers during transients upon changing the Fermi level. For instance, and according to Figure 3, if starting from the negative charge state in the k configuration in a reversebiased n-type diode, the ionization into the most stable neutral state, which adopts the CAV structure, involves a net change in the energy of ∆E(−/0) = E c − E(−/0) = 1.51 eV. If the temperature of the sample is high enough for this step to proceed within an observable timescale, an additional ionization into charge state q = +1 will follow at an even faster rate because it involves an electron emission with ∆E(0/+) = 1.04 eV only. Of course, the experimental verification of this effect depends on the intermediate transformation barriers, and as we will see further ahead, on the competition with the migration and annealing out of the defect itself.
The thermodynamic drive that controls the V Si ↔ V Si (CAV) reaction direction is the Fermi energy. From Figure 2, we find that in n-type 4H-SiC and up to T ∼ 600 • C, the vacancy is negatively charged (q = −3, −2, or −1) in the site configuration. Above this temperature, the positive CAV state is more favorable. Supposing that the Fermi level is located at mid-gap and the vacancy is in the 3/2 V −1 Si (non-equilibrium) state, Figure 3 tells us that the CAV structure is only reachable if a jump over a barrier of about 4.3 eV is performed. Alternatively, it could first emit an electron into the conduction band (with an activation energy of about 2.2 eV), followed by a jump in the neutral charge state over a 3.6 eV barrier. Neither of these alternatives is likely to occur below T ∼ 1000 • C in a timescale of seconds (compare with the migration of the carbon vacancy, which is observed above 1200 • C and has a measured barrier of 3.6 eV [47]).
On the other hand, if a V + Si (CAV) defect is found in n-type 4H-SiC without an external bias, it will be "eager" for electrons. First of all, it will quickly capture two electrons and become negatively charged, 1/2 V −1 Si (CAV). Further energy lowering is only possible by overcoming a 3.3 eV barrier via 1/2 V −1 Si (CAV) → 3/2 V −1 Si . After that, additional electrons can be captured depending on the Fermi level.
Finally, it is noted that previous semi-local calculations gave a qualitatively similar picture for the silicon vacancy, including the acceptor character and the preference for the CAV configuration in p-type material [59]. However, and despite its merits, the semi-local approach used to describe the exchange-correlation functional lacks accuracy, making it difficult to draw quantitative conclusions. For instance, the calculated electronic transition levels reported in [59] suffer from the same insufficiency that underestimates the bandgap width by about 40%. By employing the experimental bandgap in the construction of the formation energy diagram, the calculated acceptor levels came out too deep (with respect to E c ) and actually, a (−4/ − 3) transition was artificially predicted. For the same reason, the calculated V Si ↔ V Si (CAV) conversion barriers were about 30-50% smaller than those reported here.

Migration of the Silicon Vacancy
As reported by Bockstedte et al. [63], migration of a silicon vacancy in SiC must occur via exchange with Si second neighbors. The right-hand side of Figure 3 depicts the calculated migration barriers of the vacancy in 4H-SiC within the cubic sub-lattice. All charge states are included. Solid and dashed lines represent potential energy surfaces for migration along the basal and axial directions, respectively. We note that each axial mechanism involves a total of four consecutive jumps. Values in Figure 3 correspond to the energy of the highest saddle point along the way with respect to the initial and final states. Two features are readily noticeable-(i) the migration barrier increases as the defect emits electrons and becomes less negative and (ii) basal migration barriers are invariably lower than their axial counterparts (for the same charge state). This suggests that, like the carbon vacancy, diffusion of silicon vacancies could be highly anisotropic.
Before jumping to the conclusions, we should bear in mind several effects that compete with the jumps of the vacancy, namely, electron emission and capture, in addition to conversion into CAV and subsequent dissociation into a remote C Si -V C pair via migration of V C . We know that V C migrates only at T 1200 • C, at which point the Fermi level is close to the middle of the gap, and for which the ground state of the Si vacancy is actually V + Si (CAV). From here, it could capture a free-electron (releasing 1 eV) and convert it into the metastable V 0 Si species by overcoming a 4.2 eV high barrier and finally perform a migration jump (over a 4 eV barrier). Considering V + Si (CAV) as the zero-energy reference, the overall energy barrier for migration could be as low as 3.7 eV if the energy that results from the electron capture was not dissipated via phonon emission, but rather converted into kinetic energy and assist the Si nuclei to perform the migration jump. However, if the energy drop from the capture is dissipated, the migration barrier increases to 4.7 eV. It is also noted that at least two kinetic effects are expected to hinder V Si migration. Firstly, the capture of a carrier depends on their availability and cross section, and secondly, the V 0 Si → V 0 Si (CAV) return barrier is considerably lower than the migration barrier.
Alternatively, and again assuming that the Fermi energy is close to mid-gap (high temperatures), V + Si (CAV) can convert into V + Si (CAV2) (show in Figure 1), and from there V C could migrate in the positive or neutral charge state. Considering that the energy of uncorrelated neutral C Si and V C defects are about 1 eV higher than V Si (CAV), and adding this to the 3.6 eV migration barrier of V C [47], we arrive at a dissociation barrier for the CAV configuration of about 4.6 eV.

Silicon Interstitials
Now we turn to the silicon self-interstitials (Si i ). Several structures were investigated, among which we have Si atoms at empty sites (cage structures), Si-Si and Si-C split interstitials. The cage structures are denoted by the sublattice layer where the interstitial Si atom is located. Four of these sites are represented graphically in Figure 4 (left). Two types of Si-Si split structures were found, namely those with the Si-Si dimer direction possessing a component along the hexagonal axis of the crystal (referred to as axial or A-structures), and those whose Si-Si dimer is parallel to the basal plane of the crystal (referred to as basal or B-structures). Examples of these are depicted in the middle and right-hand side of Figure 4, respectively. Si-C split structures were found to be unstable and relaxed to cage configurations. These findings were already reported in previous studies [46,63,64].  Split-interstitial structures A and B are also divided into two groups, namely those occupying and ℎ sublattice Si sites. Figure 4 shows Si i (A, ℎ) and Si i (B, ) configurations for the sake of exemplification. Table 3 reports the relative stability of the lower energy Si i defects in 4H-SiC for several charge states. As it will be shown below, acceptor states are not stable, and therefore only neutral and positive charge states are reported. The energy of the most stable configuration in each charge state is set to zero. Energies for ′ and ℎ′ configurations are not included in the table. These structures are either unstable (in = 0 and +1 charge states) or their relative energies are way too high ( > 2.5 eV for = +2, +3, and +4 states) with respect to other competing structures. All cage structures retain the maximum site symmetry ( ), except the Si i (ℎ) state, which undergoes a distortion upon symmetry-unconstrained relaxation, possibly due to a pseudo-Jahn-Teller effect. For this reason, it is referred as ℎ * , and its relative energy in Table 3 is also labeled with a star. Split-interstitial structures A and B are also divided into two groups, namely those occupying k and h sublattice Si sites. Figure 4 shows Si i (A, h) and Si i (B, k) configurations for the sake of exemplification. Table 3 reports the relative stability of the lower energy Si i defects in 4H-SiC for several charge states. As it will be shown below, acceptor states are not stable, and therefore only neutral and positive charge states are reported. The energy of the most stable configuration in each charge state is set to zero. Energies for k and h configurations are not included in the table. These structures are either unstable (in q = 0 and +1 charge states) or their relative energies are way too high (E > 2.5 eV for q = +2, +3, and +4 states) with respect to other competing structures. All cage structures retain the maximum site symmetry (C 3v ), except the Si 0 i (h) state, which undergoes a distortion upon symmetry-unconstrained relaxation, possibly due to a pseudo-Jahn-Teller effect. For this reason, it is referred as h * , and its relative energy in Table 3 is also labeled with a star. Table 3. Energies (in eV) of different silicon self-interstitials in 4H-SiC, within the same charge state, with respect to the corresponding ground state (zero energy). Cage-like ({k, k , h}), splitaxial (A,{k, h}) and split-basal (B,{k, h}) structures were considered. Empty cells indicate that the respective states were unstable. All states have a low-spin configuration (S = 0 and S = 1/2 for even and odd charge states, respectively). We can readily draw a few conclusions from Table 3. Overall, cage structures are more stable in positive charge states, while split-interstitial forms are more favorable in the neutral state. The Si i defect almost exclusively inhabits the h layer of 4H-SiC. The exception is the Si +4 i state, which stabilizes in the k configuration in p-type material. As the defect becomes more positively charged (moving from left to right in Table 3), the Si i (A, h) state becomes less stable in favor of Si i (h), and the latter becomes less stable in favor of Si i (k). We can rationalize these findings in the following way: the fewer electrons bound to the Si-Si split defect, the less stable the bonds with ligands become, and with further ionization, the internal crystal field across the large h-cage tends to stabilize interstitial Si cations towards the more electronegative carbon anion in the 4H-SiC lattice (upper left region of Figure 4).
Blank cells in Table 3 mean that the respective states are not stable. Hence, cage structures can emit up to four electrons into the conduction band bottom (or capture four holes from the valence band top). On the other hand, split interstitial structures can only emit up to two electrons. The Si 2+ i (A, h) state is not stable and relaxes into the Si 2+ i (h) cage form.

Bistability of the Silicon Interstitial
We can get a better perspective of the above results by looking at the right-hand side of Figure 2. Here we depict a diagram with formation energies of the cage (blue lines), split-axial (red lines), and split-basal (green lines) configurations as a function of the Fermi level (bottom horizontal axis). Solid lines and dashed lines represent the formation energies of defects located at h and k sublattice sites. The black line labeled h * stands for the neutral distorted Si 0 i (h * ) as described above. The energy of the non-distorted trigonal Si 0 i (h) configuration is 0.9 eV higher than that of Si 0 i (h * ). The top horizontal axis represents the Fermi level location at the temperatures indicated in a sample with a free-electron density of n = 1 × 10 14 cm −3 .
In general, the formation energy of the Si vacancies is considerably lower than that of Si self-interstitials. This is particularly noticeable in n-type material and especially relevant under intrinsic conditions (high temperatures) as it has an important impact on self-diffusion. With the Fermi level at mid-gap, Si vacancies and interstitials are most stable in the positive V + Si (CAV) and double-positive Si +2 i (h) states, respectively. Hence, mutual repulsion will hider their annihilation should any of them become mobile. However, in n-type material in the range of T ≈ 400-600 • C, the vacancy adopts the site structure V − Si and could be an efficient trap for mobile Si +2 i defects. In n-type 4H-SiC the Si interstitial adopts the split Si 0 i (A, h) state (solid red horizontal line in Figure 2). Among competing configurations, at about 0.3 eV higher in the energy scale we find the split-basal structures and Si 0 i (h * ). On the other hand, in p-type 4H-SiC, Si 4+ i (h) is the most favorable state. Hence, while Si vacancies compensate n-type material, the deep donors of the Si interstitial are expected to strongly compensate p-type 4H-SiC as well.
The carbon self-interstitial has two negative-U metastable charge states, respectively; q = +1 and q = +3. These are only expected to form under non-equilibrium conditions, e.g., under illumination or thermal excitation. The effect is clearly illustrated in Figure 5, where a configuration coordinate diagram is presented for the Si i defect in n-type 4H-SiC. Only states that are relevant for transformations and migration of Si i are included in the diagram. Split-basal configurations are metastable in the q = 0, +1, and +2 charge states.
configurations can differ by a basal lattice vector. Along the hexagonal direction, they differ by half lattice vector only. However, the remaining half path is identical to the first half by symmetry if we consider a rotation of the Si-Si dimer in Si i (A, ℎ) by 2 /3 (passing by Si i (B, ℎ) ), which has a barrier of 1.3 eV . The transition-state energies indicated over the forward arrows in Reaction (7) are all relative to the Si i (A, ℎ) ground state. Overall, the migration barrier is a,m = 1.54 eV (limited by the first step), and this is summarized on the left-hand side of the configuration coordinate diagram of Figure 5.
In intrinsic material (with the Fermi level close to mid-gap) the Si i (ℎ) state is the most favorable. In this case, we found that "cage-jumping" leads to the lowest transition state energies. In this case, the mechanism is which accounts for axial and basal migration as well. The end-configurations differ by one basal lattice vector and half axial vector. The saddle point is actually a structure close to the Si i ( ′) state. Again, energies over reaction barriers are relative to the initial (final) ground state. A simplified potential energy surface for the mechanism is depicted in the upper left part of Figure 5.  In equilibrated n-type material, the defect is found in the Si 0 i (A, h) state. If the defect is provided with at least 1.45 eV (thermally or optically) so that it becomes ionized as Si 0 i (A, h) → Si + i (A, h) + e − , a subsequent transformation into a slightly lower energy state, Si + i (h), followed by further emission Si + i (h) → Si 2+ i (h) + e − will take place immediately. This is because the transformation barrier is only 0.2 eV high and the second ionization costs 0.5 eV, which must less than the first ionization. Likewise, hole emission from Si 4+ i (k) → Si 3+ i (k) + h + involves the absorption of about 1.6 eV (which is the position of the E(+3/ + 4) − E v transition involving the cage interstitial in the k site). However, the ionized Si 3+ i (k) species is not the ground state anymore, and an eventual transformation into the lower energy Si 3+ i (h) state leaves the defect in a state with a hole ionization energy of only 0.75 eV (which is the position of the E(+2/ + 3) − E v transition involving the cage interstitial in the h site). It is noted that in this case the structural transformation involves overcoming a large barrier of 2.1 eV, and as we will see, that is higher than the activation energy for migration in this particular charge state.

Migration of the Silicon Interstitial
The reporting of migration mechanisms and barriers of silicon interstitials is now carried out for n-type and intrinsic conditions. Results for p-type conditions are presented at the end of this section. Similar to the vacancy, we assume that the interstitials do not exchange carriers with the crystalline states in the course of elementary jumps. This may only occur when the defect "lands" at intermediate (meta-)stable states along the mechanisms, whose lifetime is much longer than the charge carrier emission/capture periods. Figure 2 (right) suggests that up to about 400 • C in n-type 4H-SiC, the Si i defect is expected to adopt the Si 0 i (A, h) state. The lowest energy path for migration of this species was found to be which accounts for both axial and basal migration. This is because the first and final configurations can differ by a basal lattice vector. Along the hexagonal direction, they differ by half lattice vector only. However, the remaining half path is identical to the first half by symmetry if we consider a rotation of the Si-Si dimer in Si 0 i (A, h) by 2π/3 (passing by Si 0 i (B, h)), which has a barrier of 1.3 eV. The transition-state energies indicated over the forward arrows in Reaction (7) are all relative to the Si 0 i (A, h) ground state. Overall, the migration barrier is E a,m = 1.54 eV (limited by the first step), and this is summarized on the left-hand side of the configuration coordinate diagram of Figure 5.
In intrinsic material (with the Fermi level close to mid-gap) the Si +2 i (h) state is the most favorable. In this case, we found that "cage-jumping" leads to the lowest transition state energies. In this case, the mechanism is which accounts for axial and basal migration as well. The end-configurations differ by one basal lattice vector and half axial vector. The saddle point is actually a structure close to the Si 2+ i (k ) state. Again, energies over reaction barriers are relative to the initial (final) ground state. A simplified potential energy surface for the mechanism is depicted in the upper left part of Figure 5.
As shown in Figure 2, in p-type material the ground state is Si +4 i (k). For these conditions, the migration mechanism is still of cage-jumping type and it is now Overall the migration barrier is E a,m = 3.71 eV high. The first step involves a purely axial jump of the Si 4+ i ion and it is the limiting one. The second step, in addition to an axial component, involves the displacement of Si 4+ i by a basal lattice vector. Hence, the mechanism accounts for both basal and axial migration.

Discussion
The calculated E(−3/ − 2) and E(= /−) levels of V Si in the site configuration at about E c − 0.33 eV and E c − 0.74 eV, are rather close to the depth of S 1 and S 2 traps measured at E c − 0.41 eV and E c − 0.71 eV, respectively. The large barriers that were found for the migration of V Si in different charge states (>3 eV) and those for the conversion into the CAV form (>3.5 eV) are not compatible with the observed first-order annealing kinetics with activation energy E a = 1.8 eV for S 1/2 . Such low activation energy suggests an annihilation mechanism due to the capture of mobile self-interstitials by the vacancies. The figure is not far from the calculated migration barrier of neutral Si i defects (1.54 eV). However, in order to account for the first-order kinetics, the concentration of interstitials has to be much larger than that of vacancies. One possibility is that carbon interstitials (which are easier to be displaced from the lattice in irradiated material) are the main culprits for the annihilation of V Si during anneals at above 300 • C.
In n-type material, at higher temperatures (T 600-800 • C), any available V Si defect that escaped annihilation at lower temperatures, is now more stable in the positive CAV form. This suggests that an annealing stage driven by a Fermi level effect should occur in this temperature range. According to Figure 3, migration of V Si could occur via the following mechanism, which is strongly hampered by the need of capturing an electron (or emitting a hole) in step 1, plus the fact that the barrier for the back-reaction of step 2 is lower than the forward barrier (3.58 eV against 4.22 Ev; c.f. Figure 3). The overall barrier of Reaction 10 is estimated as high as E a,m = 4.7 eV. Hence, considering the binding energy of C Si and V C of 1 eV with respect to CAV, plus the migration barrier of V C (3.6 eV), it is likely that the high-temperature stage for the annealing of the Si vacancy occurs via dissociation of CAV with an activation barrier of about E a,d = 4.6 eV. Annealing via capture of a mobile species (e.g., self-interstitials from dissolving interstitial aggregates) is also possible.
The CAV defect has calculated donor levels at E(0/+) = E c − 1.03 eV and E(+/ + +) = E c − 1.98 eV. It also has a calculated acceptor transition at E(−/0) = E c − 0.44 eV. As recently suggested by Karsthof et al. [19], the first donor level is rather close to the EH 4 and EH 5 electron traps. Other potential DLTS signals that could be assigned to CAV(+/ + +) and CAV(0/+) transitions are HK4 and EM1 [65]. These are hole and electron traps found in p-type 4H-SiC with levels at E v + 1.56 eV (HK4) and E v + 2.26 eV (EM1), respectively. Considering the calculated band gap of 3.3 eV, we place the calculated CAV donor levels at E(+/ + +) = E v + 1.32 eV and E(0/+) = E v + 2.27 eV. HK4 and EM1 can actually be found in as-grown p-type 4H-SiC with a concentration of the order of 10 12 cm −3 [65], thus being consistent with the low formation energy of CAV in the p-type material. Considering their relative locations in the gap, it is possible that EM1 (minority carrier trap in p-type material) and EH 4 (majority carrier trap in n-type material) could originate from the same defect transition. In any case, we are still left with an additional CAV level of acceptor type, predicted at E c − 0.44 eV (which should be easy to detect DLTS), but for which there is no experimental match with compatible intensity and thermal stability to the other peaks.
Interestingly, in n-type 4H-SiC, silicon vacancies have several acceptor states, whereas interstitials have several donor states. Their multiple opposite charges are expected to drive a strong Coulomb attraction, and therefore to promote annihilation. On the other hand, in p-type and intrinsic material, both vacancies (with the CAV structure) and interstitials are positively charged, and annihilation should be hindered.
During epitaxial growth and higher temperatures (T 1400 • C), the material is intrinsic. Under these conditions, the formation energy of the Si vacancy is about 6.2 eV and that of the Si interstitial is about 8.7 eV. Assuming that (i) V Si (CAV) dissociates before being able to migrate in the form of a site vacancy and (ii) the migration barrier of Si +2 i is E a,m = 2.7 eV, we arrive at an estimate for the activation energy for self-diffusion of silicon, which is solely due to the motion of silicon interstitials, of E Si a,sd = 11.4 eV. Within this picture, V Si defects do not contribute to self-diffusion. This exceeds the measured analogous quantity for carbon self-diffusion by about 4 eV [66,67], and it is consistent with the observation that silicon diffusivity is more than two orders of magnitude lower than that of carbon [68][69][70].
Regarding the migration of neutral Si interstitials in n-type 4H-SiC, the proposed mechanism is in line with that found for the migration of Si i in 3C-SiC by Bockstedte et al. [63], where a 〈110〉-kick-out sequence was found to be the most favorable path. It is noted that a strictly identical sequence is not possible in 4H-SiC due to the lower symmetry of the host crystal. That is why the Si 0 i (h * ) state (which is not far from a split configuration), appears along the way. In Ref. [63], the kick-out migration mechanism of the neutral splitinterstitial in 3C-SiC was reported with an activation barrier E a,m = 1.4 eV. Our results indicate that in 4H-SiC the barrier is about 0.1-0.2 eV higher.
For p-type 3C-SiC, the authors of Ref. [63] found a migration barrier for Si 4+ i of E a,m = 3.5 eV. Again, this figure is about 0.2 eV lower than the analogous barrier found in this work for 4H-SiC. In a more recent study, Yan et al. [46] found migration barriers for the Si i defect lower than 0.9 eV for the neutral state and about 1.8 eV for the q = +4 charge state. These results, which employed methodologies identical to those in this paper, imply a considerably faster diffusivity and lower temperature for annealing of Si i . Unfortunately, the details provided in Ref. [46] regarding the mechanisms involved are rather limited, and therefore it is rather difficult to provide a comment other than suggesting that this issue needs to be revisited.
The formation of interstitial clusters in SiC is an important issue. Hornos et al. [71] found that silicon interstitials can aggregate and form stable and electrically active complexes. They suggested that at high temperatures, the Si clusters can dissolve and emit interstitials with significant consequences to the formation of other defects. Our calculated migration barriers indicate that the kinetics of these processes should be much faster in n-type than in p-type material.

Conclusions
This paper presents a theoretical study of the transformation and migration mechanisms of Si vacancies and self-interstitials in 4H-SiC using first-principles methods. Particular attention was given to the investigation of charge state effects on the potential barriers.
Considering that the vacancy has two distinct configurations (site and CAV structures) in n-type and p-type material, it is shown that the neutral state is metastable, irrespectively of the position of the Fermi level. That is, V Si has a negative-U in the neutral state, giving rise to a (−/+) transition estimated at E c − 1.25 ± 0.02 eV. This level is likely to define the threshold for the annealing of V Si into the CAV structure in n-type material due to a displacement of the Fermi level deep into the gap. V Si is also an amphoteric center that strongly compensates both p-and n-type material.
The calculated electronic transitions and thermal stability of CAV are consistent with the DLTS emission energy and annealing of EH 4 . Accordingly, EH 4 is stable up to 1200 • C, which is interpreted as the migration of V C away from C Si . The dissociation barrier of CAV is estimated as E a,d = 4.6 eV. Si vacancies are unlikely to anneal out via migration of the site structure. That route was shown to be hindered by the need for the emission or capture of a carrier. Additionally, at high temperatures (intrinsic conditions), Si vacancies and interstitials are most stable in the positive V + Si (CAV) and double-positive Si +2 i (h) states, respectively. Hence, mutual repulsion will hamper their annihilation should any of them become mobile. However, in n-type 4H-SiC in the range of T ≈ 400-600 • C, the vacancy adopts the site structure V − Si and could be an efficient trap for mobile Si +2 i defects as well as carbon interstitials.
The Si I defect shows a strong charge-state dependent migration barrier. In n-type material the defect is neutral, and the mechanism consists of a sequence of kick-out jumps involving Si-Si split structures with an overall barrier of 1.5 eV. On the other hand, in intrinsic and p-type material the stable charge states are q = +2 and q = +4, respectively. In that case, the migration mechanism involves a jump of Si i between neighboring cages with barriers of 2.7 eV and 3.7 eV, respectively.
We find that the Si i is a negative-U defect with metastable q = +1 and q = +3 states. These are the only paramagnetic states of the defect and that could explain why it escaped detection by electron paramagnetic resonance, even in p-type material where the migration barrier is at least 2.7 eV high. Of course, other techniques are in principle sensitive to the presence of Si i in 4H-SiC, including those related to the absorption or luminescence involving electronic transition within gap states.

Funding:
The author thanks the support of the i3N projects, Refs. UIDB/50025/2020 and UIDP/50025/2020, financed by the Fundação para a Ciência e a Tecnologia in Portugal. This work was supported by the Science for Peace and Security NATO Program through project SPS G5674.