Measurement and Control of Single Nitrogen-Vacancy Center Spins above 600 K

We study the spin and orbital dynamics of single nitrogen-vacancy (NV) centers in diamond between room temperature and 700 K. We find that the ability to optically address and coherently control single spins above room temperature is limited by nonradiative processes that quench the NV center's fluorescence-based spin readout between 550 and 700 K. Combined with electronic structure calculations, our measurements indicate that the energy difference between the 3E and 1A1 electronic states is approximately 0.8 eV. We also demonstrate that the inhomogeneous spin lifetime (T2*) is temperature independent up to at least 625 K, suggesting that single NV centers could be applied as nanoscale thermometers over a broad temperature range.

The negatively-charged nitrogen-vacancy (NV) center spin in diamond stands out among individuallyaddressable qubit systems because it can be initialized, coherently controlled, and read out at room temperature [1]. The defect's robust spin coherence [2] and optical addressability via spin-dependent orbital transitions [3] have enabled applications ranging from quantum information processing [4][5][6][7] to nanoscale magnetic and electric field sensing [8][9][10]. While it has been shown that NV center spins can be optically polarized up to at least 500 K [11,12], little is known about what processes limit the spin's optical addressability and coherence at higher temperatures. Understanding these processes is important to high-temperature field sensing applications and will aid the search for new defect-based spin qubits analogous to the NV center [13,14] by identifying the aspects of its orbital structure responsible for its high temperature operation.
The NV center's optical spin polarization and optical spin readout result from a spin-selective intersystem crossing (ISC). Although optical transitions between the spin-triplet ground ( 3 A 2 ) and excited ( 3 E) states (1.945 eV zero-phonon line [ZPL]) are typically spin conserving, the 3 E state can also relax to the 3 A 2 state via an indirect pathway that involves a nonradiative, triplet to singlet ISC and subsequent transitions through at least one additional singlet [ Fig. 1(a)]. The 3 E ISC is much stronger for the m s = ±1 3 E sublevels than for the m s = 0 sublevel, which facilitates spin readout through the resulting spin-dependent photoluminescence (PL) and initializes the spin into the m s = 0 3 A 2 sublevel with high probability (P ms=0 ∼ 0.8) through repeated optical excitation [15]. Despite the singlet pathway's critical role in preparing and interrogating the spin, open questions remain regarding the number and energies of the singlets involved. Recent experiments have established that it consists of at least two singlets of 1 A 1 and 1 E symmetry separated by 1.19 eV [16][17][18], and have shown that the 1 E state lifetime is strongly temperature dependent below 300 K [17,19]. The details of the 3 E ISC remain uncertain, however, given computations that suggest the singlet pathway involves three singlets [20] and the lack of direct measurements of the 3 E to singlet transition energy.
In this work, we report measurements of single NV center spins between 300 K and 700 K. We show that thermally-activated nonradiative processes diminish the spin-selectivity of the 3 E ISC and quench the optical spin readout above 550 K. Our measurements indicate the energy barrier for nonradiative relaxation from the 3 E state is ∼ 0.5 eV. We perform electronic structure calculations that suggest an orbital structure consistent with a two-singlet decay pathway, the measured energy barrier for nonradiative relaxation, and the temperaturedependence of the 1 E state lifetime. Furthermore, we demonstrate that the spin coherence remains robust even as nonradiative processes quench the optical spin readout. In addition to demonstrating the defect's potential for high-temperature magnetic field sensing applications, these results suggest single NV centers could act as nanoscale thermometers with sensitivities on the order of 100 mK/ √ Hz between room temperature and 600 K. This high sensitivity and wide range of operating temperatures make NV centers an attractive candidate for a variety of thermosensing applications such as cellular thermometry [21] and diamond-based scanning thermal microscopy [22].
To study single NV centers at elevated temperatures, we constructed devices that combine on-chip heating and thermometry elements, solid immersion lenses (SILs) for enhanced photon collection efficiency [23,24], and onchip microwave elements for spin control [ Fig. 1(b)] [25,26]. We studied NV centers in a single-crystal diamond sample ([N] 0 S < 5 ppb) using a home-built confocal microscope [27]. We etched hemispherical SILs around single defects with a focused ion beam, patterned mi- The straight lines represent optical transitions and the snaked lines represent nonradiative intersystem crossings (ISC). Bra-ket notation denotes the effective levels used to model the photoluminescence intensity (IPL) and electron spin resonance (ESR) contrast. Level |2 (|4 ) is a grouping of the 3 A2 ( 3 E) ms = ±1 sublevels and |5 is an effective singlet level. (b) Sample schematic showing the resistive heater and resistive temperature detector (RTD), the solid immersion lens (SIL), and the antenna. The scanning electron microscope image shows the SIL and antenna. (c) IPL of a single NV center versus applied microwave frequency. Lorentzian dips with normalized amplitude C are observed at the spin resonance frequencies (B ∼ 45 G). The left (right) axis shows the normalized (absolute) IPL to illustrate the temperature-dependent decrease in C and IPL. The red lines are two-Lorentzian fits. The crystal field splitting (D) is the average of the resonance frequencies. The measurements were performed at high microwave power to saturate the spin transitions.
crowave antennas around the NV centers, and patterned a resistive heater on the sample surface. We adhered a commercial resistive temperature detector (RTD) to the diamond and connected the heater and RTD to a temperature controller to achieve stability within ±50 mK up to 700 K in air.
Continuous-wave (CW) electron spin resonance (ESR) measurements revealed that the NV center's PL intensity (I PL ) and relative I PL difference between its spin states (ESR contrast) strongly decrease above 550 K.
The measurements were carried out by applying a sweptfrequency AC magnetic field while monitoring I PL under continuous 532 nm laser excitation. We observed Lorentzian dips in I PL with normalized amplitude C centered at the ground-state spin resonance frequencies [3]. The resonance frequencies are determined by the crystal field splitting (D) between the m s = ±1 and m s = 0 states, and the Zeeman shift due to an applied magnetic field [28]. We applied a ∼ 45 G magnetic field along the defect's symmetry axis to separate the m s = ±1 states by ∼ 250 MHz. Consistent with measurements performed at lower temperatures [29,30], we observed shifts in D on the order of 100 kHz/K due to lattice expansion. In addition to this shift, both C and the off-resonant I PL showed a pronounced decay above 550 K; by 700 K the Lorentzian dips were barely observable and the offresonant I PL dropped to ∼ 20% of its room-temperature value [ Fig. 1(c)].
Thermal quenching of a point defect's I PL often results from thermally activated nonradiative processes that shorten the effective optical lifetime of an emitter [31]. Previous ensemble measurements of the temperaturedependence of the NV center's optical lifetime in bulk diamond [32] and nanodiamonds [33] have shown conflicting results. We therefore have measured the spinresolved excited-state lifetimes (τ ms=0 and τ ms=±1 ) for a single NV center as a function of temperature following the method of Ref. [19] to investigate their influence on ESR contrast and I PL [Fig. 2]. The spin-dependence of the 3 E ISC preferentially shortens τ ms=±1 and thus optical excitation leads to biexponential fluorescence decay, where the decay constants correspond to τ ms=0 and τ ms=±1 and the relative amplitudes reflect the spin polarization. To probe the lifetimes, we initialized the spin into m s = 0 using a 1.4 µs pulse from a 532 nm laser. After waiting 500 ns for the singlets to depopulate, we applied a resonant microwave pulse to rotate the spin into a superposition of m s = 0 and m s = −1. We then applied a 555 nm picosecond laser pulse and measured the resulting PL with a time-correlated photon counting module. We repeated this measurement for varying spin rotation angles and performed global Bayesian inference to infer τ ms=0 , τ ms=±1 , and the finite spin polarization at each rotation angle from the measurements (see Ref. [26] and references [34][35][36][37][38] therein). Above 550 K, τ ms=0 (T ) showed a sharp reduction with increasing temperature. In spite of this, we found the ground-state spin polarization after optical pumping was temperature independent within our measurement uncertainty up to 650 K [26].
We find that τ ms=0 (T ) is accurately described by the Mott-Seitz formula for nonradiative relaxation via multiphonon emission [31], where 1/τ ms=0 (300 K) = 1/τ rad + k ms=0 . Here 1/τ rad is the radiative rate, k ms=0 is the room temperature nonradiative transition rate from the 3 E m s = 0 sublevel to the uppermost singlet state [ Fig. 1(a)], s ms=0 is the frequency factor, and ∆E is the energy barrier for the nonradiative process. From a fit to τ ms=0 (T ), we find τ ms=0 (300 K) = (13.4 ± 0.6) ns, s ms=0 < 1.64 × 10 4 , and ∆E = 0.48 +0.15 −0.13 eV. In general, ∆E is interpreted via a one-dimensional configuration coordinate diagram. These diagrams show the dependence of total energies in different electronic states as a function of a generalized configuration coordinate ∆R. The latter measures the total displacement of all atoms along the path that interpolates between equilibrium geometries in the relevant electronic states. ∆E is then classically defined as the energy difference between the intersection point of two potential energy curves and the energy minimum of the upper curve. In such a scenario, nonradiative relaxation is interpreted as a phonon-assisted process over the energy barrier.
We combine our measurement of ∆E with electronic structure calculations to construct a configuration coordinate diagram (CCD) for the NV − center in its 3 A 2 , 3 E, 1 A 1 , and 1 E electronic states. The theoretical methodology employed relies on the use of supercells to describe defects and hybrid density functionals for an accurate treatment of the electronic structure [13,39]. Additional details regarding the computational methods are provided in the Supplementary Material [26] and references [40][41][42][43][44][45][46][47] therein. The calculated ZPL between the triplet states, 1.972 eV, agrees favorably with the experimental value of 1.945 eV. In Fig. 3, we placed the minima of 3 A 2 and 3 E states according to the experimental ZPL with zero defined by the 3 A 2 state's equilibrium energy. We mapped the potential energy surfaces for these two ∆R measures the total displacement of all atoms along the path that interpolates between equilibrium geometries in the 3 A2 and 3 E electronic states. All energies are relative to the 3 A2 state's equilibrium energy. We placed the 1 A1 state based on the assumption that its equilibrium atomic configuration and vibrational frequencies were similar to those of the 3 A2 state and using our measurement of ∆E. The 1 E state was placed 1.19 eV below the 1 A1 state based on their known energy difference. The uncertainties in the 1 A1 and 1 E energies reflect the uncertainty in ∆E (68% interval, dashed red lines). The triplet (1.945 eV) and singlet (1.19 eV) zero phonon lines are also indicated.
states along the path that linearly interpolates between their equilibrium geometries. The resulting Stokes and anti-Stokes shifts were 0.27 eV and 0.23 eV, respectively.
In constructing the CCD for the singlet states, we also assumed that the 1 A 1 state's equilibrium atomic configuration and vibrational frequencies were equal to those of the 3 A 2 state. Indeed, the 1 A 1 state's a 2 1 e 2 electron configuration [48,49] implies the two states have a similar electron density. In turn, this implies that the ionic forces are similar for these two states away from their equilibrium geometries. Therefore, we conclude the Huang-Rhys factor (S), the average number of nonequilibrium phonons emitted when the electronic state changes [31], for the 3 E to 1 A 1 transition is the same as that for the 3 E to 3 A 2 optical transition. For the latter our calculations give a value of S ≈ 3.2, close to the value of S ≈ 3.65 inferred from the relative spectral weight of the ZPL with respect to the phonon sideband [50]. The value S > 1 for the 3 E to 1 A 1 transition confirms our hypothesis, implied in Eq. (1), that the ISC is accompanied by multiphonon emission.
We interpret ∆E as the energy difference between the intersection point of the 3 E and 1 A 1 potential energy curves and the energy minimum of the 3 E state. This places the 1 A 1 state (0.76 ± 0.07) eV below the 3 E state. By fixing ∆E based on our measurements, our analysis avoids difficulties in computing the absolute energies of the singlet states resulting from their multi-determinant electronic wavefunctions [48,49]. We placed the 1 E state 1.19 eV below the 1 A 1 state [16][17][18] and find that the 1 E singlet is close in energy to the 3 A 2 triplet. Since the 3 A 2 state is the ground state, the 1 E state should be higher in energy. Within the uncertainty in ∆E, the configuration coordinate diagram fulfills this constraint and indicates that only small energy barriers for nonradiative relaxation from the 1 E state are possible in qualitative agreement with energy barriers inferred from the temperature-dependence of the 1 E state lifetime [17,19].
To demonstrate that the neutral charge state of the NV center (NV 0 ) does not interfere in the temperatureactivated non-radiative processes discussed above, in Fig.  3 we also placed the potential energy curve for the NV 0 center in its 2 E ground state based on the charge-state transition levels given in Ref. [13] and the calculated total-quadratic distortion ∆R. The resulting high ionization potential for the NV − center in the 3 E state (> 0.8 eV) shows that the escape of the electron to the conduction band minimum is negligible for the temperatures studied.
Our inference of the 1 A 1 and 1 E energies is influenced by our classical interpretation of ∆E, which is based on the consideration of a single effective phonon mode of a 1 symmetry that strongly couples to the distortion of the defect [31]. In fact, the temperature-activated ISC from the 3 E m s = 0 state to the 1 A 1 state also requires coupling to phonons of e symmetry, while this is unnecessary for the ISC from the 3 E m s = ±1 states [18]. Nevertheless, the presented model is a good approximation that provides a suitable description of the temperaturedependent ISC. The generalization of this model that takes into account quantum effects as well as coupling to phonons of different symmetry during radiative and non-radiative processes will be presented elsewhere [51].
Having interpreted ∆E, we employ a model of the ESR experiments using the density matrix formalism that includes dissipation to reveal that the observed decreases in C and I PL are predominantly caused by the temperature dependence of τ ms=0 and τ ms=±1 [52]. The model consists of irreversible transitions between five effective energy levels [ Fig. 1(a)] under optical pumping with an offdiagonal term to capture the coherent microwave driving [26]. The only temperature-dependent components are τ ms=0 , τ ms=±1 , and the known temperature-dependence of the effective singlet lifetime [19]. The model shows agreement with our measurements of C(T )/C(300 K) and I PL (T )/I PL (300 K) [ Fig. 4]. It reveals that the small decrease in C below 550 K is due to shortening of the singlet lifetime while the decreases in both C and I PL above 550 K are primarily caused by the thermally-activated nonradiative processes that shorten τ ms=0 and τ ms=±1 .
We established from pulsed ESR measurements that the NV center remains spin coherent even as nonradiative processes diminish the ESR contrast and I PL . The spin's robust coherence is due to the equal population of the sublevels of the 13 C nuclear spin impurities which decohere the spin at all but cryogenic temperatures [53]. We performed these measurements by polarizing the spin with a 532 nm CW laser pulse, coherently manipulating the spin with resonant microwaves, and reading out the spin state by monitoring the PL [1]. Both Rabi [ Fig.  5(a)] and Ramsey pulse sequences were employed (see Ref. [26] and references [54,55] therein), demonstrating that the spin can be coherently controlled at temperatures exceeding 600 K and showing I PL and ESR contrast trends consistent with the fluorescence quenching observed in our CW ESR measurements. Moreover, fits to the Ramsey oscillation decay envelope showed T * 2 is temperature independent up to at least 625 K [ Fig. 5 We also performed longitudinal relaxation measurements by optically polarizing the spin into m s = 0 and measuring the relaxation into m s = ±1 through the timedependent I PL . The measurements demonstrated that T 1 at 600 K (340 ± 50 µs) was still much longer than the period of the pulsed ESR measurements, confirming that spin relaxation does not explain the observed decrease in the ESR contrast and I PL .
The persistence of the spin coherence and the strong temperature dependence of the crystal field splitting, D, suggest the possibility of using the spin resonances for thermometry. Based on our CW ESR measurements, we find that D(T ) [Fig. 5(c)] is accurately described by a third-order polynomial, D(T ) = a 0 + a 1 T + a 2 T 2 + a 3 T 3 , between 300 K and 700 K, with a 0 = (2.8697 ± 0.0009) GHz, a 1 = (9.7 ± 0.6) × 10 −5 GHz/K, a 2 = (−3.7 ± 0.1)×10 −7 GHz/K 2 , and a 3 = (1.7 ± 0.1)× 2 versus temperature inferred from the decay envelope of fits to Ramsey oscillations. These measurements were performed at ∼500 G to polarize the NV center's 14 N nuclear spin to simplify its hyperfine spectrum. All other measurements were performed at ∼ 45 G. Error bars represent 68% intervals. (c) The crystal field splitting, D (T ), inferred from Lorentzian fits to the CW ESR data. The solid red line is the third-order polynomial fit to the data given in the main text. (d) Single-spin thermal sensitivity (η). η was calculated using Eq. (2) with all parameters from the same NV center.
10 −10 GHz/K 3 . The derivative of this polynomial yields thermal shifts ranging from 80 kHz/K at 300 K to 170 kHz/K at 700 K. The possibility of spin-based thermometry is further motivated by the established methods for measuring small static shifts in the NV center spin resonance frequencies [8,9,56,57]. Perhaps the best known method for sensing such shifts is through a Ramsey measurement, in which a superposition of two spin states is created and evolves under the influence of static dephasing fields. After a time t, a microwave control pulse is applied to convert the accumulated phase into a population difference of the two spin states that can be read out optically.
We quantify the thermal sensitivity (η) based on our measurements of the defects's temperature-dependent optical and spin properties using an expression analogous to those derived for NV center DC magnetic field sensing with a Ramsey pulse sequence [58]: A(T ) accounts for the finite photon count rate and ESR contrast and is given by where α 0 = τ F I PL (T ) and α 1 = τ F I PL (T ) (1 − C P (T )).
Here τ F is the fluorescence readout duration (350 ns) and C P (T ) is the contrast of a pulsed ESR measurement (0.24 for our Ramsey measurements). Based on the data presented in Fig. 4 and Fig. 5(b,c), we find that A (300 K) ∼ 0.02 and η ∼ 100 mK/ √ Hz between room temperature and 600 K [ Fig. 5(d)]. At higher temperatures, the sensitivity becomes limited by the quenching of the optical spin readout. This sensitivity could be enhanced for NV centers in isotopically pure diamond, which can exhibit T * 2 times greater than 100 µs [59], corresponding to η better than 10 mK/ √ Hz. These results demonstrate that NV centers offer operating temperatures and thermal sensitivities relevant to studies of thermal gradients in a variety of microscale and nanoscale systems. Recent demonstrations of the incorporation and quantum control of nanocrystal diamonds containing NV centers in living cells [60] suggest that NV centers could provide an alternative to quantum dot nanostructures for cellular thermometry [21]. In addition, as the ideal thermal and mechanical properties of diamond have motivated its use in scanning thermal microscopy tips [22], advances in the incorporation of NV centers in scanning probe tips [61,62] indicate the possibility of spin-based thermal microscopy for investigating materials such as nanostructured thermoelectrics [63,64]. Finally, given the large thermal shifts in D and the ability to optically polarize the spin at high temperatures, highdensity ensembles of NV centers [65] could be used as an in-situ thermometer for electron paramagnetic resonance experiments in analogy with chemical-shift thermometers used in nuclear magnetic resonance [66].
The results presented here provide fundamental insights on how thermally-driven nonradiative processes can influence the fluorescence properties of defect-based spin qubits. We show that the NV center's fluorescencebased spin readout is ultimately limited above room temperature by nonradiative processes that diminish the spin-selectivity of its 3 E to 1 A 1 orbital transition and upset the balance between the defect's radiative and nonradiative relaxation rates. Coupled with the persistence of its spin coherence, these measurements suggest the NV center could find application as a nanoscale thermal sensor over a broad temperature range.
This work was funded by the AFOSR, ARO, and DARPA. We acknowledge computational resources at XSEDE. A portion of this work was done in the UCSB nanofabrication facility, part of the NSF funded NNIN network. AA was supported by a grant from the Swiss NSF.

Sample fabrication
The sample used in this work was a 4 mm x 4 mm x 0.5 mm (100) Electronic Grade CVD diamond from Element Six. The sample was irradiated with 2 × 10 14 /cm 2 , 2 MeV electrons and annealed at 800 • C in forming gas to increase the nitrogen-vacancy (NV) center density. The solid immersion lenses were etched using a Ga + focused ion beam system (FEI Helios 600) at an accelerating voltage of 30 kV and an etch current of 2.8 nA. After milling, we used an Ar/Cl 2 inductively-coupled plasma etch to remove approximately 80 nm from the surface to remove damage from the Ga + ion implantation. A thick Ti/Pt metallization was used for the heaters and ring antennas (10 nm Ti, 400 nm Pt) so that the devices would operate well in air at elevated temperatures. The area of the diamond where the resistive temperature detector (Omega Engineering, RTD-3-F3105-72-G-B) was bonded with ceramic adhesive (Cotronics, Durabond 950) was first roughened with an oxygen plasma to promote adhesion. The samples were mounted with ceramic adhesive on glass to thermally isolate the diamond and were packaged in brass holders with coplanar waveguide transmission lines for microwave connections. The temperature of the sample was controlled with a PID feedback loop using a Lakeshore 340 temperature controller. The temperature stability of the sample was measured using the Lakeshore controller to be within ±50 mK of the temperature setpoint at all temperatures investigated. We also independently calibrated our thermometry by heating the sample with the controller feedback loop and separately measuring the sample temperature with an infrared microscope and found good agreement between the two measurements.

Experimental setup
The confocal microscope setup used for this work was similar to the microscope used for angleresolved magneto-photoluminescence studies described in Ref. [1]. Here we will describe additional experimental details relevant to the work described in the main text.
A 0.7 NA microscope objective (Mitutoyo Plan Apo) with a 6.5 mm working distance was used in order to provide clearance for the RTD on the sample surface. The overall timing of the pulsed measurements was controlled by both a pulse pattern generator (Agilent, 81110A) and a digital delay generator (Stanford Research Systems, DG645). Pulsed microwaves for ESR measurements were generated with a signal generator, I/Q modulator, and arbitrary waveform generator from National Instruments (PXIe-5673E) in series with a microwave amplifier (Amplifier Research, 5S1G4). An acousto-optic modulator (AOM) from Isomet (1250C) was used to modulate the 532 nm CW laser (Oxxius, 532S-150-COL-OE) for pulsed ESR measurements. The voltage output from the AOM driver (Isomet, 525C-2) was optimized with a microwave switch, attenuators, and amplifier to provide an optical extinction of ∼ 50 dB in a single pass configuration. Photon counting was performed with a silicon avalanche photodiode (Perkin Elmer, SPCM-AQRH-15-FC). For all measurements discussed in the main text the NV center's fluorescence was collected between the wavelengths of 630 nm and 800 nm. It has been shown previously that the NV center's phonon sideband does not significantly broaden over the range of temperatures measured in this work [2].
Measurements of the optical lifetime were carried out using a mode-locked Ti:Sapphire laser

NV center optical lifetime fits
We applied a Bayesian approach to infer the spin-resolved fluorescence lifetime parameters and their associated intervals for the uncertainties shown in the main text. The difficulty of fitting sums of exponentials to experimental data is a known problem and originates from the high correlation between amplitude and lifetime parameters that tends to confound them [3]. Our approach forgoes the use of asymptotic approximations in favor of calculating the actual joint posterior density to more accurately infer the parameters and their associated errors.
Before we analyze the data, we remove approximately 2.5 ns of histogram data after the initial onset of fluorescence to eliminate fast-decaying background photoluminescence detected in the sample. Because of the finite extinction of our electro-optical modulator, very small ∼2 ns fluorescence bursts with a 76 MHz spacing appear in our raw data. Of these, we removed the first two and truncated the dataset by removing points including and after the third unwanted pulse.
At each temperature, we found about six spin rotation angles were needed for good resolution of the separate lifetimes. At the driving strength used for these measurements, we could perform a π rotation in 80 ns.
For our analysis, we assume that for each microwave burst length, the fluorescence decay follows biexponential decay with a constant offset subject to random error, where is normally distributed with mean 0 and variance η 2 σ 2 Poisson . Here σ Poisson is the square root of the noiseless value of f ij (t k ) and η is a dimensionless scaling parameter strictly above unity that accounts for other random noise in the experiment. This value is inferred from analysis of experimental data to be very close to unity, serving as a measure of goodness of fit that indicates our measurements are nearly shot-noise limited, and that any systematic effects from unextinguished pulses are within the random error of the measurement. For each microwave burst length, A j , B j , and C j are free parameters while parameters τ ms=0,i and k i are shared between burst lengths at a given temperature. τ ms=±1,i is related to τ ms=0,i through the rate relation, where k i is the difference in excited state intersystem crossing rates between the two spin states.
All parameters in the analysis use non-informative and proper (bounded) prior distributions and all posterior distributions were found to be relatively insensitive to the particular form of the prior.
Because the lifetimes are shared globally between data sets at the same temperature, the posterior density typically exceeds 30 dimensions, making the use of standard numerical integration techniques infeasible. To compute marginal posterior densities for the lifetimes plotted in the main text, we employ the DREAM ZS algorithm written in MATLAB [4,5], a Markov Chain Monte Carlo (MCMC) approach that uses adaptive proposal distribution tuning, sampling from the past, and snooker update on parallel chains to rapidly explore high-dimensional posterior distributions.
In MCMC, each of the N chains executes a random walk through the parameter space following a modified Metropolis-Hastings rule to control whether a proposed d-dimensional move is accepted or rejected. Because the algorithm maintains detailed balance at each step, the target distribution after a burn-in period is the desired posterior distribution. We find good results using the recommended settings along with N = 8 parallel chains. Convergence to the target distribution was assessed both graphically and with the Gelman-Rubin statistic,R < 1.05 [6]. The point estimates for the lifetimes and other parameters are computed from the respective sample empirical means, and the highest posterior density intervals are computed using the method of Chen and Shao [7].
To infer the parameters in the Mott-Seitz equation for non-radiative decay, we apply a similar technique of bounding the base lifetime, frequency factor, and energy barrier and computing the respective marginal posterior densities. As an approximation, we assume the likelihoods for the choice of this prior since we were unable to measure above 675 K. The analysis also allows for a linear systematic temperature offset up to a maximum of 10 K at 700 K, although this correction is essentially negligible. We compute the marginal densities for the Mott-Seitz parameters and quote summary statistics in the main text as previously described.

Measurements of multiple NV centers
All measurements discussed in the main text were taken on single NV centers. The data shown in Fig. 1, Fig. 2, and Fig. 5(b,c,d) was taken on an NV center (NV1) in a solid immersion lens while the data in Fig. 5(a) was taken on another NV center (NV2) in bulk diamond. The two NV centers showed nearly identical electron spin resonance (ESR) contrast temperature dependencies [ Fig. S1]. NV1 and NV2 showed similar spin-dependent excited-state lifetime temperature dependencies up to 625 K, although the lower photon collection efficiency on NV2 precluded spin-resolved lifetime measurements on this NV center at temperatures above 625 K.

Rate equation modeling
A simplified model of the NV center with five energy levels is used to model the optically detected ESR signal observed in experiment [8]. We augment the classical rate equation model of Ref. [8] with the off-diagonal terms from the density matrix that give rise to coherent driving in the orbital ground state with dissipation. We numerically integrate this expression, a master equation of the Lindblad form, to simulate the time evolution of the system under both optical and microwave excitation that correspond to experimental conditions [9]. We use the optical decay rates measured in our experiment and the non-spin-conserving optical transition rates from Ref. [8], and assume that they are temperature-independent. To capture the temperature dependence of the transition rates, we use the best fit equation versus temperature for the singlet rates from Ref. [10] and which is consistent with the polarizations inferred in our analysis shown in Fig. S2 to within the statistical uncertainty. We summarize the rates used in Table I.

Ramsey measurements
The T * 2 values in the main text were inferred from fits to the decay of Ramsey oscillation measurements. The Ramsey measurements were carried out by first optically initializing the spin using a 2 µs 532 nm laser pulse. After waiting ∼500 ns for the singlet pathway to depopulate, we rotated the spin into a superposition of its m s = 0 and m s = −1 spin states with a π/2 pulse and allowed the spin to evolve for a variable time τ R . A second π/2 pulse was used to convert the phase of the spin acquired during free evolution into a population difference and the spin was then read out by collecting its fluorescence for the first ∼ 350 ns of the following laser pulse. The primary NV center studied in this work showed a hyperfine coupling of ∼2.5 MHz to a nearby 13 C (I = 1 2 ) in addition to its hyperfine coupling to the 14 N nuclear spin (I = 1). To simplify the hyperfine spectrum to allow for more accurate fits to T * 2 we performed the Ramsey measurements at 500 G to polarize the nitrogen nuclear spin into its m I = +1 state [11]. This resulted in a hyperfine spectrum with two visible ESR transitions showing a relative polarization consistent with Ref. [12]. For the Ramsey measurements, we detuned the π/2 pulses by -2 MHz from the lower frequency ESR transition to induce oscillations in the measured PL intensity as a function of τ R . We fit the resulting data [ Fig. S3] with the function to account for the two visible ESR transitions and to extract T * 2 . The Rabi measurements shown in the main text were carried out with a similar sampling scheme except that resonant microwave pulses of varying duration (burst length) were applied instead of the π/2τ Rπ/2 portion of the Ramsey pulse sequence.

ELECTRONIC STRUCTURE CALCULATIONS
The configuration coordinate diagram (CCD) presented in Fig. 3 of the main text depends on the values of Stokes (S) and anti-Stokes (AS) shifts pertaining to the optical transition between triplet states (i.e. relaxation energies in the excited and ground states). While these quantities influence greatly the shape of optical absorption and luminescence bands, it is not straightforward to extract them from measurements, such as those presented in Ref. [13]. Therefore, to determine the values of these shifts, we have utilized electronic structure calculations [14,15].
The calculations in this work were performed employing the hybrid density functional of Heyd, Scuseria, and Ernzerhnof (HSE06) to describe the NV center's electronic structure [16]. Atomic cores were treated via the projector-augmented-wave (PAW) scheme [17], and the plane wave cutoff for valence wave functions was set to 400 eV. To simulate the NV − center, we used a cubic supercell nominally containing 216 atoms. The Γ-point was used for the Brillouin zone integration.
To simulate the excited state triplet, we used a method outlined in Ref. [14], whereby an electron in the spin-minority channel is promoted to the e state thus leaving a hole in the a 1 state. We employed the code VASP [18]. A similar computational methodology was also used in Refs. [14,15].
The equilibrium geometry of the defect in the 3 A 2 state has been determined by fully relaxing the geometry. The resulting point symmetry was C 3v , and the electron configuration a 2 1 e 2 . In the 3 E state the C 3v symmetry is broken, but the energy associated with symmetry breaking is small (∼0.025 eV) [19]. In the main text it is assumed that the C 3v symmetry is maintained for equilibrium geometries of all electronic states of the NV − center. In particular, it implies that the symmetry is maintained along the path that interpolates between equilibrium geometries of the defect in 3 A 2 and 3 E electronic states [cf. Fig. 3]. This reasoning neglects the occurrence of the (dynamic) Jahn-Teller effect in the 3 E (and possibly 1 E) state (see Refs. [19,20] ).
The calculated parameters of the defect depend somewhat sensitively on the adopted lattice constant. The theoretical lattice constant of diamond is 3.539Å, and is thus slightly underestimated with respect to the experimental value of 3.567Å [21]. Employing the theoretical lattice constant, we obtain the value of 2.035 eV for the zero phonon line (ZPL) of the optical transition between the triplet 3 E and 3 A 2 states. The resulting S and AS shifts are 0.291 and 0.245 eV. If the experimental lattice constant is adopted, the ZPL reduces to 1.972 eV (thus in better agreement with the experimental value of 1.945 eV). Similarly, S and AS shifts decrease, to 0.274 and 0.230 eV, respectively. Thus, for a ≤1 % change in the lattice constant, the resulting change in S and AS is more than 5%. This is a likely consequence of the high rigidity of the diamond lattice. We have used the parameters from the second set of calculations to construct the CCD presented in The construction of the configuration coordinate diagram in Fig. 3 is based on the assumption that the equilibrium electron density of the 1 A 1 state is very similar to that of the 3 A 2 state along the generalized configuration coordinate. In the molecular orbital picture, the wavefunction of the The wavefunction of the 1 A 1 state is a sum of two Slater determinants [22,23]: The explicit form of this wavefunction, focusing on the two electrons in the e states, is then: [e x (1)e x (2) + e y (1)e y (2)] · 1 √ 2 [α(1)β(2) − β(1)α(2)] , with a symmetric real-space part and an anti-symmetric spin part. The real-space part of the wavefunction has thus a bosonic symmetry, and this causes obvious problems when treating such wavefunctions within density-functional theory based approaches. However, similar to the 3 A 2 state, the electron density associated with e electrons in the 1 A 1 state can be shown to be e 2 x (1) + e 2 y (1): x (1)e 2 x (2) + 2e x (1)e y (1)e x (2)e y (2) + e 2 y (1)e 2 ye (2) d2 = e 2 x (1) + e 2 y (1), (S9) Assuming that the single-particle orbitals do not change much in the two electronic states, this implies that electron densities and thus forces on atoms are very similar in the two electronic states as long as the C 3v symmetry is maintained. Furthermore, in this consideration the mixing of the 1 A 1 singlet with higher-lying singlets (with the nominal electron configuration e 4 ) is ignored.
Since the latter has a very high energy [24], this assumption seems justified.
The conclusion about very similar electron densities in the singlet and the triplet state seems to be in contrast to knowledge acquired considering Heitler-London-type wavefunctions. These occur, e.g. in studying the hydrogen molecule or coupled quantum dots. In these situations electron densities of the triplet and the singlet wavefunctions are very different due to the appearance of mixed densities of the form e x (1)e y (1) [25]. This is due to the fact that single-electron orbitals in the Heitler-London wavefunction are typically not orthogonal.
A similar reasoning applies to the electron density of the 1 E state. In this case the C 3v symmetry can be broken, but the energy related to symmetry breaking is small (see above). However, 1 E level can interact with 1 E level (electron configuration a 1 e 3 ), and thus the resulting electron density is slightly different from that of the 1 A 1 state. There are some indications that the admixing coefficient of the a 1 e 3 configuration is about 0.3 [26]. In broad terms, this indicates that the equilibrium geometry of the defect in the 1 E state is (approximately) a weighted average of the geometries in 3 A 2 and 3 E states. This information was used in placing the total energy curve pertaining to the 1 E state in Fig. 3.
In the main text we mention the experimental results of Robledo et al. [10] and Acosta et al. [27] as indicating that a small energy barrier exists for nonradiative relaxation from the 1 E state. Under the assumption that the ISC is dominated by single-phonon emission, energies of 17 meV [10] and 15 meV [27] have been inferred from the temperature dependence of the 1 E state