Strong spin–orbit quenching via the product Jahn– Teller effect in neutral group IV qubits in diamond

Arti ﬁ cial atom qubits in diamond have emerged as leading candidates for a range of solid-state quantum systems, from quantum sensors to repeater nodes in memory-enhanced quantum communication. Inversion-symmetric group IV vacancy centers, comprised of Si, Ge, Sn, and Pb dopants, hold particular promise as their neutrally charged electronic con ﬁ guration results in a ground-state spin triplet, enabling long spin coherence above cryogenic temperatures. However, despite the tremendous interest in these defects, a theoretical understanding of the electronic and spin structure of these centers remains elusive. In this context, we predict the ground-state and excited-state properties of the neutral group IV color centers from ﬁ rst principles. We capture the product Jahn – Teller effect found in the excited state manifold to second order in electron – phonon coupling, and present a nonperturbative treatment of the effect of spin – orbit coupling. Importantly, we ﬁ nd that spin – orbit splitting is strongly quenched due to the dominant Jahn – Teller effect, with the lowest optically-active 3 E u state weakly split into m s -resolved states. The predicted complex vibronic spectra of the neutral group IV color centers are essential for their experimental identi ﬁ cation and have key implications for use of these systems in quantum information science. e g pJT 2nd order in Here we E ð i JT for the result of constructive ( i = 1) and destructive ( i = 2) interference of the two orbital branches. The axial asymmetry arises from a second order effect denoted similarly by the parameter i JT . The distortion space, which allows us to fully i


INTRODUCTION
Artificial atoms in diamond are promising candidates for a wide variety of quantum technologies [1][2][3][4][5] , including as quantum repeaters for long-range quantum networks 6,7 . Many milestones have been reached using the nitrogen-vacancy (NV − ) center 8,9 and more recently the SiV −10-14 . Further exploration of novel defect candidates has included the GeV −15-18 , SnV −19-22 , PbV −23,24 , and SiV 025-27 , all of which have been observed experimentally and described theoretically 28,29 . The neutrallycharged SiV 0 has symmetry analogous to the SiV − , but its missing electron gives rise to a triplet ground state as found in the NV − , with the corresponding potential for both long spin coherence times and symmetry-protected optical transitions. Theoretical work has postulated the remaining group IV neutral (IV 0 ) centers 29 (GeV 0 , SnV 0 , and PbV 0 ) and described the negatively-charged group III defect centers 30 as isoelectronic to the SiV 0 . Calculations suggest that all of these defect candidates are thermodynamically more likely to exist in intrinsic diamond than the SiV 0 , which requires p-type doping 27 . Within this growing space of candidate artificial atom qubits, an ab initio understanding of the level structure is required to harness the advantages of each emitter in quantum science 31 .
Accurate descriptions of artificial atoms in diamond can be particularly challenging because of the dominant Jahn-Teller (JT) distortions 32 present. In such systems, the total energy of a JTunstable electronic configuration is lowered as a result of the coupling of the electronic structure to nuclear motion, introducing electron-phonon interactions. In the case of group IV 0 defects, the excited state exhibits a product Jahn-Teller (pJT) effect which results from simultaneous Jahn-Teller instabilities in two orbitals 29,[33][34][35] . The pJT interaction leads to either a dynamical or static JT effect, or a mixture of both. In the case of a dynamical JT distortion, the system is best described as a collective electron-vibration (vibronic) system. This strong coupling of electronic and vibrational states may modify electronic observables, for example a quenching of spin-orbit (SO) coupling (SOC).
Including the pJT effect is therefore critical for predictions of the zero-phonon line (ZPL) transition energies and the excited-state level structure. Previous work has found that describing pJT interactions to first order in coupling explains the observed energy splitting 25 between the optically-bright E u and dark A 2u states for SiV 029 . An important effect to consider, particularly for the heavier group IV 0 defects, is the role of spin-orbit interactions, as these defects can have coupling constants on the order of 100s of meV 28 . The interplay of SOC interactions and JT physics in the excited-state of group IV 0 centers has significant impact on the expected SO behavior if the JT effect couples the electrons and phonons strongly, as we find.
In this article, we describe the combined impact of spin-orbit and Jahn-Teller interactions in the neutral group IV centers in diamond from first principles. We describe the product Jahn-Teller effect to second order in electron-phonon coupling and find a large second order energy shift. Importantly, the effects of spin-orbit coupling are included nonperturbatively and splittings are found to be an order of magnitude smaller than expected for a purely electronic system as a result of the JT interaction. These fine structure details reveal additional physics of color center qubits in diamond and present a pathway to identify GeV 0 , SnV 0 , and PbV 0 experimentally.

Electronic structure
The group IV centers in diamond adopt a split-vacancy configuration within the diamond lattice, where the dopant group IV atom sits between two vacant carbon sites, as shown in Fig. 1a   1 and denoted by the point group D 3d . The defect introduces localized electronic orbitals comprised of the dangling bonds of the nearby carbon atoms and the dopant atom, which can be captured using density functional theory (DFT) 36 (see "Methods" section for computational details) and are labeled by their symmetry. The energetically-relevant orbitals are of e u and e g character and exist near and above the valence band of intrinsic diamond, respectively, shown schematically in Fig. 1b. Both the e u and e g orbitals are doubly-degenerate and can be further labeled by their spatial orientation, i.e., {e u } = {e ux , e uy } and similarly {e g } = {e gx , e gy }. Including spin, these levels combined can host up to eight electrons. For group IV 0 centers, six electrons are present in the (e u e g ) manifold. Equivalently, we can describe these electronic states in the basis of two defect-bound holes. We choose to adopt this convention for the remainder of this article.
The ground state has the hole configuration e 2 g (e 1 gx e 1 gy ), and prefers the triplet S = 1 spin configuration. The total defect wavefunction is of 3 A 2g symmetry, and is directly obtained from electronic structure calculations. In constructing the total wavefunction, given the symmetric triplet spin component, we ensure that the orbital wavefunction is antisymmetrized; this is given by the A symbol. The ground state orbital wavefunction can be written as A e gx e gy ¼ 1= ffiffi ffi 2 p e gx ðr 1 Þe gy ðr 2 Þ À e gx ðr 2 Þe gy ðr 1 Þ À Á . In the excited electronic configuration, one hole moves from an e g to an e u orbital. Unlike in the ground state, there exist four distinct hole occupations with this e 1 g e 1 u configuration. The antisymmetrized orbital wavefunctions are given by A e ux e gx , A e uy e gx , A e ux e gy , and A e uy e gy . We can construct the irreducible representations of the triplet subspace as linear combinations of these orbital states, as has been done previously 29 .
Jahn-Teller interaction Each of these antisymmetrized states obtained from our ab initio calculations are Jahn-Teller unstable, in that they energetically prefer a configuration with the lower symmetry C 2h point group to that with the higher symmetry D 3d point group. We note that despite the lower symmetry, the C 2h point group is still inversionsymmetric. The nuclear motion associated with these distortions is a result of interactions with phonon modes of symmetry E g . In contrast with the single JT system (E g ⊗ e), the JT distortion found in the excited state of group IV 0 systems is due to simultaneous JT interactions in both the e u and e g orbitals. This collective product Jahn-Teller behavior is denoted by E g ⊗ e u ⊗ e g and shown schematically in the right panel of Fig. 1b. Previous work has covered the single JT to second order as well as the pJT 32,33,35 to first order in electron-phonon coupling. Here, we describe the coupling of the two electronic states with the E g -type vibrational mode to second order in vibrational coupling. The Hamiltonian for Fig. 1 Product Jahn-Teller effect in the excited state of group IV vacancy centers in diamond. a Lattice configuration of the group IV 0 defects, in which the impurity atom (blue) sits between two vacant carbon sites (gray). b Simplified energy level diagram showing the energy location of the doubly-degenerate e u and e g orbitals relative to the band gap of bulk diamond. The ground state is a spin triplet and the corresponding excited state undergoes a symmetry-breaking pJT distortion (right) as a result of orbital instabilities in both the e u and e g orbitals. c Potential energy surfaces computed for the pJT system including effects up to 2nd order in coupling. Here we label the energy instability by E ðiÞ JT for the result of constructive (i = 1) and destructive (i = 2) interference of the two orbital branches. The axial asymmetry arises from a second order effect denoted similarly by the parameter δ ðiÞ JT . The black curves indicate 1D cuts through the 2D (Q x , Q y ) distortion space, which allows us to fully parameterize the system. d DFT-obtained potential energy surfaces along these 1D cuts for the SnV 0 defect. The D 3d high-symmetry point (Q x = 0 Å) is found to be unstable in two surfaces, consistent with the pJT picture. We also label the displacement amplitudes ρ ðiÞ 0 from the D 3d to the C 2h minima. The splitting Λ is a result of static electronic correlation. All values are tabulated in Table 1. this interaction can be written as: The first two lines represent linear coupling with coupling constants F u/g while the latter two represent quadratic coupling terms with coupling constants G u/g for both the e g and e u orbital branches. The nuclear component of the Hamiltonian is written withX andŶ representing bosonic operators for the phonons given by fX;Ŷg ¼ ðâ y fx;yg þâ fx;yg Þ= ffiffi ffi 2 p and the electronic component in terms ofσ i which are the standard Pauli and unit matrices acting on the e u ⊗ e g subspace. The Hamiltonian in Eq. (1) is defined within the single-excitation two-particle hole manifold, therefore the basis states are A e ux e gx , A e uy e gx , A e ux e gy , and A e uy e gy , which are captured from electronic structure calculations. For Eq. (1), we note that there can exist different sign conventions; in this work we adopt the sign convention commonly applied to Jahn-Teller Hamiltonians 28,32,37,38 .
In the pJT case, two independent solutions which are unstable at the high-symmetry point can exist. One corresponds to the constructive interference of the two JT distortions ð$ðF g þ F u Þ 2 Þ and the other to the destructive interference ð$ðF g À F u Þ 2 Þ, as shown in Fig. 1c. To find the coupling constants and solve for the coupled vibronic states, we obtain displacement ρ ðiÞ 0 and energy E ðiÞ JT , δ ðiÞ JT parameters from the defect potential energy surfaces (PES) computed from first principles electronic structure, where i = 1, 2 for the constructive and destructive pJT, respectively. For the SnV 0 color center we show the resulting adiabatic PES as a onedimensional cut along Q y = 0 in Fig. 1d (see Supplementary Fig. 1 for similar plot including all group IV 0 defects). In principle the PES are two-dimensional, with the minima being threefold degenerate (see Fig. 1c). However, due to the symmetry of the PES, this 1D cut completely parameterizes the pJT Hamiltonian. For additional details on connecting the coupling constants in Eq. (1) to our calculations, refer to Supplementary Note 2 and Supplementary  Table 1.
In these defect systems electronic correlationŴ plays a role in splitting the electronic states for reasons distinct from the Jahn-Teller physics. This correlation can be incorporated along the lines of previous work 29 , leading to the following total Hamiltonian for the system: pJT þŴ: pJT . The termŴ represents the effects of static electronic correlation, which causes the energy splitting Λ between the two branches of the potential energy surface shown in Fig. 1d.
Spin-orbit interaction Next we describe spin-orbit interactions in the pJT system. In the presence of a dynamical JT effect, expectation values of purely electronic operators can be quenched because of the coupled vibronic nature of the system, as first shown by Ham 39 . Thus it is important to analyze the effects of SO interactions with caution, as has already been demonstrated for the group IV − defects 28 . In these group IV 0 centers, the SOC Hamiltonian can be written as a product of the single-hole interactions 37 , since the spin-orbit coupling does not mix the e u and e g orbitals 40 . The SOC Hamiltonian is written aŝ Here, we introduce SO splittings λ 0 u=g for both the e u and e g orbitals, which can be obtained from ab initio calculations. The variable m s corresponds to eigenvalues ofŜ z and for the S = 1 triplet system can take on values of m s ∈ [1, 0, −1]. While SOC in its generalL ÁŜ form (with angular momentum operatorL and spin operatorŜ) also contains transverse terms, these transverse terms only couple e g /e u orbitals to a 2u orbitals which are outside the (e g , e u ) manifold of interest 40 . This consideration allows us to effectively writeĤ SOC solely in terms proportional toL zŜz , yielding Eq. (3). TheL zŜz interactions can couple the excited-state singlet manifold with the m s = 0 excited-state triplets, however we choose to consider only the triplet subspace as the ðe 1 u e 1 g Þ singlet excited states are expected to be higher in energy due to Coulomb repulsion 26 . Ultimately intersystem crossing (ISC) rates between these triplet and singlet levels will likely depend on phonon overlaps of the full diamond + defect system, however they require nonzero spin-orbit coupling and thus our analysis is important for further understanding ISC.
To capture the spin-orbit interaction in addition to the pJT physics, we find that including SOC perturbatively is insufficient, even for the SiV 0 system. Thus, we invoke a complete spinresolved orbital basis including all spin sublevels of Eq. (3). From this we perform direct diagonalization of the combined spin-orbit and Jahn-Teller system (see Supplementary Note 5), where we take all terms in Eq. (2) to be spin-independent. The solutions of this coupled Hamiltonian allows us to extract both the absolute energy shifts of our vibronic eigenstates with SO effects and the effective SO splittings between spin sublevels nonperturbatively. Table 1 summarizes the results of our work. In each of the defect centers studied, we find a significant pJT effect, with the constructive interference yielding instabilities of over 200 meV. We find the second order effects are also relatively large, with δ ð1Þ JT $ 0:3E ð1Þ JT for each of the defects studied. These second order shifts are important, as they represent the energy barrier between the three energy minima present in the 2D vibrational (Q x , Q y ) space. This energy barrier helps to determine if the system will prefer a static or dynamic JT distortion, the latter of which means the electron and phonon degrees of freedom cannot be decoupled and instead a coupled vibronic solution is required. Indeed, the system can be parameterized as strongly-coupled as given by the parameter λ = E JT / hω E , which is >2 for all cases studied here. After calculation of the parameters in Table 1, we can solve for the coupled electron-vibrational system as defined in Eqs. (1,2). plotted as a function of the expectation value of displacement from the high-symmetry D 3d minima. In both cases, we can project the solutions onto the irreducible states of the D 3d excited-state manifold, with the color legend given in panel c. We find that the lowest energy states are comprised of roughly equal contributions from the undistorted j 3 E u i and j 3 A 2u i electronic states. This is true for the quadratic coupling as well. In Fig. 2d we specifically focus on the lowest-energy vibronic solutions. The lowest vibronic state has total symmetry A 2u which is optically dark, whereas the next eigenstate is an optically-active, doubly-degenerate E u level. In first order pJT, the splitting γ (1) between these two states for SnV 0 is 8.96 meV, while including second-order coupling decreases splitting γ (2) to just 6.22 meV. Even at second order, the 3 E u state remains degenerate, however overall the eigenstates of the system shift upwards in energy by roughly 20 meV.

DISCUSSION
It is interesting to note that in general including second-order terms in the pJT Hamiltonian decreases the splitting γ between the lowest vibronic states (see Table 1). This splitting was measured experimentally for SiV 025 to be 6.8 meV; here we find a larger discrepancy to experiment in the case of quadratic coupling (γ (2) = 3.2 meV) than we do for linear coupling (γ (1) = 7.2 meV). We emphasize, however, that an inclusion of second order electron-phonon coupling more closely resembles the ab initio data, as can be seen in Fig. 1d due to the nonvanishing δ ðiÞ JT . The origin of this disagreement is unknown and beyond the scope of this work. We suggest that it may represent an energy-resolution limitation in the approach employed. Possible routes to improve the agreement with experiment, which will be analyzed in future work could be an improved model of the static correlation or the use of a different hybrid functional than HSE06. We note that inclusion of higher-order terms 41 up to fourth order in electron-phonon interactions is found to negligibly change our results, as the corresponding coupling coefficients for these terms are found to be orders of magnitude smaller than the first and second order terms.
The coupled spin-vibronic results are shown in the final panel of Fig. 2d and are found after including the SOC Hamiltonian directly. We find that the m s = 0 and m s = ±1 sublevels of the A 2u vibronic states are split, in the case of SnV 0 by 5.9 meV. The m s = ±1 sublevels of the E u states also split (here we distinguish the E u states by labels + and −). These E ± u states have a Kramers degeneracy, very much analogous to the lowest E g vibronic states of the group IV − , where j 3 E þ u i "" j i and j 3 E À u i ## j i are the degenerate, lowest energy E u states. These are split by an energy of λ u + λ g from the degenerate j 3 E À u i "" j i and j 3 E þ u i ## j i states, as shown in Fig. 2d for SnV 0 . In the absence of JT interactions this splitting λ g + λ u would be over 100 meV, however here it is only~3 meV, a direct consequence of the strong electron-phonon coupling present in the pJT system. Additional interactions such as effects of strain and spin-spin coupling could split and shift these levels further.
For all cases, the reduction factors denoted by p u/g are smaller than 0.05, indicative of a very strong quenching of the SO interaction, even more so than the group IV − color centers. This can be attributed in part to the scaling of the Jahn-Teller instability vs. the spin-orbit splitting in the two-hole case. While to first order the JT energy scales as the square of the coupling (i.e., $ ðF u þ F g Þ 2 ), the SO splitting scales linearly (i.e., λ g + λ u ). Such a scaling and the resulting JT energies intuitively explains the significant SO quenching we find in this work. We note that shifts in the absolute energies of the E u states are found to be most significant in the case of PbV 0 , where we find a redshift in the predicted ZPL of roughly 0.05 eV. All lighter defects have much weaker absolute energy shifts due to their reduced SO coupling constants. More details and discussion of the SOC approach used here can be found in Supplementary Note 5, including Supplementary Fig. 6 and Supplementary Table 2.
We note that additional theoretical work remains to be done towards understanding these emitters, including the nature of the photoluminescence spectra 42 , intersystem crossings, and the effects of nonradiative decay channels.
In conclusion, we present first principles calculations of group IV neutral artificial atoms in diamond, where we capture the product Jahn-Teller effect to second order in electron-phonon coupling and nonperturbatively describe the effects of spin-orbit interactions. Our results find significant reduction in the spin-orbit splitting due to the strong pJT. However, we also find that the spin-orbit interactions would split the lowest optically-active states into m s -resolved levels split by up to a few meV in the heavier candidates. These results provide qualitative insight into the physics of artificial atom qubits in diamond, and are of quantitative importance in experimental identification and manipulation of these centers in quantum information science.

First principles calculations
We employ constrained Kohn-Sham density functional theory calculations performed within the VASP code 36 (version 5.4.4). The Kohn-Sham wavefunctions are described using a plane wave basis set with projector-augmented wave (PAW) pseudopotentials. The basis set used corresponds to an energy cutoff of 400 eV. We also performed calculations at stringent 800 eV energy cutoffs and found no appreciable changes in the results. We employ the hybrid HSE06 exchange-correlation functional 43,44 to accurately describe the energetics of these defect systems. All calculations were performed with spin polarization. The defects are modeled in a cubic supercell equivalent to 512 carbon atoms of the diamond lattice, with a lattice constant of 6.70 Bohrs (3.545 Å). For such a supercell we sample only the Γ point of the corresponding Brillouin zone.
Ionic relaxation was performed until forces on all atoms fell below 10 −2 eV/Å. The excited-state geometries were obtained for the high-symmetry D 3d geometry by enforcing symmetry during relaxation. We also relaxed with a C 2h symmetry constraint to capture the low-symmetry structures and energy. We also verified these were indeed the lowest-energy structures by removing symmetry constraints and finding no change in ionic positions or total energy.

Solving of Hamiltonians
All subsequent calculations and solving of Hamiltonians are done numerically. We include up to 40 phonons in the expansion of the product Jahn-Teller Hamiltonian (Eq. (1) of the main text). Eigenvalues and eigenvectors are found by diagonalizing the corresponding Hamiltonian. We determine the parameters ρ ðiÞ 0 , E ðiÞ JT , δ ðiÞ JT and Λ directly from the DFT potential energy surface (e.g., Fig. 1d). The effective vibrational energy hω E can be found from these parameters similarly to the case of the single Jahn-Teller (see Supplementary Note 2). The vibronic splitting between the lowest levels to first and second order are given by γ (1) and γ (2) , respectively. SO effects are included nonperturbatively and we find significant quenching of the pure electronic SO splitting (p u,g ≪ 1), a consequence of the strong electron-phonon coupling induced by the pJT. The energy λ u + λ g corresponds to the energy splitting between the m s = ±1 levels of the lowest E u vibronic eigenstates.  x þ Q 2 y q relative to the D 3d minima. The solutions for both a and b are projected onto the D 3d symmetry-adapted electronic states and the resulting composition is represented by the color shown in c. d The effects of 2nd order JT and explicit inclusion of SOC are detailed for the lowest-energy eigenstates of the system. In 1st and 2nd order JT, the A 2u state is nondegenerate and the E u state is twice degenerate. The inclusion of second order decreases the splitting γ between these levels, while also introducing an absolute energy shift of around 20 meV. The inclusion of SOC splits the E u levels into E þ u and E À u , each with corresponding m s sublevels. The splitting between the m s = ±1 levels is given by λ g + λ u , which is strongly attenuated. The m s = 0 (labeled by jSj "#i) levels are unaffected by SOC.