Engineering the formation of spin-defects from first principles

The full realization of spin qubits for quantum technologies relies on the ability to control and design the formation processes of spin defects in semiconductors and insulators. We present a computational protocol to investigate the synthesis of point-defects at the atomistic level, and we apply it to the study of a promising spin-qubit in silicon carbide, the divacancy (VV). Our strategy combines electronic structure calculations based on density functional theory and enhanced sampling techniques coupled with first principles molecular dynamics. We predict the optimal annealing temperatures for the formation of VVs at high temperature and show how to engineer the Fermi level of the material to optimize the defect’s yield for several polytypes of silicon carbide. Our results are in excellent agreement with available experimental data and provide novel atomistic insights into point defect formation and annihilation processes as a function of temperature.


I. INTRODUCTION
Spin defects in wide bandgap semiconductors are promising platforms for several quantum technologies, including quantum photonics, and quantum sensing and communication [1,2].Among spin qubit hosts, in recent years silicon carbide (SiC) has emerged as an ideal material, due to mature growth, doping and fabrication techniques [3,4], with qubits realized with silicon vacancies (V Si ) and nitrogen-vacancy pairs (NV), divacancies (V C V Si ), and carbon antisite vacancies (CAV).The V C V Si in SiC (which we denote as VV) has attracted particular interest, due to its optical addressability [5], a near-infrared spin-photon interface [6], long coherence times [7] and high-fidelity readout via spin-to-charge conversion [8].While numerous studies of defects in SiC have focused on their physical properties, much less is known about their formation processes, whose control is critical for the integration of semiconductors hosting spin qubits within electronic and optical devices [1,3,4,[9][10][11][12].
Defects in SiC are usually generated via implantation, irradiation or pulse laser, and by subsequent thermal annealing at high temperature [1,4].Several experimental methods have been used to monitor defect formation, including electron paramagnetic resonance (EPR), photoluminescence (PL), and deep-level transient spectroscopy [13][14][15][16][17][18].Recent progress has been reported in achieving spatial localization of defects [11], as well as in controlling their charge state [19], performance and yield [14,20].In the case of the VV, one of the most studied defects in SiC, it is well established that n-doping conditions are beneficial to its formation [16,[21][22][23], and a lower bound for the annealing temperature (T Ann ) required to generate VVs, ∼ 1,000 K, has been estimated experimentally [16,[24][25][26][27].However, different experiments have reported different temperatures [8,15,24,[26][27][28][29][30][31][32][33], with an optimal T Ann often quoted around 1,150 K [24,25,34,35].The experimental determination of activation and optimal annealing temperatures remains a challenging task, because these quantities are usually inferred from the intensity of EPR/PL signals which are affected by several factors, including the charge state and concentration of defects [36], Fermi-level position (E F ) [15,27], and specific synthesis conditions [37,38].Recently, the pairing of V C and V Si into neutral VVs has been investigated theoretically, providing the first atomistic insight into the formation process [39].However, as is the case for most point defects in semiconductors, our understanding of the VV formation mechanism at the atomistic level remains preliminary and qualitative.In particular, a relation between the host E F and T Ann has not yet been established, which is of great importance to control defects' formation, and an upper bound to T Ann is yet unknown.
Moreover, the dynamics of VV [23,39], the conditions for the defect immobilization in the lattice and the effect of temperature on formation processes are only partially understood.
Addressing these open problems is difficult from an experimental standpoint, especially in the presence of limited microscopic resolution, and atomistic simulations are key tools to gain detailed insights.
Here we present a general computational protocol, based on first principles calculations, to study the formation of point defects in covalently bonded materials; in addition we provide specific predictions on optimal conditions for the formation of double vacancies in SiC.We focus on the cubic phase (3C-SiC) for its simplicity with only one type of lattice site, and we discuss implications of our results for hexagonal polytypes.We determine the preferred pathways leading to the VV formation and optimal values of T Ann and E F , and we elucidate the interdependence of these parameters.Our results point at the importance of considering multiple charge states of defects, as well as of configurations that are not thermodynamically stable, for accurate predictions of formation pathways.On the other hand, the sampling of paths with different spin states has a negligible impact on our predictions.

A. Computational strategy
Pathways -We studied defect dynamics and transformations during the thermal annealing process following defect generation by e.g.particle irradiation.We considered several possible processes relevant to the formation of VVs in 3C-SiC, based on previous studies [21][22][23]39], and on our chemical intuition; they are summarized in Fig. 1.In addition to formation, we also considered dissociation processes, specifically CAV → C Si + V C (where C Si is an isolated carbon antisite) and VV → V C + V Si , which involve multi-step migrations of mono-vacancies (MV) (not explicitly shown in Fig. 1).
We did not consider interstitial (e.g.Si or C interstitial), and substitutional (e.g.N substitution or C Si ) defects; the former are expected to be annealed out once the paths described in Fig. 1 occur [22,23], and the latter are immobile at ∼ 1,000 K [40].Simulation protocol -The simulation protocol used in our work is presented in Fig. 2. We studied the processes displayed in Fig. 1 using density functional theory (DFT) calculations with both the Perdew-Burke-Ernzerhof (PBE) and dielectric-dependent hybrid (DDH) functionals, and we considered several charge (q) and spin (s) states (see Methods).Specifically, we considered different s states to determine the minimum energy path and energy barrier E b as a function of q, and for a given pathway, we obtained an effective barrier, E b, EFF , as a function of the Fermi level E F : where  of E F , while E b exhibits steps at charge transition levels.The expression of E b, EFF in Eq.
1 assumes that charge state equilibration processes are faster than the transformation of defects into different configurations.We verified the validity of this assumption at high T (∼ 1,000 K; see SI).We emphasize that thermodynamically unstable q states may participate and play an important role in defect transformation processes, since exploring those states may lead to lower effective barriers.
As mentioned above, the Fermi level is a parameter in Eq. 1, and we estimated the experimental conditions that may lead to specific, desired values of E F based on charge neutrality conditions and the electronic properties of the system (see SI).
We estimated the entropy change ∆S from the initial to the transition state by computing the difference in free energy barriers ∆G b between 0 and 1,500 K, where G at 0 K is: and ξ is a collective variable (the choice of collective variables is described in the Methods and SI); U is the potential energy and x are atomic coordinates.We calculated G at high temperature, specifically T = 1, 500 K, using first-principles molecular dynamics (MD) and the adaptive biasing force method (see Methods), and we estimated ∆S as ∆G b /T (see SI).
Due to the computational cost, we obtained ∆S for only three paths (see SI).
Once we obtained E b, EFF and ∆S, we could compute the temperature T a , above which a given process is thermally activated, and for which we used the harmonic transition state theory: where Γ 0 denotes an attempt frequency and Γ a jump frequency.The values chosen for Γ 0 , Γ and ∆S are given in the SI.A simple sensitivity analysis, also in the SI, shows that in Eq. 3 the prefactor is relatively insensitive to the choice of these values.In addition, we systematically investigated the effect of thermal expansion and that of entropy on computed activation temperatures, amounting to variations in T a of less than 10 % (see SI).

B. Theoretical predictions
We start by presenting our results for 3C-SiC and we report our predictions for the activation temperature T a for various processes.
Activation Temperature -In Fig. 3 we show T a as a function of E F , where lines indicate the values above which a given process can occur.We find that for all values of the Fermi level, the onset of V Si migration occurs at temperatures lower than those activating V C diffusion, consistent with the results of previous studies [16, 17, 21-23, 25, 26, 34, 35, 41, 42].Above 1,000 K, V Si can diffuse, and hence when migrating it may lead to the formation of VV.
Interestingly, our calculations show that the pairing of mono-vacancies is facilitated by the Coulomb interaction between V − Si and V +1/+2 C . Indeed, we find that for 1.46 are respectively the most stable charge states of the two mono-vacancies (whether V C is in charge state +1 or +2 depends again on the Fermi level).It is important to note that for E F < 1.85 eV, a simple consideration based on energy barriers E b would yield T a ∼ 1,500 (1,700) K as the temperature required for a carbon vacancy to migrate in a stable charge state q = +1 (+2).However, upon computing effective barriers, we find a process with lower T a (as low as ∼ 1,300 K); such process involves intermediate charge states that are not thermodynamically stable but nevertheless allows for paths with lower barriers.Specifically, we find that thermal vibrations and changes in carrier density at high charge states, and that the latter migrate through a path with a barrier lower than that of V

+1/+2 C
, before returning to the original charge state (see SI).
The migration of V Si discussed above is a necessary but not a sufficient condition for the formation of VV.We find, in agreement with previous studies [16,[21][22][23], that it is important to realize, at the same time, n-type conditions.In particular E F should be above 1.46 eV (see Fig. 3).Indeed, under p-type conditions (E F < 1.46 eV), the V Si → CAV process is energetically favored over monovacancies diffusion (region B in Fig. 3) and V Si is trapped into the CAV complex and becomes immobile [21,22].Once formed, CAV remains stable as both the back conversion to V Si (CAV → V Si ) and its dissociation (CAV → C Si + V C ) are unlikely below 1,500 K due to high free energy barriers.Instead, under n-type conditions (E F > 1.46 eV) the migration of V Si is an energetically favored process and VV creation may occur (region A in Fig. 3).We note that in general, the higher E F , the more favorable the conditions for VV formation for several reasons.Increasing E F leads to a lower MV migration barrier and increased mobility for V Si , and to a higher (reduced) barrier for , leading to a lower probability of CAV formation and a higher probability of CAV conversion to V Si .
An additional condition for the formation of VVs is an annealing T below 1,300 K.We find that VV can migrate for T > ∼ 1,300 K and will likely either form large complexes (e.g. VV + V C ) [43] or diffuse and eventually move to the surface of the sample (region C in Fig. 3).These processes undermine the stability and abundance of double vacancies.In addition we find that V C is immobile below 1,300 K, which is overall a favorable condition for VV formation.Indeed we expect V C to be abundant in experimental samples, due e.g. to a formation energy lower than that of V Si and other point defects, and it can be incorporated in a VV + V C cluster if it migrates.Therefore, we suggest that T Ann should be < 1,300 K for optimal yield, stability and localization of VV defects.
Note that larger vacancy clusters can be formed also by incorporating V Si [21] into VV, and this undesirable process can only be mitigated by reducing the concentration of V Si .
Unfortunately, once VV is formed, charge-state engineering [12,14,20] is not an effective tool to hinder the formation of larger vacancy clusters because the most stable state of VV is neutral for E F above mid-gap.
Other defects of potential concern for the stability and formation of the VV are single interstitials (C i and Si i ).For example, C i could be re-emitted from C i clusters at high T , and subsequently aggregate with VV.Fortunately, C i emission is unlikely to occur below 1,300 K, according to previous DFT calculations (barrier > ∼ 4 eV) [43][44][45].Nonetheless, VV could be annihilated by the presence of C i if some weakly bonded C i clusters turn out to be present in the sample.
To estimate the optimal annealing T in the range of (1,000, 1,300) K, we consider the dependence of the Fermi level on doping densities.In 3C-SiC, maintaining E F > 2 eV requires a rather large doping density > ∼ 10 18 cm −3 (see Fig. 4).Therefore, it is conceivable that a desirable Fermi level range is 1.46 < E F < 2 eV, over which T a for V Si migration is a constant, roughly equal to 1,200 K.For T Ann between 1,200-1,300 K the Fermi level would be lower, for a fixed doping density, than in the range 1,000-1,200 K, hence possibly leading to V Si trapping at CAV defects.Therefore we conclude that the optimal T Ann ∼ 1,200 K.
So far, we have identified a suitable range of T Ann under n-type conditions, (1,000, 1,300) K, with an optimal value of 1,200 K.However, as shown in Fig. 3 (region D), there exist conditions at which VV may form also below 1,000 K, as long as there are V-V defects present in the sample.Surprisingly, barriers are lower for the pairing of third than second neighbors' vacancies.These results may help understand the conditions required for VV formation in small SiC nanoparticles (with diameter less than 10 nm), observed at lower T Ann , e.g.∼ 440 K [46], than in the bulk, since in nanoparticles MV separation distances are usually smaller.
Fermi level and defect density-We now turn to exploring how conditions identified above for VV formation may be achieved experimentally, by controlling for example the Fermi level and density of defects.In addition to the electronic properties of the system, E F depends on T , initial sample doping and of course defect density (see SI) which, at each given time of the annealing process is the most elusive parameter.The spatial distribution of defect density may be non-uniform and it depends on specific dose and energy of particles used during the bombardment of the sample [1,4,12,38].In spite of these uncertanties, it is interesting to obtain a qualitative estimate of the doping conditions necessary to achieve the desirable Fermi level values for the formation of VVs.In Fig. 4, E F is calculated at several T and for various doping and defect densities.We find that the presence of V C and/or CAV would induce n-doping while the presence of V Si would induce p-doping in the sample.Hence, the required condition to reach E F > 1.46 eV within (1,000, 1,300) K, is that at least one of the following concentrations-n-doping (e.g.amounts (see Fig. 4).Further, the V Si to CAV conversion process, although it renders V Si less mobile, helps increasing E F which in turns facilitates VV formation.These results emphasize the complex, interdependent role of multiple defects in tuning E F and ultimately leading to the formation of VVs.
We conclude this section by discussing VV formation properties in hexagonal (hex) polytypes, e.g., 4H-SiC.The extension of our results for 3C-SiC (where only k-sites are present) to hexagonal lattices (where both h-and k-sites are present) should be considered as a qualitative prediction.In hexagonal samples, the variation in stability and barriers of defects occupying different lattice sites is negligible, compared to the energy scale of ∼ several eV of most barriers computed in our calculations [16,18,21,23,42].The position of the valence band maximum (VBM) is nearly the same in cubic and hexagonal SiC, while that of the conduction band minimum (CBM) is higher in hex-SiC [47] (see Fig. 3).We find that the creation of VVs is more facile in hex-SiC; indeed the conditions of regions A (corresponding to formation and stability of VV) and D (corresponding to pairing of nearby vacancies) can be obtained in a slightly wider range of temperatures than in 3C-SiC (for values of the Fermi level attainable in hex-SiC) and, importantly, for lower doping densities.For example, in 4H-SiC, under intrinsic condition, E F ∼ 1.6 eV is larger than 1.46 eV, and with a moderate n-doping > 10 15 cm −3 , it may be increased above 2.0 eV at 1,200 K (see SI).We predict that an appropriate T Ann for hex-SiC is in the range of (900, 1,300) K, with an optimal value around 1,200 K.
Our predictions are in excellent agreement with several experimental observations.To synthesize VV, most experiments adopted T Ann in a range of ( Further, T Ann should be lower than ∼ 1,300 K to suppress VV diffusion, thus ensuring its stability and immobilization, with the optimal T Ann estimated to be ∼ 1,200 K.However, VV can also be created at lower T from neighboring V C -V Si pairs; these may be present after irradiation or implantation, and may be prominent in SiC nanostructures, suggesting that the formation of VVs in small nanoparticles should occur at lower T than in the bulk.
Our findings also suggest that VV signals may be detected at low annealing temperatures, which however should not be interpreted as lower bounds for V Si diffusion.Moreover, we predict that VV formation in hex-SiC can be more facile than in 3C, due to a larger band gap and higher CBM position, which allow for the use of lower doping densities and lead to a slightly broader range of favorable annealing T .Our results are in excellent agreement with experiments, while providing new and improved understanding of formation mechanisms at the atomistic level.The knowledge obtained here may benefit the controlled fabrication and device integration of VV, assisting its applications for quantum technologies.
Importantly, the computational protocol and strategies developed here, based on first principles calculations, are general and can be readily extended to investigate defects in other covalently bonded materials.Multiple paths with different charge states should be considered to understand point defect formation processes, taking into account thermodynamically unstable ones, which may facilitate the exploration of low barrier paths at high T .Our findings show that it is key to conduct calculations of effective barriers as a function of the Fermi level, which itself depends on T , and not only of barriers between thermodynamically stable states.In addition, it is critical to consider not only formation but also annihilation pathways to obtain faithful predictions of formation processes.Unexpectedly, we found that although important for accurate quantitative predictions, thermal expansion and entropic contributions are not critical to determine general trends of activation temperatures for different paths.
One important problem that remains to be addressed is the influence, on defects' formation, of the specific synthesis procedures, e.g., by irradiation of the sample.Using our computed energy barriers as input, one can simulate real-time defect evolution, e.g., via kinetic Monte Carlo methods, which could then provide information about optimal annealing times.These possible directions are worthy of future explorations.

A. Density functional theory calculations
We performed DFT calculations using the Qbox [48] and the Quantum Espresso [49] codes.We used the PBE [50] and DDH [51] (15% exact exchange) functionals, optimized norm-conserving Vanderbilt pseudo-potentials [52], a plane-wave kinetic energy cutoff of 60 Ry.We conducted calculations in 216 atom supercells with lattice constant 4.416 Å, and with either the Γ point or a 2×2×2 Monkhorst-Pack grid to sample the Brillouin zone.The lattice constant was determined by first-principles MD (FPMD) simulations in the NPT ensemble at 1,500 K at the PBE level of theory.We considered structural relaxations as converged when residual forces on atoms were < 0.01 eV/ Å.We considered charge state q from -2 to 2 for all defects, expect for V C where q = 0, 1 or 2; spin state for even (odd) number of electrons, where S: singlet, D: doublet, T: triplet, Q: quartet.We chose not to employ empirical force-fields, which would have allowed for the use of larger supercells, as they are not appropriate to simulate q and s degrees of freedom; in addition we found that in several cases many of the popular force-fields used for SiC cannot reproduce DFT results.

B. Nudged elastic band calculations
We carried out climb image nudged elastic band (CI-NEB) simulations at the PBE level with 2×2×2 k-point grids, by coupling Qbox with the PASTA [53] code.We used spring constants of 2 eV/ Å2 and force tolerance of 0.02 eV/ Å.We determined the most stable spin state among [S, T] or [D, Q] for each NEB image at a given charge state q; the corresponding total energies and atomic forces were then used to update NEB images to determine the minimum energy path and energy barriers E b .In this way, E b is only a function of q.For most pathways studied here, the most stable spin state remain the same along the whole path; for those paths for which we observed a change of spin states, we found that the energy splitting between different spin states at the transition-state is generally small, i.e. less than 10 % of E b .
We then computed total energies for converged images, at the PBE and DDH level of theory using only the Γ point [18].We denote the barriers obtained in this way as

C. Formation energy calculations
The formation energy of defect X in charge state q, E f (X q ), was computed as: where E tot (X q ) is the total energy of a SiC supercell with X q ; E tot (SiC)is the total energy of the pristine SiC supercell; µ C and µ Si are chemical potential of C and Si; n C and n Si are number of added (+) or removed (−) C or Si atoms to form X, respectively; E F is the Fermi energy referred to the VBM; E corr (X q ) is the energy correction for spurious electrostatic interactions present in supercell calculations.
Using relaxed configurations at the PBE level of theory, we computed the total energy and electrostatic potential using Quantum Espresso and the DDH functional.We obtained E corr using the method developed by Freysoldt, Neugebauer, and Van de Walle [54].We used a dielectric constant equal to 9.72.The chemical potential µ C was calculated as the energy per atom in diamond; µ Si was calculated as µ SiC − µ C , where µ SiC is the energy per formula unit in bulk SiC.The results are shown in Fig. S1.
Finally, binding energies between defects are required to compute the barriers of the CAV/VV dissociation processes.We estimated the C Si and V C binding energies as ∼ 1 eV from previous studies [21,22].We directly computed the V C and V Si binding energies, which are ∼ 3 eV for E F near the mid-gap of 3C-SiC.

D. Enhanced sampling calculations
We computed free energies of defect transformations by coupling the Qbox and SSAGES [55] codes.We used Qbox to perform FPMD in the NVT ensemble and the adaptive biasing force method [56] in SSAGES to calculate free energy gradients.We utilized the collective variable (CV) ξ C/Si [22]: where R C/Si are the coordinates of moving C/Si atoms; m i is the mass of the i th gate atom; R i is the coordinate of the i th gate atom; M is the total mass of the gate atoms; e projection is the unit projection vector (see Fig. S2).
We carried out free energy calculations for the three processes presented in Fig. S2, where the definition of CVs and gate atoms is specified.For each path, we performed FPMD simulation at 1,500 K using a time step of 1 fs, for ∼ 370 ps.For computational efficiency, we used the PBE functional, 40 Ry kinetic energy cutoff and the Γ point; we considered defects only at q = 0 and s = T in our MD simulations.
To elucidate the effect of T on barriers, we also computed free energy profiles at 0 K.
We used the same functional and cutoff as in FPMD simulations for consistency; we used the Sequential Least SQuares Programming method [57] in the SciPy package to carry out constrained optimizations along one-dimensional CVs.   a Only the first step in VV migration path is simulated, as shown in Fig. S2C.
b F refers to forward direction from the initial to the transition state; B refers to backward direction from the final to the transition state.

NOTE 3: CALCULATION OF EFFECTIVE BARRIERS
During transformations occurring at high temperature (T ), point defects may be in several charge (q) and spin (s) states different from the thermodynamically stable ones.Hence we should consider different energy barriers E b (q, s).Moreover, the transition between different q and s states may occur at elevated T due, e.g., to vibrational effects.For these reasons, we used effective barriers E b, EFF [1,2] to describe atomic processes occurring at high T , instead of simple barriers E b .
In our NEB calculations, we first determined the most stable s state for each image at a given q; the corresponding total energies and atomic forces were then used to determine the minimum energy path and E b .Hence the final E b obtained in this way is only a function of q, with the effect of the s degree of freedom (DOF) included implicitly.Our treatment of the spin DOFs assumes an instantaneous equilibration of spin states during defects' transformations.Although we could not estimate the timescale of s transitions via spin-orbital-coupling or spin-phonon interactions at ∼ 1,000 K, we found that in general considering different spin states affects only slightly the computed E b (see Methods in the manuscript).
We computed E b, EFF based on E b and defect formation energies, considering only the charge DOF (see Computational strategy in the manuscript).We assumed the charge state q to be preserved during defect transformations, due to the short lifetime of barrier crossing over the transition state, and we considered q transitions at the initial and final states of a given path.
In the case of the dissociation of complex defects, involving multiple steps, E b, EFF was estimated from the binding energy and diffusion barriers.For instance, E b, EFF for the CAV dissociation process was obtained as the sum of the binding energy ( C Si & V C ) and the For E b, EFF to be accurate, the transition between different q states needs to be fast relative to the transformation of defects into different configurations.We estimated the q transition, via carrier capture or emission, is indeed fast at high T , as indicated below.
Under equilibrium condition, the q transition rate (g) [2,3] can be estimated as: where σ is the capture cross section; v is the average thermal velocity of carriers; N is the effective density of states (DOS); γ is the degeneracy factor; ∆E is the energy difference between defect levels and the closest band edge; k B is the Boltzmann constant.

NOTE 4: CALCULATIONS OF THE ACTIVATION TEMPERATURE
According to the harmonic transition state theory [10,11], the jump frequency Γ can be calculated as: where Γ 0 is the attempt frequency and G b the free energy barrier of a given process.
During defect transformations, we assumed the system volume to be constant, hence: where ∆U is the change in internal energy from the initial to the transition state of a given path; ∆S is the change in entropy from the initial to the transition state of the path.Note that we computed ∆S for three paths only, due to the computational cost, and found values varying within ∼ (0.85, 2.90) k B (see Note 2).Based on the harmonic approximation and classical statistics, the change in kinetic energy is 0 eV, and the change in potential energy is E b (barrier at 0 K), enforced by the equi-partition theorem.In addition, ∆S is constant, as determined by calculations of phonon frequencies [10,11].Therefore, we obtain: The activation temperature T a is defined as the T above which a process is thermally activated.Based on Eq. 4, T a can be written as: Similar to previous studies, we approximated Γ 0 to be 1.6 × 10 13 Hz [12]; jump frequency Γ to be 0.1 Hz [13,14]; ∆S to be 2.3 k B .This value corresponds to our estimate of the ∆S from the stable VV configuration to the transition state (see Note 2).We obtained a prefactor (inverse of the quantity within square brackets in Eq. 5 ) of 331 K/eV.A simple sensitivity analysis shows that such prefactor is relatively insensitive to the choice of Γ and ∆S.For example, by varying Γ in the range of (0.01, 1) Hz (with all other parameters fixed), the prefactor changes by < ∼ 7 %; by varying ∆S in the range of (0.85, 2.90) k B , the prefactor changes by < ∼ 4 %.
We systematically investigated the thermal expansion and entropic effects on computed energy barriers.We found that considering lattice expansion at 1,500 K leads to a minor change of E b (which is on the order of several eV) of approximately ∼ ± 0.1 eV for most processes, with the exception of small carbon-clusters' formation for which differences in energy barriers are ∼ ± 0.3 eV.We found that due to entropic effects, our computed free energy barriers are lowered by ∼ (0.11, 0.38) eV at 1,500 K, relative to those obtained at 0 K (see Note 2), consistent with estimates based on the harmonic approximation [12].As a result, we estimate that the variation of T a due to thermal expansion and entropic effects is less than 10 % .to the VV creation processes.We determined E F using the following equations [5,15,16]:

n-& p-doping
b) n 0 + N a = p 0 + N d + X∈(V C ,V Si ,CAV) q q × N (X q ) c) N (X q ) ∝ N X g X q exp(−E f (X q )/k B T ) where n 0 is the electron density; p 0 is the hole density; N C is the effective conduction band DOS; N V is the effective valence band DOS; E g is the band gap; N a is the acceptor density; N d is the donor density; X q stands for defect X in charge state q; N (X q ) is the density of X q ; N X is the total density of defect X; g is the degeneracy factor; taken as 1; E f is the defect formation energy (see Note 1); E V is the valence band maximum (VBM) energy.
Eq. 6-b expresses the charge neutrality condition, incorporating the effects of defects.
These equations need to be solved self-consistently.Here, we performed a line-search by step-wisely increasing E F from VBM to the conduction band minimum; we determined E F as the value, which makes Eq. 6-b satisfied with a minimal error.The electronic properties of SiC were obtained from the Appendix C of the book [4].For simplicity, we ignored the T dependence of these parameters: 1) E g at 300 K was used; 2) E f diagram at 0 K was used (Fig. S1).We deduced the DOS effective masses based on the measured N C and N V at

Fig. 1 .
Fig. 1.Atomic pathways investigated in this work.A: Monovacancy dynamics, including carbon (V C ) and silicon (V Si ) vacancy migration, and V Si and carbon antisite vacancy complex (CAV) inter-conversion.B: Pairing of second (V-V 2 ) and third (V-V 3 ) neighbors V C and V Si vacancies to form a double vacancy VV.Only V-V up to third neighbors were considered, due to the size of our supercells.C: VV migration path with the lowest barrier, where several steps are illustrated.V C C Si V C complex in 3 is denoted as VCV.

Fig. 2 .
Fig. 2. Computational Protocol.The calculations highlighted in blue (red) were carried out at zero (finite) temperature (T ).We specify the functionals used in the calculations (PBE and DDH, see text) and the computed quantities: defect formation energies (E f ), energy and effective energy barriers (E b and E b, EFF ), Gibbs free energy bariers (G b ), entropy differences (∆S) and activation temperature (T a ).We obtained E PBE b

Fig. 3 .
Fig. 3. Computed activation temperature T a as a function of the Fermi level E F (referred to the top of the valence band).The processes considered in this work are shown on the right hand side of the figure and they are summarized in Fig.1; the notation V-V 2/3 → VV refers to V-V 2/3 → VV @ V Si (Fig.1B), a path with lower barrier than V-V 2/3 → VV @ V C .The lines indicate the temperature, as a function of E F , above which the process indicated on the right hand side is activated.Regions A, B, C and D where specific processes occur are described in the text.The arrows indicate the conduction band minimum of the various polytypes, which was computed by aligning their respective valence band maxima; they are shown in increasing order of energy for 3C, 6H and 4H of SiC.

Fig. 4 .
Fig. 4. Fermi level (E F ) as a function of defects density.The Fermi level is referred to the top of the valence band.We show results for different temperatures (1,000 K (A), 1,200 K (B) and 1,400 K (C)); we consider, separately, initial n-doping of the sample, and specific concentrations of carbon (V C ) and silicon (V Si ) vacancies, and carbon antisite vacancy complexes (CAV); we also consider two additional cases: (i) same concentration of V C and V Si (V C = V Si ); (ii) same concentration of V C and CAV (V C = CAV).The dashed line indicates the value of mid-gap for 3C-SiC; the green-region for E F > 1.46 eV indicates favorable conditions for the formation of the double vacancy in the range of temperatures between 1,000 and 1,300 K (see text).
We computed the correction to apply to PBE results in order to estimate DDH barriers as: [E DDH b @ Γ − E PBE b @ Γ].We added such correction to E PBE b @ 222 (barriers computed with the 2×2×2 k-point grid) to obtain E DDH b @ 222.Here, we assumed that the minimum energy paths at the PBE and DDH level of theory are similar; energy difference calculated with the Γ point differed only slightly from those obtained with the 2×2×2 k-mesh.The results reported in the main text were obtained with E DDH b @ 222.

Fig. S1 .
Fig. S1.Formation energy of defects in 3C-SiC obtained under C-rich conditions, obtained using DFT and the DDH functional.The geometry of these defects can be found in Fig. 1.See Methods in the manuscript.

NOTE 5 :
CALCULATIONS OF THE FERMI LEVELAfter irradiation or implantation of SiC samples, multiple defects can be created including interstitials, antisites, substitutionals and vacancies.In order to obtain an accurate value of E F we need to consider both external doping and the charge state of the defects created in the sample.

Fig. S3 .
Fig. S3.Fermi level as a function of T in 3C-SiC.The Fermi level is referred to the top of the valence band.We consider the doping or defects density in the range of (10 14 , 10 19 ) cm −3 .A: presence of n-or p-doping only.B: presence of carbon vacancy (V C ) only.C: presence of silicon vacancy (V Si ) only.D: presence of carbon antisite vacancy (CAV) only.E: presence of V C and V Si of the same amount.F: presence of V C and CAV of the same amount.The value of mid-gap for 3C-SiC is indicated by a grey dashed line; the green-region for E F > 1.46 eV indicates the suitable conditions for the VV creation.

Fig. S4 .
Fig. S4.Fermi level as a function of T in 4H-SiC.The Fermi level is referred to the top of the valence band.We consider the n-or p-doping density in the range of (10 14 , 10 19 ) cm −3 .The value of mid-gap for 4H-SiC is indicated by a grey dashed line; the green-region for E F > 1.46 eV indicates the suitable conditions for the VV creation.
F .We then identified favorable conditions for the formation of VVs and discussed how suitable values of E F can be obtained via careful tuning of doping or defects densities.Our calculations show that one should use n-doped samples with E F > 1.46 eV during annealing, to ensure the stability of single V Si , and T Ann > ∼ 1,000 K to activate V Si migration for aggregation with V C .
Ann 17] determined by PL or EPR maximum intensities and found to be ∼ 1,150 K, in agreement with our calculations of ∼ 1,200 K.We emphasize that, depending on the experimental setup, the decrease in signal above ∼ 1,150 K is not necessarily related only to changes in VV concentration.We predict VV can be stable up to 1,300 K, above which its density decreases due to diffusion.This is consistent with the significant drop of VV signals in experiments as T > 1,300 K[15], and with the highest PL and EPR intensities detected at 1,300 K[16,17].III.DISCUSSIONBy combining DFT calculations with semilocal and hybrid functionals, nudged elastic band and first principles MD simulations, we obtained a detailed, atomistic description of the VV formation process in 3C-SiC.We computed energy barriers and activation temperatures for multiple defects and pathways as a function of the Fermi level E