Ab initio description of highly correlated states in defects for realizing quantum bits

Coupled localized electron spins hosted by defects in semiconductors implement quantum bits with the potential to revolutionize nanoscale sensors and quantum information processing. The present understanding of optical means of spin state manipulation and read-out calls for quantitative theoretical description of the active states, built-up from correlated electrons in a bath of extended electron states. Hitherto we propose a first-principles scheme based on many body perturbation theory and configuration interaction and address two room temperature point defect qubits, the nitrogen vacancy in diamond and the divacancy in silicon carbide. We provide a complete quantitative description of the electronic structure and analyze the crossings and local minima of the energy surface of triplet and singlet states. Our numerical results not only extend the knowledge of the spin-dependent optical cycle of these defects, but also demonstrate the potential of our method for quantitative theoretical studies of point defect qubits. The electronic states of point defects in semiconductors can be studied via first principles, suitable for large systems with thousands of electrons. The nitrogen-vacancy centre in diamond and the divacancy complex in silicon carbide are promising candidates for quantum applications. The non-radiative decay from their optically allowed excited states is pivotal to initializing and reading their spin. Michel Bockstedte and colleagues now model the electronic states of these defects using a method that relies on many-body theory and is free of empirical parameters. They map the optical transitions, as well as the spin relaxation dynamics and the role of the spin-orbit and electron-phonon coupling interactions involved, in good agreement with experimental data. The method is suitable for large systems, and can therefore be used to model qubits in other semiconductors.


INTRODUCTION
Calculation of highly correlated electron states in contact with a bath of extended states is a long-standing problem 1 that is of immediate importance in such divers fields as the description of exotic phases in transition metal oxides [2][3][4] or the mangeto-optical properties of point defects in semiconductors [5][6][7][8] that comprise transitions between correlated localized in-gap levels and delocalized defect resonances placed in the range of extended states. To this end, we have developed a practical first principles method that can be employed for large systems including~2300 electrons. Here, we demonstrate on the most prominent point defects acting as solid state quantum bits, i.e., the nitrogenvacancy center in diamond 9 and the divacancy in silicon carbide, 10,11 that highly correlated electron states occur in these quantum bits that significantly interact with the electron bath of the solids. Their electronic states also strongly couple to the coordinates of the nuclei constituting the defects. Such electronic correlation and electron-vibration coupling play an essential role in the spin-selective non-radiative decay from the optically allowed excited state of these point defects that is the key of quantum bit initialization and readout.
The highly correlated electron states are built-up from localized orbitals with strong mutual Coulomb interaction, which is responsible for the high correlation. Such localized orbitals appear in the dangling bond orbitals of vacancies in covalent solids, for instance, nitrogen-vacancy (NV) center of diamond 12 and divacancy (V C V Si ) in silicon carbide (SiC). 13 In the former, the vacancy sits near a nitrogen atom that substitutes one carbon atom in the diamond lattice with an extra electron donated by other defects; in the latter two nearby vacancies form the neutral divacancy complex in SiC. In both cases, the optical and magnetic activity of the defects are related to the deep levels of the defects, particularly, caused by the dangling bonds of the carbon vacancy (see Fig. 1). These electrons form high-spin S = 1 states and highly correlated low-spin S = 0 multiplets (see Fig. 2) where the S = 0 1 A 1 multiplet is known to be responsible for the spin selective non-radiative decay from the S = 1 3 E optically allowed excited state of NV center. [14][15][16] The non-radiative decay is mediated by intersystem crossing (ISC) involving spin-orbit and electron-phonon couplings between the triplet and singlet states. Thanks to the effective ISC in these defects, 9,14,15,17 their electron spin can be initialized and read out optically, which makes it possible to apply single NV centers in diamond and divacancies in SiC as qubits. These solid state qubits are strong candidates to realize nanoscale sensors [18][19][20][21] and quantum communication devices. 22,23 It is an immediate quest to calculate these processes ab initio, in order to predict how to optimize the qubit operation for known qubits, and find new candidates. [24][25][26] The inevitable first step in determining the ISC routes and processes is to accurately calculate the highly correlated states and their energies as a function of the nuclei coordinates of the point defects.
A state-of-the-art route to calculate the excitation energies and states of defects in semiconductors 27,28 at ab initio level is the many-body perturbation theory within the GW approach and solving the Bethe Salpeter equation (BSE). 29,30 However, GW+BSE method yields results for the NV center in diamond 5 that are at variance with the models derived from experimental data. 16 Configurational interaction (CI) quantum chemistry calculations 6,31 on tiny diamond clusters and parameterized semi-empirical Hubbard-model calculations indicate 7,8 that BSE is unable to treat the highly correlated states arising in NV center. Nevertheless, unrealistically small defect models or methods with inherent semiempirical parameters do not allow systematic calculations of highly correlated states. An alternative route to defects in solids 3 is indicated by dynamic mean field theory where the physics of strong local correlation in the material is mapped to an effective local Anderson impurity problem with a few strongly interacting localized (d or f) orbitals of the host atoms (or the impurity atom) and an interaction with a bath of interant electrons. 2 The effective interaction can be constructed from ab initio for instance within the constrained RPA 32 and a solution of the impurity problem can be obtained via a truncated configuration interaction approach. 33 In case of the defects of interest in our work the strong correlation arises from the interaction of the defect orbitals with extended defect resonances and band states. The proper choice of the impurity problem thus is a task.
Motivated by the dynamic meanfield theory we have developed a methodology based on an screened effective interaction derived from many body theory and a configuration interaction approach that is free from empirical parameters and is applicable to necessarily large clusters, including periodic boundary conditions, to calculate the states and energies of highly correlated states. With this method we explore the magneto-optical properties and spin-relaxation dynamics of the two leading solid state qubits.

RESULTS AND DISCUSSION
CI-CRPA approach Description of the multiplet levels of the defect is equivalent to treat static correlation among the localized defect states within the band gap and extended defect resonances via the manyelectron hamiltonian whereâ y ,â are the creation and annihilation operators, respectively, acting on the single electron states i, j, k, l, and t ij accounts for a single electron part and the second term describes the electron-electron interaction via the bare Coulomb term (V) or as implemented here by an effective screened interaction (V eff ) as described below. Inclusion of the latter is essential for a tractable basis of single particle states allowing for a minimal subset of physically well motivated states, namely the defect states within the band gap and ideally only a few defect resonances, valence or conduction band states (see Fig. S1 in the Supplementary  Information). The effective interaction V eff is treated within the constrained random phase approximation (CRPA). 32 It is the screened Coulomb interaction using the CRPA dielectric function ε CRPA , i.e. V eff ¼ ε À1 CRPA V. ε CRPA is calculated like the RPA one but excluding transitions among those Kohn-Sham states that form the CI basis. V eff employed together with the RPA screening of the subspace ε subspace RPA reproduces the self-energy in the single-shot GW approximation of the full system, 32 i.e., the self-energy P ¼ ε À1 RPA V¼ε subspaceÀ1 RPA V eff . The CI basis is constructed from Kohn-Sham orbitals of the defect ground state employing hybrid density functional theory (DFT) as outlined in Methods. For intrinsic defects in diamond and SiC, this approach yields band gaps and defect levels in good agreement with experimental findings and sophisticated GW calculations. 34 Spin-symmetry in the many body hamiltonian is   7 The CI with bare Coulomb interaction is executed with a basis including the a and e levels as well as valence band states with The results for this basis are not fully converged regarding the 1 A 1 state but match a full CI calculation for a small cluster of ref. 6 Ab initio description of highly correlated states M Bockstedte et al.
enforced by utilizing a spin-restricted DFT calculation with equal occupation of both spin channels. The matrix elements t ij include the kinetic energy and interaction with the ion cores in a full CI approach. We here employ a limited basis set and therefore also include the interaction with the remaining itinerant electrons via the DFT mean field potential. As we outline in the Supplementary Information, we express the matrix by t ij ¼ ε HSE i δ ij À h dc ij using the Kohn-Sham eigenvalues ε HSE i and a double counting correction for the electron-electron interaction between the basis states.
Multiplet states of NV − and V C V Si We first demonstrate our ab initio methodology on the NV center in diamond that was modeled in a 512-atom supercell containing 2046 valence electrons in the system (see also Methods section). We are able to set the level of approximations and compare the results obtained by previous methods (see Fig. 2). This comparison highlights the complex physics of highly correlated states of defects for quantum applications and the need of parameter free treatment in realistically large cluster to treat these states. Particularly, our findings indicate that the screened interaction is a pivotal ingredient. We calculated the multiplet spectrum for the bare Coulomb interaction (denoted CI in Fig. 2) and effective interaction screened with the CRPA or RPA dielectric functions (CI-CRPA and CI-RPA, respectively). For the 3 E and 1 E states results from CI (2.2 and 0.7 eV) and CI-RPA (1.9 and 0.4 eV) are found within 0.2 eV of the CI-CRPA value. The most striking differences are found for the 1 A 1 , which is relevant for the ISC, and 1 E′ states. CI-RPA suffers overscreening effect because of double counting (see Supplementary Information). The molecular cluster calculations using full CI calculations 6 or cluster Hubbard-model 7 predict unreliable energy for the critical 1 A 1 state. In such small clusters bulk-like states are not fully developed both with respect to their extent and the energy separation between occupied and empty states, which yields a dynamical screening distinct from that of bulk diamond. Our CI-CRPA results for 1 A 1 and 1 E qualitatively agree with the GW + BSE periodic cluster approach (refs. 5,8 see "ext. Hub. BSE column" in Fig. 2 for the latter), however, 3 E and 1 E′ are found at different energies, mainly as the BSE basis set is limited to electron-hole pairs that is insufficient for highly correlated states. Our CI-CRPA hamiltonian reproduces the results obtained for an extended Hubbard hamiltonian. 8 There the hamiltonian parameters were fitted such as to reproduce GW quasiparticle energies of the localized defect levels and a resonance. However, this fitting procedure may be ambiguous as the number of resonant states needed for the correct description of the system cannot be defined in a consistent manner. Our method does not contain any fitting procedure and the number of states in the CI basis can be chosen to achieve convergent results (see Fig. S1 in Supplementary Information). In the particular case of the NV center, by taking only the in-gap defect states in the CI basis and the corresponding V eff already produce near convergent multiplet energies within 80 meV. This explains the success of the extended Hubbard hamiltonian approach 8 for this particular defect. Nevertheless, the defect orbitals (dangling bonds) can principally hybridize with the valence band electrons into several relevant defect resonances, thus any methodology with simple fit to a single resonant state does not provide consistent results in general as we show below.
We also employed our CI-CRPA method for the divacancies in 4H SiC. Among the four distinct configurations of the defect we here focus on the axial divacancy at hexagonal sites (see Supplementary Information for the divacancy configurations in 4H SiC and results for all configurations) that we briefly call here divacancy in the context. It is modeled in a 576-atom periodic cluster involving 2296 valence electrons. The CI-CRPA basis has to include valence bands and resonant states down to E V − 0.9 eV, specifically the a C resonance that is subject to a substantial exchange splitting in spinpolarized DFT (cf. Figure 1), in order to yield consistent and convergent results for this defect in SiC. This finding underlines the necessity for such methodology that does not contain any fitting parameters in the description of highly correlated electron states. The calculated low energy spectrum, including the highly correlated multiplets, is indeed very similar for diamond NV and SiC divacancy, just the energies of divacancy's states are scaled down. The divacancy's e Si -related many electron states are well decoupled from this range of the spectrum. This result fully demonstrates from ab initio that the NV center in diamond and the divacancy in SiC are indeed similar systems. This has been confirmed in recent resonant photoluminescence excitation experiments. 17 We find that the m s = ±1 3 E states are built-up from a single Slater-determinant, thus the constraint DFT approach (see Methods) is justified to explore the adiabatic potential energy surface (PES) of the 3 E excited state. 3 E and 1 E states are subject to Jahn-Teller distortion 35,36 where the energy of the multiplets is calculated as a function of the coordinates of the nuclei along the PES to minimize the energy of the 3 E state with the stable minimum ES 2 − cf. Figure 3. Furthermore a second metastable configuration ES 1 arises from non-linear Jahn-Teller coupling and many body effects as this was first discussed for the negative vacancy in silicon 37 (cf. Fig. S2 and the Supplementary Information for a detailed discussion of all Jahn-Teller distortions). Here we first focus on the minimum ES 2 .
For the NV center, the calculated zero-phonon line (ZPL) energies of the transitions 3 E → 3 A 2 and 1 A 1 → 1 E are about 0.2 eV lower in energy than the experimental values but the relaxation energy upon optical excitation of about 0.255 eV is well reproduced. The calculated energy difference of 0.5 eV between the 3 E and 1 A 1 states in the groundstate geometry of NV is in the right energy region estimated from experimental data. 16 The 1 E′ state lies high in energy and corroborates the expectation 38 of a strong 1 E → 1 E′ transition that could not be observed below a photon energy at and below 2 eV. It is, however, located above the ionization threshold to NV 0 at 2.2 eV. In the case of divacancies in 4H SiC experimental data is available for the ZPL energies that are reproduced within 0.1 eV by CI-CRPA (see Supplementary  Information). The small and consistent discrepancy between the theoretical data and CI-CRPA results on ZPL energies might partially come from the employed approximation in the calculation of electron correlation but also from the strong electron-phonon renormalization effect known in diamond 39 where the latter is beyond the scope of this study.

Spin relaxation dynamics
The spin selective non-radiative decay is a key mechanism for spin state initialization and readout of NV center and divacancy in qubit applications and has not been fully understood. The dynamics of the decay process depends on the coupling and ISCs of the triplet and singlet states and the interplay with phonon modes. The most important phonons are those that govern the Jahn-Teller distortion of the triplet and, as we discuss, give rise to a crucial pseudo Jahn-Teller coupling 40 among the singlet states. In the group theory approach to the spin-relaxation dynamics [14][15][16] multiplet states are identified with fundamental multiplet states of the high symmetry ground state, so far neglecting further coupling among states of the same representation via the Coulomb interaction or Jahn-Teller coupling. Here we advance the understanding of the relaxation dynamics using the adiabatic picture and based on the coupled fundamental multiplet states that are eigenstates of our CI-CRPA hamiltonian. The picture that arises from this analysis is depicted in Fig. 3c and d. The multiplet wavefunction of the singlet states in the ground and excited state configuration ES 1 and ES 2 of the NV center and divacancy can be seen in Supplementary Tables S5 and S6. The possible decay processes in the low energy branch in the NV center have not been completely understood so far, thus we consider these first. We find that the 1 E state, which splits to 1 A and 1 A′ in low symmetry Jahn-Teller configurations, already contains a contribution from the ae x y ð Þ À ae x y ð Þ multiplets (within the hole notation, e x (e y ) as indicated in the insets of Fig. S2) in the ground state geometry, thus the 1 E state is connected to the m S = ±1 state of the ground state 3 A 2 triplet via the perpendicular component of the spin-orbit coupling. Furthermore, the dominant contribution to 1 A 1 , the e x e x j iþ e y e y determinant, dynamically couples to the 1 E singlet state through a pseudo Jahn-Teller effect 40 driven by e symmetric localized phonon modes (compare 1 A 1 and 1 A states for the undistorted GS and distorted configurations ES 1 and ES 2 in Table S5). At the distorted configurations we observe a contribution of ≈13%. Since, e x e x j iþ e y e y and the m S = 0 spin state of 3 A 2 are coupled through the parallel component of the spin-orbit coupling, our results demonstrate a link between 1 E and the m S = 0 spin state of the triplet ground state. These results explain the experimental observations of decay processes between 1 E and all the spin states of 3 A 2 state 14 and extend the theory of spin state initialization and readout of the NV center. The spin selective decay in the high energy branch of the NV center has been recently discussed in the literature, 16,41 however, the contribution of the aa j i multiplet has been neglected so far. Our results show non-negligible (3.5%) contribution from this multiplet to the 1 A 1 state, which decreases the calculated coupling strength of 1 A 1 to the m S = ±1 states of 3 E, thus it reduces the rate of the spin selective decay as well.
From the comparison of the singlet wavefunctions of the NV center in diamond and divacancy in 4H SiC, one can see that both defects exhibit similar shelving states. The only notable difference is the stronger mixing of 1 A 1 state and the 1 E state in divacancy, which may result in faster decay from the 1 E state toward the m S = 0 spin state of the 3 A 2 ground state. From the calculated energy gap Δ of the 3 E and 1 E levels of the NV center in diamond and divacancy in 4H SiC, one can see that in the latter case the gap is much smaller, only 0.7 eV. The phonon spectral function of the divacancy for Δ = 0.7 eV is non-negligible, thus there is non-zero overlap of 3 E state and the phonon sideband of 1 E that give rise to an additional decay path between these states. The additional decay path explains the deduced rate of 15 MHz for the 3 E m S = 0 to singlets and 21 MHz 3 E m S = ±1 to singlets transitions in divacancy of 3C SiC 17 which has a similar electronic structure to that of the considered divacancy configuration in 4H SiC. 42 Note that, in the case of NV center in diamond, the 3 E m S = 0 to singlets transition rate is two orders of magnitude smaller than the 3 E m S = ±1 to singlets transition rate.
The above arguments show that decay rates depend not only on the matrix elements provided here but also on the spin-orbit coupling strength, non-radiative matrix elements, and the phonon spectral density. The complex PESs of the 3 E and 1 E with two distinct minima that arise from (pseudo) Jahn-Teller coupling (see Supplementary Information) as well as the coupling of multiplet states with two sets of E modes require further attention. In case of the spin-dependent relaxation rates of the divacancy, only the full set of E modes enables the proper inclusion of the ISCs that occurs between 1 E and 3 A 2 along the path GS-ES 1 (cf. Fig. S2) and between 3 E and 1 A 1 at ES 1 and ES 2 . Neglect of the former would lead to underestimated transition rates. The ISC between 1 E and 3 A 2 along GS-ES 1 is absent in case of the NV center and given the much smaller Jahn-Teller coupling in the 1 E mode a description with two E modes 35,41 should effectively cover the main physical mechanisms. Nevertheless, these findings require further investigations, in particular, regarding the calculation of spin-selective decay rates and non-adiabatic coupling that are beyond the scope of the present work.

Optical transitions
Optical transitions are at the heart of the photon-spin interface of qubit defects. These include radiative lifetimes and give access to spin-manipulation and charge state control. In the following and the Supplementary Information we employ ab initio multiplet states for the investigation of such optical properties.
In the case of NV center we calculated optical transitions from in-gap states to the diamond bands (see Supplementary Information) that are relevant for the optical excitation 43 and electrical readout 44 of the center. From these data one can calculate the radiative lifetimes τ NV of the excited states (see Fig. S4 and the Supplementary Information). The calculated 3 E τ NV = 8 ns is close to the experimental value at ≈12 ns. 16,45 For the singlet we obtain 1 A1 τ NV = 1878 ns that is two orders of magnitude longer than that for the triplet. This ab initio result sheds light on the experimental fact that not just the competition with non-radiative decay but the inherently tiny optical transition dipole moment is responsible for the very weak infrared PL signal of NV center (cf. the absorption spectrum in Fig. S4b). This result clearly indicates that sensing protocols based on infrared absorption of NV center 46 have serious limitation. The calculated absorption from the shelving 1 E state to the conduction band is considerable, and its transition energy is close to the energy of the 532-nm laser that is often applied in NV experiments. Taking into account the consistent underestimation of the optical transition energies, it can be concluded that 532-nm optical excitation cannot ionize NV center from the shelving state which explains the negative contrast in the electrically detected magnetic resonance spectrum. 44 Optical ionization or the strong transition 1 E → 1 E′ require larger photon energies. Since the 1 E → 1 E′ transition is embedded in the continuum of bound to free transitions just below an even Ab initio description of highly correlated states M Bockstedte et al.
stronger transition to the conduction bands at 2.62 eV (cf. Fig.  S4b), excitation protocols utilizing this transition have to deal with competing ionization via non-radiative processes. For the divacancy in 4 H SiC we obtain the lifetimes 3 E τ VCVSi = 63 ns and 1 A1 τV C V Si = 3109 ns, respectively, where the known PL lifetime of 3 E state is 15 ± 3 ns 19 that inherits direct radiative and non-radiative processes between the triplet excited state and ground state, and the non-radiative ISC processes via the singlets. It has been recently found in experiments 17 that the direct rate consisting of both the radiative and non-radiative ones and the ISC rate is about the same in divacancy. This results in ≈30 ns lifetime associated with the direct decay processes which is still about a factor of 2 shorter than the calculated direct radiative decay. We argue that the phonon-mediated direct non-radiative decay from 3 E to 3 A 2 is much more competitive in the SiC divacancy than that in NV diamond which naturally explains this discrepancy. In NV diamond, the optical transition energy (≈1.94 eV) is about twice larger than that of SiC divacancy (≈1.1 eV) but with similar optical transition dipole moments and phonon energies. This implies that the direct non-radiative decay by phonons are strongly suppressed in NV diamond while it is relatively efficient in SiC divacancy. We further note that the ≈0.7 eV 1 A 1 → 1 E optical transition does not likely occur in SiC divacancy because of the very slow radiative rate.
The present ab initio method employed here to analyze and predict properties of prototypical examples of qubits in solid can be applied to a broad class of materials with using realistically large models that exhibit highly correlated electron states interacting with extended electrons.

METHODS
The defect centers are analyzed in four stages: (i) The high spin ground state and the relaxed geometry of defects is obtained via spin polarized hybrid-DFT calculations. (ii) The lowest optically excited electronic configurations are determined by our DFT-based CI-method. (iii) Within the constraint DFT (CDFT) approach we occupy the defect states according to these excited states and relax the geometry of the defect in these states. This approach already proved to be successful. 25,26 (iv) The ground state and excited state geometries are connected by a linear path. For a few (i.e., 4-6) distinct geometries along this path we calculate the excitation spectrum using the CI-method (cf. main text and Supplementary Information). In the step (i) to (iv) all DFT-calculations were performed within the adiabatic approximation. Thus the potential energy surfaces obtained in step (iv) also rely on the adiabatic approximation.
With our CI-method we thus obtain high-spin and intermediate low-spin excited states as well as the important ISCs. The defect centers and their environment are represented in large super cells containing 512 atoms for the NV center and 576 atoms for the defects in 4 H SiC. The projectoraugmented plane-wave method along with the screened hybrid exchange-correlation-functional HSE06 47,48 are employed as implemented in VASP. [49][50][51] VASP was extended to calculate the CRPA response function and the matrix elements of the electron-electron interaction. We applied 420 eV plane wave cutoff and the usual carbon, nitrogen, and silicon PAW potentials provided by the VASP package. In the geometry optimization procedure the threshold on the forces acting on the atoms was set to 0.005 eVÅ. These parameters provided convergent results.

Data availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request.