{\emph{Ab initio} magneto-optical spectrum of Group-IV -- Vacancy color centers in diamond

Group-IV -- Vacancy color centers in diamond are fast emerging qubits that can be harnessed in quantum communication and sensor applications. There is an immediate quest for understanding their magneto-optical properties, in order to select the appropriate qubits for varying needs of particular quantum applications. Here we present results from cutting edge \emph{ab initio} calculations about the charge state stability, zero-phonon-line energies, spin-orbit and electron-phonon couplings for Group-IV -- Vacancy color centers. Based on the analysis of our results, we develop a novel spin Hamiltonian for these qubits which incorporates the interaction of the electron spin and orbit coupled with phonons beyond perturbation theory. Our results are in good agreement with previous data and predict a new defect for qubit applications with thermally initialized ground state spin and long spin coherence time.


I. INTRODUCTION
In recent years vacancy-impurity defects in diamond have become of high interest and are important because they show great potential in various quantum technology applications. In particular, the spin properties of the negatively charged silicon-vacancy [SiV(−)] color center [1][2][3][4] with a zero-phonon line (ZPL) energy at 1.682 eV and S = 1/2 spin has been recently studied for qubit applications [5][6][7][8][9][10] . As this defect has inversion symmetry (D 3d point group) it does not directly couple to external electric field, and as a consequence, SiV(−) possesses narrow 11 inhomogeneous linewidth and negligible spectral diffusion 7,8 . In addition, ∼70% of the total emission occurs in ZPL emission 5 , with a corresponding Huang-Rhys factor of 0.3. These properties are promising for realizing solid-state sources of indistinguishable single photons for quantum communication applications 8 . The fine splitting in the ground and excited levels caused by spin-orbit coupling (SOC) is harnessed to realize Λ scheme qubit operation 12 , however, the spin coherence time of 35 ns (see Refs. 7 and 13) is short because of the fast scattering of the electrons between the sublevels in the ground state mediated by the dynamic Jahn-Teller (DJT) effect even at T=4.5 K. Improvement on the spin coherence time can be achieved by cooling down the system below T=0.5 K 14,15 that suppresses the density states of the phonons that can mediate this process. Let us mention here that SiV(−) has also been proposed for optical thermometry at the nanoscale using the temperature dependent shift of the ZPL 16 . We further note that the neutral SiV, SiV(0), with a ZPL at 1.31 eV and S = 1 spin 17 , exhibits spin coherence time almost up to a second 18 and relaxation time nearly a minute 18 . SiV(0) is associated with the KUL1 electron paramagnetic resonance (EPR) center 4,17,19,20 . We note that SiV(0) can be found in special diamond samples where boron and silicon doping should be simultaneously realized during the diamond chemical vapor deposition (CVD) growth 18 .
Inspired by the success of SiV(−) color center, germanium, another Group-IV element in the periodic table, was introduced into the diamond lattice, either by CVD growth [21][22][23] or high pressure high temperature (HPHT) synthesis [24][25][26][27] . In all samples, a new ZPL line at 602 nm (2.06 eV) was observed in the PL spectrum, and unambiguously identified as a Ge related center because of the isotopic shift in the ZPL line 24 and its first vibronic peak 25 , and the anharmonicity of this peak 28 in the PL spectrum. Theory predicted 21,29 that this defect possesses the same D 3d symmetry like SiV does. The observed optical transitions support this conclusion 21,24 , thus the defect is indeed GeV(−). The 2.06-eV line has narrow (5 nm) ZPL linewidth even at room temperature, short excited state lifetime 30 (6-7 ns). The majority of the emission is concentrated in the ZPL emission 24 , with a Huang-Rhys factor of 0.5. Very recently, the spin relaxation and coherence times have been observed in Gedoped diamond samples at T 1 ∼ 0.34−25 µs and T 2 ∼ 19 ns (see Ref. 27), respectively, at T=2 K. We note that the signatures of GeV(0) with S = 1 spin have been observed in EPR spectrum in HPHT diamond 31,32 . This defect might have improved spin properties similar to those of SiV(0). The optical signature of GeV(0) has not yet been identified.
Simultaneously with the preparation of this paper, Snrelated PL centers have been reported in Sn implanted diamond 33,34 , where Sn is another Group-IV element next to germanium in the periodic table. In particular, the 620.3-nm PL signal showed the same fine level structure in PL like that of SiV(−) and GeV(−) 34 . By the use of this analog, they concluded that 620.3-nm center is associated with SnV(−). The spin properties and other charge states of this defect have not been reported. Next to tin, lead is the next Group-IV element in the periodic table. To the best of our knowledge, no Pb-related color centers have been reported in diamond so far.
Understanding the magneto-optical properties and spin coherence time of Group-IV -Vacancy color centers is of immediate interest and high importance in the fast emerging field of solid state qubits. Here we present a systematic study on the magneto-optical properties of the Group-IV -Vacancy defects, including Si, Ge, Sn and Pb impurities, by means of cutting edge first principles meth-ods. In Sec. II, we describe the first principles methodology for calculating the electronic structure, spin-orbit interaction and electron-phonon coupling of the systems. We then continue with detailed description of the results in Sec. III where we discuss the photostability and spin Hamiltonian of the qubits. We find that PbV color center exhibit superior spin properties over the other Group-IV -Vacancy color centers. Finally, we conclude the results in Sec. IV. We give additional data and derivation on the developed spin Hamiltonian in the Appendices.

II. METHODS
We characterize point defects embedded in diamond within spin-polarized density functional theory (DFT) as implemented in the vasp 5.4.1 code 35 . Our DFT method is within Born-Oppenheimer approximation, as ions are treated as classical particles. By varying the positions of ions one can achieve an adiabatic potential energy surface (APES) map of the system. The global minimum of APES defines the optimized geometry of the system. We reach this minimum upon relaxing the atomic positions till the force acting on every ion falls below 10 −2 eV/Å. We embed the point defects into a 512-atom diamond supercell, which suffices to sample the Brillouin-zone only at the Γ-point for converged charged density. A relatively low energy cutoff (370 eV) for the expansion of the plane waves within the applied projector-augmentation-wavemethod (PAW) 36,37 yields converged results. We calculated the excited states with the constrained-occupation DFT method (CDFT) 38 . We relaxed the atomic positions by minimizing the forces acting on them in the excited electronic state within CDFT method. The electronic structure is calculated using HSE06 hybrid functional 39,40 within DFT. This technique reproduces the experimental band gap and the charge transition levels in Group-IV semiconductors within 0.1 eV accuracy 41 . This procedure also yields excellent results for the zerophonon-line energy of SiV(−) center in diamond 4 . We determine the adiabatic charge transition levels or photoionization energy thresholds as where the ∆E q is the total energy correction for the point defect with q charge by following the procedure of Lany and Zunger 42,43 , while the E q tot the total energy of the system including the ions and electrons. We will provide these charge transition levels graphically with respect to the valence band maximum in Fig. 1.
For the calculation of the phonon sideband of the PL spectrum, we calculated the vibration modes of the defects in a quasiharmonic approximation in the groundstate at high symmetry, with equal occupation of the degenerate orbitals in the band gap. We used the numerical derivatives of the forces to generate the Hessian matrix that we diagonalized to obtain the phonon frequencies and normal modes. The geometry is preoptimized with the very strict force criterion of 10 −4 eV/Åfor the vibration calculations. We consistently applied the computationally powerful Perdew-Burke-Ernzerhof (PBE) functional 44 in this procedure that reproduces the calculated HSE06 phonon spectrum of the perfect diamond within 5 cm −1 and the quasilocal modes of SiV(−) within 2 cm −1 .
We determined the spin-orbit coupling (SOC) in noncollinear approach as implemented in VASP 5.4.1 for the negatively charged Group-IV -Vacancy centers, or briefly, XV centers. We set the quantization axis of the spin along 111 direction, the C 3 rotation axis of the XV point defects. We determined the SOC parameters by HSE06 DFT functional, which has generally provided accurate results for the spin-orbit splitting of nitrogenvacancy center in diamond 45 . SOC is a small perturbation to the electronic structure of the system, thus we fixed the atomic positions in the high symmetry D 3d configuration as obtained from the previous spin-polarized DFT geometry relaxations. As was reported in our previous study 4 for SiV center and is also shown in Sec. II for other XV centers, a double degenerate e g{x,y} level appears in the gap, which will be occupied by three electrons in the negatively charged state, which can be treated as a single hole on this state. After applying the SOC on the system, e g± = 1 √ 2 (e gx ± ie gy ) states split by λ 0 (see Fig. 1 for the electronic structure) which comes from the z component of SOC. By using CDFT procedure it is feasible to introduce the hole either on the e g+ or e g− state, and the calculated total energy difference is the strength of spin-orbit coupling. The Hamiltonian of SOC coupling is the following: whereL z is the effective orbital moment operator (L = 1) of the electron, whileŜ z is the electronic spin. We find that half of this total energy difference is equal within 10 −7 eV to the split of e g+ and e g− Kohn-Sham levels when these two states are occupied by half-half electrons.
Thus, the strength of SOC can be calculated by the halfhalf occupation of the e states with following the SOC splitting of these e g states. The negative sign of λ 0 accounts for the fact that e g particle is a hole and not an electron. The SOC Hamiltonian (2) can be expressed with the Pauli matricesĤ SOC = λ0 2 σ y where that represent the e g± electron states. We note that the dynamic Jahn-Teller (DJT) effect quenches the orbital moment, at least partially, known as the Ham effect [46][47][48] . The intrinsic λ 0 that we calculate directly from Kohn-Sham orbitals, is severely reduced by the Ham reduction factor p, thus λ Ham = pλ 0 reduced value is observed in the experiments.
We now discuss two cases: (i) the electron-phonon coupling manifested as DJT effect is significantly larger than spin-orbit coupling, so one can solve first the electronphonon system and then calculate the spin-orbit ener-gies as a first perturbation acting on the resultant vibronic wavefunctions; (ii) the electron-phonon coupling and spin-orbit coupling are in the same order of magnitude, so the orbital, spin, and phonon degrees of freedom of the wavefunction are strongly coupled which requires the exact diagonalization of the sum of spin-orbit and electron-phonon Hamiltonians, respectively. We call the first case as Ham reduction factor solution whereas the second case the exact diagonalization procedure.
We first describe the case (i) where we define the electron-phonon Hamiltonian caused by DJT effect. The DJT effect is the interaction of an E g type quasilocalized vibration mode with the e g electron orbital, known as E × e DJT system. The Hamiltonian of such system is the following: where the a † x,y , a x,y are the creation and annihilation operators of the E g vibration mode, respectively. The system is a twodimensional harmonic oscillator with frequency ω e , where the terms labeled by F and G parameters are the corresponding linear and quadratic part of DJT. The σ i operators are twodimensional Pauli matrices representing the e gx and e gy electrons, respectively. F and G parameters are the values describing the Jahn-Teller distortion through (x,ŷ) = 1 √ 2 a † (x,y) + a (x,y) operators. F and G parameters can be easily derived after the APES of the DJT system is determined. The energy gain from the symmetry distortion (see Chapter 3.2. in Ref. 47) to one of the three global minima is E JT = F 2 2(hωe−2G) , and the barrier energy separating the these three minima is δ JT = 4EJTḠ hωe+2G . For SiV(−), GeV(−), and with less extent, for SnV(−) in the ground state, the energy gain from DJT, i.e., E JT is orders of magnitude larger than the energy of SOC coupling. Thus, SOC can be evaluated as a perturbation on the DJT groundstate wavefunction, where the electrons and vibrations are entangled. We determine the eigenvalues of the Hamiltonian in Eq. (4) numerically with the following series of expansion, [c nm |e g± |n, m + d nm |e g∓ |n, m ] . (6) We then express the SOC splitting within perturbation theory as λ Ham = 2 Ψ ± Ĥ SOC Ψ ± = pλ 0 . We limit the series expansion of the twodimensional harmonic oscillator |n, m up to ten quanta, thus n + m ≤ 10. From this expansion, the p reduction factor can be expressed as p = n,m c 2 nm − d 2 nm with c nm and d nm expansion coefficients (see Appendix B for derivation) .
We next describe case (ii) where DJT and SOC energies are comparable, so the corresponding Hamiltonians in Eqs. (4) and (2) should be added and solved simultaneously with coupled orbital, spin and phonon wavefunctions. Here we expanded the polaronic wavefunctionΨ Γ with the spin degrees of freedom as nm |e g± |n, m |χ + d χ nm |e g∓ |n, m |χ ] (7) and then directly diagonalize the Hamiltonian including the SOC and DJT effects simultaneously, where χ can be either ↑ or ↓ spin state. This solution represents a coupling between spins and phonons and goes beyond the perturbation theory of SOC acting on the polaronic wavefunctions. Here, the subindex Γ refers to the total angular momentum of the wavefunction which is either 3 /2 or 1 /2 (see also Appendix C). We note that the simultaneous treatment of SOC and DJT in E × e JT systems has been only considered for very small molecules at ab initio level in the literature so far 49,50 . However, extensive ab initio study on point defects in solids, with a model consisting of hundreds of atoms, has not yet been performed, to the best of our knowledge. As we will show below our methodology is able to reproduce previous experimental data.

III. RESULTS
All the Group-IV impurities reside in the symmetric split-vacancy configuration in diamond according to our calculations that results in a D 3d symmetry. The split-vacancy configuration may be labeled as V-X-V which implies that the 'X' impurity atom lies at an interstitial position, more precisely, at the inversion point of the diamond lattice halfway between two adjacent vacancies, or briefly, in divacancy. Nevertheless, the quantum optics groups labeled these defects by XV in the literature, thus we consequently used this notation in the Introduction and the rest of the paper. According to this description, this type of defects exhibits an inversion symmetry which adds a parity to the wavefunctions, odd (ungerade) and even (gerade), labeled by u and g, respectively. One can construct the possible single electron orbitals from the defect-molecule model. There are 6 carbon dangling bonds pointing towards the impurity atom, from which 6 orbitals emerge in D 3d symmetry: a 1g +a 2u +e u +e g . The corresponding orbitals are filled with 10 valence electrons where 6 electrons comes from the dangling bonds and 4 electrons from the impurity atom. As discussed previously for SiV defect 4 , only the e g orbital appears in the gap which is filled by two electrons with parallel spins in the neutral charge state (see also Fig. 1) whereas the a 1g and a 2u levels fall inside the valence band and e u level is resonant with it. The same e g orbital occurs in the gap for the other XV defects. The in-gap optical excitation involves an e u orbital resonant with the valence band (VB) that pops up in the gap in the excited state electronic configuration. Special bound exciton states might occur between in-gap ↔ band edge optical excitation that might also lead to photoionization. The electronic structure of the neutral XV complexes implies that charge states from (2+) to (2−) can exist depending on the position of the (quasi) Fermi-level in diamond. Here, we focus our study on the (−) charge state (coherent dark states in the prototype SiV(−) qubit 10 ) and the neutral charge state (optical spin polarization with long spin coherence time of prototype SiV(0) 18 ), i.e., their charge state stability. We provide trends for the magneto-optical properties of the negatively charged XV complexes.

A. Charge state stability of the XV defects
A general trend in the electronic structure of XV defects is that the e g level shifts up with increasing atomic number of Group-IV impurity atom (see Fig. 1) which has a consequence on the formation energies and the corresponding adiabatic charge transition levels too. In order to readily see the trends, the calculated (−)→(0) and (0)→(−) transition energies are also plotted in Fig. 2, where the first and second transition is associated with promoting an electron from the defect level to the conduction band and from the valence band to defect level, respectively, but the plot depicts the charge transition levels with respect to the valence band maximum. The stability window of (−) state shifts up in the gap with increasing atomic number of Group-IV impurity atom. SnV(−) and PbV(−) can only be stable by providing substitutional nitrogen donors in the diamond sample.
The photostability of PbV(−) requires special attention as the (−) state can be converted to (0) by illuminating the sample in the visible region, by ∼2.6 eV, that it is close to its ZPL energy of about 2.4 eV (see Table I).
On the other hand, photoexcitation in the visible region may convert (−) state to (2−) charge state for SiV, GeV and SnV defects. For isolated SiV(−), ultraviolet (UV) light would be needed to reionize (2−) to (−) by single photon absorption that is difficult to realize in the experiments. On the other hand, violet and blue illumination can reionize GeV(2−) and SnV(2−) to GeV(−) and SnV(−), respectively. The (2−) charge state is a spin singlet, nevertheless, a shelving triplet bound exciton state might exist that can be accessed by optical pumping at ∼2.4 eV and ∼1.8 eV for SnV(2−) and PbV(2−), respectively. The stability window of S = 1 (0) state of GeV, SnV and PbV also shifts up in the gap with increasing the atomic number of the Group-IV impurity atom.
One can conclude that SnV(−) defect should be very photostable whereas isolated SiV is trapped in the (2−) charge state once SiV(−) has been photoionized into that state. On the other hand, there is a small energy margin between the calculated neutral excitation energy (∼2.43 eV) and (−)→(0) transition energy (∼2.6 eV) for PbV(−), thus this defect can be photoionized into the (0) charge state by blue illumination into the phonon sideband.
We note that the neutral XV(0) with S = 1 ground state may also act like a qubit, with presumably long electron spin coherence time in good quality of diamond. Our calculations show that SnV(0) and PbV(0) defects can be engineered into typical diamond samples where nitrogen contamination occurs, in contrast to the case of SiV(0) defect, which requires boron doping of diamond. Formation of GeV(0) requires very low nitrogen concentration or compensation of nitrogen donors by acceptors.

B. Photoluminescence of XV(−) defects
The photoluminescence can be described by spontaneous emission from the optically lowest energy excited state to the groundstate. The lowest energy excited state can be understood as promoting an electron from the e u orbital to the e g orbital in XV(−) defects. As a consequence, the excited state is 2 E u and the ground state is 2 E g in the negatively charged state (see Ref. 4 and Fig. 1). Both 2 E u and 2 E g states are dynamic Jahn-Teller systems 4,6,51 . We find a general trend in the calculated Jahn-Teller energy, E JT , as a function of the atomic number of the Group-IV impurity atom, where E JT is defined as the total energy difference in the high symmetry D 3d geometry and the lowest energy C 2h geometry in the adiabatic potential energy surface (APES) (see Fig. 3). The ZPL energies are calculated by taking the lowest APES energy in C 2h symmetry both in the ground and excited state (see Table I) that we call here an "average" method. Interestingly, the calculations do FIG. 1. Kohn-Sham single particle levels and charge transition levels of XV defects for X=Si, Ge, Sn, and Pb. We note that the effect of spin-orbit interaction is not included here. We represent the ↑ (↓) spin channel with triangles pointing upwards (downwards). The filled (empty) triangles depicts occupied one electron (hole) orbitals. In the (+) and (−) charge states, the symmetry is lowered from D 3d to C 2h within adiabatic potential energy surface in the DFT calculations. We show the charge stability window for the given charge state referenced to the valence band maximum of diamond at the left panel in each figure. The excitation processes are also shown for XV(−) color centers. The eu orbital is resonant with the valence band that pops up in the excited electron configuration. not predict linearly increasing trend in the ZPL energies by increasing the atomic number of the Group-IV impurity but the ZPL of SnV(−) should be smaller than that of GeV(−). This agrees well with the very recent experimental data 34 .
Regarding XV(−) qubits or qubit candidates, we show the calculated average ZPL energies compared to the experimental ones when available (see Table I). The calculated average ZPL energies, that do not contain the zeropoint energies and spin-orbit couplings in the ground and excited state, somewhat overestimate the experimental ones. By an accurate calculation of the zero-point energies within our DJT treatment together with the spinorbit coupling, brings the computed ZPL energies closer to the experimental ones, with respect to the average ZPL energies. Nevertheless, both methods are accurate within 0.1 eV (see Fig. 2). This gives us confidence that the calculated ZPL energy of PbV(−) is well predicted. As can be read out from Table I the nature of ZPL transition involves polaronic states together with spin-orbit effects (see Appendix C). Thus, perturbation effects, e.g. strain or temperature, on ZPL energies should involve the complex analysis of this coupled spin-orbit-phonon system. The calculation of spin-orbit coupling and electronphonon coupling will be described in the next sections. Because the excited states can be well described by promoting an electron from the e u orbital to the e g orbital according to our CDFT calculations, therefore, the optical transition dipole moment for XV(−) defects can be well approximated by calculating the optical dipole moment µ between these Kohn-Sham states. The radiative lifetime of these color centers then can be calculated 56 as (e 4 u e 3 g )~λ In the ground state (g) the global APES minimum can be found at Rg distance (C 2h distortion) from the high symmetry D 3d geometry while this is Ru for the u excited state. EJT is the Jahn-Teller energy is the energy difference between the D 3d and C 2h configurations. δJT energy barrier between the three local minima of C 2h configurations (few meVs) is not shown here for the sake of clarity. On the right panel (c), we schematically depict the phonon sideband of the ground and excited state where blue (red) lines represent the excitation (luminescence) of the system into the phonon sideband, whereas the black line defines the zero-phonon line transition (ZPL). The spin-orbit coupling (λ) split 2 Eg and 2 Eu ZPL states apart into E g1/2 ,E g3/2 ,E u1/2 ,E u3/2 Kramers doublets where 3 /2 and 1 /2 refers to an effective total angular momentum of the electronic states. The doublets may further split into individual e ↑/↓ g,u± states under external static magnetic field.
where n is the refractive index andhω is the excitation energy. The calculated radiative lifetimes are listed and compared to the observed PL lifetime of the XV(−) defects in Table II. A general trend is that the computed radiative lifetime somewhat decreases with increasing atomic number of the impurity atom but basically they are all in the same order of magnitude. The predicted short radiative lifetime (≈3 ns) of PbV(−) center is favorable for quantum emitter applications. We note that the observed PL lifetime of SiV(−) is significantly shorter than its computed radiative lifetime. We attribute this effect to the strong non-radiative processes which we tentatively assign to the ionization process that competes with the neutral excitation process (see also a recent photoluminescence excitation observation in Ref. 57). The calculated (2 − |−) charge transition level of SiV(−) at ≈2.05 eV is very close to the energy of the phonon sidebands of neutral excitation whereas the energy between these two increases with increasing atomic number of the impurity atom that should suppress this type of nonradiative processes. In addition, the optical gap also significantly larger for GeV(−), SnV(−), and PbV(−) than that of SiV(−), which also suppresses the direct nonradiative decay induced by phonons in the former color centers.
The phonon sideband in the PL spectrum is determined within the Huang-Rhys (HR) theory (see the original theory in Ref. 59 and our implementation in Ref. 60   TABLE II. The calculated radiative lifetimes (τ rad ) versus the observed photoluminescence lifetimes (τPL) for XV(−) color centers at cryogenic temperatures. We used the experimental ZPL energy where available in the calculation of τ rad . We note that τPL involves both radiative and non-radiative processes. that is based on Ref. 61). We calculate the HR spectra between two statically distorted Jahn-Teller structures as depicted in Fig. 3. By this way, we can take into account the contribution of the e g phonons in the PL spectrum that are responsible for the most intense phonon sidebands for SiV(−) and GeV(−) defects. These e g phonons are bulk-like in nature and do not localize on the defect. The participation of these e g phonons in the phonon sideband is a consequence of the dynamic Jahn-Teller nature of the groundstate and excited state. We note that quasilocal vibration modes are also visible as relatively sharp features in the PL spectrum of SiV(−) and GeV(−) at about 62 meV and 43 meV, respectively, that are not reproduced by our method. Based on the calculated vibration spectrum of SiV(−) defect in our previous study 62 , we associate this feature with the e u quasilocal vibration modes of the defects that involve the (x, y) motions of the impurity atom. Principally, the usual Franck-Condon approximation on the luminescence of polyatomic systems does not allow the participation of ungerade modes in the PL process, thus the observations might be explained by invoking the Herzberg-Teller effect that goes beyond the Franck-Condon approximation (see Ref. 62 for details). This issue is beyond the scope of the present study, and we rather focus on the general trends in the PL spectra. As can be seen in Fig. 4 the contribution of phonons to the PL spectrum increases with increasing atomic number of Group-IV impurity. In particular, the contribution of a 1g phonons significantly increases because the geometry change between the ground and excited state's geometries increases, as measured by the Huang-Rhys factor S. This can be understood by the size of the heavy atoms that cannot readily be accommodated by the divacancy of diamond and they start to substantially distort the diamond lattice. We note that this trend is disadvantageous for creating very efficient single photon sources emitting light dominantly in the zero-phonon emission, nevertheless, the calculated S = 1.60 factor for PbV(−) is still much smaller than that S ≈ 3.5 for NV(−) center in diamond (see Ref. 45 for detailed analysis).

C. Effective spin-orbit coupling in XV(−) defects
A splitting occurs in both the 2 E u and 2 E g levels due to spin-orbit coupling (SOC) between the S = 1/2 electron spin and the double degenerate orbital forming two Kramers doublets that we call zero-field-splitting (ZFS) where zero-field refers to zero external magnetic field. In DJT systems, the spin-orbit coupling can reduce the effective SOC by the p Ham reduction factor 46,63 . We use here the same ab initio theory to calculate the p factor as we demonstrated for the 3 E excited state of the negatively charged nitrogen-vacancy center in diamond 45 . This requires to calculate the full APES in the corresponding electronic state as depicted in Fig. 3. Here we note that E JT decreases with increasing atomic number of Group-IV impurity (see Table III). This results in smaller damping of SOC with increasing atomic number of Group-IV impurity. The intrinsic spin-orbit coupling is determined as described in Sec. II. Accurate calculation of spin-orbit coupling requires scaling method in giant supercells (see Ref. 45 and Appendix A). The results are summarized in Table III.
The general trend is that the λ 0 rapidly grows with increasing atomic number of Group-IV impurity. In the 2 E g ground state, there is a turning point for PbV(−) defect where λ 0 is greater than E JT , thus SOC is not a small perturbation w.r.t. electron-phonon coupling but the electron-phonon Hamiltonian has to be parallel diagonalized together with the spin-orbit Hamiltonian (see Sec. II and Appendix C). As a consequence, the estimated ZFS between the sublevels of 2 E g ground state is around 18.7 meV with the Ham reduction scheme (λ Ham ), and 18.1 meV (λ) with the exact diagonalization (see Tab. III for details). The p Ham reduction parameter also increases with increasing atomic number of Group-IV impurity because E JT decreases whereas the vacancy related ω e vibration modes are relatively insensitive to the type of Group-IV impurity atom.
We calculated the p reduction factors for the 2 E u excited state of the XV(−) defects too (see Table III). In the excited state, the Ham reduction factors are significantly smaller than those in the groundstate because the E JT energies are larger in the excited state than those in the groundstate (c.f. Table III). We think that the larger electron density in the interstitial region around the impurity atom in the 2 E u state contributes to form long bonds between the carbon atoms around the impurity atom and thus makes the Jahn-Teller distorted structure more favorable in that state than in the 2 E g state. Surprisingly the Ham factor p scheme provides reasonable results even for the optically excited state of PbV(−) for which λ 0 > E JT , nevertheless, exact diagonalization of the adjoint DJT and SOC Hamiltonians is needed for accurate results. We show a graphical interpretation for the spin-orbit and electron-phonon coupled systems in Fig. 5. We expand the series expansion with the spin degrees of freedom in Eq. (6) then directly diagonalize the Hamiltonian including the SOC and DJT effects simultaneously. As long as λ 0 E JT SOC is only a small perturbation over the polaronic DJT groundstate, and the 2 E g,u 4× degenerate level is split into double degenerate E g,u 3 /2 and E g,u 1 /2 levels. We label this splitting with λ in Fig. 5 that is directly observed in the fine structure of the ZPL optical emission. The perturbative approach of SOC is valid mostly for the SiV(−) system, especially on its optically ground state (a), thus approximation of λ with λ Ham is valid. If SOC energy is higher than E JT then Ham reduction scheme still provides surprisingly good results when compared to those from exact diagonalization. Considerable deviations only begin to appear for PbV(−) (see also Table III). For SiV(−), the plotted data is symmetric respect to the x = 0 axis as the E u 1 /2 is increased by pλ 0 /2 energy and E u 3 /2 is decreased by the same amount. However, for the excited state of PbV(−), this is not symmetric. The SOC systematically shifts the eigenvalues to the left in x-axis. In this case, the E u 1 /2 state contains larger contribution from spin-orbit favored e ↓ u− and e ↑ u+ states rather than from the unfavorable e ↑ u− and e ↓ u+ states, indicating that the DJT and SOC Hamiltonians should be solved simultaneously.

D. Spin Hamiltonian for Group-IV -Vacancy qubits and its implications
By applying an external constant magnetic field, the spin double degenerate levels of the XV(−) color centers  may split. Previously, a spin Hamiltonian was deduced for SiV(−) qubit 6 to describe this feature that we further develop based on the coupled spin-orbit-phonon Hamiltonian for the groundstate (g) and excited state (u) as followŝ where g S = 2.0023 is the g-factor of the electron, µ B is the Bohr-magneton of the electron, B is the external homogeneous magnetic field, B z is its z component with z-axis parallel to the symmetry axis of the defect. We note that the hypefine interaction between the electron spin and nuclear spins in the diamond lattice or with the impurity atom is not considered here. In the braces we merge the different effects into single effective parameters where we follow the nomenclature of Ref. 6 for the two common parameters λ and f , whereas parameter δ f appears as a new parameter according to our derivation. Our derivation reveals the microscopic origin of the spin Hamiltonian merged parameters that will be discussed below. The operatorL z acts on the E g,u orbitals as ±1, where g and u refers again to the parity of the wavefunctions with g groundstate and u optically allowed excited state.
The derivation of the terms can be found in Appendix C and D. Here we discuss all the resultant terms in our spin Hamiltonian in details. The first term in Eq. (9) contains an effective spin-orbit splitting. We first note the negative sign which originates from the threeelectron many-body E g and E u states. As a consequence, E g,u 3 /2 is lower in energy than the E g,u 1 /2 , in contrast to TABLE III. The calculated basic parameters of the APES such as EJT Jahn-Teller energy, δJT barrier energy and thehωe energy of the effective eg phonon driving the Jahn-Teller effect are shown as well as the calculated λ0 intrinsic spin-orbit coupling and p Ham reduction factor, and the deduced λHam = pλ0 effective spin-orbit coupling for the 2 Eg optically ground state and 2 Eu excited state of XV(−) defects. The calculated zero-field splitting (λ) is for comparison to the experimental value (λexp). We determine the calculated λ values beyond the simple Ham reduction theory where we treat the SOC and DJT Hamiltonians simultaneously (see Fig. 5 for graphical interpretation). the previous assignments. Here, we label the states with m j quantum numbers which is the sum of m l = ±1 orbital angular momentum and the m s = ± 1 /2 spin quantum numbers. Since the orbitals are coupled to phonons the Ham reduction factors p can be different forẼ g,u 3 /2 andẼ g,u 1 /2 polaronic wavefunctions (see Eqs. (C8)), and the final p g,u will be the average of the two (see Eq. (C5)).
In addition, the vibronic zero-point energy of these states can also differ, in principle, that will change the energy gap between these two states that we label by K g,u as defined in Eq. (C10). We note that K g,u can be neglected for SiV(−) and GeV(−) but it becomes substantial for SnV(−) and PbV(−). In the second term two reduction factors appear. The p g,u is already introduced above and caused by electronphonon coupling. The g g,u L orbital reduction factor was previously discussed by Stevens 64 , thus we call it Stevens' orbital reduction factor. This originates from the fact that orbital angular momentL z is only an effective operator as the D 3d point group of the XV(−) systems does not respect the full O(3) rotational symmetry as illustrated in Appendix D. We find that the g g L is in particular substantially smaller than one, and significantly reduces the effectiveL z (see Appendix D).
The third term is the usual Zeeman term for electron spins. The fourth term provides a correction to the g S constant, thus modifies it to a tensor where the corrected zz component is caused by the polaronic nature of thẽ E g,u 3 /2 andẼ g,u 1 /2 states where the δ g,u p scales with the difference of the corresponding Ham reduction factors (see Eqs. (C9) and (C5) in Appendix C). This is a new correction which is negligible for SiV(−) but already appears for the excited state of GeV(−), and both for the ground and excited state of SnV(−) and PbV(−).
We propose that theΥ JT operator in a previous study 6 that was associated with the Jahn-Teller effect originates from the residual strain around the individual SiV(−) centers, thus we propose to label as aΥ strain operator instead. We note that the strain was also considered in that study which can be very strong (in the order of 100 GHz) in nanodiamond samples and much smaller in the bulk diamond samples (few GHz) that is fairly described in Chapters 4.1 and 4.3 in Ref. 65. We find that the strong electron-phonon coupling occurs for the phonons with the energy in the order of 10 meV that is clearly manifested in the PL spectrum of SiV(−). This energy region is significantly larger than the spin-orbit energy, thus the strong electron-phonon coupling, i.e., the DJT effect is manifested via Ham reduction of the spin-orbit coupling. We note that the very low energy acoustic E g phonons might be considered as static strain which would treat these distortions by static Jahn-Teller effect instead of DJT. In any case, the final form of that Hamiltonian is the strain Hamiltonian, thus the two effects cannot be distinguished in experiments.
Based on these considerations, we simulate the Zeeman splitting of the corresponding states with a magnetic field aligned to 100 direction in Fig. 6 where we assume that no strain acts on the XV(−) color centers. In these simulations, we used ab initio electron-phonon deduced parameters but the intrinsic spin-orbit energies (λ 0 ) were scaled in order to reproduce the experimental zero-field-splittings (λ). We followed this procedure in order to directly compare our results to the experimental spectrum for SiV(−), GeV(−) and SnV(−). The g g,u L parameters were fitted to obtain the experimental spectrum of SiV(−) and we applied these Stevens' orbital reduction factors for the other XV(−) qubits. This is a simple approximation that could lead to a slightly underestimated values for the ground state of GeV(−), SnV(−) and PbV(−), where the d orbitals of the impurity atom may contribute to g g L . Nevertheless, our Calculated eigenvalues of the adjoint DJT (Eq. (4)) and SOC (Eq. (2)) interaction for the XV(−) color centers. The two lowest eigenvalues for each figure correspond to the 2 Eg,u vibrionic ground states for g ground (a-d) and the optically allowed u excited (e-h) states. All the considered eigenvalues are doubly degenerate in spin dimension because of the Kramers degeneracy. Each state consists of pure ↑ or ↓ spin state, thus the fourfold degeneracy of 2 Eg,u is fulfilled. We label the energy difference of the two lowest energy states by λ that is directly observed in the fine structure of the ZPL in the PL spectrum known as zero-field-splitting. Along x axis we depict the eigenvalues with respect to their partially quenched spin-orbit coupling strength LzSz , thus one can directly read out the p factors from this figure. δp shows deviation of the Ham reduction factors on E g,u3/2 and E g,u1/2 states. The larger the δp the less accurate is the treatment of SOC as a perturbation over DJT. We note that a mirror symmetry at x = 0 shows up for the ground state of SiV center in the entire vibronic spectrum that demonstrates that SOC can be treated as a perturbation over JT effect. The systematic left-shift at x axis for the vibronic spectrum of SnV and PbV defects implies that SOC is comparable with the JT coupling.
simulations should result in relatively accurate spectra. For SiV(−), our procedure resulted in g g L = 0.328 and g u L = 0.782. The strong reduction in the ground state can be understood by the shape of the corresponding orbital that we depict in Fig. 8 in Appendix D. For instance, E g+ state not only transforms as m l = +1 wavefunction but also as m l = −2 wavefunction. The linear combinations of the two leads to a significant reduction in the effective interaction with the external magnetic field, i.e., relatively small g L . In the E u excited state, the impurity p orbitals contribute to the interaction with the magnetic field unlike in the E g ground state, thus the reduction parameter is significantly larger for the excited state (see Appendix D). The list of all parameters can be found in Table IV in Appendix C. For SiV(−), our spin Hamiltonian can well reproduce the curvatures of the experimental Zeeman spectrum. We apply the same Hamiltonian for the other XV(−) qubits. Beside the obvious growing ZFS going from smaller to larger atomic number of impurity atoms, the general trend is that the curves are steeper for larger atomic number of impurity atoms because of the enhanced g zz values which is caused by the complex spin-orbit-phonon coupling (the fourth term in Eq. (9)).
For PbV(−) the calculated zero-field-splitting λ in the 2 E g ground state is about 4.4 THz or 18.2 meV which means that the 2 E g 3 /2 groundstate will be thermally filled with 100% and ≈96% occupation at cryogenic and liquid nitrogen temperature, respectively, that can be useful for quantum optics protocols. By applying the theory developed for the estimation of the coherence time of the SiV(−) (see Ref. 13 and the Supplementary Material in Ref. 34), we find that the decoherence process caused by the acoustic phonons is completely quenched at cryogenic temperatures for PbV(−) because of the large λ between the two branches of the 2 E g ground state, and the coherence time of the electron spin should be only limited by the nuclear spins or electron spins in the diamond crystal. That is much more practical for quantum communication applications compared to the millikelvin cooling needed for similar coherence time in SiV(−) qubit 14 . We note that the energy gap in SnV(−) should result in microsecond to millisecond coherence time for the electron spin going from 4 K to 1 K measurement temperature 34 .

IV. SUMMARY AND CONCLUSION
We performed a systematic study on the magnetooptical properties of Group-IV -Vacancy color centers in diamond by means of ab initio density functional theory calculations. We identified the photostability of these centers that can act as solid state qubits. We developed a novel spin Hamiltonian for these qubits in which the electron angular momentum and spin as well as the phonons are strongly coupled and identified such terms that have not been considered so far but are important in understanding their magneto-optical properties. We solved ab initio this complex problem for the model of these color centers consisting of up to 1000-atom supercells, and were able to reproduce previous experimental data. Furthermore, we identified SnV(−) and PbV(−) qubits with long spin coherence time at cryogenic temperatures where the spin state of PbV(−) can also be thermally initialized at these temperatures. Our ab initio toolkit and spin Hamiltonian analysis serve as a template for similar studies in 3D materials such as silicon carbide or 2D materials such as hexagonal boron nitride or transitional metal dichalgonides (TMD) or dioxides (TMO) which are fast emerging materials hosting qubits or single photon sources. In particular, the TMD and TMO materials exhibit strong spin-orbit couplings induced by the transition metal ions in the crystal in which strong mixing of spin-orbit and electron-phonon coupling are expected in the defects acting as qubits, and they should be treated at equal footing.
The vibronic wavefunctionΨ ± caused by electronphonon interaction can be expanded as written in Eq. (6) that mixes the electronic orbitals and |n, m e g phonons. The spin-orbit coupling then should be calculated for thẽ Ψ ± wavefunction, By expandingΨ ± in Eq. (B2) we arrive at FIG. 6. Magnetic field dependence of the ZPL splitting where the magnetic field is aligned along 100 direction. Zero energy is aligned to no spin-orbit and no-phonon coupling solution. The λ0 SOC coupling is fitted to obtain the experimental zerofield-splitting for SiV, GeV, and SnV in order to have a direct comparison to the experimental data whereas full ab initio result is shown for PbV. The actual parameters are listed in Table IV   This term gives negligible correction to the ZFS of GeV(−) and SiV(−) but becomes significant for PbV(−) as shown in Table IV.

Appendix D: Origin of Stevens' orbital reduction factor
The orbital angular momentum of an atomic orbital may be reduced in the potential created by surrounding ions that reduce the spherical symmetry 64 . We show that similar quenching can take place for vacancy-type defects in diamond. In the particular case of XV(−) color centers, the carbon dangling bonds in the vacancies form double degenerate E g± and E u± states. The E g,u± states were considered as m l = ±1 states 6,65 which led to the assumption of g g,u L = ±1. However, one should notice that the E g,u state will also transform as m l = ∓2 under D 3d symmetry. We illustrate this by plotting the e g,u± orbitals in comparison to the m l = ±1 and m l = ∓2 wavefunctions in Fig. 8. It can be observed that the e g,u± orbitals can be rather described as α|m l = ±1 + √ 1 − α 2 |m l = ∓2 where α is a coefficient of the m l = ±1 contribution. This results in a |g L | that is smaller than 1. Since the p orbitals of the Group-IV atom can contribute to m l = ±1 in the E u excited state, in contrast to the case of E g ground state, g g L is smaller than g u L in XV(−) qubits. We note that the full ab initio calculation of the expectation value ofL z is not straightforward. The implementation of spin-orbit coupling in vasp is based on the assumption that the spin-orbit interaction is predominant close to the core of the atoms and it is negligible in the interstitial regions, which is considered to be a well-justified approximation. In that case, the spin-orbit integrals can be distributed to the spherical volumes (PAW spheres) around the ions with using the atomic wavefunction projectors and potential. However, as Fig. 8 demonstrates, the interstitial regions can significantly contribute to the expectation value ofL z , which gives the interaction between the electron wavefunction and the external magnetic field. Therefore, we rather did not estimate this quantity from the values in the PAW sphere of the ions.   (9) for XV(−) systems. λ0 is scaled to reproduce the experimental λ for SiV(−), GeV(−) and SnV(−). The p, δp reduction parameters are calculated ab initio from the SOC and DJT coupled Hamiltonian (see Fig. 5 for visual interpretation). The orbital gL reduction factors were fit to reproduce the experimental curves for SiV(−), and the same gL factors were applied to the other XV(−) systems. The definition of the parameters shown here can be found in Eqs. (C8)-(C10