Exhaustive characterization of modified Si vacancies in 4H-SiC

The negatively charged silicon vacancy ($\mathrm{V_{Si}^-}$) in silicon carbide is a well-studied point defect for quantum applications. At the same time, a closer inspection of ensemble photoluminescence and electron paramagnetic resonance measurements reveals an abundance of related but so far unidentified signals. In this study, we search for defects in 4H-SiC that explain the above magneto-optical signals in a defect database generated by Automatic Defect Analysis and Qualification (ADAQ) workflows. This search reveals only one class of atomic structures that exhibit silicon-vacancy-like properties in the data: a carbon antisite ($\mathrm{C_{Si}}$) within sub-nanometer distances from the silicon vacancy only slightly alters the latter without affecting the charge or spin state. Such a perturbation is energetically bound. We consider the formation of $\mathrm{V_{Si}^-+C_{Si}}$ up to 2 nm distance and report their zero phonon lines and zero field splitting values. In addition, we perform high-resolution photoluminescence experiments in the silicon vacancy region and find an abundance of lines. Comparing our computational and experimental results, several configurations show great agreement. Our work demonstrates the effectiveness of a database with high-throughput results in the search for defects in quantum applications.


I. INTRODUCTION
Silicon carbide (SiC) has great potential for quantum information and sensing technologies due to the large-scale high-quality manufacturing and ease of integration into existing semiconductor devices. The silicon vacancy in SiC [1] is a well-studied defect with applications such as qubits, sensors, and single photon emitters [2][3][4][5][6][7]. Several of these applications work at roomtemperature [4,6]. Many efforts are dedicated toward the controlled fabrication of such defects with various techniques. Current endeavors focus on the fabrication of silicon vacancy arrays with both laser writing [8] and ion implantation [9][10][11]. To guide the ion implantation, several molecular dynamics studies of implanting hydrogen, helium, and silicon in SiC have been carried out [12,13]. All these efforts show that the silicon vacancy in SiC is a viable candidate for various quantum information and sensing technologies.
The negative charge state of the silicon vacancy has C 3v symmetry, a ground state with spin-3/2, and a rich many-body structure [14]. In 4H-SiC, there are only two non-equivalent positions for the silicon vacancy, denoted h for hexagonal-like layer and k for cubic-like layer. However, electron paramagnetic resonance (EPR) experiments show at least six different signals related to a spin-3/2 defect. Two of these signals (T V 1a and T V 2a ) have been accredited to the h and k configurations, respectively [7]. The remaining four additional related signals (T V 1b and T V 2b ) [15,16] and (R 1 and R 2 ) [16] remain unidentified.
In addition to EPR, several photoluminescence (PL) measurements have reported siliconvacancy-like signals. Five of these signals have been measured within a 19 MHz range, and there are additional small peaks whose origins are unknown [17]. Using high-resolution photoluminescence excitation, nine additional signals have been measured in 0.15 meV range without spectral diffusion [18]. Banks et al. state that the strain and electric field perturbation are low, and these signals most likely correspond to other defects. This is echoed by Ramsey et al., who state that the most likely source of such signals is nearby defects perturbing the silicon vacancy [19]. Furthermore, in an array of silicon vacancies created with ion implantation, 30% of the measured spots showed even larger spectral drift, as large as 3 nm (5 meV), towards smaller energies [9]. In addition, PL signals more than 30 meV from the silicon vacancy, have been seen in other experiments [15,[20][21][22][23]. The most prominent line is around 1.412 eV (878 nm). This line is also found in connection to measurements done for the L-lines-additional lines that span a 15 meV range next to the h silicon vacancy and are tentatively suggested to be vibronic replicas of the silicon vacancy [24].
Before assigning the above signals to a novel defect, it is necessary to consider the role of other perturbative effects such as thermal vibration, strain, and surface termination. Thermal vibrations can be eliminated as a potential source since many reported measurements are performed below 4.9 K. Similarly, the signals cannot be attributed to phonon replicas since they appear at higher temperatures (15 K) at a larger energy difference (37 meV) from the zero phonon line (ZPL) [25].
They could emerge from strain, where a large shift in ZPL (26 meV) has been reported for 6H-SiC nanoparticles [26], which corresponded to a 2.3% basal strain [27]. This extreme shift observed in nanoparticles is probably not reasonable for bulk. Udvarhelyi et al. found that the shift of the silicon vacancy is larger for axial strain than for basal strain and concluded that a shift of several meV would be possible in bulk. It is also possible that the surface effect may shift the ZPL.
However, all these effects are small. Some of the PL shifts appear to be quite large, beyond the range of 30 meV. Thus, these signals most likely correspond to an unknown defect.
Previous efforts to explore specific defect realizations that explain the EPR measurements have tested modifying the silicon vacancy with carbon vacancies along the c-axis [28,29]. This approach was extended to include all possible vacancies and antisites along the c-axis [30]. Moreover, carbon antisites as second nearest neighbors slightly off the c-axis have also been considered.
As reported, in Ref. 30, one of the tested configurations could be responsible for the R 1 EPR signal. However, no experimental agreement with PL has been found. The above manually tested point defects models are limited in scope since they do not include impurities or defect clusters containing interstitials. No thorough large-scale search for point defects that can explain the experimental observations has been conducted yet.
In this paper, we exhaustively show how a silicon vacancy modified by a carbon antisite is the only candidate among thousands of considered defects and characterize its different configurations. Sec. II outlines how we search for point defects that explain the observations related to the silicon vacancy in data produced in high-throughput calculations [31][32][33]  To narrow down the number of possible defects that fit the observed experimental data, we turn to the high-throughput data [31][32][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.

B. 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 8 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 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.

III. MODIFIED VACANCY-SILICON VACANCY MODIFIED BY CARBON ANTISITE
Around each silicon vacancy configuration (h and k) in 4H-SiC, there are three non-equivalent second nearest neighbor sites where a carbon antisite can be placed. Hence, six different configurations exist for the closest modified vacancies, which are shown in Figure 1a.  ADAQ generated the closest configurations. In addition, we also manually place the carbon antisite farther away from the silicon vacancy to study how the distance and orientation between the two defects affect the properties. Figure 1b and c show these configurations. Here, the local site symmetry around the silicon vacancy is used to find the non-equivalent modified silicon vacancies. The following nomenclature is used to keep track of the different configurations: conf iguration layer distance index. Conf iguration is either h or k and refers to the silicon vacancy. Distance refers to the distance between the silicon vacancy and the carbon antisite inÅ. Layer refers to the position of the carbon antisite with respect to the silicon vacancy, either inplane (p), below (b), and above (a) the plane of the silicon vacancy. For example, the six closest configu-rations are denoted h a 3.1 , h p 3.1 , h b 3.1 , k a 3.1 , k p 3.1 , and k b 3.1 . These labels do not always correspond to a unique configuration, hence, in cases, where there are multiple non-equivalent atoms with the same distance and layer, these are separated with an index of I, II, or III; see Figure 1b.
Carbon antisite positions along the c-axis are labeled with (layer prime) to indicate the layer distance from the silicon vacancy plane. From our compact nomenclature, the carbon antisite shares the same local environment (h or k) with the silicon vacancy for even primes (such as p, a ), whereas for odd primes (such as a, a ) the antisite is placed in different layer than the vacancy. The ordering of h and k planes in 4H-SiC results in unique antisite positions for h and k vacancies. The closest modified vacancy with C 3v symmetry for the k silicon vacancy is a carbon antisite at 5Å whereas the nearest modified vacancy with C 3v symmetry for the h configuration is at a distance of 10Å. As we explain later in Sec. V B, the k site modified vacancies with C 3v symmetry are the likely candidates for the unidentified defects in EPR experiments.

IV. RESULTS
This section presents our high-accuracy computational results for the modified vacancy and experimental PL measurements. See Ref. VII for method details. The theoretical results include formation energy, optical properties like ZPL and transition dipole moment (TDM), electronic structure, and zero field splitting (ZFS). All the presented data is found in Table I and Table II. Formation energy (eV) 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 b 5.0 .
Next, we discuss the optical properties of modified vacancies. Figure

B. 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. VII B 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 ad- Energy (eV) Energy (eV) ditional lines which are listed in Table III and Table IV with the polarization given for the most prominent lines. 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 Table I and Table II, the 3rd nearest neighbor site k b 4.4 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 a  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.
configuration has the smallest D value (−13.34, see Table II or Figure 5c) of all the tested configurations and could be the defect responsible for the R 1 EPR signal [16]. The same configuration 30) also has the lowest D value. We draw the same conclusion as the authors in Ref. 30, the k a 5.0 configuration is responsible for R 1 EPR signal. With the additional configurations considered in this paper, it is clear that no other configuration has a smaller ZFS value than k a 5.0 . The other modified vacancy with C 3v symmetry, k b 5.0 , was not considered in Ref. 30 [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.

B. 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 does 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 [45].
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 V, and also marked in Table III and Table IV. 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 1 EPR signal, k a 5.0 has ZPL around 40 meV from the isolated h silicon vacancy with a parallel polarization, see Table II 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 1 signal [16].
Hence, the k a 5.0 configuration is identified as the source of the R 1 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 V shows that the h b 3.1 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 1 and e states compared to the isolated case, mentioned in Sec. IV. 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 contradicts this hypothesis. Alternatively, these lines could originate from strain, surface effects, or another defect. Several possible modified vacancy candidates from Table I match the experimental spectra.
By comparing the HSE results of the modified vacancies (Table I) with experiment (Table III), the best matches for both ZPL energy and polarization are: 862.2 nm is h a 6.9 and 863.2 nm is h b 7.8 . The 864.1 nm is matched to k p 6.2 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 p 5.4 I is either 866.6 nm, 869.6 nm (assigned), or 873.4 nm, and h b 5.4 could be either of 864.3 nm (assigned), 865.6 nm, 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.

C. 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 2+ , 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. Fig. 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 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 [47]. 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 [48,49].

VI. CONCLUSION
In this paper, we have searched for a defect that explains the experimental signals in the vicinity of silicon vacancy in a defect database generated by ADAQ containing 8355 single and double intrinsic defects. Based on this search, the modified silicon vacancy with a carbon antisite is the only defect candidate. The modified vacancy has similar properties as the silicon vacancy, and the carbon antisite is a small enough perturbation that does not change the spin or charge state of the silicon vacancy. This defect has several configurations which have been identified with the new PL experiments performed in this paper. Foremost, the k a 5.0 configuration has been identified to be the source of the R 1 EPR signal and 880.0 nm PL line. Moreover, the h b 3.1 configuration, which is one of the most stable and closest configurations, is most likely responsible for the 875.9 nm PL line. Other configurations explain lines closer to the isolated vacancy and might explain line broadening as well. However, for these lines, further studies and experiments are needed to rule out other effects, such as strain, before any conclusion can be drawn. Finally, our work demonstrates a high-throughput approach to search for point defects that explain experimental results is a promising direction forward.

A. 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) [50,51] 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) [52,53] 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 −6 eV for PBE, and 10 −4 eV for HSE. The structural minimization criterion is set to 5 × 10 −5 eV for PBE, and 10 −2 eV/Å for HSE. The pseudopotentials for C is labeled 05jan2001 and for Si is labeled 08april2002. Two supercell sizes have been employed with 576 (6x6x2) and 2304 (12x12x2) 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: • Formation energy [54] for the charged defect: 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 [55].
• The binding energy H b is defined as: 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: 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.
• ZPL is calculated as the total energy difference between the ground and excited state [56].
The excited state is found by using constrained DFT (∆-SCF method) [57,58]. The relative ZPLs can be compared within an accuracy of 100 meV [45]. The lowest energy optical excitation is from the occupied a 1 state to the unoccupied a 1 state in one spin channel [59].
• 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].
• ZFS is approximated using the dipole-dipole interaction of the spins [59,61] and is cal-