State-dependent phonon-limited spin relaxation of nitrogen-vacancy centers

Understanding the limits to the spin-coherence of the nitrogen-vacancy (NV) center in diamond is vital to realizing the full potential of this quantum system. We show that relaxation on the $|m_{s}=-1\rangle \leftrightarrow |m_{s}=+1\rangle$ transition occurs approximately twice as fast as relaxation on the $|m_{s}=0\rangle \leftrightarrow |m_{s}=\pm 1\rangle$ transitions under ambient conditions in native NVs in high-purity bulk diamond. The rates we observe are independent of NV concentration over four orders of magnitude, indicating they are limited by spin-phonon interactions. We find that the maximum theoretically achievable coherence time for an NV at 295 K is limited to 6.8(2) ms. Finally, we present a theoretical analysis of our results that suggests Orbach-like relaxation from quasilocalized phonons or contributions due to higher-order terms in the spin-phonon Hamiltonian are the dominant mechanism behind $|m_{s}=-1\rangle \leftrightarrow |m_{s}=+1\rangle$ relaxation, motivating future measurements of the temperature dependence of this relaxation rate.


INTRODUCTION
The nitrogen-vacancy (NV) center in diamond is a solid-state quantum system that has attracted interest for a range of potential applications, such as quantum information processing [1,2] and sensing of the local magnetic [3], electric [4,5], strain [6,7], and thermal [8,9] environment. Critical to the performance of the NV center in many of these applications is its long electronic spin coherence time at room temperature. Coherence times of approximately 2-3 ms have previously been demonstrated in bulk diamond at room temperature [10,11]. However, despite the use of sophisticated dynamical decoupling sequences and high-purity diamond samples, these values remain far short of the measured lifetime of the |m s = 0 state (often called the T 1 time), prompting consideration of the validity of treating the NV as a two-level system [11]. Furthermore, the phononlimited lifetimes of the |m s = ±1 states have not previously been investigated, leaving an incomplete picture of the NV's spin-lattice relaxation dynamics. As a result, the ultimate limit to the NV electronic spin coherence time at room temperature has remained unclear.
Numerous experimental measurements of the lifetime of the |m s = 0 state in a range of samples with varying defect concentrations have indicated that spin-lattice interactions are the primary mechanism for spin relaxation at room temperature in high-purity bulk diamond samples [12]. Specifically, two-phonon Raman transitions with major contributions from acoustic and quasilocalized phonons are understood to account for the experimentally observed temperature dependence of the rate of relaxation out of |m s = 0 near room temperature [12][13][14][15]. Prior work on phonon-induced relaxation has focused on the lifetime of the |m s = 0 state, however, neglecting possible relaxation on the |m s = −1 ↔ |m s = +1 transition. Recent work with near-surface NVs has demonstrated fast relaxation on the |m s = −1 ↔ |m s = +1 transition, with relaxation rates of up to approximately 2 × 10 3 s −1 and 2 × 10 5 s −1 observed in shallow NVs [16] and NVs in nanodiamond [17] respectively. In both contexts this relaxation is attributed to electric field noise emanating from the diamond surface. Deep in bulk diamond, far less electric field noise can be expected, leaving open the question of the phonon-induced relaxation rate of the |m s = −1 ↔ |m s = +1 transition as well as the role of other impurities in the diamond lattice.
In this Letter we present experimental results indicating that phonon-induced relaxation on the |m s = −1 ↔ |m s = +1 transition occurs at roughly twice the rate of relaxation on the |m s = 0 ↔ |m s = ±1 transition in ambient conditions, setting an ultimate limit on the maximum achievable NV electronic spin coherence time of 6.8(2) ms for a superposition of |m s = 0 and |m s = ±1 at room temperature (295 K) in high-purity bulk diamond. We discuss the importance of considering the full picture of the NV qutrit, rather than the simplified qubit model. In addition, we show that the standard theoretical model of two-acoustic-phonon Raman processes invoked to explain NV spin-lattice relaxation, which is thought to significantly contribute to relaxation on the |m s = 0 ↔ |m s = ±1 transition, predicts a rate of zero for the |m s = −1 ↔ |m s = +1 transition as a consequence of time-reversal symmetry. This suggests that two-quasilocalized-phonon Orbach-like processes or processes involving a spin-phonon Hamiltonian to second order in the atomic displacements are a dominant source of spin-phonon relaxation on the |m s = −1 ↔ |m s = +1 transition at room temperature. Future studies of re-laxation on the |m s = −1 ↔ |m s = +1 transition as a function of temperature could resolve the relative contributions of these different mechanisms.
As a matter of terminology, we note that the existing literature refers to the same transitions with various names. The |m s = −1 ↔ |m s = +1 and |m s = 0 ↔ |m s = ±1 transitions have been referred to as doublequantum and single-quantum transitions respectively, but for mixed eigenstates these terms may be misleading. The names of the transitions were generalized to account for this in Ref. [17], and we will adopt that naming scheme here. Specifically, with |H; i denoting the energy eigenstate with majority component |m s = i , we call the |H; −1 ↔ |H; +1 and |H; 0 ↔ |H; ±1 transitions the qutrit and qubit transitions respectively. In referring to the eigenstates going forward, we will omit the H identifier from kets, referring to |H; i simply as |i .

METHODS
Experiments were conducted under ambient conditions using a homebuilt confocal microscope with a temperature stable to 295 K within ±1 K. Applied magnetic fields were typically < 100 G, with few measurements conducted at higher fields. (Table S1 of the Supplemental Material [18] contains a complete summary of the data collected for this paper, including applied magnetic fields.) The components of the applied magnetic field are defined as on-and off-axis with respect to the NV's spatial axis of symmetry. Spin polarization and readout were accomplished with approximately 1 mW of 532-nm light.
To extract the qubit and qutrit relaxation rates, defined as Ω and γ respectively as in Fig. 1(a), we use stateselective π-pulses to measure the population decay into and out of each spin state individually, as described in [16,17]. Under a classical three-level population model, subtraction of these curves yields single-exponential decays from which we isolate the relaxation rates for each transition. In particular, with P i,j (τ ) denoting the population in |j at time τ after initialization in |i , we determine Ω and γ using the differences defined where a is the fluorescence contrast between the two subtracted data sets. Fig. 1(c) shows a representative result of the single-exponential decay curves F Ω and F γ . In ensemble measurements, where multiple NVs with different orientations contribute to the measured signal, the applied magnetic field and the subtraction of population curves isolates the dynamics of the single NV orientation resonantly addressed by the microwaves from among the four orientations allowed by the diamond's tetrahedral lattice.
Relaxation measurements were performed on native NVs in three chemical-vapor-deposition (CVD) grown bulk diamond samples from different growers with different native NV concentrations. Confocal microscope images of the three samples can be seen in Fig. 2(a-c) ) is described by a single exponential with decay rate 3Ω (cyan curve). Relaxation out of |+1 (yellow diamonds) is described by a biexponential decay that depends on both Ω and γ (yellow curve). Solid lines are fits to population dynamics equations of the ground state [17]. (The deviation of the yellow curve from 1 at τ = 0 is due to π-pulse infidelity.) For data shown, Ω = 56(3) s −1 , γ = 132(11) s −1 . (c) Semilog plot of subtracted population curves, FΩ and Fγ in Eqs. 1 and 2, used to extract Ω and γ respectively. Displayed data is from the same measurement series as that shown in (b).
concentration of approximately 10 −3 ppb. The NV concentrations of samples A and B are low enough to allow for measurements to be performed on single NVs using confocal microscopy. The NV concentration of Sample C, from Chenguang Machinery & Electric Equipment Company, is high enough to prohibit measurement on single NVs in our confocal microscope, at approximately 0.1 ppb. For this sample we performed measurements on the ensemble of NVs within the confocal volume of the microscope. In the data, NVs are indexed first by sample and then numerically to distinguish single NVs where appropriate. Using the ratio of approximately 300 substitutional nitrogen defects per NV − determined in [19] for native defects in CVD-grown (100)-oriented samples, we estimate the substitutional nitrogen concentrations to be approximately 3 × 10 −3 ppb, 3 × 10 −1 ppb, and 3 × 10 ppb in samples A, B, and C respectively.

RESULTS
The measured relaxation rates γ and Ω under small off-axis magnetic fields (quantitatively, B ⊥ < 1 G) for native NVs deep in the bulk in the three samples are shown in Fig. 2(d). We find that neither γ nor Ω depends on the NV concentration over a range covering four orders of magnitude from approximately 10 −5 -10 −1 ppb. This indicates that the relaxation we observe is not the result of interactions between defects. In the small off-axis magnetic field regime, we find that relaxation on the qutrit transition is approximately twice as fast as relaxation on the qubit transition. Specifically, we measure the weighted averages γ = 117(5) s −1 and Ω = 59(2) s −1 . Importantly, this indicates that significant dynamics are missed if the NV is treated as a qubit by considering only one of the |0 and |±1 subspaces; population leakage to the third state actually occurs faster than relaxation between the states intended to comprise the qubit. The standard T 1 relaxation experiment that is most commonly used with NVs consists of measuring the lifetime of the |0 state. The T 1 times quoted in the literature based on these measurements can be related to our values of Ω by defining T (0) 1 = 1/3Ω [16]. As two examples, the T (0) 1 times reported from both Naydenov et al. [10] and Bar-Gill et al. [11] in isotopically pure diamond samples at room temperature corresponds to Ω = 56(4) s −1 , in agreement with our measurements. To our knowledge, the qutrit relaxation rate γ has not been systematically studied for NVs deep within bulk diamond samples.
We next examine the limit to coherence time T 2 imposed by the qubit and qutrit relaxation rates. We introduce, T 2,max , the maximum theoretically achievable coherence time due to incoherent relaxation between states. For a qubit basis formed of |0 and either |±1 [16,17], For the values of Ω and γ measured at small off-axis magnetic fields, we find T 2,max = 6.8 (2) ms. This value represents the maximum achievable coherence time at room temperature for NVs deep in high-purity bulk diamond samples if all pure dephasing sources are eliminated. As a caveat, we note that T 2,max could in principle be extended by engineering the vibrational modes of the lattice to suppress relaxation due to spin-lattice interactions.
As an example of the implications of this limit, T 2,max determines the minimum theoretically achievable sensitivity of a room temperature single NV electronic spin to magnetic fields. If we assume 100% charge state initialization and spin state polarization as well as zero dead time between measurements, the ultimate limit is given by [20,21]: where T is the total averaging time, is the reduced Planck constant, g ≈ 2 is the electronic Landé g-factor, µ B is the Bohr magneton, and C ≤ 1 describes the measurement efficiency. In our experiment C ≈ 0.02. Using the value of the qubit T 2,max determined by our measurement of γ and Ω, we calculate δB min = 69(1) pT/ √ Hz for the idealized case C = 1. As a practical example, this corresponds to the ability to detect the magnetic field from a single nuclear spin 30 nm away from the NV after 1 s of measurement time.
There has been considerable interest in using coherent superpositions of the |+1 and |−1 states for magnetic field sensing. This "double-quantum" basis eliminates noise from temperature drift and fluctuations of the electric and scaled-strain field along the NV axis while also doubling the frequency shift due to the magnetic field along the NV axis [22][23][24][25]. However, a sensor in this basis is also more susceptible to the faster qutrit decay γ. We next investigate the sensitivity of an ideal NV sensor making use of this basis, where the maximum theoretically achievable coherence time is [16] Using the same values of Ω and γ as before, we find the maximum theoretically achievable coherence time for the double-quantum basis is T 2,max = 5.7(2) ms, which is approximately 1 ms shorter than that in the "singlequantum" basis. This corresponds to a magnetic field sensitivity of δB min = 38(1) pT/ √ Hz, indicating that despite the shorter coherence time, an overall sensitivity advantage is conferred due to the factor of two increase in signal. We note that if we set γ = 2Ω in accordance with our measurements for small off-axis magnetic fields, we find T 2,max = 2/(5Ω) in the single-quantum basis and T 2,max = 1/(3Ω) in the double-quantum basis. Figure 3 shows measurements of Ω and γ for varying applied magnetic field magnitudes and orientations plotted as functions of B and B ⊥ , the on-axis and offaxis components of the magnetic field respectively. Linear regressions are shown to demonstrate correlations between the rates and the field components. As functions of B , Ω and γ display slopes of −0.01(2) s −1 /G and 0.02(5) s −1 /G respectively. As functions of B ⊥ , Ω and γ display slopes of 0.07(5) s −1 /G and 0.81(14) s −1 /G respectively. Interestingly, the slope of γ as a function of B ⊥ (Fig. 3(d)) indicates there is a significant but small correlation between the two variables, with γ accruing an approximate 1% increase per G of applied field. The analysis of directional correlations shown here provides additional information which is not typically captured in studies of phonon-limited relaxation that focus on temperature dependence. Consideration of such properties may help uncover and explain rich spin-lattice dynamics in NVs. In the discussion of maximum coherence times above and for the remainder of this work we focus on measurements conducted under small off-axis magnetic fields because superior phonon-limited coherence times are offered in this limit and because NV applications are commonly conducted with magnetic fields applied along the NV axis.

DISCUSSION
We now turn to a discussion of the origins of relaxation on the qutrit transition. In high-purity bulk diamond at room temperature, NV spin relaxation is attributed to spin-lattice interactions [12][13][14][15]. The analysis by Jarmola et al. [12] found that a sum of the form provides a good fit to experimental measurements of the rate of relaxation out of |0 as a function of tempera-ture in a range of samples. The fit parameter A 1 (S) is a sample-dependent constant, while A 2 = 2.1(6) × 10 3 s −1 and A 3 = 2.2(5)×10 −11 s −1 are coefficients weighting the relative contributions of the exponential and T 5 terms. The parameter ∆ = 73(4) meV is an empirical activation energy. Only the latter two terms in Eq. 6 provide significant contributions at temperatures T 50 K, in agreement with the earlier analysis of Ref. [14] and the work of Ref. [13], which considers a wider variety of lowtemperature relaxation phenomena. While past results uniformly emphasize the significance of the exponential and T 5 terms at high temperatures, claims regarding which of these terms is dominant at room temperature have been mixed. Using the numbers from Ref. [12], we calculate that the exponential term contributes approximately 120 s −1 and the T 5 term contributes approximately 50 The mechanism behind the exponential term has been identified as a two-quasilocalized-phonon Orbach-like process. The empirical activation energy associated with this term is roughly consistent with ab initio calculations that have revealed a band of quasilocalized phonon modes centered at 65 meV [26][27][28]. The T 5 term is attributed to a two-acoustic-phonon Raman process. The theoretical basis for this identification comes from Ref. [29], in which Walker demonstrates a two-acoustic-phonon Raman process leading to transition rates with a thermal scaling of T 5 .
The Supplemental Material [18] contains derivations of the two-acoustic-phonon Raman process in Ref. [29] and a phenomenological model of a two-quasilocalizedphonon Orbach-like process. Both derivations are based on a spin-lattice Hamiltonian taken to first order in the atomic displacements and applied to second-order in time-dependent perturbation theory. The processes may proceed by either absorption-followed-by-emission or emission-followed-by-absorption, as shown in Fig. 4. Due to the time-reversal symmetry of the spin-lattice Hamiltonian, a generic two-phonon Raman process to first-order in the atomic displacements and second order in perturbation theory only provides a nonzero contribution to γ by the inequality of these two variants of the process. Because this restriction is inconsistent with the assumptions in Ref. [29] leading to a T 5 thermal scaling of a two-acoustic-phonon Raman process (which has previously been invoked to explain the T 5 term in Eq. 6), we conclude that this theoretical model cannot explain our observation of significant relaxation on the qutrit transition.
In contrast, a phenomenological model of the effect of two quasilocalized phonons results in an Orbach-like process (the exponential term in Eq. 6) that contributes nonzero γ and Ω relaxation rates. These observations indicate that two-quasilocalized-phonon Orbach-like processes may be the dominant mechanism for relaxation on the qutrit transition at room temperature. Along with

Energy
Energy FIG. 4. Two examples of two-phonon processes driving the qubit transition. In both processes, a phonon of energy ωem is emitted and a phonon of energy ω abs is absorbed. Importantly, the emission-first process (a) is associated with a different suppressive factor than the absorption-first process (b) in perturbation theory.
the observation that γ > Ω, this suggests that contributions from quasilocalized phonons may also be dominant for relaxation on the qubit transitions. While the simplicity of the ratio γ/Ω ≈ 2 is suggestive, we are unable to say conclusively whether this is a coincidence arising from various contributions to the two rates, or whether it is because a single mechanism dominates both the qutrit and qubit relaxation rates and naturally gives rise to the ratio between the rates.
Other possible contributions from two-phonon processes arise due to the action of a spin-lattice Hamiltonian taken to second order in the atomic displacements and applied to first order in perturbation theory [30]. These contributions are not typically considered for the NV; where they have been considered, symmetry arguments such as the one discussed here have not been applied [13]. We include a brief discussion of these contributions in the Supplemental Material [18]. We leave for future work a complete investigation of the effect of the secondorder terms in the spin-lattice Hamiltonian in light of symmetry considerations. Ultimately, the difference in expected temperature scalings of the relaxation rates induced by the various processes provides a way to potentially disentangle their contributions, motivating future measurements of γ, and the ratio γ/Ω, as a function of temperature. Ab initio calculations of the coupling constants describing the spin-lattice Hamiltonian may also be critical in explaining the ratio γ/Ω ≈ 2 [31].

CONCLUSIONS
We have presented observations of spin statedependent relaxation rates under ambient conditions in high-purity bulk diamond samples. We have measured the rate of phonon-limited relaxation on the qutrit transition, γ, at room temperature. Our measurements provide a complete experimental picture of the three-level relax-ation dynamics in high-purity diamond at room temperature, indicating T 2,max = 6.8(2) ms as the limit on the maximum theoretically achievable coherence time of an NV electronic spin at 295 K. Furthermore, we have shown that the standard model of two-acoustic-phonon Raman processes predicts no relaxation on the qutrit transition, motivating measurement of the temperature dependence of γ and providing evidence that quasilocalized phonons or contributions from higher order terms in the spinphonon Hamiltonian provide dominant contributions to spin-lattice relaxation on the qutrit transition at room temperature.  Here we derive expressions for the relaxation rates due to two-phonon processes consistent with the exponential and T 5 terms in Eq. 6 of the main text. The calculation considers processes involving acoustic phonons in the bulk and quasilocalized phonons, which have together been shown to provide significant contributions to the experimentally observed temperature dependence of the lifetime of |0 at T 50 K [1,2]. In addition, we demonstrate that due to time-reversal symmetry that two-acoustic-phonon Raman processes do not contribute to γ in the standard theory while two-quasilocalized-phonon Orbach-like processes can in general contribute to both the qutrit and qubit relaxation rates.
The calculation proceeds by application of Fermi's golden rule to second order, treating the interaction between the NV center and the diamond lattice as a perturbation. The eigenbasis is simultaneously diagonal in the Hamiltonian governing the NV and the Hamiltonian governing the lattice.

NV Hamiltonian
We consider an NV ground-state triplet Hamiltonian of the form where h = 2π is the Planck constant, D gs = 2.87 GHz is the ground-state zero-field splitting, and gµ B = 2.8 MHz/G is the NV electronic spin gyromagnetic ratio. The spin-1 operators are denoted S = (S x , S y , S z ) and the magnetic field is B.

Lattice Hamiltonian
The lattice Hamiltonian for a perfect lattice (i.e. a lattice without any defects) is where i indexes the atoms comprising the lattice, and α spans the Cartesian components (x, y, z). The mass of a carbon atom, which comprises the lattice, is denoted m. The atomic momenta and equilibrium displacements are p iα and u iα . The interaction between atoms is described to second order in the displacements with isotropic coupling Φ(∆ ii ′ ), where ∆ ii ′ = R i − R i ′ is the distance between the equilibrium positions of the ith and i ′ th atoms. Though this Hamiltonian is insufficient on its own to describe quasilocalized modes, which arise due to the presence of defects in the lattice, the quasilocalized modes may be considered perturbations of existing modes within the bulk of the crystal since their frequencies fall in the band of bulk vibrational modes [3][4][5]. Quasilocalized modes will be treated phenomenologically in a subsequent section. Examination of the perfect lattice Hamiltonian reveals some important properties of acoustic phonons and demonstrates the general method by which lattice vibrations can be quantized.
We diagonalize the lattice Hamiltonian with discrete Fourier transforms over the N atoms that comprise the lattice: where k indexes the vibrational mode with wavevector k and polarization unit vector e k . The αth component of e k is denoted e kα . We may then abbreviate the wavevector space momentum and position as With −k indicating a vibrational mode with the same polarization but opposite wavevector as k, the transformation of Eqs. 3

and 4 results in the Hamiltonian
where ω k describes the energy associated with one quantum of the vibrational mode k. We can rewrite the Hamiltonian using the creation and annihilation operators b † and b, which obey Using Eqs. 8 and 9, we obtain which gives the energy of the lattice in terms of the occupation numbers of the vibrational modes. The eigenstates of the Hamiltonian are Fock, or number, states. A single quantum of a mode can be treated as a bosonic quasiparticle, the phonon. In general, the quantization of lattice vibrations described by Eq. 10 can be used as the Hamiltonian for a more complicated lattice by extending the range of k to cover vibrational modes beyond those of a perfect lattice.

Interaction Hamiltonian
A phenomenological model of the spin-lattice interaction to second order in the atomic displacements u iα is described by: where |m and |m ′ are eigenstates of the spin Hamiltonian. The λ mm ′ iα are coupling constants to first order in the atomic displacements and the λ mm ′ ii ′ αα ′ are coupling constants to second order in the atomic displacements. To avoid an unnecessary proliferation of terms, we will consider only the first-order contributions as we demonstrate how this interaction Hamiltonian may be expressed in terms of the phonon creation and annihilation operators. Excluding the second-order terms, Eq. 11 can also be written where the F β are linear combinations of the operators |m m ′ |, β indexes the combinations, and the λ βiα are coupling constants. For the NV ground-state triplet, consideration of time-reversal and spatial symmetries imposes the restrictions which may be demonstrated by either linear algebra or group theory [6]. Working from Eq. 12, the Fourier transform given by Eq. 4 yields Using Eq. 9 to expand in terms of the creation and annihilation operators, We define and obtain the spin-phonon Hamiltonian It can easily be shown that the spin-phonon Hamiltonian including second-order terms is analogously where We note that if we assume the most significant contributions to λ βk come from atoms close to the defect, then for long wavelength phonons in the context of the spin-phonon Hamiltonian to first order in the atomic displacements we can write Because the system is invariant under translation, and Eq. 21 simplifies to wherek is the unit vector along k, v s is the speed of sound in diamond, and we have used the dispersion relation for acoustic phonons ω k = v s |k|. Eq. 23 displays the proportionality λ βk ∝ √ ω k , a property which will be invoked to obtain the T 5 scaling of two-acoustic-phonon Raman processes.

Evaluation of Fermi's golden rule for Raman processes
From time-dependent perturbation theory, Fermi's golden rule allows us to calculate Γ i→f , the transition rate from an initial state |i to a final continuum state |f . To second order, where V αβ = α| V |β are the matrix elements of the perturbation, E α is the energy of |α , and |m is an intermediate state.
Treating the spin-phonon Hamiltonian as a perturbation, the eigenstates of the system are described by the tensor product of the NV Hamiltonian eigenstates and the lattice Fock states. We assume the initial Fock state occupation numbers n k are described by the Bose-Einstein distribution for temperature T , where k B is the Boltzmann constant. We consider all final states with the NV component of interest, regardless of the vibrational mode occupation numbers. That is, the phonon-induced NV transition rate is given by ms,n k ,n k ′ ,... is the transition between an initial state |m s ⊗ |..., n k , ..., n k ′ , ... to a final state |m ′ s ⊗ |..., n ′ k , ..., n ′ k ′ , ... . Because the majority of phonons at room temperature have energies much greater than the splitting between the levels of the NV ground-state triplet and the NV has no low-lying excited states, relatively few phonons satisfy the energy requirements for direct two-phonon processes. Therefore, we consider only Raman processes, in which one phonon is absorbed and another is emitted. The NV transition rate due solely to Raman processes is Evaluating Eq. 27 using the spin-phonon Hamiltonian of Eq. 19, we obtain a general expression for the transition rate associated with two-phonon Raman processes with intermediate states restricted to the multiplet containing the initial and final states: where the matrix element H msm ′ s kk ′ and H msm ′ s k are defined and is the energy difference between the NV states |m s and |m ′ s . Evaluating the squared magnitude in Eq. 28, We see that two-phonon Raman processes arise in three ways according to the three major terms in the large parenthesis in Eq. 28. The first term, corresponds to the second-order terms in the spin-phonon Hamiltonian taken to first order in perturbation theory. The mixed second term, results from the interference of the first-and second-order terms in the spin-phonon Hamiltonian. The third term, corresponds to the first-order terms in the spin-phonon Hamiltonian taken to second order in perturbation theory. For the remainder of this supplement, we explicitly consider only the third term, since this is the premise of both the exponential and T 5 terms in Eq. 6 of the main text. We will comment on the possible contributions of the other terms in a subsequent section. With this simplification, we can write This generic second-order two-phonon Raman process may be interpreted as follows. For each intermediate NV state m ′′ s , the terms in the sum in Eq. 33 represent an absorption-followed-by-emission variant (first term in the sum, depicted in Fig. 4(a) of the main text) and an emission-followed-by-absorption variant (second term in the sum, depicted in Fig. 4(b) of the main text) of the two-phonon process. As we will see, it is important that these processes are not identical, but are associated with different denominators and may interfere with one another. Eq. 33 will be the starting point for our analyses of transitions involving acoustic phonons and quasilocalized phonons, which will be treated separately.
Properties of the generic second-order two-phonon Raman process under time-reversal Before considering the effects of different types of phonons, it is helpful to examine the properties of the generic two-phonon process under time-reversal. The restriction that the spin-lattice Hamiltonian is symmetric under timereversal excludes odd powers of the spin operators from appearing in the Hamiltonian. The operator β λ βk F β is thus also symmetric under time reversal. With Θ denoting the time-reversal operator, the S z eigenstates obey Θ |m s = (−1) ms |−m s .
Thus the matrix elements H msm ′ s k follow where we have used the antiunitarity of time-reversal and taken the Hermitian conjugate. More concisely, Inspecting Eq. 33, it may seem appropriate to make the approximation where we have assumed the delta function and used the fact that the energy scale of the NV is small in comparison to the energies of typical phonons at room temperature. In this approximation, Γ ms↔m ′ s is only nonzero if the quantity is nonzero in general. We will now show that A msm ′ s kk ′ is exactly zero for the transition |m s ↔ |−m s by using the time-reversal symmetry described by Eq. 37. The demonstration will in fact apply broadly to |m s ↔ |−m s transitions within a multiplet of an integer-spin system. The terms of the sum over m ′′ s in A ms−mskk ′ may be divided into two simple cases.
Case 1, m ′′ s = 0: If m ′′ s = 0, then The opposite of this term enters A ms−mskk ′ via the intermediate −m ′′ s . Therefore, for each term associated with a particular m ′′ s , there is a cancellation that occurs with a term from −m ′′ s and the total contribution from intermediates m ′′ s = 0 is 0. Case 2, m ′′ s = 0: This is just a special case of Eq. 41. We have for which the opposite term comes from the remaining part of the m ′′ s = 0 contribution. Because A ms−mskk ′ = 0, it is clear that the approximation of Eq. 38 is insufficient to explain relaxation on |m s ↔ |−m s transitions. More specifically, nonzero relaxation rates for a two-phonon process on |m s ↔ |−m s transitions are contingent on the difference between the denominators for the absorption-followed-by-emission and emission-followed-by-absorption variants of the process. The finding A ms−mskk ′ = 0 can also be easily verified for the specific case of the NV's qutrit transition by evaluating the matrix elements explicitly. It should be noted that the same symmetry argument clearly does not hold if m s = −m ′ s , as in the case of the NV's qubit transition. A similar symmetry argument to the one presented here has historically been applied to half-integer-spin systems, where it is referred to as Van Vleck cancellation [7].

Two-acoustic-phonon Raman processes
The two-acoustic-phonon Raman process with T 5 thermal scaling presented by Walker in [8] depends on the approximation of Eq. 38, which we have just shown to be insufficient to describe relaxation on the qutrit transition. For completeness, we will demonstrate the rest of the derivation here, showing how the temperature dependence of the process arises. We restrict consideration to acoustic vibrational modes (denoted aco) and make the approximation of Eq. 38. This yields From here, our derivation will diverge slightly from the original derivation in [8]. In particular, while the original derivation was made with arguments specific to a simple cubic lattice, we will generalize those arguments for an arbitrary lattice. Transforming from wavevector space back to position space and taking the long wavelength approximation of Eq. 23 allows the matrix elements H msm ′ s k to be expressed Substitution of Eq. 44 into Eq. 43 yields We observe that B msm ′ s kk ′ depends only on the spatial characteristics of the vibrational modes k and k ′ . In the continuum limit, we may therefore reasonably replace B msm ′ s kk ′ by a spatial average B msm ′ s defined where the integrals each cover the unit sphere and γ and γ ′ each index the one longitudinal and two transverse polarizations. The unit vector n is the normal associated with the surface element dS, and the function P γ (n) returns the γth polarization unit vector associated with n. We take the continuum limit with the transformations The sums over vibrational modes become integrals over frequencies with measure provided by the density of states of acoustic phonons. We assume the Debye model density of states given by where ω D is the Debye frequency and V is the volume of the crystal. We have defined the rectangle function We obtain Evaluating the integral over ω ′ , where d = N/V is the atomic number density.
With the change of variable Eq. 55 can be rewritten Because the integral depends only weakly on temperature for T 300 K, Eq. 59 exhibits a thermal scaling of approximately T 5 at room temperature. In the limit ω D → ∞, Eq. 59 exhibits an exact T 5 scaling.
We now enumerate three different schemes by which acoustic phonons may contribute to two-phonon relaxation as alternatives to the derivation shown in this section.
Alternative 1, restrict intermediate NV states to higher-lying states beyond the ground-state triplet: Because the splittings associated with these states are much greater than typical phonon energies at room temperature, we approximate ∆ msm ′′ s ≫ ω k for any m s , m ′′ s , and ω k . This results in a scaling of T 7 . Alternative 2, neglect the spatial characteristics of the acoustic phonon coupling constants λ βk : If we do not make the approximation of Eq. 38 and instead neglect the spatial characteristics of the coupling constants λ βk (which is necessary to make the calculation tractable in this case), then H m ′ s m ′′ s k H m ′′ s msk ′ for phonons of the same energy. Now for any m s and m ′ s the rate Γ ms↔m ′ s is only nonzero by the inequality of the denominators for the absorption-followed-by-emission and emission-followed-by-absorption variants in the generic two-phonon process of Eq. 33. This assumption results in a scaling of T 3 .
Alternative 3, consider the spin-lattice Hamiltonian up to second order in the atomic displacements (see Eq. 32), which was originally proposed to explain quadrupolar nuclear spin-lattice relaxation [9]. The second-order contributions to the spin-phonon Hamiltonian taken to first order in perturbation theory exhibit a T 7 thermal scaling for acoustic phonons, which has been shown to contribute negligibly to Ω at room temperature [2]. The time-reversal symmetry argument developed in the previous section does not apply to this mechanism since there is no interference of terms by which cancellation can occur. We leave exploration of the mixed term in Eq. 32 for future work.