Exhaustive characterization of modified Si vacancies in 4H-SiC

2.1 ADAQ software and data

To narrow down the number of possible defects that fit the observed experimental data, we turn to the high-throughput data [31–33] produced by ADAQ [34] that is implemented with the high-throughput toolkit httk [35]. ADAQ is a collection of automatic workflows designed to speed up the search for point defects. It generates defects and calculates the most important properties such as total energy and ZPL for one excitation for different charge and spin states in a screening workflow. For detailed description of ADAQ, see Ref. [34]. In brief, ADAQ runs density functional theory [36, 37] (DFT) calculations using the Vienna Ab initio Simulation Package (VASP) [38, 39] (v.5.4.4) with the semi-local exchange-correlation functional of Perdew, Burke, and Erzenerhof (PBE) [40]. Due to the many VASP invocation for the different charge, spin, and excitation for point defects; ADAQ is necessary to handle the vast amount of computations.

Previously, ADAQ was employed on 4H-SiC and screened 8355 single and double intrinsic defects in 4H-SiC [32]. The detailed report on the results collected in a database will be presented elsewhere [31, 33]. The defects were generated with the settings for ADAQ to include double defects with a maximum distance of 3.5 Å. These settings roughly correspond to point defect clusters with second nearest neighbors. To keep track of the most stable defects, the concept of the defect hull is introduced [32] which consists of the point defects with the lowest formation energy for a given stoichiometry. The defect hull is analogous to the convex hull of stability used to discuss the thermodynamical stability of bulk materials.

2.2 Search for silicon-vacancy-like signals

Here, we present two different ways of searching through the defect database generated by ADAQ for a defect that explains the experimental measurements. First, the EPR measurements show that the silicon-vacancy-like defect has spin-3/2. By searching through the ground state results for the 8355 single and double intrinsic defects in 4H-SiC, 39 defect configurations with spin-3/2 are found. Of these configurations, 24 contain a silicon vacancy. These can further be reduced by only including defects that are a maximum of 1 eV above the defect hull and have a positive binding energy. The remaining eight final entries consist of the 2 isolated silicon vacancy configurations and 6 configurations of a cluster containing a silicon vacancy with a carbon antisite at the second nearest neighbor. Given these 4 search criteria that (i) the defect spin is limited to S = 3/2, (ii) the defect included a silicon vacancy, (iii) the defect has a positive binding energy, and (vi) the defect is a maximum of 1 eV energy above the defect hull. The only point defects that fit all of them are the silicon vacancy and a cluster consisting of a silicon vacancy with a carbon antisite at the second nearest neighbor. Hence, we will refer to any combination of silicon vacancy and carbon antisite as modified vacancy in the rest of the paper. The modified vacancy is on the defect hull for the stoichiometry of two missing silicons and one extra carbon. Since ADAQ is not limited to only defects along the c-axis, it allowed us to find additional configurations compared to previous searches [30].

Second, the PL results show similar ZPLs in regions next to the silicon vacancies which are at 1.352 and 1.438 eV [7]. Note that due to the use of the PBE functional, the ZPLs are shifted down by 0.2 eV which is discussed in detail in Ref. [34]. Since all defects are calculated with the same level of theory, the search criteria for ZPL are larger than 1 eV and smaller than 1.3 eV. Combined with a maximum of 1 eV above the defect hull, these criteria give 9 final entries. Here, 2 divacancy configurations emerge due to the wide ZPL search range. Disregarding them, we are left with 7 configurations, 1 silicon vacancy and 6 modified vacancies. A ZPL is missing for a silicon vacancy due to the settings in the screening workflow which is also discussed in detail in Ref. [34]. The estimated ZPLs for the modified vacancies are 1.09, 1.13, 1.20, 1.22, 1.24, and 1.25 eV.

Either way, when one searches for a defect with silicon-vacancy-like properties with results from EPR or PL, the modified vacancies are found in both cases. With these search criteria, we exhaustively examine all defects in the database. The modified vacancy is the only defect of the considered single and double intrinsic defects for 4H-SiC that can fit the experimental observations and is worth additional study. It is important to underline that ADAQ workflows are implemented with the accuracy settings needed for high-throughput calculations. This means significantly more accurate calculations are needed to characterize the modified vacancy properly. Below we present these calculations and results.

4.1 Theoretical results

Figure 2a shows the formation energy of the six closest configurations of the modified silicon vacancies obtained by the HSE functional. The formation energy trend for modified vacancies resembles silicon vacancy formation energy since the carbon antisite is neutral across all Fermi energies. The only difference is that the modified vacancy has a stable positive charge state in contrast to the isolated silicon vacancy. In the negative charge state, the binding energy for these defects is around 0.3 eV. The positive binding energy and the fact that the modified vacancy is on the defect hull make it a stable defect. Figure 2b shows the binding energy for all the considered negatively charged configurations in the 576 and 2304 atom supercells calculated with the PBE functional. Two different supercells are used to ensure that the defect self-interaction is low for the configurations with large separation. The perturbative effect of carbon antisite decreases as the antisite is placed farther away from the silicon vacancy, and the binding energy approaches zero beyond a separation of 10 Å. However, there are peaks in the binding energy at periodic intervals of about 3 Å (see the added guideline in Figure 2b). The peaks appear when the carbon antisite is placed at a silicon site with a large spin density overlap, which is plotted in Figure 2c. For an isolated silicon vacancy, the spin density is localized on carbon sites (Figure 1b), and the placement of carbon antisite adjacent to spin density results in attractive interaction. Two configurations also have negative binding energy, most notably k 5.0 b .

Figure 2:

Energetics and the site dependent behavior of modified silicon vacancies. (a) The formation energy (HSE results) for the closest modified vacancy configurations (colored lines) compared to the isolated silicon vacancy and carbon antisite (black and grey lines); (b) the binding energy as a function of distance between the silicon vacancy and the carbon antisite. Notice the guideline to show the oscillating behavior of the binding energy. Here, the PBE results for the two different supercell sizes are shown. (c) The spin density overlap (PBE results) for the carbon antisite on the isolated silicon vacancy defect states. Sites with large spin density overlap in (c) correspond to the local maxima in binding energy in (b).

Next, we discuss the optical properties of modified vacancies. Figure 3a shows the ZPL and TDM for the modified silicon vacancies compared with the isolated silicon vacancies. Here, one can see that the ZPL for the six closest modified silicon vacancies decreases by ≈75 meV, and the TDM increases by ≈1 Debye compared to the isolated counterparts. As the antisite moves further away, the ZPL and TDM get closer to the isolated silicon vacancy, as shown in Figure 3b and c.

Figure 3:

The optical data for the modified silicon vacancies. (a) The ZPL (both PBE and HSE results, the PBE results are shifted so that the isolated vacancies are align with the HSE results) and TDM for all the configurations of the modified silicon vacancies where the isolated silicon vacancies are marked with vertical dashed lines. The six closest configurations are marked with red arrows. (b) And (c) zoomed in versions close to the k and h isolated silicon vacancy, respectively.

Figure 4 shows the eigenvalues for the modified vacancies. As the carbon antisite is moved farther away from the silicon vacancy, the eigenvalues approach the isolated case. The isolated silicon vacancy has three unoccupied states in one spin channel (marked in blue). The lowest is a single degenerate a₁ state, and above it, there is a double degenerate e state. For the six closest configurations, the eigenvalues show the largest displacement with a prominent upward shift for the occupied states. The unoccupied states have the largest splitting between the previous degenerate e states. Except for the below configurations ( h 3.1 b and k 3.1 b ) that barely split the degenerate e states but reorders the a₁ and e states compared to the isolated case. This outcome explains why the below configurations have a much lower ZPL change, see Figure 3a, compared to the above and planar configurations at the same distance. At around 6 Å, the eigenvalues are indistinguishable from the isolated case. This agrees with similar trends for ZPL and binding energy, thus highlighting the extent of local site symmetry away from Si vacancy. This eigenvalue shift is similar to the shift due to compressive strain of the isolated silicon vacancy [27].

Figure 4:

The Kohn–Sham eigenvalues (HSE results) for the modified vacancy for the (h) and (k) configurations compared with the isolated silicon defects. On the x-axis are the different configurations, where the subscript is the distance between the carbon antisite and the silicon vacancy. The dashed lines denote the eigenvalues for the occupied defect states of isolated silicon vacancy for both h and k configurations. The red and blue filled dots represent occupied states in each spin channel. Unfilled dots represent unoccupied states.

Figure 5 shows the ZFS for the modified vacancy compared to the isolated silicon vacancy. The six closest configurations have much larger E and D values than the isolated case. However, as the distance between the defects increases, the values approach the isolated silicon values due to the exponential decay of the overlapping spin densities. This change corresponds well with the difference in the eigenvalues as the distance increases.

Figure 5:

ZFS for the modified silicon vacancies (HSE results). (a) The D value and (b) the E value for the six closest configurations of the modified silicon vacancies compared with the isolated silicon vacancy; (c) the D-tensor and (d) the E-tensor for the remaining configurations. The data is grouped by distinct colors to easily identify the carbon antisite distance from the silicon vacancy.

4.2 Experimental PL measurements

In this section, we present results from three high-purity semi-insulating (HPSI) 4H-SiC samples, which all exhibit strong PL from silicon vacancy. See Section. VIIB for experimental details. The PL spectra are displayed in Figure 6 which shows several additional lines in a ≈20 nm range in the vicinity of each the isolated the h and k silicon vacancies. We have counted up to 63 additional lines which are listed in Tables 3 and 4 with the polarization given for the most prominent lines.

Figure 6:

PL for the high-purity semi-insulating 4H-SiC samples. (a) Show the range of 860–885 nm containing the V1 and V1′ signals and (b) show the range of 910–935 nm containing the V2 signal. The figures show there are several small additional lines next to the V1 and V2 lines. These additional lines appear in as-irradiated samples and only some faint signals remain after annealing at 800 °C. The polarization is also shown for the 400 °C annealed sample.

5.1 General trends

The magneto-optical properties of the modified vacancy depend on the position of the carbon antisite. As the carbon antisite moves farther away from the silicon vacancy, one would expect the ZPL to approach the isolated silicon vacancy. Surprisingly, at around 6 Å, the ZPL shift is as significant as for the closest configurations, see Figure 3a. This trend is present for both supercell sizes and functional used, and is explained by the overlap of the silicon vacancy spin densities with the position of the antisites. At 3.1 and 6.2 Å away from the vacancy, the overlap is large. When carbon atoms are replaced on these silicon sites, they significantly affect the total energy, Kohn–Sham eigenvalues, ZFS, and ZPL. For the ZFS, the closest configurations have order of magnitude larger D values compared to the isolated silicon vacancy. However, the low symmetry makes it hard to detect.

As the carbon antisite moves farther away, the defect orbitals around the silicon vacancy are modified less, and the values resemble the isolated silicon vacancy. We argue that even if the global symmetry is low (which makes it difficult to measure EPR signal [30]), the local site symmetry around the silicon vacancy still remains high (C_3v) for certain configurations. One can find an example of this in Tables 1 and 2, the 3rd nearest neighbor site k 4.4 b has C_1h symmetry but a low E value indicates that the electronic structure is close to C_3v symmetry, see Figure 4. Thus, the presence of electronically benign defects such as antisites near spin defects should be studied with regards to the orientation of spin density.

Of the tested modified vacancy configurations, two have C_3v symmetry within 10 Å. The k 5.0 a ″ configuration has the smallest D value (−13.34, see Table 2 or Figure 5c) of all the tested configurations and could be the defect responsible for the R₁ EPR signal [16]. The same configuration (named V S i − -C_Si(1)-k in Ref. [30]) also has the lowest D value. We draw the same conclusion as the authors in Ref. [30], the k 5.0 a ″ configuration is responsible for R₁ EPR signal. With the additional configurations considered in this paper, it is clear that no other configuration has a smaller ZFS value than k 5.0 a ″ . The other modified vacancy with C_3v symmetry, k 5.0 b ″ , was not considered in Ref. [30]. It has a large ZFS which could be related to T _V1b . However, it has negative binding energy. Several other candidates could fit this measurement; like k 5.4 a , h 5.4 p II, and h 6.2 a ″ II for example; but they have larger E values. Further experiments are needed before any conclusion can be drawn. The C_3v candidates for h vacancy are ≈10 Å away and are unlikely to have distinct EPR signals.

The modified vacancies are created under irradiation conditions, and it is crucial to understand the behavior of these defect complexes at high temperatures. Importantly, at what temperature would the modified vacancies anneal out? Even if their binding energies are small, the defects must be mobile to separate. The energy barriers for the different single defects depend on the charge state of the defect. But the barrier ranges are around 1–2 eV for the interstitials, 2–5 eV for the carbon vacancy, 3–4 eV for the silicon vacancy, and ≈10 eV for both antisites in 4H-SiC [41]. The large energy barrier makes the antisites stable regardless of the annealing temperature. The silicon vacancies would be mobile at high temperatures and either encounter a silicon antisite and transform to a carbon vacancy or bind to a carbon antisite creating the modified silicon vacancies discussed in this paper. However, the carbon vacancies become mobile before the silicon vacancies due to the lower energy barrier. They could interact with either the silicon vacancies to form divacancies or the carbon antisites to form silicon vacancies. Both these cases would affect the modified vacancies. The carbon vacancies become mobile at around 800 K (found with transition state theory [42, 43]). In Figure 6, the additional lines disappear between 400 and 800 °C annealing. Hence, we put forward that the annealing behavior of the additional lines is due to the mobile carbon vacancies that interact with the modified vacancies. The carbon vacancy either interacts with the silicon vacancy which transforms the modified vacancies into modified divacancies [44] or interacts with the carbon antisite, turning the modified vacancy into an isolated silicon vacancy.

Are there other defects that can explain the experimental measurements? If these signals correspond to an unknown point defect containing a silicon vacancy, the modified vacancy with a carbon antisite is the only valid candidate considering double defects at close distances. The number of spin-3/2 point defect clusters with silicon vacancies in the database of 8355 defects is low. This extensive search shows that the modified vacancy is the most likely defect. It has stable stoichiometry. At larger distances, one could argue that any perturbing defect could give rise to the results presented in this paper. However, there are a few factors that limit such scenarios. Let us first consider an antisite that modifies the silicon vacancies. Even though the carbon and silicon antisite have similar formation energies (carbon is slightly lower), the clusters formed by a silicon vacancy and silicon antisite have the same stoichiometry as a carbon vacancy and most likely will anneal to that defect. The modified vacancy has a stoichiometry that cannot transform into any other single or double defect. Similar reasoning can be applied for the interstitials: a carbon interstitial and silicon vacancy would transform into a carbon antisite, and a silicon interstitial and silicon vacancy (a Frenkel pair) would cancel each other out. In the case of silicon vacancy modified with a carbon vacancy, it would lead to the formation of a divacancy. A cluster with two silicon vacancies does not bind together [33]. This only leaves the modified vacancy with a carbon antisite as a probable candidate of the considered defects.

5.2 Identification

Here, we compare the theoretical and experimental data on magneto-optical properties of the investigated defect complexes to look for matches. Due to their larger TDM, the modified vacancies would show stronger PL signals than the isolated silicon vacancies if the concentrations were the same. However, since the modified vacancy is a double defect consisting of a carbon antisite and silicon vacancy, the concentration of modified vacancies will expectedly be lower than isolated silicon vacancy in irradiated and annealed samples. Hence, we expect the modified vacancies to appear with lower intensity than the silicon vacancies, which is the case for the new lines in Figure 6. The concentration differences do not affect the identification.

The identification is mainly based on the PL values and supported with EPR values when available. For the PL lines without ZFS, the assignment is done by comparing the ZPL values and polarization. The theoretical HSE ZPL values are shifted to align the isolated silicon vacancies and remove systematic errors. The relative errors between configurations are in the meV range [58]. Then the least squared distance (lsd) between theoretical ZPL (eV), para (fraction), perp (fraction) and the experimental values is found with the following equation:

When one theoretical configuration matches with multiple experimental values, the smallest lsd value is used to assign the line. The assigned configurations are summarized in Table 5, and also marked in Tables 3 and 4.

Many of the predicted ZPLs from the modified vacancies fall close to the isolated value. This makes it difficult to exclude other effects like strain [27], L-lines [24], or surface effects. Hence, we first focus on the lines farthest away from the silicon vacancy, at least 15 meV from the isolated line. The primary candidates for the R₁ EPR signal, k 5.0 a ″ has ZPL around 40 meV from the isolated h silicon vacancy with a parallel polarization, see Table 2 and Figure 5a. This is in excellent agreement with the 880.0 nm line observed in PL experiments with an identical polarization (Figure 6 a). Furthermore, this line is also seen in the same sample used to measure the R₁ signal [16]. Hence, the k 5.0 a ″ configuration is identified as the source of the R₁ signal and 880.0 nm line.

Next, we consider the five lines with perpendicular polarization (863.6, 864.8, 874.9, 875.9, and 882.5 nm) in Figure 6a. Table 5 shows that the h 3.1 b configuration is the ideal candidate for the 882.5 nm line (also suggested for 874.9 and 875.9 nm). This complete change of polarization is due to the reordering of the a₁ and e states compared to the isolated case, mentioned in Section 4. Furthermore, the 863.6 nm line is best matched to the h 5.4 a configuration and 864.8 nm to h 9.3 p . The latter identification uses theoretical PBE ZPL data from the 2304 atoms to avoid nonphysical defect-defect interactions of the smaller supercell. Both the h 3.1 b and h 9.3 p configurations lie on local maxima in the binding energy, see Figure 2b.

There are also several additional lines around the isolated k silicon vacancy, see Figure 6b. In contrast to the h configurations, the PL signals are observed on both sides of the isolated vacancy line V2. Again, focusing at least 15 meV away from the isolated line, the 929.8 nm line is best described by the k 8.9 b ′ ′ ′ configuration (also suggested for 926.1 nm). Closer to the isolated vacancy, there are two bright peaks on either side: 910.7 and 921.7 nm, which both are missing good candidates. Since the peaks disappear upon annealing at high temperatures, they may correspond to structural defects. Coincidentally, both peaks are almost equally spaced from the V2 line by 8.1 meV. Furthermore, the 915.3 nm line is best described by the k 4.4 a configuration. The k 7.8 b ′ ′ ′ configuration could be any of the 910.7, 911.2, 913.0, 913.4, or 914.6 nm lines since they are close and have the same polarization. It is tentatively assigned to the 913.4 nm line, which has the smallest lsd.

Within the 15 meV range next to the isolated h silicon vacancy, we observe a set of lines in the same region as L-lines [24], see Figure 6a. However, there is no one-to-one correspondence due to disagreement in the observed intensities and energy spacing from earlier study [24]. The L-lines are tentatively suggested to be vibronic replicas of the silicon vacancy [24]. However, the L-lines disappear upon annealing at 800 K, whereas the silicon vacancy remains [24] which contradict this hypothesis. Alternatively, these lines could originate from strain, surface effects, or another defect. Several possible modified vacancy candidates from Table 1 match the experimental spectra. By comparing the HSE results of the modified vacancies (Table 1) with experiment (Table 3), the best matches for both ZPL energy and polarization are: 862.2 nm is h 6.9 a ″ and 863.2 nm is h 7.8 b ′ ′ ′ . The 864.1 nm is matched to k 6.2 p with excellent polarization but the ZPL match is terrible, which is reflected in the high lsd. This result is a product of the automatic matching. h 5.4 p I is either 866.6, 869.6 nm (assigned), or 873.4 nm, and h 5.4 b ″ could be either of 864.3 nm (assigned), 865.6, or 865.9 nm. The 862.7 and 862.9 nm are too similar to the isolated silicon vacancy for a definite assignment. Note that we did not find any good matches between theory and experiment for the 867.1 and 877.7 nm lines with mixed polarization. Further studies of these lines are needed.

5.3 Outlook

For the modified vacancies with the most distant carbon antisites, the ZPLs are close to the isolated silicon vacancy. Here, it is difficult to identify the individual lines. Hence, one can consider the linewidth broadening in addition to the ZPL position. There is a difference in linewidth of the V2 centers depending if the implantation is done with He⁺ or Si²⁺, 0.3 nm compared with 1 nm [10]. However, it is difficult to know which effect (surface effect, strain, or other defects) contributes toward broadening. When comparing these experimental results to the molecular dynamics simulations done in Ref. [13], one sees that He goes deeper than Si (cf. Figure 4 in Ref. [13]). If one disregards surface effects and strain, more modified vacancies could be created when the implantation is carried out with silicon ions compared to helium. Many modified vacancy configurations (see Figure 3b) may contribute to increased linewidth. However, they should anneal out at around 800 °C, and one would expect to see decreased linewidth after such annealing. Such a decrease in linewidth is not clearly observed. A distant carbon antisite could account for the defects seen in Ref. [18]. One can see carbon antisites in the molecular dynamics simulations done in Ref. [13]. Even after annealing, there are carbon antisites left (cf. Figure 8 in Ref. [13]). A possible way to verify if the modified vacancies contribute to the broadening of the isolated vacancy ZPL is to map the positions of the carbon-13 around defects with, for example, the method presented in Ref. [45]. Using this method for an enriched sample, one could look for carbon-13 antisites that could affect the silicon vacancy.

Like modifying silicon vacancy with a carbon antisite, one can imagine modifying other defects, such as divacancy. When the modified silicon vacancy anneals out, one possibility is that they transform to modified divacancy by interacting with a mobile carbon vacancy. There are already additional lines (PL5-6) next to the divacancy which have been assigned to stacking fault [46]. If there are more lines, modified divacancy is a good starting guess.

The silicon vacancy has been suggested as a qubit candidate. It has unique properties compared to the NV center in diamond [2], for example, the Si vacancy has a much smaller ZFS, in the MHz range. Due to this low ZFS, an external magnetic field is required to split the ground state [3]. The closest modified vacancies have the same properties as the isolated silicon vacancy and should also be suitable qubit candidates. In addition, they have an order of magnitude larger ZFS that allows operation at lower external magnetic fields. The lower symmetry of the modified vacancies does not have to be seen as a drawback. The low symmetry configuration of the divacancy has shown promise for quantum technologies [47, 48].

7.1 Computational details

The computations were carried out with Vienna Ab initio Simulation Package (VASP) [38, 39] (v. 5.4.4), which, in turn, uses a plane wave basis set and the projector augmented wave (PAW) [49, 50] method. For the manual calculations, we employ the semi-local exchange-correlation functional of Perdew, Burke, and Erzenerhof (PBE) [40] and the screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE06) [51, 52] with the mixing parameter α set to the standard value (25%). For the PBE calculations, the plane wave energy and kinetic energy cutoff are set to 600 and 900 eV, respectively, which are the settings used in ADAQ. Whereas, for the HSE calculations, these are reduced to 420 and 840 eV, respectively. The total energy criterion is set to 10⁻⁶ eV for PBE, and 10⁻⁴ eV for HSE. The structural minimization criterion is set to 5 × 10⁻⁵ eV for PBE, and 10⁻² eV/Å for HSE. The pseudopotentials for C is labeled 05jan2001 and for Si is labeled 08 april 2002. Two supercell sizes have been employed with 576 (6 × 6 × 2) and 2304 (12 × 12 × 2) atoms, respectively. In this paper, only Γ-point sampling of the Brillouin zone with Gaussian smearing is used, and Ψ k = Ψ − k * is the only symmetry used.

The following defect properties [34] are calculated in this paper:

1. Formation energy [53] for the charged defect:

   (2) Δ H D q ( E f , μ ) = [ E D , q − E H ] − ∑ i n i μ i + q E f + E corr ( q ) .

   where E_D,q and E _H are the total energies of the charged defect supercell and neutral host supercell, μ _i are the chemical potentials (only the rich limits are used: stoichiometric condition), n _i are the added or missing elements, q is the charge, and E _f is the Fermi energy. For the charge correction term E_corr(q), the Lany-Zunger correction is used [54].

2. The binding energy H ^b is defined as:

   (3) H b = Δ H D , q ( E f , μ ) [ A ] + Δ H D , q ( E f , μ ) [ B ] − Δ H D , q ( E f , μ ) [ A B ] .

   here A and B are the separate components and AB is the cluster. Defects with a positive binding energy are considered to be stable. For the modified vacancy, this equation becomes:

   (4) H b = Δ H D , − ( E f , μ ) [ V S i ] + Δ H D , 0 ( E f , μ ) [ C S i ] − Δ H D , − ( E f , μ ) [ V S i − C S i ] .

   for a given modified vacancy (AB), the reference values for components A and B are taken from the corresponding isolated defects. For example, a modified vacancy with a silicon vacancy at h and carbon antisite at k uses the corresponding energies for these defects in the isolated case.

3. ZPL is calculated as the total energy difference between the ground and excited state [55]. The excited state is found by using constrained DFT (Δ-SCF method) [56, 57]. The relative ZPLs can be compared within an accuracy of 100 meV [58]. The lowest energy optical excitation is from the occupied a₁ state to the unoccupied a₁ state in one spin channel [59].

4. TDM between the ground and excited state is calculated using the wave functions from the relaxed ground and excited state, which provide an accurate polarization [60].

5. ZFS is approximated using the dipole–dipole interaction of the spins [59, 61] and is calculated with the implementation in VASP. From the D-tensor, the D value is expressed as D = 3D _zz /2 and the E splitting as E = (D _yy − D _xx )/2

7.2 Samples and experimental details

All the samples are cut from the same wafer and irradiated to 17 cm-2 with electrons of energy 2 MeV. The samples are considered to have negligible strain due to the good PL agreement with the isolated silicon vacancy and divacancy. Two of the samples were annealed at different temperatures for half an hour from 200 to 800 °C for investigating the behavior in the range of the silicon vacancy zero phonon lines at different temperatures (Due to insignificant differences only 400 and 800 °C samples are shown in Figure 6). All the PL measurements were performed with a Jobin Yvon HR460 monochromator. The monochromator is equipped with 1200 grooves/mm grating blazed at 750 nm and an InGaAs multichannel detector, ensuring resolution ≈1 Å and optimum sensitivity in the silicon-vacancy range (wavelengths ≈858–980 nm). All the samples are mounted in a variable temperature closed-cycle cryostat and cooled down to 3.5 K. We use a 730 nm diode laser for exciting V1 and 890 nm from tuneable Ti-sapphire laser as an excitation laser for V2 to eliminate the background luminescence from V1.