First principles predictions of magneto-optical data for semiconductor defects: the case of divacancy defects in 4H-SiC

Study and design of magneto-optically active single point defects in semiconductors are rapidly growing fields due to their potential in quantum bit and single photon emitter applications. Detailed understanding of the properties of candidate defects is essential for these applications, and requires the identification of the defects microscopic configuration and electronic structure. Multi-component semiconductors often host two or more non-equivalent configurations of point defects. These configurations generally exhibit similar electronic structure and basic functionalities, however, they differ in details that are of great importance whenever single defect applications are considered. Identification of non-equivalent configurations of point defects is thus essential for successful single defect manipulation and application. A promising way to identify defects is via comparison of experimental measurements and results of first-principle calculations. We investigate a possibility to produce accurate ab initio data for zero-phonon lines and hyperfine coupling parameters that are required for systematic quantum bit search. We focus on properties relevant for the possible use of the divacancy defect in quantum bits in 4H-SiC. We provide a decisive identification of divacancy configurations in 4H-SiC and clarify differences in prior predictions of 4H-SiC divacancy zero-phonon photoluminescence lines.

To manipulate these centers on single defect level and to reconstruct their Hamiltonian, it is essential to identify the microscopic structure, electronic structure, and spin-configuration of the center. State-of-the-art experimental teqchniques used in experimental point defect investigation are for instance, photoluminescence (PL) or absorption spectroscopy, electron spin resonance (ESR), deep-level transient spectroscopy (DLTS), and Raman spectroscopy, that probe different characteristics of the centers. Gathering all the available information about a considered center can provide an appropriate working model. However, there are numerous unidentified defect centers in most of the commercially available semiconductors 31 .
In semiconductors where multiple non-equivalent sites exist in the primitive cells, each point defect can have several different configurations. These distinguishable configurations exhibit different properties and thus different applicability in qubit and single photon emitter applications.
The identification of such non-equivalent configurations is particularly challenging. For the nonequivalent configurations of divacancy related qubits in 4H-SiC two contradictory identifications have been presented, which rely on either the calculated zero-phonon photoluminescence (ZPL) lines 32 or the zero-field splitting parameter (ZFS) 33 . Furthermore, recently more divacancy related centers were reported than the possible number of non-equivalent divacancy configurations in SiC 18,19 , which makes the identification even more puzzling.
Identification and characterization of point defects are greatly facilitated by first principles theory. In supercell or cluster models, a small part of the material that embeds a single point defect is directly modeled in electronic structure calculations. This way many properties of the defects, such as spectral properties, charge transition levels, ESR parameters can be obtained, etc. In the literature, one can find different strategies how these quantities can be obtained in first principles calculations [34][35][36][37][38] . However, when comparing computational results to experiments, special care must be taken to account for numerical uncertainty and limitations in the theoretical methods 38 .
So far there has been limited discussion about how to generally achieve an accuracy sufficient for identification of non-equivalent defects. In the present paper, we address this issue.
We assess the accuracy of first principle calculations of ZPL and hyperfine interaction parameters to create guidelines for theoretical point defect calculations that allow non-equivalent defect configurations to be identified. In particular, we consider the divacancy defect in 4H-SiC and consistently identify PL1-PL4 room temperature qubits by comparing convergent magneto-optical date with the experiment. This defect has been studied with ab initio calculations before, but the present study helps clarifying previous results that (as discussed above) have not been fully consistent 32,33 . However, a main purpose of the present investigation is also to identify a scheme capable of reliably generating data via high-throughput calculations 39,40 useful for identification of, essentially, any previously unknown point defect. For this intended use, it is imperative to identify methods that reliably produce sufficiently accurate results, but also take minimal computational effort.
The rest of this paper is organized as follows. Section II describes the basic properties of SiC and divacancy point defect in 4H polytype of SiC. In Section III gives details on the first principle methods used in this work. Section IV presents the results and discussion of our first principles point defect calculations. In section V, we demonstrate how to use our results to identify divacancy configurations in 4H-SiC. Finally, section VI summarizes our findings.

II. DIVACANCY IN SIC
SiC is a polytypic semiconductor with more than 250 polytypes synthesized. The most commonly used forms are 3C, 4H, and 6H-SiC. The 3C polytype, shown in Fig. 1(a), has cubic symmetry with a single C and Si atom in the primitive cell. The 4H polytype, in Fig. 1(b), has hexagonal symmetry and 8 atoms in the primitive cell of which 2 are non-equivalent for both Si and C (see Fig. 1(a)). The 6H polytype is also of hexagonal symmetry, has 12 atoms in the primitive cell, and 3 are non-equivalent for both Si and C (see Fig. 1(c)). Hence, a single site defect in 4H has 2 distinguishable configurations, and 3 in 6H. A pair defect then has 4 and 6 configurations, re-3 spectively. The non-equivalent sites in 4H and 6H-SiC are refered to as h and k (4H), and h, k 1 , k 2 (6H). Here, h refers to a site in an hexagonal-like environment, and k to a cubic-like environment. In this paper, we focus on the four possible configurations of divacancy in 4H-SiC; hh, kk, hk, and kh, where the V Si −V C notation is used. For two of these configurations, hh and kk, the V Si −V C axis of the defect is parallel to the hexagonal axis of 4H-SiC and possess C 3v point group symmetry. The other two, hk and kh, have lower C 1h symmetry. hh and kk configurations are often called as axial configurations, while hk and kh as basal configurations.  Due to the C 3v symmetry in hh and kk, the dangling bonds of silicon and carbon vacancies form two fully symmetric a 1 and two double degenerate e states, which are occupied by 6 electrons 41 .

ES
The a 1 states are fully occupied, with the one localized on the silicon dangling bonds falling into the valence band, while the other one is localized on the carbon atoms and appears in the band gap of 4H-SiC near the conduction band edge. The e state is localized on the carbon dangling bonds and is located in the middle of the band gap and occupied by two electrons with parallel spin in the spin-1 ground state the divacancy. The other empty e state falls into the conduction band, as shown in Fig. 1(e). In the case of basal configurations, the low symmetry crystal field splits the e states into a and a and transforms a 1 to a .
Due to the spin-1 ground state and localized nature of the defect states, a strong dipole-dipole interaction can be observed between the unpaired electrons which causes a splitting of the spin sublevels even at zero magnetic field. For the divacancy defect, this zero-field-splitting is approximately 1.3 GHz. In SiC there are two intrinsic paramagnetic nuclei, the spin-1/2 13 C with 1.07% natural abundance and the spin-1/2 29 Si with 4.68% natural abundance, which can interact with the spin of divacancy and cause hyperfine structure in ESR spectrum. The spin density and important nuclei sites which yields resolvable hyperfine splitting of 10 − 100 Mhz are shown in Fig. 1(d).
In the single particle picture, the optically excited state of lowest energy can be constructed by a spin conserving promotion of the electron from the higher a 1 state to the e state, see Fig. 1(e).
Due to the partial occupancy of the e C state, the excited state is Jahn-Teller unstable, which causes spontaneous distortion of the atomic configurations of axial divacancies. In the many particle picture, six multiplets form in the triplet excited state that split according to the spin-spin and spin-orbit interactions. The divacancy in 4H and 6H-SiC has an electronic configuration similar to that of the divacancy and NV-center in diamond 11,18,41,42 , thus the many particle picture derived for the NV-center [42][43][44] can be applied for the divacancy too.

III. METHODOLOGY
The ZPL line is the energy difference between the ground state and excited state. These states can be seen in Fig. 1(e). The energy is obtained by using Kohn-Sham (KS) density functional theory 47,48 (DFT). In the excited state calculation, constrained occupation DFT scheme 49 is ap-plied, and accordingly a KS particle is promoted from the a 1 state to the empty e state in the minority spin channel and self-consistent energy minimization, including geometry relaxation, is carried out. In absolute values, one cannot expect better than 100 meV accuracy from this scheme, due to the single Slater determinant description of the excited state. This uncertainty of the theoretical method is an order of magnitude larger than the accuracy requirement of non-equivalent configuration identification (10 meV, which is the typical difference of ZPL energies for the divacancy defect). On the other hand, as the crystal field potential differs qualitatively only from the second neighbor shell and the defect state localization is decaying exponentially in this region, the electronic structure of the non-equivalent configurations can be considered as nearly identical. Hence, ZPL energy differences follow from the potential perturbations acting on the defect orbitals. Therefore, to identify the non-equivalent configurations, our DFT calculations must capture those effects that are caused by a perturbing crystal field potential. As the potential has a direct effect on the density and energy, such identification is likely to be possible through DFT ZPL energy calculations. In other words, for relative differences of the ZPL energies, one may expect better than 100 meV in constrained occupation DFT calculations. In the following, this topic is investigated in details by assessing technical and theoretical limitations of ZPL energy calculations.
We apply three exchange-correlation functionals in our calculations; the semi-local functionals of Perdew, Erzenerhof, and Burke (PBE) 50 and of Armiento and Mattsson (AM05) 51 ; and the screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE06) 52,53 . The hybrid functional is computationally much more expensive than the semi-local functionals, however, the band gap of semiconductors are closer to experiment 54 and accurate results in hyperfine field 55 as well as in zero-phonon line calculations 37,49 . All functionals are computed using the PBE pseudopotential labeled 05jan2001 for C and 08april2002 for Si.
The recommended procedure for the HSE06 hybrid functional is to start from a semi-local density, hence the following scheme is introduced. In practice, we employ the Vienna Ab initio Simulation Package (VASP) 56,57 , which uses the plane wave basis set and the projector augmented wave (PAW) 58,59 method to describe the KS 6 states and handle the effects of the core electrons. Since we need highly accurate results, we use comparatively high settings for those convergence parameters that are not further discussed in this study. The stopping criterion for the self-consistent field calculations and for the structural minimization are 10 −6 eV respectively 10 −5 eV (energy difference) for the PBE functional. For HSE06 functional, the settings are instead 10 −4 eV and 10 −2 eV/A (force difference). The grid for the Fast Fourier transformation (FFT) is set to twice the largest wave vector in order to avoid wrap around errors. For the HSE06 functional, the FFT grid for the exact exchange is set to the largest wave vector. This produces some noises in the forces but good energies. The abovedescribed settings ensure a numerical accuracy in the order of 1 meV for the calculated total energies. A Monkhorst-Pack 60 k-point grid is used for those calculations that use more k-points than the gamma point. In the rest this work, unit cell atom counts always refer to the number of atoms in a pristine supercell, if not otherwise specfied.
The zero-point energy shift of the ZPL energies due to the different vibrational properties of the ground and excited states is assumed to be small enough to not interfere with our conclusions.
It is neglected in the present study.
For hyperfine field calculations 55 , we use the implementation included in VASP, which gives the hyperfine tensor that describes the interaction between nuclear spin and electronic spin. This interaction produces a small splitting in energy levels which can be measured in experiments.
For zero-field-splitting, we follow the method of Ref. 61

IV. ACCURATE POINT DEFECT CALCULATIONS
In this section, results from zero-phonon line energy, hyperfine field, and stress calculations are presented.

Supercell size
First, we demonstrate the convergence of ZPL energies on supercell size. We begin by using supercells that retain the hexagonal symmetry of the primitive cell, the PBE functional, and Γ-point  the ZPL energies appear to converge to within 1 meV. Using this converged distance ≈ 30Å in the c-direction as well gives a supercell (10 × 10 × 3) consisting of 2400 atoms. We will use the converged values at this size as a benchmark for with which to compare other methods.

ZPL energy [eV]
Basal supercell size [Å] No. of atoms Note that, the observed finite size dependence is a combination of several effects, including the convergence of the charge density of 4H-SiC, exponential decrease of the self-interaction wavefunction of the defect, and the relaxation of the strain cased by the defect as the supercell size is increased. In the following subsections, we investigate these effects separately.

Brillouin zone sampling
In contrast to the straightforward study above of the convergence of ZPL on supercell size in a Γ-only k-point calculation, we now turn to convergence in smaller cells with higher Brillouin zone sampling. While a smaller supercell will increase errors due to vacancy-vacancy interaction, the aim is to investigate if one can still reach results that are accurate enough with significantly less computational overhead. Fig. 3 shows the PBE ZPL energies for different supercells of hexagonal and rectangular symmetry. As can be seen, the level of BZ sampling is important for all the considered supercells. The order of the ZPL energies of the non-equivalent divacancy configurations largely depends on the k-point convergence. The 576 atom supercell requires a 2 × 2 × 2 BZ sampling to provide the convergent order of PL lines, even though the absolute values are slightly smaller than for those convergent 2400 atom supercell. These results indicate that careful convergence in k-point density can produce the same order

Functional and geometry dependence
Next, we investigate how the choice of the functional and the way how the geometry optimization is performed affect the calculated ZPL energies. Fig. 6

A. Hyperfine field calculations
In this section, we discuss the methodological requirements needed for accurate hyperfine tensor calculations. As it was established previously, the most accurate values can be obtained by HSE06 functional including core state polarization effects 55 . Previously, however, no supercell size and BZ sampling tests were carried out.
First, we investigate the k-point convergence of hyperfine tensor elements. We use the PBE functional and study different supercell sizes. In the tests we consider two hyperfine parameters, the isotropic Fermi-contact contribution, A Fc = (A xx + A yy + A zz ) /3, and the dipole-dipole splitting, A dd = |(A xx + A yy ) /2 − A zz |. The obtained convergence curves for a set of 29 Si and 13 C sites close to a hh divacancy are depicted in Fig. 8. As one can see different sites exhibit different convergence behavior, which can be explained by the fact that the defect has different extension in different directions. In the basal plane, the extension is larger thus the defect states may require denser k-point grid in this direction. Furthermore, one can see that the criteria of complete convergence is similar to the criteria obtained for the ZPL energies.
Next, we investigate the supercell size dependence of A Fc and A dd hyperfine parameters. To do so we calculate the deviation of hyperfine parameters for numerous sites from the values obtained in a 2400 atom supercell calculation. In Fig. 9, relative errors of A Fc and A dd are shown for three smaller supercells. For all considered sites, the distance from the silicon vacancy site of the divacancy is also provided in the figure. As one can see, the relative error in the Fermi contact term increases dramatically with increasing distance of the divacancy and nuclear spins. Similar    Fig. 1(d). The right axis shows the hyperfine parameters for the 2400 supercell.
tendency can be observed for the dipole -dipole hyperfine term. Furthermore, relative errors drastically decrease with increasing supercell size. The mean relative errors for different supercells are provided in Table II. Note that substantial errors were obtained even in the 128 atom supercell, while the 576 atom supercell results are nearly identical with the absolutely convergent 2400 atom supercell results. These observations can be explained by the overlap of the state from the defect and its replicas. As further hyperfine interaction of the nuclei sensitively depends on the spin density localization in the farther neighbor shells, where defect state overlap can occur, they exhibit an enhanced finite size effect.  Calculations were carried out in 72, 128, and 576 atom supercell by using the PBE functional.
In summary, our results indicate that accurate hyperfine field calculations require supercells as large as 576 atoms or ≈ 20Å lateral size. In smaller supercells, such as 128 atom supercell, only closest hyperfine field of the nuclei can be determined with reasonable accuracy.

B. Defect induced stress in supercell calculations
Point defects induce a distortion in their host crystal, which relaxes with increasing distances from the point defects. In finite supercell models, the size of the model is usually not sufficient to As can be seen in Fig. 10, the stress indeed relaxes with increasing supercell size, however, it does not reach zero even in the largest 2400 atom supercell. Furthermore, the small difference between the 576 and 2400 atom supercell results may suggest that 576 atom supercell is convergent in terms of stress.
In order to make a rough estimate how the observed stress affects the ZPL energies, we imagine that the NV center is subject to 10 kbar external pressure, which is approximately the difference of the stress observed in the smallest and largest supercells. By using the experimental pressure dependence of the ZPL energy, 0.575 meV/kbar 62 , we obtain 5.75 meV. In SiC the effect could be larger, due to the smaller bulk modulus of SiC, however, presumably still remains in the order of 10 meV. Therefore, we conclude that the stress has a minor effect on the calculated ZPL energies.

V. IDENTIFICATION OF DIVACANCY RELATED ZERO-PHONON LINES IN 4H-SIC
In this section, we present an example of how accurate ZPL, ZFS, and hyperfine field splitting calculations can be used for identification of the divacancy configurations in 4H-SiC. We now consider a hypothetical case where one has been given the experimental ZPL lines in Table III for the four different defect configurations. Our aim is to identify the type of defect (vacancy, divacancy, interstitial, etc.), and match each line to a corresponding configuration (hh, hk, kh, and kk.) First, the variation of ZPL energies between different defect types are usually on a scale > 100 meV. Hence, it should be straightforward to make the identification of defect type with access to k-point-converged results for HSE06 for the 96 atom unit cell with the geometry converged using PBE or AM05. Even the PBE results for 96 atom supercell could work to identify the defect type, even though the absolute error is larger, the relative error is on the same scale (cf. Fig. 3(b) and Fig. 5(a)). Next, we turn to the identification of the defect configuration corresponding to each line. The conclusion from the present work is that the accuracy of the ZPL energies with the methods described here are not of sufficient quality to identify the configurations from them alone (cf. Fig. 7). As has been discussed above, while PBE predicts the correct order for the ZPL, HSE06 does not. However, if we in addition to the ZPL also have experimental results for the ZFS and the hyperfine tensor, they can be used to aid the identification. These quantities require calculations using HSE06 with one k-point and a supercell converged in size (as described in Sec. IV A and Ref. 55,61 61 . This property is assumed to have the same convergence as hyperfine field. The results of these calculations are summarized in Table III. Comparing with prior work, our identification agrees with that by Falk et al. 33 , which uses the ZFS parameter and ab initio simulations. The ZFS parameter was calculated using the data from a 1200 atom supercell, with Γ-point sampling, using the PBE functional. Our computational results predicts the wrong order and as discussed in the present work more accurate calculations do not resolve this but additional properties is needed for a correct identification.

VI. SUMMARY
This work discusses how to appropriately use ab initio calculations to facilitate identification of defect types and configurationsin semiconductors. Specifically, we have shown how to correctly identify the different non-equivalent divacancy configurations in 4H-SiC using ZPL, ZFS, and hyperfine field splitting calculations.
The value and order of the zero-phonon lines are dependent on the choice of functional, the way the geometry is optimized, supercell size, and k-points density. A comparably small supercell that has been converged with respect to the number of k-points produces sufficient accurate results at lower computational cost than the typical setup of a large supercell with one k-point. The absolute value of the zero-phonon line energy depends strongly on the lattice constant. A smaller lattice constant gives larger zero-phonon line energy and vice versa. The functional affects both the order and absolute value. It turns out that using only ZPL data is not enough to successfully identify the different non-equivalent configurations, but ZFS and hyperfine field can provide additional information.
For hyperfine field, the size of the supercell is the most important factor. Close to the defect, the values do not vary as the supercell size changes. But as the supercell size decreases, the hyperfine field gets less accurate the further away from the defect one calculates it.

21
The approach that produces both most accurate zero-phonon lines and hyperfine field values, at an affordable cost, is to use a large enough supercell that only Γ-point sampling is needed, relax it with PBE or AM05 functional, and run a single self-consistent HSE06 calculation. Using this proposed algorithm, we have shown how to correctly identify the different non-equivalent divacancy configurations in 4H-SiC.