Orbital and Spin Dynamics of Single Neutrally-Charged Nitrogen-Vacancy Centers in Diamond

The neutral charge state plays an important role in quantum information and sensing applications based on nitrogen-vacancy centers. However, the orbital and spin dynamics remain unexplored. Here, we use resonant excitation of single centers to directly reveal the fine structure, enabling selective addressing of spin-orbit states. Through pump-probe experiments, we find the orbital relaxation time (430ns at 4.7K) and measure its temperature-dependence up to 11.8K. Finally we reveal the spin relaxation time (1.5s), and realize projective high-fidelity single-shot readout of the spin state ($\geq98\%$).

Defect centers in solids are a promising class of systems for quantum science and technology [1,2]. They combine bright optical transitions, access to long-lived electronicand nuclear-spin registers and compatibility with solidstate device engineering. Of particular prominence is the negatively-charged nitrogen-vacancy center (NV − ) in diamond, which has enabled recent advances in quantum information science [3,4] and quantum sensing [5][6][7].
Alongside NV − , the nitrogen-vacancy defect can exist in both the neutral-(NV 0 ) and -with sufficient Fermilevel engineering -positive-(NV + ) charge states. These additional charge states can be used as a resource in a number of applications, such as spin-to-charge conversion for improved spin-state read-out [8,9], classical data storage in NV ensembles [10], and deliberate charge-state switching for improved nuclear-spin coherence under ambient conditions [11,12].
Conversely, for experiments based upon NV − , undesired conversion to NV 0 can be a hindrance: active charge-state initialization protocols have been used to counter this [13,14]. For quantum networks, stochastic conversion from NV − to NV 0 is an important decoherence mechanism for nuclear-spin quantum memories [15].
Despite the importance of NV 0 , understanding of many of its properties remains elusive. In particular, the orbital-and spin-dynamic timescales are unknown. Also, while recent magnetic circular dichroism (MCD) measurements on ensembles [16,17] give insight into the NV 0 fine structure, no direct observation has been reported. Building an understanding of the system and its associated dynamic processes is important for improving control in NV quantum devices. Moreover, the knowledge gained may offer new insights into the physics of other impurities in solids [18]. Finally, NV 0 may prove to be a powerful quantum system in its own right.
Here, we develop protocols combining resonant excitation of both NV 0 and NV − . We apply these novel pro- * These authors contributed equally to this work. tocols to reveal the orbital and spin dynamics of single NV 0 centers in diamond, as well as to realize initialization and single-shot readout of the NV 0 spin state. We perform our measurements on single NV centers at cryogenic temperatures [19], see Fig. 1(a). The NV center is adressed with microwave (mw) pulses (NV − ground-state spin transitions) as well as with polarization-controlled λ red = 637 nm (NV − zero-phonon line (ZPL)) and λ yellow = 575 nm (NV 0 ZPL) laser light. We apply an axial magnetic field of B z = 1890(5) G to induce significant Zeeman splitting.
The ZPL of the NV 0 center has been conclusively attributed to this defect [20][21][22][23][24][25]. A combination of abinitio calculations and symmetry arguments led to the proposal of ground states of 2 E symmetry, which can be optically excited to a 2 A 2 manifold [26,27]. An additional metastable 4 A 2 quartet state was also predicted, and has been observed by electron paramagnetic resonance (EPR) measurements under excitation of the NV 0 ZPL [28]. A splitting of the transitions of the two orbital states E x and E y has been measured [27,29]. However, the associated fine structure has not been observed in PL or EPR measurements.
We start by performing spectroscopy using the experimental procedure sketched in Fig. 1(b). For each frequency step, we (1) probabilistically prepare the emitter in NV 0 by applying strong laser excitation resonant with the NV − ZPL, in combination with weak mw driving [19] to induce the conversion NV − → NV 0 . We then (2) apply polarized yellow light, during which time all single-photons above 650 nm are integrated.
The measured spectra ( Fig. 1(c)) show four transitions -the first direct spectroscopic observation of the NV 0 fine structure. These observations validate the model of Barson et al. [16], and we hence follow their theoretical description below. Under the secular approximation, the ground-state Hamiltonian of NV 0 can be described by H = gµ BŜz B z + lµ BLz B z + 2λL zŜz + ⊥ (L − +L + ). g is the spin g-factor, µ B is the Bohr magneton, l is the orbital g-factor, λ is the spin-orbit interaction parameter and ⊥ is the perpendicular strain parameter.L z = σ z andŜ z = 1 2 σ z are the orbital and spin operators defined in terms of the Pauli matrix σ z , whileL ± = |± ∓| with |± = ∓(1/ √ 2(|X ± i |Y )) are the orbital operators defined within the basis of the strain eigenstates {|X , |Y }. The z-axis is defined parallel to the NV axis.
The resulting level structure is presented in Fig. 1(d). The 2 E ground state is composed of a pair of doublet states with opposite spin-orbit parity (lower spinorbit branch: {|+, ↓ , |−, ↑ }; upper spin-orbit branch: {|−, ↓ , |+, ↑ }). The degeneracy of each doublet is lifted by orbital-and spin-Zeeman contributions under the applied magnetic field. Conversely, the 2 A 2 excited state exhibits no spin-orbit structure, but is rather split by the spin-Zeeman effect alone. These contributions lead to four spin-conserving transitions. The contributing ground state for each observed transition is indicated in Fig. 1(c).
We find that the luminescence of the transitions depend significantly on the polarization of the excitation light (see Fig. 1(c)). Differing transition amplitudes for orthogonal polarizations can be attributed to optical selection rules that are strongly dependent on ⊥ [16,19]. Based upon these observations, we develop a method to extract ⊥ and simultaneously the fine structure param- . By fitting spectra from three individual NV centers against our theoretical model, we find l = 0.039 (11) and λ = 4.9(4) GHz. These values are roughly a factor of 2 larger than those found previously using NV-ensemble MCD measurements [30]. Crucially, the data in Fig. 1(c) shows that resonant optical excitation in this magnetic field regime allows for state-resolved addressing, enabling the heralded preparation of specific states and investigation of the system dynamics. To date, only the excited-state lifetime, τ exc , of 21 ns has been reported [31]. Here, we investigate the orbital-and spin-relaxation timescales of the ground state, τ orbit and τ spin , see Fig. 1(d).
In order to unambiguously measure the dynamics of NV 0 , we design and implement a charge-resonance (CR) protocol that realizes high-fidelity heralded preparation into NV 0 , with the λ = 575 nm laser resonant with a chosen optical transition, see Fig. 2(a). The CR protocol (1) can be broken down as follows. First, a heralding signal confirms preparation in NV − , with the λ = 637 nm lasers on resonance with the NV − transitions. Next, a strong red optical pulse induces charge state conversion, after which a chosen NV 0 transition is excited with yellow light. If the photon counts obtained during the 'NV 0 check' exceed a pre-set threshold, the protocol is completed. Further details are given in the supplementary information [19]. After the CR protocol, we perform the experimental sequence on NV 0 (2). Finally, we detect whether undesired conversion to NV − occurred during the experimental sequence, and then perform read-out of the NV 0 state (3). Note that the CR protocol prepares a specific spin state of the NV 0 center. For circular polarisation we typically start the experiment by heralding the |↓ spin state. For linear polarisation, however, due to their close spectral vicinity, the CR check heralds either the |↓ or |↑ spin state.
In Fig. 2(b) we show time-resolved pump measurements. Here, the yellow laser is gated by an acoustooptic modulator (AOM), with measured rise-/fall-time of 30(5)/7(1) ns. Upon opening the AOM, we observe a rapid increase in fluorescence due to optical cycling, which is then damped as population is pumped out of the driven state. By fitting the steady-state fluorescence counts for L/H polarization we extract a saturation power of 2.5(2)/1.8(1) nW and saturation counts of 105(2)/103(2) kcts/s, see Fig. 2(c). As the optical power is increased, coherent optical Rabi oscillations are observed. In Fig. 2(d), we plot the fitted frequency of these oscillations, revealing the expected P yellow dependence. When the AOM is closed the fluorescence decays with τ exc = 22(1) ns (inset Fig. 2(b)), which is consistent with literature [31].
To uncover the recovery timescale after pumping we turn to pump-probe spectroscopy. Example time-traces are shown in Fig. 3(a). The resulting data is well described by an exponential recovery with a single timescale associated with how fast the system relaxes [19] once illumination is turned off. At base temperature of our cryostat (T = 4.65(3) K), we extract τ recovery = 0.43(6) µs. We attribute these fast dynamics to orbital relaxation processes, i.e. |+ ↔ |− and τ orbit = τ recovery .
We repeat the pump-probe measurements across a range of temperatures. The fitted recovery times are shown as rates R recovery = 1/t recovery in Fig. 3(b). After an initial linear increase a rapid increase is observed at higher temperatures. At these higher temperatures, the required time resolution exceeds the AOM switching time constants, which we take into account in the fitting procedure [19].
The initial linear increase (∝ T ) can be attributed to single-phonon processes, while high-order processes appear to govern the recovery rate at higher temperatures [32,33]. Here, we fit individually to a two-phonon Raman process (∝ T n ) and a two-phonon Orbach process (∝ exp[−∆/k B T ]), with k B being the Boltzmann constant. For the Raman process the fit returns n = 13(2) (14(3)) for the lower (upper) spin-orbit branch; a physical explanation for such values is currently lacking. For the Orbach process we find a characteristic energy scale of ∆ = 12(2) meV (∆ = 13(4) meV) extracted from a fit to the lower (upper) spin-orbit branch. ∆ is associated to the energy splitting to the first vibronic level of the NV 0 ground state, predicted to be a Jahn-Teller system [20,34]. The value found here agrees with the bulk absorption measurements of Davies [20] (13.6(7) meV), and with recent density-functional theory calculations (21.4 meV) [34], suggesting that the measured increase of R recovery is predominantly due to two-phonon Orbach processes. While a detailed model is beyond the scope of this work, we expect that our findings will aid in the further understanding of the vibronic structure of NV 0 . Now we turn to the spin dynamics of NV 0 . Here, we exploit polarization control to selectively prepare, address, and readout the NV 0 spin state. These measurements are all performed on timescales τ orbit = 0.43(6) µs and thus average over the orbital basis; we will therefore only refer to the spin states. In all experiments below, we use L polarisation, addressing the |↓ state. We herald the preparation of |↓ by applying 25 nW for 250 µs, and proceed when more than 25 photons are detected. After a delay of 0.1 ms, we perform a charge-state check with red excitation, followed by a second yellow readout, see Fig. 2(a)(3). We then repeat this experiment, but with a delay of 10 s between the yellow readouts, allowing for relaxation processes to occur. The resulting histograms are shown in Fig. 4(a).
In the first case (dark colors), we observe a single dominant population which can be modelled by a Poissonian distribution with mean photon count 25.2 (2), and that we attribute to |↓ . In the second case (light colors), we additionally observe a second distribution, again modelled as a Poissonian distribution with mean pho- ton count 0.171 (4). A charge-state measurement of NV − performed before each read-out shows that only a small fraction of the population (P NV − ∼ 1%) is found in the unwanted charge state -which we discard from the histograms -and that the majority of low-count events can be attributed to a dark state of NV 0 . As the populations evolve without laser excitation, the dark state must be part of the ground state manifold; we therefore assign this state to the second spin state |↑ . A read-out threshold of 5 photons (solid line, Fig. 4(a)) discriminates the two spin states. We now sweep the delay time between initialisation and read-out. The measured populations of |↓ (P ↓ ) and |↑ (P ↑ ) are plotted in Fig. 4(b), showing relaxation to a mean population of 0.494 (6). The data is consistent with a spin-1/2 T 1 process of characteristic timescale τ spin = 1.51(1) s. Note that the observed value is a lower bound of the intrinsic spin relaxation, as it may be limited by leakage of resonant laser light. By setting the initial and long-time population in |↓ to be 1 and 0.5 respectively, we obtain a lower-bound for the single-shot read-out fidelity, F RO = 1 2 (F |↓ +F |↑ ) ≥ 98.2(9)%, where F |s is the probability to assign |s after preparing |s [19].
To investigate the cycling nature of the driven optical transition we now repeat the measurement under 5 nW of resonant yellow excitation, see Fig. 4(c). We find that P ↓ decreases on a timescale faster than can be explained by spin-relaxation alone, showing that the optical excitation induces spin pumping. Possible spin-mixing channels are given either in the 2 A 2 excited state or via an intersystem crossing, which might be offered by the 4 A 2 state. We also find a significant increase of P NV − due to opticallyinduced charge conversion [19,35]. However, this slows once |↓ is depleted as |↑ is a dark state for optical excitation. Beyond this, P ↑ reduces with τ spin , and charge conversion continues. We find a high state preparation fidelity for |↑ of 99 +1 −10 % after 600 ms, but with an absolute population in the NV 0 |↑ state of only 22(2) %.
To reveal the respective rates we develop a three-levelrate equation model that we fit to our data, using the measured spin-relaxation time as a fixed input (solid lines, Fig. 4(c)) [19]. For the applied power of 5 nW, we extract characteristic timescales of 27(1) ms (90(4) ms) for the charge conversion (spin pumping) process. From this we can estimate the cyclicity of the |↓ state within this regime to be 0.98(8) × 10 5 cycles, mainly limited by In a second experiment the 5 nW yellow excitation is stroboscopically interleaved with strong NV − → NV 0 ionisation pulses [19], see Fig. 4(d). Again we observe a gradual decrease of P ↓ , and an increase of both P ↑ and P NV − , but then P NV − growth stops and even inverts. This observation can be explained via the picture that the removal of an electron from NV − prepares a random spin-state in NV 0 , eventually populating the dark state |↑ . Competing rates between this spin-selection process and spin relaxation lead to the observed steady state populations. We again fit a three-level rate equation model, using the previously obtained parameters as fixed inputs [19], and extract a timescale for ionisation of 18(4) ms. The rate equation model does not accurately describe the behaviour at long timescales, which is likely due to a reduction of the NV 0 spin-relaxation time under red excitation and strong NV − microwave driving [19].
As a final step, we develop a master equation simulation to capture the full dynamics of the NV 0 center [19]. In Fig. 2(b) we plot the simulated excited state population (solid line), using the uncovered NV 0 timescales and spectral properties. We match the Rabi frequency to the measured optical power and further include a spectral average over a Gaussian distribution of detuning values with FWHM = 2π × 20 MHz. We find excellent agreement with our experimental fluorescence data, emphasizing a consistent understanding of the NV 0 dynamics.
In conclusion, we have developed a novel toolbox for the study and control of single neutrally-charged NV centers in diamond. We have uncovered the dynamic timescales and demonstrated single-shot readout and initialization-by-measurement of the NV 0 spin, each with high fidelity. In future investigations, coherent control of the spin states may be obtained. Detailed modelling of the defect may give new insights into the observed temperature-dependence of the orbital dynamics. On the application side, protection of nuclear spin quantum memories from dephasing by NV 0 may be achieved by microwave spin locking in both orbitals, or by feedback based upon the NV 0 spin read-out demonstrated here. Finally, at reduced temperatures that suppress the orbital dynamics, NV 0 may prove to be a powerful system for quantum technologies in its own right.
We thank Michael Barson Our experiments are performed on single nitrogenvacancy (NV) centers in type-IIa bulk diamond (Element Six, CVD grown, <111> oriented), using a cryogenic (Montana Cryostation, 4 K) home-built confocal microscope setup. Enhanced photon-collection efficiency is achieved by fabrication of solid immersion lenses [37] and an anti-reflection coating [38]. For phonon-sideband (PSB) detection, a dichroic mirror (Semrock, pass above 650 nm) and an additional long-pass filter are used to block reflections of the excitation lasers. Photon emission is detected via an avalanche photo-diode (APD, Laser components, quantum efficiency ∼ 80%), with a total collection efficiency of ∼ 3% (∼ 10%) of the NV 0 (NV − ) PSB.
We apply a magnetic field of B z = 1890(5) G along the symmetry axis z of the NV center via a permanent magnet. A slight misalignment of the field leaves a small perpendicular magnetic field component of B ⊥ = 10(5) G. As the ratio B ⊥ /B z is small, we neglect the effect of the perpendicular magnetic field. In addition to the magnetic field, local strain and electric fields can alter the NV center level-structure [39]. For the NV center used in the main text we observe the level structures as depicted in Fig. S1.
To address the NV 0 charge state we apply resonant optical excitation of the NV 0 zero-phonon-line (ZPL) (λ = 575.17 nm, ω = 2π × 521.22 THz). The laser frequency can be manually tuned to each 2 E to 2 A 2 transition. In the NV − charge state, selective excitation of ZPL transitions (λ = 637.25 nm, ω = 2π×470.45 THz) enables optical readout (RO, can be coherently adressed with microwave (mw) pulses delivered via gold strip lines on top of the diamond surface. For the NV used in the main text (NV A), we extract an NV − perpendicular strain of NV − = 4.2(1) GHz, from the observed optical transition frequencies, see Fig. S1. We note that it is unclear how this relates to the strain in NV 0 , as the susceptibility of the NV 0 states to electric fields is currently unknown.
FIG. S1. Level structure for the two charge states NV 0 and NV − of a single NV center. The optical transitions used within this work are indicated by yellow (red) solid arrows for NV 0 (NV − ). Charge-state switching between the two charge states is achieved via two-photon absorption of the respective ZPL laser [35].
Further, charge state conversion may result in differing local charge environments for the two charge states [40].

Photo-luminescence measurements
For the photo-luminescence (PL) measurements presented in Fig. 1 of the main text, the following methodology was used. First, the NV center is prepared in the neutral charge state by the application of a 50 µs optical SP pulse (1 µW) in combination with weak mw driving (ν mw Rabi ∼ 150 kHz) of the NV − m s = 0 − → ← − m s = −1 ground-state spin transition. This method prevents optical spin-pumping into a NV − dark state (m s = 0). Second, 500 pW of yellow light (P yellow ) is applied for 15 µs at the NV 0 ZPL transition, during which all single-photon detection is integrated. The red-yellow procedure is repeated N = 150 times before the frequency is stepped by 1 MHz.
A total of 20 full scans were made for each polarization setting, which were collated to produce the final PL data. For each PL scan, the fluorescence maxima are found via a peak finding routine (python, scipy.signal.find peaks cwt). Further, all data sets are shifted to the mean frequency of all fluorescence maxima and summed. We typically observe shifts of the maxima by up to 200 MHz due to spectral diffusion. PL scans of the lower and upper spin-orbit branches were done in two separate measurements. To avoid systematic shifts of the splitting between the two spin-orbit branches, the NV center is reset by strong green (λ = 515 nm, 10 µW) illumination in between these two measurements, cancelling potentially accumulating effects of spectral diffusion from the red-yellow scans.
In a second set of experiments, we study the linewidth of the observed optical transitions for various values of P yellow , see Fig. S2(a). We observe a broadening of the lines with increasing power. As a result, in the high power regime the fine structure is no longer resolved. The extracted full width at half maximum (FWHM) is plotted in Fig. S2(b) as a function of P yellow (i.e. ∝ optical Rabi frequency). For power broadening, a linear dependence on the Rabi frequency is expected, while at low powers the FWHM is limited by the intrinsic linewidth of the defect. In Fig. S2(b) we fit the dependence under the assumption that the linewidth can be described by a convolution of a lifetime-limited Lorentzian profile (f L = 1/(2π × τ exc ) = 7.6 MHz) with Gaussian broadening terms. The resulting Voigt FWHM can be approximated by f ≈ 0.5446f L + 0.2166f 2 L + f 2 G [41]. The Gaussian component is given by the convolution of a power-dependent term from power-broadening, and a power-independent term arising from spectral diffusion: f G = a 2 P yellow + b 2 . The fit shows good agreement with the observed behaviour. We find a spectraldiffusion-limited linewidth of 30.3(3) MHz, a factor of 4 above the transform limit.
Extraction of the NV 0 fine structure parameters The obtained PL measurements carry information of the parameters of the ground and excited state Hamiltonians. Based on these measurements we develop a methodology to extract the NV 0 orbital g-factor l and spin-orbit interaction parameter λ.
Beside the fine-structure constants, the PL spectrum of the NV 0 center depends on both magnetic field and strain (stress within the crystal and electric fields). These dependencies are captured within the NV 0 Hamiltonian of the 2 E ground state, see Ref. [16]: Here, g is the spin g-factor, µ B is the Bohr magneton,L z = σ z andŜ z = 1 2 σ z are the orbital and spin operators defined in terms of the Pauli matrix σ z . The last term of the Hamiltonian shows the influence of perpendicular strain ⊥ , whereL ± = |± ∓| with |± = ∓(1/ √ 2(|X ± i |Y )) and {|X , |Y } is the basis for the strain eigenstates. Note that parallel strain is not included as it does not affect the relative energy of the ground state levels. For the excited state 2 A 2 , the Hamilonian reads as and does not show a dependency on strain. These two Hamiltonians lead to the energy level structure as presented in Fig. S1. The four resulting transition frequencies are shown in the PL spectra in Fig. S3(a). The corresponding eigenstates in the ground state are indicated. Depending on the polarisation (circular right (R) or left (L); linear horizontal (H) or vertical (V)), the amplitude of the observed PL varies. From Voigt fits to the individual PL lines we extract the transition frequencies and PL amplitudes. From the transition frequencies, we determine the energy splitting, ∆ spin , between the two spin-states, |↑ and |↓ , associated with each spin-orbit branch, and the energy splitting ∆ spin-orbit between the two spin-orbit branches.
Further, the PL amplitudes can be directly related to the transition strength for a given polarization and transition. In the absence of strain, circularly-polarized transitions are expected, as one quantum of orbital angular momentum has to be transferred upon excitation. Accordingly, under such excitation, full PL contrast would be expected between the transitions within each spinorbit branch ('orbit contrast'). However, under large strain, the ground state is better described within the strain eigenbasis [16], with associated linearly-polarised transitions. In this scenario, full contrast would instead be expected between the spin-orbit branches ('spin-orbit contrast'), whilst no 'orbit contrast' would be expected within each branch. As a consequence the observed contrasts can be used to determine the strain.
To extract the contrast, we repeat PL scans of a spinorbit branch for several angles of the quarter-wave plate (QWP) and half-wave plate (HWP). In the case of circular polarisation, we normalize each PL scan by the respective integrated total counts. The fitted amplitudes A circ  for all measured angles, A circ,norm For linear polarisation, for each angle we take the mean amplitude of the transitions within a spin-orbit branch, A lin = (A lin 1 + A lin 2 )/2. We then normalize each angle by the mean for all measured angles, A lin,norm = A lin /A lin . Figure S3(b,c) shows example plots for the resulting normalized amplitudes for varying circular and linear polarizations. Measurements on the respective other spin-orbit branch give similar results.
As a next step we extract the transition contrast for circular polarizations (Fig. S3(b), 'orbit contrast') and linear polarizations (Fig. S3(c), 'spin-orbit contrast') from the amplitudes of fits with a sine function. These contrasts are plotted in Fig. S3(d) for three independently measured NV centers, NV A, NV B, and NV C. Further, we plot the mean transition energy splittings ∆ spin-orbit and ∆ spin in Fig. S3(e) and (f).
These three data sets are now fitted against the Hamiltonian with l, λ, and ⊥ for each NV as free fitting parameters. Our results are summarized in Tab. I (method 1). The fitted strain results are used to place the experimental data points in Fig. S3(d-f), see dashed lines in Fig. S3(d). Solid lines show the result of our fit. The transition energy splittings ∆ spin and ∆ spin-orbit are well described by our theoretical model. While a good match of the transition contrast for NV A is found, NV B and C show a discrepancy. A possible explanation could be non ideal polarization settings at the position of the NV for these data sets, resulting in less clean rotation around the Poincaré sphere, i.e. mixed circular and linear polarizations when rotating the QWP/HWP. A mixed polarization leads to both reduced 'orbit contrast' and reduced 'spin-orbit contrast'.  As a second method we repeat our fitting procedure, but this time fixing the individual strain values to NV − of each NV, obtained from NV − , see Fig. S1. While it is not known how the NV − strain translates to strain in NV 0 , a correlation is expected. From the fit we obtain l and λ, see Tab. I (method 2). Within error the two methods give the same values. In the main text we report the mean of the two values: l mean = 0.039 (11), λ mean = 4.9(4) GHz.
We now compare our data to the transition frequencies calculated using the fine structure parameters of Barson et al. [16], l lit = 0.0186(5), λ lit = 2.24(5) GHz, (dashed lines, Fig. S3(e,f)). Strikingly, with these parameters, our data for ∆ spin cannot be reproduced for any strain value, strongly indicating that the discrepancy in fine structure parameters cannot be explained by systematic errors in our method to extract the strain of the NV.
We note that Barson et al. have used NV ensemble magnetic-circular dichroism measurements while we here observe PL of single NV centers.

Charge-resonance check
As described in the main text, one of the key components of the experiments in this work is the introduction of a charge-resonance (CR) check for NV 0 . This check allows heralded preparation of the NV in the neutral charge state, while also preparing the red lasers on resonance with the NV − optical transitions, and the yellow laser on resonance with one of the four NV 0 transitions. The full procedure is shown in Fig. S4 and outlined in the following.
We first prepare the negatively-charged state. A strong green pulse (12 µW) is applied for 300 µs in order to prepare the system ('reset'). Next, we simultaneously apply a combination of the red RO (1 nW) and SP (3 nW) light for a duration of 70 µs, during which time we integrate all single-photon counts incident on an avalanche photo-diode (APD) ('check NV − '). If a set photon-count threshold is exceeded, we have high confidence that the red lasers are well on resonance with their associated transitions and that the NV is in NV − , and proceed to the next step. If the count is below the threshold, but above zero, then it is assumed that the NV is in the negative charge state and close to resonance, but not yet in a satisfactory regime. In this case, the red check is repeated until the threshold is passed. In the case that any red check produces a zero photon-count, the green pulse ('reset') is reapplied to reset the charge state or to induce significant spectral diffusion bringing the NV − transitions back in resonance with the red lasers. We note that the low powers used for resonant excitation itself cause minimal spectral diffusion.
After the NV − check, an ionization pulse is applied to prepare NV 0 ('ionise'). Here, we apply 5 nW (10 nW) of RO (SP) light respectively, for a total duration of 1 ms. While the ionisation probability is low (∼ 2%), the chosen powers ensure that spectral diffusion is minimal.
To verify that the NV has been successfully transferred to the neutral charge state and to confirm that the yellow light is on resonance with a single transition we apply 25 nW of yellow light for 250 µs ('check NV 0 '). In this check, we either exceed the threshold, in which case we proceed to the main experiment, or we return to the 'check NV − ' step. The counts of each successful 'check NV 0 ' step are saved.
In our experiments, we use a single yellow laser only. We note that polarization of this laser affects the spinstate prepared after the CR check. When exciting with linear polarization, one of the two spin states is prepared in each experimental repetition. Which spin state is prepared may vary due to spectral diffusion between repetitions in combination with the close spectral vicinity between the spin states (see Fig. S3(a)). However, when exciting with circular polarization, a single NV 0 spinstate is selectively addressed and heralded throughout: the probability to false-herald a non-targeted transition is negligible. While the laser frequency corresponds to a transition associated with a specific spin-orbit state, the check heralds a mixed orbital state as it takes significantly longer than the orbital relaxation time. We note that a general (not spin-selective) resonance check with higher efficiency could be achieved by adding a second laser to address a NV 0 transition corresponding to the opposite spin-state.
Finally, after the experiment has been completed, we perform the 'check NV − after' step, which can be used to detect transfer from NV 0 to NV − during the experiment. This is enabled by the fact that transfer between the NV charge states induces minimal spectral diffusion, as witnessed by our recharging data (see corresponding supplementary section).
We then return to the 'check NV 0 ' step due to two reasons. First, the detected counts can be used to readout the NV 0 spin state in the case recharging has not happened. Second, for most experimental repetitions we remain in NV 0 , leaving a high probability that the 'check NV 0 ' step is passed again. In such cases, it is not necessary to repeat the 'check NV − ' step, reducing measurement overhead time.
We now show an example case for dynamics under optical pumping (see Fig. 2 of the main text). In the experiment, we herald a specific spin-orbit state via the 'check NV0' step, simultaneously bringing the yellow laser on resonance with the respective transition. As an example, we can herald the |+, ↓ state. However, due to orbital relaxation dynamics being much faster then the duration of the state check we effectively herald an orbitally mixed |↓ state, i.e.
Under optical pumping, we can describe the (rotating frame) system Hamiltonian, H(t), and the total collapse operator, n C n as: The diagonal elements in the Hamiltonian of Eq. 7 correspond to the detuning of each transition with respect to the laser frequency, while the off-diagonal elements enable Rabi driving between respective levels. In our example, the laser is on resonance with the |+, ↓ − → ← − |0, ↓ transition. The transition |−, ↑ − → ← − |0, ↑ has a detun-ing of ∆ = 160 MHz. An additional detuning, δ, is randomly sampled from a Gaussian distribution with FWHM = 2π × 20 MHz, to account for the effects of imperfect laser resonance checks or small spectral diffusion. The Rabi frequency is given by Ω(t) = α P (t), with α being a proportionality factor and P (t) the time dependent laser power. Note that we neglect driving of the (far-detuned) upper spin-orbit states and thus omit the corresponding terms. In our simulations we implement P (t) with rise/fall times as independently measured in our experiment, see Fig. S5(a). The elements in Eq. 8 correspond to relaxation processes between certain levels, see Fig. S1. Here we use the respective timescales as extracted within the main text. In Fig. S5(b-c) we plot the simulated expectations values for two different temperature scenarios, both for a driving power of 5 nW. In both cases we observe that the |+, ↓ state is depopulated, while population in the opposite orbital state, |−, ↓ , grows via the excited state, |0, ↓ . Initially, a damped Rabi oscillation between |+, ↓ and |0, ↓ is observed, before a steady state population is reached. In the case of lower temperature and hence a slower orbital relaxation constant, τ orbit , a stronger orbital pumping is observed. After the laser light is switched off, the excited state decays with τ exc , while the two orbital states relax back to an equal population with the time constant τ orbit .
In Fig. 2 of the main text we plot the expectation value of the excited state |0, ↓ , and find an excellent agreement with the experimental fluorescence counts. This confirms the accuracy of our theoretical model, and shows that the observed NV 0 dynamics are well understood. We note that the timescales for spin relaxation, τ spin , and spin pumping, τ exc,spin , are much longer than the pumping duration, and hence are not observed in this set of experiments.

Recharging dynamics
An important feature of our experiments is the ability to switch between the NV − and NV 0 charge states through resonant excitation. Ionisation (NV − − → NV 0 ) under resonant excitation has previously been experimentally studied [35,42], and an ionisation mechanism was proposed that combines a two-photon and an Auger process [35]. While a mechanism for the recharging process (NV 0 − → NV − ) was also proposed, the power dependence under resonant excitation has not been measured [35]. Here we present such measurements.
We first herald the defect in the NV 0 state, see 'check NV 0 ' step in the previous section, with the yellow laser resonant to a single transition in the lower spin-orbit branch. We then apply yellow recharging light for a certain time, before we read-out the NV − population via the 'check NV − after' step. In Fig. S6(a) (Fig. S6(b)) we plot the NV − population as a function of recharging time for linear (circular) polarization, for a few selected powers. For all powers we observe a growth in NV − population that eventually approaches unity both under high-and low-power excitation. This important observation shows both that high-fidelity switching can be achieved, and that the red read-out is robust: spectral diffusion remains minimal even for second long experiments. Additionally, the ability to reach near-unity NV − population within a few-hundred µs suggests that significant population is not trapped in the 4 A 2 level of NV 0 . Possible explanations for this observation are that the inter-system crossing rate is small, the 4 A 2 lifetime is short, or the recharging process can also occur from the 4 A 2 level itself.
For a range of recharging powers, we perform sweeps as exemplified in Figure S6 growth functions. In the simulations the NV − state is implemented as a dark state that can be populated via an excitation of the 2 A 2 NV 0 state to the conduction band. Orbital dynamics are neglected as their corresponding timescales are much faster then the observed timescales for recharging. Further, in order to perform the simulations in reasonable computational time, the recharging rate, spin-pumping rate, and spin-relaxation rates are re-scaled by a fixed factor of four orders of magnitude. Remarkably, our simulations show both qualitative and quantitative agreement with the experimental data. At very short times small deviations arise as the (unmodified) excited-state lifetime becomes comparable to the re-scaled recharging process.
We generally anticipate two recharging timescales. A fast timescale, τ fast , is associated with resonant recharging. A slow timescale, τ slow , arises from NV 0 spinpumping and spin-relaxation causing the driven transition to become dark, i.e. via a spin flip |↓ − → ← − |↑ . Recharging from the dark state occurs due to a combination of off-resonant excitation and relaxation back to the resonantly driven spin-state. Note that a slow timescale is expected for both linear and circular polarisation. The initial CR check prepares a single spinstate, and at low powers, both polarisations perform spin-selective driving as power-broadening is significantly less than the detuning between the |↑ and |↓ optical transitions. The behaviour can thus be described by a double-exponential growth function: f (t recharge ) = Ae −t recharge /τ fast + (1 − A)e −t recharge /τ slow . Where it is not possible to fit a second timescale due to the dominance of the fast timescale (high powers), a single exponential growth function is used, f (t recharge ) = Ae −t recharge /τ fast .
In Fig. S6(c), we plot all of the fitted rates. Two key features are apparent. First, for both linear and circular polarisation, the fast timescales broadly overlap and follow a linear power-dependence, as expected for a twophoton process for which the first step is resonant and in saturation. From a linear fit we extract a recharging rate of 9.3(6) Hz/nW. Second, the slow timescale is seen to be significantly faster for linear than for circular polarisation. This difference arises from the fact that off-resonant excitation is strongly suppressed by circular polarisation, preventing both spin-pumping back to the resonant transition and off-resonant recharging. The extracted rates from the Master equation simulations are in good agreement with our experimental data both for the fast and the slow timescales.
Importantly, this characterisation can be used to determine the frequency of charge-resonance checks required in each experiment to ensure that the NV 0 population remains high.

Pump-probe spectroscopy
Here, we describe the procedures followed for the pump-probe spectroscopy measurements, Fig. 3 in the main text.
For each set temperature of the cryostat, we measure a series of time traces as exemplified in Fig. 3(a) of the main text. For each time trace, we integrate the total photon counts during the first 40 ns after opening the AOM for the probe and pump measurements, respectively, and take their ratio. We then fit the recovery behaviour against the delay time, t delay , following the function f (t delay ) = a + A(1 − e −(t delay −t0)/T ), where a is the steady-state fluorescence under pumping, A is the peak fluorescence of the mixed state, t 0 is the time at which the pump is turned off, and T is the characteristic recovery time.
At higher temperatures, a challenge for these measurements arises as the required time resolution exceeds that of the AOM rise-and fall-times, see Fig. S5(a), and hence the optical pulse is not turned off for the entirety of the set t delay . We correct for this systematic error by calibrating the delay times as those for which the optical pulse falls below 90% of the full amplitude, measured using a fast photodiode. However, we note that some corrections may remain, which may lead to an underestimate of the orbital relaxation rates in those measurements. In future work, such measurements could be improved using a fast electro-optic modulator to gate the pulses.

Rate Equations
In order to extract the timescales for recharging and spin-pumping, as shown in Fig. 4 of the main text, we develop a three-level model for which we can derive analytic solutions for use in a fitting routine.
We start with the scenario of Fig. 4(c). Given the relatively slow timescales under consideration (for 5 nW of yellow excitation), we choose to neglect the orbital basis (which is effectively mixed, see Eq. 6). This leaves three levels, |↑ , |↓ and NV − , which we denote as U (t), D(t), N (t). We consider the processes of resonant recharging (r, from D − → N ), resonant spin-pumping (p, from D − → U ) and spin-relaxation (s, from D − → ← − U ), and neglect off-resonant recharging (from U − → N ), offresonant spin-pumping (from U − → D), and ionisation (from N − → D, U ), which are all expected to have negligible rates in this parameter regime. The dynamics between these levels are thus described by the following coupled equations: We impose initial conditions and derive analytic solutions to these equations (solutions available in code form upon request), which we then use as fitting functions for the measured populations. The spin relaxation time is fixed to the independently measured time of 1.51(1) s. All other parameters are free in the optimisation, leading to the following best-fit values: 1.51 0.027(1) 0.090(4) 0.960(6) 0.012 (5) Here, τ spin = 1/s, τ recharge = 1/r, τ pump = 1/p, and c 1 and c 2 are the respective populations in |↓ and |↑ after initialisation.
We now move to the scenario shown in Fig. 4(d). In this case, we stroboscopically interleave periods of 5 nW yellow excitation with periods of 500 nW of red NV − spin-pumping excitation coupled with hard microwave πpulses regularly spaced by 1.25 µs. To simplify this into a rate-equation model, we make the following assumptions (along with those already made for the previous case). First, we combine the stroboscopic driving into continuous driving with averaged resonant recharging (r, from D − → N ) and ionisation (i, from N − → D, U ) rates. We assume that charge conversion from NV − to NV 0 results in a completely mixed NV 0 spin state: that is, the rate, i, equally couples to both D and U . We note that analytic solutions incorporating an asymmetry in these couplings were derived, but that any asymmetry could not be constrained by the fitting procedure. Finally, we assume that the microwaves only couple the (not-described) spin levels of NV − -preventing spin-pumping in that charge state -and so only influence the rate i. Under this model, we arrive at the following set of coupled equations: We again impose initial conditions and derive analytic solutions to these equations (solutions available in code form upon request), which we then use as fitting functions for the measured populations. For this scenario, we fix r, p, c 1 and c 2 to the values obtained from the fit to the previous case, as the initialisation step and yellow excitation parameters are unchanged. The spin-relaxation rate, s, is fixed to twice its previous value (t spin = 1.51/2), as the time-axis on this plot is the total yellow excitation time, which is half of the total sequence time. This leaves only the ionisation rate unfixed, for which we obtain a fitted value of τ ion = 1/i = 0.018(4) s. While the behaviour at short times is well described by the model, the long-time behaviour is not well captured. If all fit parameters are unconstrained, we find that the long-time behaviour is well described, and still obtain agreement to within 50% for all values, aside from for τ spin , for which we now fit 0.14(1) s. This could indicate that the presence of red excitation and/or strong microwave driving induces an additional NV 0 spin-relaxation mechanism.

Spin cyclicity
In order to extract the cyclicity of the |↓ transition under resonant excitation as reported in the main text we use the following procedure. Here, we are interested on how many photons are scattered before the |↓ spin state is left, i.e. before leaving the three-level system {|+, ↓ , |−, ↓ , |0, ↓ }. The experiment starts when |+, ↓ is heralded, which can then be excited to |0, ↓ (for P yellow = 5 nW a Rabi flop takes about 42 ns). From here the electron decays (τ exc = 22(1) ns) either to |+, ↓ or to |−, ↓ , each with a 50% probability. If it decays to |+, ↓ it can immediately be excited again, while |−, ↓ corresponds to a dark state that can only relax back to |+, ↓ via orbital relaxation (τ orbit = 0.43(6) µs). This process continues until either a charge-or spin-flip event happens. From Fig. 4(c) we have extracted characteristic timescales of 27(1) ms (90(4) ms) for the charge conversion (spin pumping) process. Within this regime the cyclicity is limited by charge conversion, for which we calculate 0.98(7) × 10 5 scattered photons before leaving the |↓ manifold.
We note that at lower powers the recharging process, which requires a two-photon excitation, is expected to be suppressed. Here, the cyclicity will then be limited by a spin-flip process, enabling 3.2(2) × 10 5 cycles.

Readout fidelity
In the main text, we report single-shot read-out (RO) of the NV 0 spin state with fidelity, F RO ≥ 98.2(9)%. Here we outline the characterisation procedure.
The readout fidelity is defined by the relation F RO = 1 2 (F ↓|↓ + F ↑|↑ ), where F i|j is the probability to assign the spin state |i after preparing |j . Note that this assumes that each spin state can be initialised perfectly. Imperfect initialisation will lead to a decrease in the achievable fidelity: the calculated RO fidelities are thus a lower bound.
In the presented experiments, we use a single yellow laser for both state initialisation and readout. As this means that we are only able to herald the |↓ state with high fidelity, we use the intrinsic spin relaxation process of NV 0 to prepare a mixed state, which can also be used to calculate the fidelity for |↑ read-out. The histograms presented in Fig. 4(a) of the main text are used for this calculation.
First, to calculate F ↓|↓ and F ↑|↓ , we herald |↓ using a threshold of 25 photons for the 'check NV 0 ' step of Fig. S4. After a brief delay (0.1 ms), we check for any residual population in NV − ('check NV − after'), which is discarded. We then perform the 'check NV 0 ' step again to read out the spin population. From 3000 experimental shots, we discard 52 cases, corresponding to an NV − population of 1.7(2)%. Of the remainder, 98.4(2)% of cases match or exceed the chosen threshold of 5 photons, see main text. We thus have: F ↓|↓ = 98.4(2)%, and F ↑|↓ = 1.6(2)%. We note that the initialisation fidelity of the |↓ heralding step (25 nW for 250 µs) is likely limited by spin pumping to |↑ . From independent spin-pumping measurements we expect a reduction of that population by ∼ 0.8% over the 250 µs, though this is partially mitigated by the initialisation threshold of 25 photons.
To calculate F ↑|↑ and F ↓|↑ , we repeat the procedure, but now wait for 10 s between the first 'check NV 0 ' step and the 'check NV − after' step, preparing the mixed state. We anticipate a state preparation infidelity of < 0.2% arising from the finite waiting time.
Combining the results we arrive at a single-shot readout fidelity, F RO ≥ 98.2(9)%.