Spin-relaxation times exceeding seconds for color centers with strong spin-orbit coupling in SiC

Spin-active color centers in solids show good performance for quantum technologies. Several transition-metal defects in SiC offer compatibility with telecom and semiconductor industries. However, whether their strong spin-orbit coupling degrades their spin lifetimes is not clear. We show that a combination of a crystal-field with axial symmetry and spin-orbit coupling leads to a suppression of spin-lattice and spin-spin interactions, resulting in remarkably slow spin relaxation. Our optical measurements on an ensemble of Mo impurities in SiC show a spin lifetime $T_1$ of $2,4$ s at $2$ K.

Spin-active color centers in semiconductors have attracted significant interest for the implementation of quantum technologies, since several of these systems combine long-lived spin states with a bright optical interface [1][2][3][4]. Long distance spin entanglement has been achieved for a variety of defects as stationary nodes [5][6][7][8]. However, finding suitable emitters that combine long-lived spins, short excited-state lifetimes and optical transitions compatible with telecommunication fiberoptics infrastructure in an industrially established material has remained elusive. Silicon carbide, a wide bandgap semiconductor with mature fabrication technology, hosts a range of defect centers with optical transitions near or at the telecom range [9,10], including several defects containing transition metal (TM) impurities [11][12][13][14][15][16][17]. The electronic and spin properties of these defects derive largely from the character of the d-orbitals of the TM under the action of a crystal field determined by the lattice site [18][19][20]. Furthermore, the presence of a heavy atom in the defect implies that spin-orbit coupling (SOC) plays a significant role in the electronic structure of these color centers. Generally, the presence of strong SOC is expected to degrade spin properties by coupling spin and orbital degrees of freedom [17]. However, the extent to which the spin lifetimes of TM defects are limited by the strong influence of SOC is not clear.
We report here on slow spin relaxation with T 1 exceeding seconds below 4 K for a Mo defect ensemble in 6H-SiC, indicating that the defect spin is surprisingly robust with respect to spin relaxation despite the presence of strong SOC. In order to understand this, we measure the spin-relaxation time of the Mo defect in SiC between 2 and 7 K and identify the main processes leading to spin relaxation in this temperature range. We analyze the manifestation strength of these processes, while considering the electronic structure of the defect, and find that SOC is in fact responsible for suppressing several * These authors contributed equally to this work. spin-relaxation mechanisms in this system, leading to unexpectedly long T 1 .
We focus on the Mo defect associated with an optical transition line at 1121.3 nm in 6H-SiC, which consists of a Mo impurity in a Si substitutional site of quasi-hexagonal symmetry ( Fig. 1(a), h site) [13,15]. The defect is positively ionized, such that after binding to four neighboring carbons, the TM is left with one active unpaired electron in its 4-d shell. In this lattice site and only considering the rotational symmetry of the defect, both ground and excited states are two-fold orbitally degenerate. Due to this degeneracy, the orbital angular momentum is not quenched [20]. In the presence of spin-orbit coupling, the orbital degeneracy is broken, giving rise to two Kramer's doublets (KD) in the ground and two KD's in the excited state ( Fig. 1(b)) [13,15], resembling what is observed for the silicon-vacancy in diamond [21]. Concerning its electronic structure, this defect shares several similarities with the V defect in SiC, which has optical transitions fully compatible with telecommunications infrastructure [13,16,18] Each KD is a doublet composed of a time-reversal pair: the doublet splits as an effective spin-1/2 system in the presence of a magnetic field, but its degeneracy is otherwise protected by time-reversal symmetry (see [22] for a more extensive summary). The energy difference between the two spin-orbit split KD's in the ground state, ∆ orb in Fig. 1(b), is expected to be approximately 1 meV [15]. This is in accordance with the appearance of a second zero-phonon line (ZPL) in resonant photoluminescence excitation (PLE) experiments at approximately 8 K (Fig. 1(c)). Finally, the presence of sharp phonon replicas of the ZPL in the photoluminescence spectrum indicates that the defect center couples strongly to localized vibrational modes ( Fig. 1(b,d)).
We measure the spin relaxation time of this defect by means of a pump-probe experiment as shown in Fig. 2. Experiments are performed on an esemble of defects in a 6-H SiC sample previously investigated in Ref. [13], and further described in the SI [22]. We create pulses out of a CW laser beam by using a combination of an electrooptical phase modulator (EOM) and a Fabry-Pérot (FP) cavity. The EOM generates sidebands from our CW laser at frequency steps determined by an RF input signal; by tuning the FP cavity to transmit this sideband only, we create pulses that turn on/off as the RF generator turns on/off. These pulses resonantly drive optical transitions between ground and excited states, and we measure the photon emission into the phonon sideband with a singlephoton counter after filtering out the laser line. We apply a magnetic field non-collinear with the c-axis of the sample such that spin-flipping transitions between ground and optically excited states are allowed ( Fig. 2(a-b)) [13]. In order to counteract slow ionization of the defects (see [22, Sec. III] for further details), we apply a repump laser in between measurements [13,[23][24][25].
In the pump-probe experiment (Fig. 2(c)), the initial response of the sample to pulse P1 provides a measure of the population in the bright spin state |G1 ↓ at thermal equilibrium (see the caption of Fig. 2 for an explanation of the terms dark and bright spin sublevels in this work). The sharp increase (decrease) of the PL signal as the pulse turns on (off) indicates that both the optical decay rate and the Rabi frequency are relatively fast (see [22] for further details). In fact, an orbital decay rate of the order of a few MHz, in agreement with previous estimates [13], would explain the absence of a PL tail as the laser pulse turns off. Optical excitation provided by P1 polarizes the spin ensemble in a dark spin sublevel (|G1 ↑ in Fig. 2(b)) within the first few microseconds, as evidenced by the decrease and subsequent saturation of the PLE signal. The leading-edge response of the ensemble to a second (probe) pulse P2 reflects the recovery of the population in the bright spin sublevel (|G1 ↓ in Fig. 2(b)) during the time τ between the two optical pulses. Between 2 K and 7 K, we observe a monoexponential recovery of h 2 towards h 1 as a function of τ ( Fig. 2(d)), which must correspond to the spin relaxation time T 1 given the considerations presented below.
We repeat this experiment at zero magnetic field in order to confirm that the PL darkening observed within the first few µs of optical excitation corresponds to optical pumping of the ensemble into a dark spin sublevel. In this case, we observe no leading-edge peak in the photoluminescence, indicating that no optical pumping occurs in this timescale if the spin sublevels are degenerate [22, Fig. S2]. The absence of PL darkening at zero magnetic field implies that we cannot trap Mo centers in state |G2 for observable times. This is the case if the optical decay rate between states |E1 and |G2 is smaller than the relaxation rate between the two ground state KD's |G1 and |G2 in the temperature range investigated. (We do expect the optical decay rate between states |E1 and |G2 to be small, since we observe no lines corresponding to this transition in PL and PLE scans.) Thus, we conclude that after optical pumping in the presence of a magnetic field, no significant population is trapped in state |G2 , and that relaxation dynamics via this level does not affect the signals used to derive T 1 . Furthermore, we investigate the timescale associated with bleaching due to ionization [24,25], which is found to be several orders of magnitude slower than the one associated with spin dynamics [22, Sec. III].
The temperature dependence of the spin relaxation time spans several orders of magnitude, going from 2.4 s at 2 K to 83 µs at 7 K (Fig. 3). Concerning phononmediated spin relaxation, the mechanisms leading to spin relaxation are well established [20,26]. One (direct) and two-phonon (Raman, Orbach) processes are relevant when transitions between the levels involved are thermally accessible. Direct spin flip processes are expected to lead to spin relaxation rates that grow linearly with temperature (T −1 1,direct ∝ T ). In contrast, two phonon processes give rise to spin relaxation rates that grow superlinearly with temperature (T −1 1,Raman ∝ T 5<n<11 and where ∆ is the energy of the relevant excited state) [20,[26][27][28]. Additionally, interaction with paramagnetic moments in the material are expected to lead to temperature independent spin-relaxation. We can fit the data presented in Fig. 3 to a combination of these processes and identify the temperature regimes where they are relevant (see [22, Sec. VII] for additional details concerning the fitting procedure).
The defect has a rich electronic structure, with orbital and vibrational degrees of freedom at energies that are thermally accessible between 2 and 7 K ( Fig. 1(b)). Thus, we expect two-phonon processes involving real (Orbach) and virtual (Raman) transitions into both vibrational and orbital excited states to contribute to the temperature dependence of T −1 1 . Surprisingly, however, we find that not all available states contribute significantly to spin relaxation, leading to unexpectedly long spin lifetimes below 4 K.
Above 4 K, we identify an exponential growth of T −1 1 as a function of temperature, indicating the prevalence of spin-relaxation via Orbach processes in this temperature range. From the fit presented in Fig. 3, we extract ∆ = 7 ± 1 meV for the energy of the excited state involved in these spin flips. This energy matches the difference observed in PL experiments between the zero-phonon line and the onset of the phonon-sideband emission ( Fig. 1(d)). From the PL spectrum, we extract ∆ vib ≈ 10 meV for the energy of the first phonon replica ( Fig. 1(c)); nonetheless, the phonon replicas are significantly broadened, such that the first available vibrational levels are lower in energy. Thus, we identify vibronic levels where the spin is coupled to localized vibrational modes of the defect as the relevant excited state for twophonon mediated spin relaxation processes. As expected, these spin-relaxation mechanisms do not depend on the magnitude of the magnetic field [22]. We note that we cannot observe Orbach processes involving states |G2 . Between 2 and 7 K, the contribution of these processes to the spin relaxation is expected to be close to saturation, giving a flat temperature dependence in this range. Nonetheless, the long spin relaxation times observed at 2 K provide an upper bound to the spin relaxation rates due to two-phonon processes mediated by |G2 , and are evidence that these processes are slow. This leads us to conclude that transitions between these two KD's are strongly spin-conserving, similarly to what is observed for SiV color centers in diamond [29,30]. The presence of the additional orbital degree of freedom does not significantly degrade the spin relaxation time of this defect since hopping between states |G1 and |G2 preserves the spin state (but does contribute to decoherence [31]).
Below 4 K, the slow spin-relaxation rates observed must result from three mechanisms: i ) processes whose temperature dependence is saturated in this range (as treated above); ii ) direct one phonon transitions between the two spin sublevels |G1 ↓ and |G1 ↑ and iii ) temperature independent processes (such as spin-spin paramagnetic coupling). The observation of T 1 above seconds below 4 K is evidence that all three processes must be relatively slow for this particular system. The first contribution has been discussed above. Below, we elaborate on how the electronic structure of this defect leads to a suppression of the latter two processes.
The two spin-sublevels pertaining to a KD are strictly time-reversal symmetric with respect to each other, such that they must be degenerate eigenstates of any external field that preserves time-reversal symmetry. For this reason, electrostatic fields arising from the presence of phonons cannot, to first order, cause a direct spin-flip in a one-phonon type of transition between pure KD pairs [20, 22, Sec. IV]. A magnetic field or the interaction with nearby nuclear spins breaks time-reversal symmetry, and leads to mixing between spin sublevels pertaining to different KD. It is only due to this mixing that direct spinflips via interactions with single phonons can occur. This mixing is inversely proportional to the energy separation between the various KD, which is in turn largely determined by the spin-orbit splitting. Thus, a large spin-orbit coupling protects the KD character of the ground state spin doublet, suppressing direct spin-flipping processes.
Additionally spin-orbit coupling leads to a highly anisotropic Zeeman splitting of the ground-state spin sublevels [13,22], which hinders its interaction with a bath of paramagnetic impurities in the SiC crystal. Firstly, the spins are insensitive to magnetic fields perpendicular to their quantization axis, such that small fluctuations in the local magnetic field are not likely to induce a spin flip. Secondly, the spin doublet always has a Zeeman splitting governed by a g factor that is at maximum 1.6 [13], such that resonant spin flip-flop interactions with neighboring paramagnetic impurities with g ≈ 2 are largely suppressed.
The same arguments presented above to explain the long spin-relaxation times observed in this defect indicate that obtaining control of the defect spin via electric or magnetic microwave fields is a challenge. In first order, the ground state electronic spin is not allowed to couple to magnetic fields perpendicular to its quantization axis, or to electric fields (see [22, Sec-III] for a detailed grouptheoretical analysis supporting this). As a consequence, one cannot directly drive spin resonances between the two ground state spin sublevels. Mixing of the ground state KD with states lying at higher energies due to the presence of magnetic, hyperfine or electric fields can, however, lift these restrictions as is observed for V defects in SiC [16]. Although this is expected to provide a contribution to spin relaxation mechanisms, it is also a necessary requirement in order to control the ground state electronic spin via microwaves.
As discussed above, in the particular case of a TM at a quasihexagonal site, such as the Mo center, a combination of rotational symmetry and strong SOC protects the effective spin from flipping via two independent mechanisms. On the one hand, it gives rise to a unique Zeeman structure that is insensitive to perpendicular magnetic fields and non-commensurate with the behavior of other paramagnetic centers in the crystal. On the other hand, SOC generates a ground state effective spin which is an isolated KD, protected by time-reversal symmetry. In this doublet, the crystal field locks the orbital angular momentum along the axis of rotational symmetry of the defect. Via the strong spin orbit coupling, the electronic spin is stabilized, giving rise to robust effective spin states with long spin-relaxation times.
Considering these processes alone, stronger SOC is thus expected to lead to slower spin relaxation rates in TM color centers with an odd number of electrons and rotational symmetry. Nonetheless, our work shows that the presence of localized vibrational modes is pivotal in generating spin-flips, their energy determining the temperature where the onset of two-phonon relaxation mechanisms happens. Since the energy of the localized vibrational modes depends non-trivially on the reduced mass of the defect, whether or not defects containing heavier TMs (where SOC is more prevalent) will exhibit longer spin-relaxation times remains a question. Thus, it would be interesting to investigate these processes in defects containing 5-d electrons in SiC, where SOC is enhanced further. Tungsten defects, for example, have been observed with an optical transition at approximately 1240 nm and odd number of electrons, although their particular lattice site and charge state are not fully determined yet [32].
More generally, SOC should not be regarded as a detrimental feature when investigating solid-state defects for quantum communication applications. Transition metal defects in SiC offer interesting opportunities, such as charge state switching [33,34], emission in the nearinfrared [14] and long spin-lifetimes [17]. The maturity of the SiC-semiconductor industry means that a wide range of defects has been identified [18]. Nonetheless, their characterization with respect to optical and spin properties is still vastly unexplored. The first pulse P1 probes the population at thermal equilibrium in the ground state spin sublevel |G1 ↓ (which we call bright state in this work), subject to resonant excitation into the optically excited state. This pulse polarizes the spin population into the spin sublevel which is not subject to resonant optical excitations (which we here call dark state). After a variable delay τ , a second pulse P2 probes the recovery of the population of the bright state. (inset) The ratio of the leading edge peaks in the PLE signal (h2/h1) recovers monoexponentially as a function of the delay between the two pulses, such that we can extract T1. with temperature can be accurately described by a combination of direct (one phonon), Raman and Orbach processes (two-phonon processes). The two-phonon processes happen due to coupling to an excited state approximately 7 meV above the ground state KD, compatible with coupling of the electronic state of the defect to localized vibrational modes (see main text).
A. Henry, and E. Janzén, Optical identification and electronic configuration of tungsten in 4H and 6H-SiC, where δ is the baseline height, h 0 is the peak height at thermal equilibrium minus the baseline δ, and t is the delay between both pulses. This allows us to create a fit function with only two free parameters: T 1 and q, which represents the fraction h 0 h 0 +δ . c. Deadtime pile-up effect For the experiments in this work we used an avalanche photodiode single photon counter (SPC) for detection. After every detection event the detector has to reset its state, for which a predefined deadtime is used (10 µs in our case).
This creates a pile-up effect in our experiments when integrating many pulse sequences. We investigated how this causes an error in determining T 1 .
The probability of having a count at time t can be described by the combined probability that a photon is present and detected at the counter P photon multiplied by the probability of not having measured a photon before for the duration of the deadtime ∆t. We get the equation The outcome is plotted in Fig. S1 for a square PLE response to a pulse where P photon (t) = 100 kHz during the pulse, which saturates the SPC. Before the pulse, the detector is always ready to receive a photon, thus the initially measured count rate relates accurately to the number of photons present. For the duration of the deadtime (10 µs) the detector will be recovering and thus the count rate drops. This is periodic, but it clearly damps out to a steady-state count rate due to a spread in the actual detection times. This will give a positive error 1 ( 2 ) for the measured value of h 1 (h 2 ). At small t the magnitude of the error 2 is low and increases to 1 as t T 1 . If we rewrite Eq. 2 to represent the measured values it becomes where q is represents the corrected fraction h 0 h 0 +δ− 1 . The final term in this equation makes the measured h 2 h 1 approach its asymptote faster than in reality. Thus, fitting h 2 h 1 to Eq. 2 will always yield shorter T 1 values compared to the true spin-flip times. Simulating the pile-up effect for our measurements yields errors in T 1 varying from 2% to 7%. We choose to not correct for this and report the lower-bound values for the T 1 spin-flip times. Background PLE and dark counts have been neglected in this analysis as we found them to be of little influence so long as they remain below 10% of the peak count rate.

II. ZERO-FIELD MEASUREMENT
To confirm that the T 1 lifetimes measured are actually from the spin from the |G1 and not influenced by some other process, we performed zero-field experiments. Both groundstate spins are then degenerate and no spin pumping is expected. The results are shown in Fig. S2. The leading-edge peak in PLE as seen in Fig. 2(c) in the main text vanishes in Fig. S2(a). We note that thermal effects from laser driving at high intensities become more prominent at zero magnetic field, since darkening of the PLE no longer occurs. Additionally, we checked again the optical polarization dependence for these conditions at B = 0 mT, but this confirmed that the driving was still only sensitive to the component parallel to the c-axis.
In Fig. S3 we show the evolution of the PLE response with magnetic field strength.

III. CHARGE-STATE SWITCHING
When resonantly addressing optical transitions in the Mo defect in SiC, the PLE response drops over time. This is ascribed this to charge-state switching of the defect. We investigate the timescales of this bleaching in order to rule out its influence on the measurements of the T 1 spin-flip times.
The experimental approach is depicted in Fig. S4(a). First, a repump beam (770 nm, pulsed) counteracts any prior bleaching [1], resetting the charge state for 60 seconds. Next, a probe beam resonant with the ZPL illuminates the sample, slowly bleaching the Mo defects.
We can track the bleaching timescale by measuring the PLE response as function of time.
After another 1000 seconds the repump beam is incident on the sample together with the probe beam, which allows us to track the recovery timescale. Finally, the repump is switched off for 60 seconds to check the initial decay of the PLE response. This sequence is repeated four times. We use a 2 mW repump beam and a 200 µW probe focused to approximately 100 µm diameter in the sample. The magnetic field strength is 100 mT at an angle of 57 • with the c-axis.
The results are shown in Fig. S5. The bleaching between 60 and 1060 seconds occurs according to two timescales, both are fit with an exponential decay for the orange and blue curve. The yellow curve is an exponential fit to the recovery by the repump laser.
All three timescales are plotted versus temperature in Fig. S4(b). For the T 1 experiments, where the repump beam was only on in between measurement runs, the two bleaching scales are most relevant. Both occur at rates that are at least one order of magnitude slower than the observed spin-flip times. Any bleaching occurring at faster timescales should have been visible in the zero-field measurements from section II. Thus, the effect of bleaching on measuring T 1 can be deemed negligible.
Note the fast decay of PLE for 4 K (Fig. S5(a)) after the repump laser is blocked at 1300 seconds. We ascribe this to fast spin flips induced by the repump laser. When performing the experiments to measure T 1 with an omnipresent repump laser, we observe that T 1 is reduced by an order of magnitude (for similar laser powers). This fast decay is not visible in Fig. S5(b-f) since the spin-flip times are already quite short at higher temperatures.  seconds). When the repump is off, but the system is being renonantly driven, the PLE signal decays exponentially with two characteristics time constants (blue and red fits). As the repump is turned on, the PLE signal recovers quickly, with a single time constant (yellow fit).

IV. ELECTRONIC STRUCTURE -GROUP THEORETICAL APPROACH
The Mo defect is characterized by a single active spin in the 4d shell of a molybdenum impurity at a Si substitutional lattice site [1]. In this configuration, the Hamiltonian of the defect has a three-fold rotational symmetry and three vertical mirror planes, such that the eigenfunctions of the electronic orbital state transform according to the symmetry group C 3v .
Despite the extensive literature available on the effect of crystal field and spin-orbit coupling on the electronic states of these defects [2,3], we are unaware of a comprehensive report on the effect of symmetry on the coupling terms between the various specific spin sublevels and external fields based on the double group representations of the defect symmetry, and we present this analysis here. We apply a group-theoretical approach to obtain the symmetry of the eigenfunctions associated with a defect center in a crystal field of C 3v symmetry in the presence of spin-orbit coupling. Furthermore, we obtain and explain the selection rules governing the interaction of the electronic spin with magnetic and electric fields.
The group theoretical rules governing the selection rules presented in the end of this work do not rely on a particular basis set for the description of the electronic wavefunction. By this, we mean that even if we consider hybridization of the wavefunction of the bare transition metal atom/ion with the nearest carbon atoms due to covalent bonding, the symmetry of the crystal field Hamiltonian is preserved such that the new, modified wavefunctions will still obey the selection rules arising from a group-theoretical analysis. Nonetheless, it is instructive to start from an analysis of the effect of the Hamiltonian on the 10 spin-orbital states arising from a single electron sitting in one of the d-orbitals of the transition metal.
The transition metal at a silicon substitutional site shows tetrahedral coordination due to bonding to the 4 nearest carbons. In this configuration, the 5 d orbitals split into an orbital doublet and an orbital triplet (which transform as the irreps E and T 2 of the symmetry group T d ), where the triplet lies highest in energy. Due to the hexagonal character of the lattice, the tetrahedral symmetry of these sites is lowered to C 3v , with the rotational axis aligned parallel to the growth axis of the crystal. Upon this symmetry reduction, the triplet further splits giving rise to an orbital doublet and a singlet, which transform respectively as E and A 2 . The effect of this symmetry lowering operation is expected to be largest in the lattice sites of quasi-hexagonal symmetry (h), and to only modestly affect the lattice sites of quasi-cubic symmetry (k) (Fig. S6(a)).
A wavefunction transforming as a non-degenerate irrep of a given point-group cannot have an effective orbital angular momentum (in other words, the orbital angular momentum is quenched) [4]. However, this requirement is lifted in the presence of degeneracies, such that the eigenfunctions of the Hamiltonian transforming as E are allowed to have a non-zero orbital angular momentum. Thus, in order to fully describe our system, we must consider the effect of spin-orbit coupling. In a group-theoretical approach, this is done by extending the group of interest to include 2π rotations which bring a spin ↑ into a spin ↓ [5]. That is, this is done by considering the eigenfunctions as basis states of the irreps of the double group associated with the C 3v group, here denoted byC 3v .
In the double group including the effect of spin-orbit coupling, three irreps describing how odd spin wavefunctions transform are added to the group. These irreps are Γ 4 , which is doubly degenerate, and Γ 5,6 , two irreps that are connected by time-reversal symmetry and must thus be degenerate in the presence of time-reversal symmetry. The orbital singlet transforming as A 2 gives rise to a Kramer's doublet (KD) transforming as Γ 4 , whereas an orbital doublet transforming as E splits into two KDs, of which one transforms as irrep Γ 4 , and the other transforms as irreps Γ 5,6 . Thus, the symmetries mentioned above split the 10 states arising from an electronic configuration 2 D into 5 Kramer's doublets, of which 2 transform as Γ 5,6 , and 3 transform as Γ 4 (Fig. S6(a)). The character table of the double groupC 3v is given in Tab. SI. Additionally, we explicitly show what are the transformation properties of the vectors x, y, z and the axial vectors R x , R y , R z , as well as how the cubic harmonics z 2 , x 2 − y 2 , xy, xz, yz transform under the operations of the group.
We can investigate the role of small magnetic and electric fields in driving transitions between different KD (coupling between different KDs), and spin resonances (coupling between the two eigenstates pertaining to a single KD) in the framework of group-theory, given that these fields are small enough that the symmetries of the Hamiltonian H 0 are preserved.
The selection rules between two wavefunctions can be obtained in a straight-forward way. If |ψ i and |φ i are two eigenstates of the Hamiltonian H 0 transforming respectively as irreps Γ i , Γ i , the selection rules with respect to a perturbative Hamiltonian H are given by the product ψ i | H |φ i . In order for this matrix element to be non-zero, it must transform as a scalar, that is, as the totally symmetric irrep A 1 [5]. Thus, the product of the represen- where the perturbation H transforms as Γ j and * denotes complex conjugation, must contain the totally symmetric irrep A 1 . Table SII gives the decomposition of the various products of Γ 4 , Γ 5,6 , in terms of irreps of the C 3v group, and translates this into the selection rules governing the coupling between various spin states.
Optical transitions between various sets of KDs are allowed due to coupling to E , E ⊥ , which belong to irreps A 1 and E, respectively. We can extract polarization selection rules from table SII. Electric field driven transitions between two KDs transforming as Γ 5,6 will be polarized along the symmetry axis of the defect; transitions between two KDs transforming as Γ 4 can be polarized along any direction; transitions between a KD transforming as Γ 5,6 and a KD transforming as Γ 4 are only allowed for light polarized perpendicular to the symmetry axis. These properties are summarized in Fig. S6(b). This means that only electric or magnetic fields of E symmetry (that is, in the xy plane) are capable of coupling and mixing states |G1 (which transforms as Γ 5,6 ) and |G2 (which transforms as Γ 4 ).
Transitions and energy splittings within each of the KDs can also be understood based on the symmetry of the defect. The anisotropic Zeeman structure observed for the ground state spin doublet, which is insensitive to magnetic fields perpendicular to the crystal symmetry axis [1] can be understood purely based on the properties of the group. A magnetic field along the symmetry axis of the defect transforms as R z , whereas a magnetic field perpendicular to this axis transforms as R x , R y . Within a doublet which transforms as Γ 5,6 , no coupling is allowed with a magnetic field perpendicular to the symmetry axis since Γ * 5,6 ⊗ E ⊗ Γ 5,6 = E ⊃ A 1 . This is not the case for a doublet transforming as Γ 4 , such that the spin sublevels that transform as Γ 4 are allowed to couple to magnetic fields in the plane, and will not have g ⊥ = 0. Thus, we conclude that the ground state doublet belongs to the irrep Γ 5,6 ( Fig. S6(c)). As long as the quantization axis of the defect spin (states |G1 , belonging to irrep Γ 5,6 ) points parallel to the symmetry axis of the defect, we cannot rotate the spin via microwave spin resonances, since these spins are insensitive to magnetic or electric fields perpendicular to this axis.
Finally, we note that if two spin sublevels are strictly connected by time-reversal symmetry (that is, they are a pure KD), they cannot be connected by operators that preserve time-reversal symmetry. This was proven by Kramer and became what is known as Kramer's theorem [5]. Thus, within a pure KD, electric fields are not capable of driving transitions between the two spin-sublevels.

V. SIMULATION OF RAW DATA
Due to the large number of available states for defect (vibrational levels, orbital state |G2 in Fig.1 of the main text, ionized states), it is not straight forward to obtain quantitative information from the shape of the raw data plots presented in the main text. Nonetheless, we can apply a rate equation model to reproduce the data and, upon carefully taken assumptions, obtain a bound for the values of the optical decay time and Rabi driving frequency in our experiments.
In order to minimize the set of free parameters and facilitate the analysis of the behavior of the system, we simulate this defect center as a three-level system ( Fig. S7(a)), where the ground (state 1) and excited (state 3) states can be coupled by an optical field with Raby frequency Ω R . From the optically excited state, the system can decay either back into the ground state with a rate Γ 31 , or into a shelving state 2 with a rate Γ 32 . Additionally, population can be transferred between states 2 and 1 at a rate Γ 21 , and between 1 and 2 at a rate Γ 12 = e −∆/kT Γ 21 , where ∆ is the energy difference between 2 and 1 and kT denotes the thermal energy of the system. We simulate the system with a simple set of rate equations for the populations of each state (P 1 , P 2 , P 3 ), without treating coherences explicitly. Finally, we consider that the photoluminescence observed is proportional to the population of the optically excited state, P 3 .
We try to reproduce the typical shape of the raw PL data obtained experimentally ( Fig. S7(b)) in order to obtain a set of reasonable values for the Rabi driving frequency and optical decay rates in our system. We assume that states 1 and 2 correspond to the ground state spin sublevels, |G1 ↓ and |G1 ↑ respectively. In this way, we can write Γ 32 in terms of Γ 31 by assuming that the branching ratios correspond to the overlap of the spin states in ground and optically excited state [1].  (Fig. 3).
We note that differences based on the exact value of Γ 31 are barely noticeable. Thus, we cannot restict our estimate for Γ 31 further. Nonetheless, we can restrict the expected values of Ω R by comparing the calculated traces presented in Fig. S7(c,d) and the raw trace presented in Fig. S7(b). We not that if Ω R is very small, of the order of a kHz, PL darkening is almost absent, unlike what is seen in experiment, where PL darkening is significant. In contrast, if Ω R is of the order of a few MHz, the defect darkens completely within the time of the driving pulse. This is also in disagreement with the experimental data. Thus, we conclude that the Rabi frequency in our experiements is of a few tens to hundreds of kHz.
Furthermore, section II shows that we do not see any PL darkening when we perform the time-resolved measurements described in the main text at zero magnetic field. In this case, state 2 in our model corresponds to the orbital state |G2 from the main text. We calculate the population in state 2 after optically driving the system for approximately 500 ms, with Rabi frequencies of the order of a few tens of kHz, and present these results in Fig. S8. We only transfer significant population into 2 (leading to PL darkening) when the optical decay rate into state 2, Γ 32 is larger than the rate at which the system leaves state 2, Γ 21 . Since we do not observe any PL darkening, we conclude that Γ 32 Γ 21 such that, within the time of our measurements, no significant population is transferred into the orbital state |G 2 , and the presence of this state does not influence our measured value for the spin T 1 .

VS TEMPERATURE
The spin-lattice relaxation of single spins of substitutional defects in solid state materials arises from a modulation of the crystal field potential in time due to the presence of phonons, which perturbs the stationary crystal field (V (0) ) and couples various eigenstates of the time-independent Hamiltonian to each other [4]. Thus, the probability of a spin flip to occur depends largely on the matrix elements of the time-dependent crystal field V (1) between the various electronic levels accessible to the defect. We thus define the terms V orb and V vib , which indicate the order of magnitude of the matrix elements of V (1) connect- Fig.1 in main text for definitions).
The direct process, a one phonon interaction driving transitions between states |G1 ↓ and |G1 ↑ directly, is expected to show a temperature dependence of the kind T −1 where ω is the Zeeman splitting between the spin sublevels |G1 ↓ and |G1 ↑ . Since the states |G1 ↓ and |G1 ↑ are each other's time-reversal pair and V (1) preserves time-reversal symmetry, the matrix elements | G1 ↓| V (1) |G1 ↑ | are identically zero (see section IV). Nonetheless, the presence of a magnetic field or hyperfine interaction perturbs states |G1 ↓ and |G1 ↑ by mixing in states higher in energy, in such a way that we expect the direct process to be present with a magnitude roughly proportional to ( ω) 4 |V (1) | 2 ∆ 2 T , where V (1) is now the matrix element of the time-dependent crystal field coupling states |G1 ↓ , |G1 ↑ to a generic excited state |E lying an energy ∆ above |G1 ↓ , |G1 ↑ . All of the excited states shown in Fig. 1(b) are expected to contribute to this process, such that mixing with both the higher KD (|G2 ↓ , |G2 ↑ ) and the vibronic states (|G1 vib ↓ , |G1 vib ↑ ) should be considered. In this way, we expect a dependence of the kind T −1 Additionally, a spin polarization in the defect can decay back to its equilibrium value via two-phonon processes comprising transitions into real (Orbach process) or virtual (Raman process) excited states. The former gives rise to an exponential temperature dependence of the type T −1 1 ∝ | G1 ↓| V (1) |E E| V (1) |G1 ↑ |∆ 3 exp(−∆/k B T ) in the limit of ∆ k B T , where ∆ is the energy difference between a generic excited state |E and the KD |G1 ↓ , |G1 ↑ . Orbach processes relative to transitions into states |G2 ↓ , |G2 ↑ are expected to give rise to a strong temperature dependence at temperatures below 1 K and saturate at higher temperatures, when ∆ orb ∼ k B T , and its exponential behavior is thus not visible in our data. In contrast, Orbach processes relative to transitions into vibronic levels are expected to contribute significantly to the temperature dependence of T −1 1 at a few K, since ∆ vib ∼ 10 meV k B T between 2 and 8 K (Fig. 1(d)) .Thus, we expect the Orbach process to give rise to a temperature dependence of the kind T −1 1 ∝ |V vib | 2 ∆ 3 vib exp(−∆ vib /k B T ). Finally, second order Raman processes give rise to a temperature dependence of the kind The parameter Γ 0 is included to account for temperature independent processes of spin relaxation. The fit quality does not improve significantly if we consider n = 9 instead of n = 5 for the Raman process involved. Thus, we are unable to determine which levels are involved in the Raman transitions responsible for spin relaxation between 3 and 4 K.
In either case, however, the exponential increase of the spin relaxation rate above 4 K is accounted for by an Orbach process where two phonons drive transitions between the ground state |G1 and the its vibrational excited state |G1 vib flipping its spin. From the fit, we extract ∆ vib ≈ 7 meV, consistent with the energies of the phonon-coupled states responsible for the PSB emission in the photoluminescence spectrum. Additionally, from the fitting parameters C R and C O , we get values of V (1) vib with a consistent order of magnitude of a few tens of meV.
Finally, we note that the data can also be fit by a power law model of the type T −1 1 = αT + βT γ , with γ ≈ 13. Spin-lattice relaxation of this type has been previously reported for heavy ions in solid state environments [6]. In that work, a Raman process is observed with a power dependence with γ ≈ 11 > 9. They justify the large power observed by noting that the spin sublevels in the KD are not exactly each other's time-reversal conjugate, in such a way that the 'Van Vleck' cancellation does not happen completely [4]. We exclude this as a relevant model in our case due to the fact that the power dependence necessary to explain our data is much higher, with γ ≈ 13. Additionally, the consistency of ∆ vib observed by fitting the data with the energies observed in the PSB of the PL spectrum indicates that the rapid increase of the relaxation rate observed above 4 K is indeed related to exponential Orbach processes involving |G1 vib ↓ , |G1 vib ↑ .