# Supplemental Material for “Temperature-dependent spin-lattice relaxation of
the nitrogen-vacancy spin triplet in diamond”

M. C. Cambria  Department of Physics, University of Wisconsin, Madison,
Wisconsin 53706, USA  A. Norambuena  Centro de Optica e Información Cuántica,
Universidad Mayor, Camino La Pirámide 5750, Huechuraba, Santiago, Chile  H. T.
Dinani  Centro de Investigación DAiTA Lab, Facultad de Estudios
Interdisciplinarios, Universidad Mayor, Santiago, Chile  Escuela Data Science,
Facultad de Ciencias, Ingeniería y Tecnología, Universidad Mayor, Santiago,
Chile  G. Thiering  Wigner Research Centre for Physics, P.O. Box 49, 1525
Budapest, Hungary  A. Gardill  I. Kemeny  Y. Li  Department of Physics,
University of Wisconsin, Madison, Wisconsin 53706, USA  V. Lordi  Lawrence
Livermore National Laboratory, Livermore, CA, 94551, USA  Á. Gali  Wigner
Research Centre for Physics, P.O. Box 49, 1525 Budapest, Hungary  Department
of Atomic Physics, Institute of Physics, Budapest University of Technology and
Economics, Műegyetem rakpart 3., 1111 Budapest, Hungary  J. R. Maze  Instituto
de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago,
Chile  Centro de Investigación en Nanotecnología y Materiales Avanzados,
Pontificia Universidad Católica de Chile, Santiago, Chile  S. Kolkowitz
<EMAIL_ADDRESS>Department of Physics, University of Wisconsin, Madison,
Wisconsin 53706, USA

††preprint: APS/123-QED

##  I Experimental details

Experiments were conducted using a homebuilt confocal microscope with dual
support for low-temperature (9 K to ambient) and high-temperature (ambient to
474 K) operation modes (Fig. S1). Switching between modes is accomplished with
a flip-mounted mirror that routes the optical path between the low- and high-
temperature setups. The low-temperature setup consists of an attocube
attoDRY800 closed-cycle cryostat with a 0.82 NA cryogenic objective (attocube
LT-APO/VISIR/0.82). The high-temperature setup consists of a homebuilt copper
hot plate with a 0.55 NA long working distance (LWD) objective (Edmund Optics
59-879). The front of the LWD objective is actively cooled by air convection
from a 12 V $50\times 50\times 10$ mm case fan. In order to minimize
vibrations, the case fan is mounted on a bench separate from the optical table
and airflow is directed to the objective via a hose attached to the fan. An
apertured aluminum plate is used to isolate the hot plate from the airflow. In
our experience, an objective which has a shorter working distance or is not
actively cooled can reach high temperatures and sustain permanent damage as a
result. The temperature of the hot plate is stabilized by a custom
proportional-integral-derivative (PID) feedback loop. In both operation modes,
the diamond sample is thermally anchored to the cold finger/hot plate using
conductive silver paint (Ted Pella Leitsilber). To accurately report the
temperature of the diamond sample in high-temperature operation mode, we use
the NV as a temperature sensor by inverting the temperature dependence of the
zero field splitting reported in Ref. [1]. For low-temperature operation mode,
we calibrated the cryostat’s integrated sample temperature sensor against an
additional temperature sensor (Lakeshore CX-1050-SD-HT-1.4L-QT) mounted in the
same manner as a diamond sample.

Figure S1:  Diagram of the experimental apparatus supporting spin-to-charge-
conversion readout and dual low- and high-temperature operation modes.

As stated in the main text, the experimental sequence used to extract the
relaxation rates is the same as that used in prior works [2, 3, 4]. Briefly,
the sequence consists of the following steps: (1) optical polarization into
$m_{s}=0$, (2) a state-selective $\pi$ pulse to transfer population into a
target initialization spin state, (3) some relaxation time $\tau$, (4) a
state-selective $\pi$ pulse to transfer population from a target readout spin
state to $m_{s}=0$, and (5) optical readout of the population in $m_{s}=0$.
This procedure is repeated for different values of $\tau$ and
initialization/readout states, and the relaxation rates are extracted from the
resultant population curves. Further details are available in the above-
mentioned works [2, 3, 4].

Optical polarization is achieved using a 1 ms pulse of approximately 1.5 mW of
532 nm light.The population in $\ket{0}$ is readout optically by collecting
fluorescence from the phonon sideband above 638 nm. Microwaves are delivered
using a coplanar waveguide (CPW) mounted on top of the diamond sample. The CPW
has an $\Omega$-shaped loop along the signal trace which is shorted to the
ground trace. The loop features an 1 mm diameter aperture in the center to
allow for optical access to the sample underneath. It was manufactured by
PCBWay and consists of a single 0.5 oz gold-plated copper layer on a 25 micron
polyimide substrate. A coaxial connection is provided by a coplanar MMCX end-
launch soldered onto pads at the end of the CPW.

Data from two different chemical-vapor-deposition grown samples is presented.
The NVs in both samples are formed naturally during diamond growth. Sample A,
from Chenguang Machinery & Electric Equipment Company, displays an NV
concentration of approximately 1 ppb. Measurements in this sample were
performed on the ensemble of NVs within a single confocal volume, which we
estimate contains about 30 (150) NVs of a given orientation for low-
temperature (high-temperature) operation mode. Ensembles in high-temperature
operation mode contain more NVs than ensembles in low-temperature operation
mode due to the difference in NA between the low- and high-temperature mode
objectives. Sample B, from Element Six, displays an inhomogeneous NV
concentration with different areas exhibiting concentrations between $10^{-5}$
and $10^{-3}$ ppb. The higher-concentration regions are localized in
irregularly-spaced planes, each approximately 5 microns thick. Measurements in
sample B were conducted with single NVs isolated within the higher-
concentration planes. The NV concentrations quoted in the main text are
estimated by NV counting or by comparing the observed fluorescence rate to
that of a single isolated NV. In both samples, measurements were conducted
with NVs several microns below the diamond surface so as to avoid the effects
of surface noise.

In several instances, we use a modified version of the standard relaxation
rate experiment sequence described in prior works [2, 3, 4] in order to
achieve more efficient measurements. For all measurements in sample A (the
higher NV concentration sample), we use Knill composite $\pi$ pulses in place
of standard $\pi$ pulses to drive spin state transitions. Knill pulses consist
of five standard $\pi$ pulses with differing phases (30°, 0°, 90°, 0°, 30°)
applied in series. Knill pulses, like other composite $\pi$ pulses, are more
robust against off-resonance errors than standard $\pi$ pulses and so help to
mitigate the effect of inhomogeneous broadening within ensembles [5, 6]. For
measurements in sample B (the lower NV concentration sample) below room
temperature, we use spin-to-charge-conversion (SCC) readout in place of
standard optical readout. In the standard readout scheme, NVs are excited with
green illumination (${\sim}1.5$ mW) and the resultant spin-state dependent
fluorescence is collected. Because the spin repolarizes during this process,
readout is limited to roughly 300 ns during which $<1$ photon is collected on
average per experiment for single NVs. Each experiment must therefore be
repeated $>10^{5}$ times for sufficient averaging, leading to prohibitively
long single-NV measurements at low temperatures. In contrast, in SCC readout a
spin-state selective optical ionization pulse (${\sim}30$ mW for ${\sim}100$
ns at 638 nm) maps the NV spin state onto its charge state, which is then
readout via fluorescence under weak (${\sim}10$ $\upmu$W) yellow illumination
[7]. The charge state is robust under the readout illumination, enabling long
readouts (${\sim}10$ ms) during which tens of photons may be collected. For
long experiments where the overhead of SCC readout is relatively small, SCC
readout results in an enhanced measurement sensitivity that allows for the
practical measurement of slow relaxation rates with single NVs. We note that
unlike in the original description of SCC readout [7], we do not use a
shelving pulse, as we found the shelving pulse to have no discernible effect
on readout fidelity during characterization. In order to verify that the
modifications discussed in this paragraph (Knill pulses and SCC readout) do
not introduce errors or artifacts into our measurements, we ran relaxation
experiments at room temperature with standard and SCC readout as well as with
standard $\pi$ pulses and Knill $\pi$ pulses and confirmed in each case that
the modified sequences yield the same relaxation rates.

The assumption of a single-valued $\Omega$ to describe phonon-limited
relaxation on both the $\ket{0}\leftrightarrow\ket{-1}$ and
$\ket{0}\leftrightarrow\ket{+1}$ transitions has been verified empirically in
previous work [4] and is consistent with all experimental data collected for
this study.

##  II Fit parameters for the proposed model

The fit parameters for the proposed model in the phonon-limited regime are
presented in Table S1. The sample-dependent constants extracted from the same
fit are presented separately in Table S2.

$A_{1}$ (s-1) |  $B_{1}$ (s-1) |  $\Delta_{1}$ (meV) |  $A_{2}$ (s-1) |  $B_{2}$ (s-1) |  $\Delta_{2}$ (meV)
---|---|---|---|---|---
$0.58(8)\text{\times}{10}^{3}$ |  $1.51(17)\text{\times}{10}^{3}$ |  68.2(17) |  $9(2)\text{\times}{10}^{3}$ |  $4.8(14)\text{\times}{10}^{3}$ |  167(12)
Table S1: Fit parameters for $\Omega$ and $\gamma$ according to the proposed
model described by Eqs. (4) and (5) of the main text. Uncertainties are
$1\sigma$.

##  III Comparison to prior measurements

Figure S2:  Comparison of prior experimental measurements of the single- and
double-quantum relaxation rates (panels a and b respectively). Data is
reproduced from Refs. [8, 9, 10, 11, 4, 12].

Several experimental studies have measured the temperature dependence of the
single-quantum relaxation rate in bulk diamond samples [8, 9, 10]. In a recent
study [12], Lin et al. report measurements of the temperature dependence of
both the single- and double-quantum relaxation rates above room temperature.
Measurements of the single- and double-quantum rates are reported in Refs.
[11] and [4], but only at room temperature. The results from these prior
studies are shown alongside the results presented in this work in Fig. S2. Our
measurements of the single-quantum relaxation rate are consistent with the
majority of prior results in low-density samples ($[\text{NV}]<1\text{ ppm}$)
to within approximately 10% [9, 10, 11, 4]. Interestingly, in the high-density
samples ($[\text{NV}]\sim 10\text{ ppm}$) measured by Redman et al. and
Jarmola et al. (sample S2), single-quantum relaxation rates were approximately
30% faster than the typical values seen in low-density samples [8, 10].
Jarmola proposed that this discrepancy is due to temperature-dependent cross-
relaxation, an effect which is not captured by the sample-dependent constant
terms in the model described by Eqs. (4) and (5) of the main text.

The single-quantum relaxation rates measured by Lin et al. in Ref. [12] are
consistent with the rates measured by Redman and Jarmola in high-density
samples, but surprisingly, the measurements in Ref. [12] were conducted using
a low-density sample with a reported NV concentration of only 2 ppb. Lin’s
observation of fast single-quantum relaxation in a low-density sample may
indicate that the temperature-dependent cross-relaxation proposed by Jarmola
involves a separate defect species; in this case NV concentration would not be
a reliable indicator of the magnitude of cross-relaxation. In addition, Lin’s
measurement of the double-quantum relaxation rate at room temperature for low-
density samples is around 70% faster than the value established by previous
work [11, 4] and replicated in this work. The temperature dependence of the
double-quantum relaxation rate reported by Lin is also qualitatively different
than that reported in this work; in Lin’s measurements, the double-quantum
relaxation rate increases with temperature more slowly than the single-quantum
relaxation, such that $\Omega>\gamma$ at $T>450\text{ K}$. In contrast we find
$\Omega<\gamma$ at all temperatures in the phonon-limited regime. It is
difficult to reconcile this observation with our findings, even allowing for a
temperature-dependent cross relaxation process involving a high concentration
of dark defects in Lin’s sample. We conclude that the results of Ref. [12] are
not entirely unprecedented, but are probably not applicable to the vast
majority of samples studied by the NV community.

Figure S3:  Extended version of Fig. 4 from the main text. Sample B (dashed
lines) exhibits shorter maximum relaxation-limited coherence times than sample
A (dotted lines) due to the high sample-dependent constant for the double-
quantum rate $\gamma$ in Sample B. As in the main text version of this figure,
the data points show the longest experimentally measured coherence times in
the single-quantum subspace at various temperatures as reported in Refs. [13,
14, 15].

##  IV Sample-dependent behavior at low temperatures

Figure S4:  Detail views of sample-dependent behavior at low temperatures for
the single- and double-quantum relaxation rates (panels a and b respectively).

At temperatures below 125 K, sample-dependent constants begin to make non-
negligible contributions to the relaxation rates. The values of these
constants as determined from the fit of the proposed model (Eqs. (4) and (5)
in the main text) to the experimental data are shown in Table S2. We note that
sample-dependent constants have been attributed to interactions between
defects [10, 16]. Accordingly, it may be expected that samples with higher
defect concentrations exhibit larger sample-dependent constants. The samples
studied for this work do not display a significant difference between their
respective constants for $\Omega$. Interestingly, we observe that the constant
for $\gamma$ is around 5 times larger in sample B ($[\text{NV}]\approx
10^{-3}\text{ ppb}$) than in sample A ($[\text{NV}]\approx 1\text{ ppb}$).
While a full investigation of this effect is out of scope for the present
work, we suggest that the unexpectedly large constant for $\gamma$ in sample B
may be associated with a dark defect species present at higher levels in
sample B than in sample A. Fig. S3 shows the effect of the sample-dependent
constants on the relaxation-limited coherence times for samples A and B. We
see that coherence times in sample A can in principle be over three times
longer than in sample B at low temperatures as a result of sample-dependent
relaxation.

$A_{3}\ (\text{Sample A})$ (s-1) |  $B_{3}\ (\text{Sample A})$ (s-1) |  $A_{3}\ (\text{Sample B})$ (s-1) |  $B_{3}\ (\text{Sample B})$ (s-1)
---|---|---|---
0.013(8) |  0.06(2) |  0.010(5) |  0.30(6)
Table S2: Sample-dependent constants from the fit of the model described by
Eqs. (4) and (5) of the main text to the experimental data. The sample-
dependent constant for $\Omega$ ($\gamma$) is denoted by $A_{3}$ ($B_{3}$).
Figure S5:  Normalized residuals for the proposed model (Eqs. (4) and (5) of
the main text) and the prior model (Eqs. (S2) and (S3) of this document). For
clarity, the extent of each plot is limited to $\pm 2.5\sigma$; some outliers
are not visible. Histograms of the residuals shown at the edges of the plots
are consistent with a normal distribution of mean 0 and variance 1 (solid
lines), indicating both models provide excellent fits to the data.

##  V Comparison between the proposed and prior analytical models

In prior work [8, 10, 16], the analytical model used to describe the
temperature dependence of the single-quantum relaxation rate $\Omega$
typically contains terms accounting for low-temperature behavior, an Orbach-
like term which scales with the occupation number
$n=[\mbox{exp}(\Delta/k_{\text{B}}T)-1]^{-1}$ at a phonon energy $\Delta$, and
a $T^{5}$ term attributed to the effect of low-energy weakly coupled acoustic
phonons. Here we develop an extended version of this prior model so that we
can compare it to the model we propose in the main text (Eqs. (4) and (5)). As
discussed in the main text, prior experimental work has tacitly assumed
$\gamma=0$ [8, 9, 10]. In these contexts the lifetime of $\ket{0}$ is denoted
$T_{1}$, which is equivalent to $1/3\Omega$. For clarity we denote this
quantity $T_{1}^{(0)}$ here. As a starting point, we consider the expression
used by Jarmola et al. in Ref. [10]:

$\displaystyle\frac{1}{T_{1}^{(0)}}=A_{1}(S)+\frac{A_{2}}{e^{\Delta/k_{\text{B}}T}-1}+A_{3}\,T^{5}.$
(S1)

We note that the occupation number $n$ is written in exponential form. As the
Orbach-like term describes scattering of phonons at a specific energy
$\Delta$, the full temperature dependence of this term is $n(n+1)$ [4].
Jarmola et al. simplify this by taking the low-temperature approximation
$n(n+1)\approx n$ for $n\ll 1$; we avoid this approximation here. With these
considerations, we extend Eq. (S1) to describe both the single- and double-
quantum relaxation rates $\Omega$ and $\gamma$:

$\displaystyle\Omega(T)=A_{1}\,n(n+1)+A_{2}\,T^{5}+A_{3}(S),$ (S2)
$\displaystyle\gamma(T)=B_{1}\,n(n+1)+B_{2}\,T^{5}+B_{3}(S).$ (S3)

Fits to the experimental data according to Eqs. (S2) and (S3) are nearly
identical to those provided by the proposed model described by Eqs. (4) and
(5) of the main text and shown in Fig. 2 of the main text. The normalized
residuals for the two models are shown in Fig. S5. We see that both the
proposed model (Eqs. (4) and (5) of the main text) and the prior model (Eqs.
(S2) and (S3)) provide excellent fits to the experimental data.
Quantitatively, the reduced chi-squared metric for the proposed model is
$\chi_{\nu}^{2}=1.25$, slightly better than that for the prior model,
$\chi_{\nu}^{2}=1.29$. As such, in the temperature range studied in this work
it is not possible to say which model is more physically accurate from the
fits alone. However, the terms $T^{5}$ and $n(n+1)$ exhibit dramatically
different scalings at higher temperatures than those accessed in this or prior
studies. Accordingly, it may be possible to demonstrate the breakdown of the
prior model at temperatures of around 700 K, where the prior model predicts
relaxation rates around 50% (20%) faster than the prior model for $\Omega$
($\gamma$).

##  VI Comparison between versions of the proposed model with different
numbers of effective phonon modes

Figure S6:  Semi-log plots of fits to the experimental data for versions of
the proposed model with one, two, and three modes (panels a, b, and c
respectively). For clarity, the sample-dependent constants are removed from
the plotted functions after the fits are performed and the plots are
restricted to the phonon-limited regime above 125 K.

The analytical model we propose (Eqs. 4 and 5 of the main text) assumes that
the NV spin-phonon spectral function is best approximated by a sum of two
strongly coupled effective phonon modes. It is sensible to ask whether a
similar model consisting of just one mode provides an adequate fit to the
experimental data, or if a model with more modes provides a superior fit. Fig.
S6 demonstrates the fits provided by models consisting of one, two, and three
effective phonon modes. We see that the the one-mode model (panel a) does not
provide an adequate fit to the data, as reflected by its poor reduced chi-
squared statistic of $\chi_{\nu}^{2}=3.88$. The one-mode model identifies an
activation energy of 80.5(6) meV, which lies between the energies identified
by the two-mode model, 68.2(17) and 167(12) meV. The two-mode model provides
an excellent fit to the data (panel b), with a reduced chi-squared of
$\chi_{\nu}^{2}=1.32$. Adding a third mode does not improve the fit in a clear
visual way (panel c), and improves the reduced chi squared only marginally to
$\chi_{\nu}^{2}=1.25$. While the activation energies in the two-mode version
of the model are consistent with our own ab initio calculations and prior
spectroscopic results, the activation energies in the three-mode version of
the model less well-motivated, at 55(14) meV, 84(13) meV, and $2.5(5)\times
10^{2}$ meV. We note that the highest activation energy here is nonphysical,
lying well above the maximum phonon energies in diamond.

Figure S7:  Ratio of $\gamma$ to $\Omega$ as a function of temperature in the
phonon-limited regime. The lines are the ratios of the fits to the measured
relaxation rates $\gamma$ and $\Omega$ according to the proposed model
described by Eqs. (4) and (5) of the main text. As in Fig. 2 of the main text,
the dotted (dashed) line includes sample-dependent constants for sample A (B).
For both samples, the ratio passes near 2 at room temperature, as observed in
Ref. [4].

##  VII Ratio of $\gamma$ to $\Omega$ as a function of temperature

In Ref. [4], we observed that $\gamma\approx 2\Omega$. We previously
speculated that this factor of 2 may be fundamental, implying that the ratio
should be independent of temperature in the phonon-limited regime. Fig. S7
demonstrates that this is not the case, as the ratio decreases nearly linearly
from approximately 2.5 at 150 K to approximately 1.5 at the highest
temperatures accessed in this work. The decreasing ratio indicates that the
single-quantum transition is more strongly coupled to the higher-energy
effective phonon mode than is the double-quantum transition.

##  VIII Ab initio calculation of the spin-lattice relaxation rates

We proceed with the ab initio calculation of the spin-lattice relaxation rates
by first computing the spin-phonon matrix elements up to second order, then
applying Fermi’s golden rule, and finally taking the continuum limit,
converting sums over matrix elements into integrals over spectral functions.
We note that we apply Fermi’s golden rule with the random phase approximation
such that interference terms are dropped.

###  VIII.1 Calculation of the spin-phonon matrix elements

In this section we describe our computational methods for evaluating the spin-
phonon matrix elements for the NV center. While the present discussion is
limited to the NV center for specificity and clarity, we emphasize that our
approach can be easily generalized to other systems with zero field
splittings. While the particular spin and symmetry of the NV center limit the
number of matrix elements that must be calculated here, more complicated
systems can also be treated using the same techniques at the expense of added
computational complexity. Similarly, our discussion assumes that the zero
field splitting arises from the spin-spin interaction, but analogous
calculations may be carried out for systems with zero field splittings that
results from the spin-orbit interaction.

We calculate the spin-phonon matrix elements using plane wave supercell
density functional theory (DFT). We apply the VASP implementation [17, 18] of
the plane wave supercell DFT calculations with the projector augmented wave
(PAW) method. The supercell consists of a 512-atom simple cubic diamond with a
single negatively charged NV defect embedded in the center. We keep the
optimized diamond lattice constant (3.567 Å) fixed during the geometry
optimization procedure of the defective supercell. The forces acting on the
atoms are lower than $10^{-4}$ eV/Å in the global energy minimum of the
adiabatic potential energy surface (APES). We apply $\Gamma$-point sampling of
the Brillouin-zone. We use a cutoff of 370 eV for the plane wave basis and a
cutoff of 740 eV for charge augmentation in all calculations. The phonons
together with the associated normal coordinates are calculated by building up
the Hessian matrix as the first numerical derivative of the forces acting on
the atoms, where the atoms are moved by $0.015$ Å in each direction and the
resultant APES is fit to a parabola around the global energy minimum. The zero
field splitting tensor (or $D$-tensor) is calculated within the PAW-method
[19] as implemented by Martijn Marsman. In the global energy minimum of the
APES, the $D$-constant is equal to $(3/2)D_{zz}$ after diagonalization of the
$D$-tensor. Spin-polarized Perdew-Burke-Ernzerhof (PBE) functionals [20] are
employed in these procedures.

The spin-phonon matrix elements are then obtained by exploiting the dependence
of the $D$-tensor on the normal coordinates:

$\overleftrightarrow{D}(\bm{R})=\overleftrightarrow{D}(\bm{R}=0)+\sum_{i}\frac{\partial\overleftrightarrow{D}}{\partial
R_{i}}\biggr{|}_{\bm{R}=0}R_{i}+\frac{1}{2}\sum_{ij}\frac{\partial\overleftrightarrow{D}}{\partial
R_{i}\partial R_{j}}\biggr{|}_{\bm{R}=0}R_{i}R_{j}\text{,}$ (S4)

where a homebuilt script is used to extract the coefficients in Eq. (S4). In
order to evaluate the second-order derivatives, we consider only the diagonal
terms which satisfy $i=j$ and distort the $C_{3v}$ symmetric atomic positions
by all degenerate $e_{x}$, $e_{y}$ phonon modes of the supercell by
$\sqrt{(\Delta\bm{R})^{2}}=0.1\;\text{\AA}\sqrt{\text{a.m.u.}}$. The second-
order spin-flipping matrix elements $V_{+0}^{ll}$ and $V_{+-}^{ll}$ then
determine the $D$-tensor according to the symmetry-adapted expression:

$\begin{split}\overleftrightarrow{D}(\bm{R})=&\overset{{\textstyle\overleftrightarrow{D}(R=0)}}{\overbrace{\left(\begin{array}[]{ccc}\\!-\frac{1}{3}D\\\
&\\!\\!\\!\\!\\!\\!-\frac{1}{3}D\\\
&&\\!\\!\\!\\!\\!+\frac{2}{3}D\end{array}\right)}}+\text{linear\;terms}+\left(\begin{array}[]{ccc}\\!-\frac{1}{3}\\\
&\\!\\!\\!\\!-\frac{1}{3}\\\
&&\\!\\!\\!+\frac{2}{3}\end{array}\right)\sum_{i}3V_{00}^{ii}{R}_{i}^{2}+\\\
+&\sum_{l}V_{+-}^{ll}\left[\left(\begin{array}[]{ccc}1\\\ &\\!\\!-1\\\
&&\end{array}\right)(X_{l}^{2}-Y_{l}^{2})+\left(\begin{array}[]{ccc}&1\\\ 1\\\
&&\end{array}\right)2X_{l}Y_{l}\right]\\\
+&\sum_{l}\sqrt{2}V_{+0}^{ll}\left[\left(\begin{array}[]{ccc}&&1\\\ &\\\
1\end{array}\right)(X_{l}^{2}-Y_{l}^{2})+\left(\begin{array}[]{ccc}\\\ &&1\\\
&1\end{array}\right)2X_{l}Y_{l}\right]+...\end{split}$ (S5)

where $R_{i}$, $X_{l}$ and $Y_{l}$ are the dimensionless coordinates (not
normal coordinates) for the phonon mode at energy $\hbar\omega_{i}$ or
$\hbar\omega_{l}$. We note that while the index $l$ only covers the $e$ modes
once, the index $i$ covers all $a_{1}$, $a_{2}$, $e_{x}$, $e_{y}$ modes and
thus runs over the $e$ modes twice. Therefore, we omit $e_{x}$’s degenerate
partner $e_{y}$ modes from our calculations because their effect is the same.
Eq. (S5) can be transformed into the spin-phonon interaction $\hat{V}$ by:

$\begin{split}\hat{V}=\overleftarrow{S}\overleftrightarrow{D}\overrightarrow{S}=D(\hat{S}_{z}^{2}-\frac{1}{3}S(S+1))+\text{linear\;terms}+\sum_{i}3V_{00}^{ii}\bigl{(}\hat{S}_{z}^{2}-\frac{1}{3}\hat{S}(\hat{S}+1)\bigr{)}\hat{R}_{i}^{2}\\\
+\sum_{l}V_{+-}^{ll}\left[\bigl{(}\hat{S}_{x}^{2}-\hat{S}_{y}^{2}\bigr{)}(\hat{X}_{l}^{2}-\hat{Y}_{l}^{2})+\bigl{(}\hat{S}_{x}\hat{S}_{y}+\hat{S}_{y}\hat{S}_{x}\bigr{)}\bigl{(}\hat{X}_{l}\hat{Y}_{l}+\hat{Y}_{l}\hat{X}_{l}\bigr{)}\right]\\\
+\sum_{l}\sqrt{2}V_{+0}^{ll}\left[\bigl{(}\hat{S}_{x}\hat{S}_{z}+\hat{S}_{z}\hat{S}_{x}\bigr{)}(\hat{X}_{l}^{2}-\hat{Y}_{l}^{2})+\bigl{(}\hat{S}_{y}\hat{S}_{z}+\hat{S}_{z}\hat{S}_{y}\bigr{)}\bigl{(}\hat{X}_{l}\hat{Y}_{l}+\hat{Y}_{l}\hat{X}_{l}\bigr{)}\right]+...\end{split}$
(S6)

where
$\overleftarrow{S}=\bigl{(}\overrightarrow{S}\bigr{)}^{\dagger}=\bigl{(}\begin{array}[]{ccc}\hat{S}_{x}&\hat{S}_{y}&\hat{S}_{z}\end{array}\bigr{)}$.
We expand the dimensionless coordinates in terms of the phonon creation and
annihilation operators: $\hat{R}_{i}=(b_{i}^{\dagger}+b_{i})/\sqrt{2}$ and
$\\{\hat{X},\hat{Y}\\}_{l}=(b_{\\{X,Y\\}l}^{\dagger}+b_{\\{X,Y\\}l})/\sqrt{2}$.

Figure S8:  Ab initio calculation of the first-order spin-phonon coupling
coefficients (thin lines) and spectral function (thick lines) for a single NV
center in a 512 atom supercell. Like in Fig. 3 of the main text, we plot the
coupling coefficients and spectral function for the magnitudes of the matrix
elements $\hat{S}_{z}\hat{S}_{+}$ (blue), $\hat{S}_{+}^{2}$ (red), and
$\hat{S}_{z}^{2}-\tfrac{1}{3}\hat{S}^{2}$ (black). The first-order spectral
function has the same double-peaked lineshape as the second-order spectral
function. The spectral function is much smaller than the energies of the
associated phonons, indicating that first-order interactions make negligible
contributions to relaxation. As in Fig. 3 of the main text, the gray lines and
confidence intervals mark the locations of the effective phonon modes as
identified by the fit of the proposed model to the experimental data.

###  VIII.2 Relaxation due to first-order spin-phonon interactions

The first-order analog of Eq. (2) in the main text describes the rate at which
first-order interactions drive Raman transitions. According to Fermi’s golden
rule taken to second order in perturbation theory,

$\Gamma_{2(m_{s}m_{s}^{\prime})}^{(1)}(T)=\frac{2\pi}{\hbar}\sum_{ll^{\prime}m_{s}^{\prime\prime}}\left[n_{l}(n_{l^{\prime}}+1)\left|\frac{V_{m_{s}m_{s}^{\prime\prime}}^{l}V_{m_{s}^{\prime\prime}m_{s}^{\prime}}^{l^{\prime}}}{(E_{m_{s}^{\prime\prime}}-\hbar\omega_{l})-E_{m_{s}}}\right|^{2}\delta(\Delta
E_{-})+n_{l^{\prime}}(n_{l}+1)\left|\frac{V_{m_{s}m_{s}^{\prime\prime}}^{l}V_{m_{s}^{\prime\prime}m_{s}^{\prime}}^{l^{\prime}}}{(E_{m_{s}^{\prime\prime}}+\hbar\omega_{l})-E_{m_{s}}}\right|^{2}\delta(\Delta
E_{+})\right]\text{,}$ (S7)

where the $V_{m_{s}m_{s}^{\prime}}^{l}$ are the first-order spin-phonon
coupling coefficients from the first-order spin-phonon interaction
$\hat{V}^{(1)}=\sum_{lm_{s}m_{s}^{\prime}}V_{m_{s}m_{s}^{\prime}}^{l}(a_{l}^{\dagger}+a_{l})$.
As in the main text, $\Delta
E_{\pm}=E_{m_{s}^{\prime}}\pm(\hbar\omega_{l}-\hbar\omega_{l^{\prime}})-E_{m_{s}}$
is the energy difference between the final and initial states of the composite
system. Making the approximation that the energy of the NV is small in
comparison to the phonon energy,

$\Gamma_{2(m_{s}m_{s}^{\prime})}^{(1)}(T)=\frac{4\pi}{\hbar}\sum_{ll^{\prime}m_{s}^{\prime\prime}}n_{l}(n_{l^{\prime}}+1)\left|\frac{V_{m_{s}m_{s}^{\prime\prime}}^{l}V_{m_{s}^{\prime\prime}m_{s}^{\prime}}^{l^{\prime}}}{\hbar\omega_{l}}\right|^{2}\delta(\hbar\omega_{l^{\prime}}-\hbar\omega_{l})\text{.}$
(S8)

In the continuum limit we replace the coupling coefficients with the first-
order spectral function
$F_{m_{s}m_{s}^{\prime}}^{(1)}(\hbar\omega)=\sum_{l}|V_{m_{s}m_{s}^{\prime}}^{l}|^{2}\delta(\hbar\omega-\hbar\omega_{l})$
and obtain

$\Gamma_{2(m_{s}m_{s}^{\prime})}^{(1)}(T)=\frac{4\pi}{\hbar}\sum_{m_{s}^{\prime\prime}}\intop_{0}^{\infty}d(\hbar\omega)n(\omega)[n(\omega)+1]\frac{F_{m_{s}m_{s}^{\prime\prime}}^{(1)}(\hbar\omega)F_{m_{s}^{\prime\prime},m_{s}^{\prime}}^{(1)}(\hbar\omega)}{(\hbar\omega)^{2}}.$
(S9)

We estimate the first-order spectral function for a macroscopic diamond by
convolving the first-order coupling coefficients with a Gaussian of standard
deviation $\sigma=7.5\text{ meV}$:

$F_{m_{s}m_{s}^{\prime}}^{(1)}(\hbar\omega)=\int
d(\hbar\omega^{\prime})\sum_{l}\absolutevalue{V_{m_{s}m_{s}^{\prime}}^{l}}^{2}\delta(\hbar\omega^{\prime}-\hbar\omega_{l})\times\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(\hbar\omega^{\prime}-\hbar\omega)^{2}}{2\sigma^{2}})$
(S10)

The ab initio first-order coupling coefficients and spectral function are
displayed in Fig. S8. Unlike for second-order interactions, Raman transitions
involving first-order interactions are not directly related to the first-order
matrix elements, but rather depend on products of the first-order matrix
elements according to Eq. S9. Nevertheless, we can estimate the first-order
interaction contributions to relaxation by comparing the magnitude of the
spectral function to the typical phonon energies. In validation of the
argument presented in the main text, the ab initio first-order matrix elements
are on the order of MHz, much smaller than the energies of the associated
phonons ($\sim$10 THz). We calculate that first-order interactions drive
relaxation at rates several orders of magnitude slower than second-order
interactions.

###  VIII.3 Relaxation due to second-order spin-phonon interactions

From Fermi’s golden rule taken to first order in perturbation theory, the
second-order interaction
$\hat{V}^{(2)}=\sum_{ll^{\prime}m_{s}m_{s}^{\prime}}V_{m_{s}m_{s}^{\prime}}^{ll^{\prime}}(a_{l}^{\dagger}+a_{l})(a_{l^{\prime}}^{\dagger}+a_{l^{\prime}})$
from Eq. (S6) drives Raman scattering at the rate shown in Eq. (2) of the main
text, which we reproduce here:

$\Gamma_{1(m_{s}m_{s}^{\prime})}^{(2)}(T)=\frac{2\pi}{\hbar}\sum_{ll^{\prime}}\left|V_{m_{s}m_{s}^{\prime}}^{ll^{\prime}}\right|^{2}\big{[}n_{l}(n_{l^{\prime}}+1)\delta(\Delta
E_{-})+n_{l^{\prime}}(n_{l}+1)\delta(\Delta E_{+})\big{]}$ (S11)

Again assuming that the NV energies are small in comparison to the phonon
energies,

$\Gamma_{1(m_{s}m_{s}^{\prime})}^{(2)}(T)=\frac{4\pi}{\hbar}\sum_{ll^{\prime}}n_{l}(n_{l^{\prime}}+1)\left|V_{m_{s}m_{s}^{\prime}}^{ll^{\prime}}\right|^{2}\delta(\hbar\omega_{l^{\prime}}-\hbar\omega_{l})$
(S12)

We transform Eq. (S12) into its integral form:

$\Gamma_{1(m_{s}m_{s}^{\prime})}^{(2)}(T)=\frac{4\pi}{\hbar}\intop_{0}^{\infty}\intop_{0}^{\infty}d(\hbar\omega)d(\hbar\omega^{\prime})n(\omega)(n(\omega)+1)F_{m_{s}m_{s}^{\prime}}^{(2)}(\hbar\omega,\hbar\omega^{\prime})\delta(\hbar\omega^{\prime}-\hbar\omega)\text{,}$
(S13)

by introducing the second-order spectral function
$F_{m_{s}m_{s}^{\prime}}^{(2)}(\hbar\omega,\hbar\omega^{\prime})=\sum_{ll^{\prime}}|V_{m_{s}m_{s}^{\prime}}^{ll^{\prime}}|^{2}\delta(\hbar\omega-\hbar\omega_{l})\delta(\hbar\omega^{\prime}-\hbar\omega_{l^{\prime}})$.
The $\delta(\hbar\omega^{\prime}-\hbar\omega)$ delta function in Eq. (S13) can
be integrated out:

$\Gamma_{1(m_{s}m_{s}^{\prime})}^{(2)}(T)=\frac{4\pi}{\hbar}\intop_{0}^{\infty}d(\hbar\omega)n(\omega)(n(\omega)+1)F_{m_{s}m_{s}^{\prime}}^{(2)}(\hbar\omega,\hbar\omega)\text{.}$
(S14)

Following Eq. S10, we estimate the second-order spectral function for a
macroscopic diamond by convolving the second-order coupling coefficients with
a Gaussian of standard deviation $\sigma=7.5\text{ meV}$:

$\sqrt{F_{m_{s}m_{s}^{\prime}}^{(2)}(\hbar\omega,\hbar\omega)}=\int
d(\hbar\omega^{\prime})\sum_{l}\absolutevalue{V_{m_{s}m_{s}^{\prime}}^{ll}}\delta(\hbar\omega^{\prime}-\hbar\omega_{l})\times\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(\hbar\omega^{\prime}-\hbar\omega)^{2}}{2\sigma^{2}}).$
(S15)

We plot the square root of ab initio spectral function in Fig. 3. of the main
text. Ignoring degeneracies, we note that the
$\hbar\omega^{\prime}=\hbar\omega$ constraint also enforces that
$l=l^{\prime}$ and so we only consider the $V_{m_{s}m_{s}^{\prime}}^{ll}$
diagonal matrix elements in our ab initio calculation.

Figure S9:  Comparison between ab initio theory and experiment. (a) Dotted
lines show relaxation rates obtained by evaluating Eq. (3) in the main text
with the ab initio second-order spin-phonon spectral function shown in Fig. 3
of the main text. Dashed lines show fit of the analytical model (Eqs. (4) and
(5) in the main text) to the experimental data with sample-dependent constants
set to 0. (b) Ratio of the ab initio relaxation rates to the analytical model
rates. In the phonon-limited regime (gray line) the ab initio theory
underestimates the experimentally measured relaxation rates by approximately
16% for $\gamma$ and a factor of 8 for $\Omega$ at room temperature.

Fig. S9a compares the relaxation rates predicted by the ab initio theory in
comparison to the fit of the analytical model to the experimental data. We
obtain the predicted relaxation rates by evaluating Eq. (3) of the main text
using the calculated second-order spin-phonon spectral function for the
single- and double-quantum transitions shown in Fig. 3 of the main text. The
ab initio theory achieves near quantitative agreement with experiment and
captures several important features (Fig. S9b). In particular, we observe that
the theory reproduces the experimental finding that $\gamma>\Omega$. The
theory curves also closely follow the lineshapes of the experimental curves,
as evidenced by the flat ratios between the ab initio and model rates in the
phonon-limited regime.

###  VIII.4 Characteristics of the most strongly coupled phonon modes

For a generic solid-state spin system that features a zero-field splitting, it
is to be expected that the most strongly coupled phonon modes are those which
induce large displacements of the specific atoms that primarily host the
electronic wavefunction. We therefore predict that the most strongly coupled
modes are also (quasi)localized modes for defects with highly localized
electronic wavefunctions. Indeed, for the NV center we find that the most
strongly coupled phonon modes are associated with large displacements of the
three carbon atoms that support the unterminated “vacancy lobes” of the NV’s
${}^{3}\text{A}_{2}$ electronic ground state triplet wavefunction. The most
localized vibrational mode ( 18% probability on the three carbon atoms) with
doubly degenerate “E” symmetry is at 62.4 meV and exhibits the largest spin-
phonon coupling coefficients at 2 MHz for $V_{+-}^{ll}$ and 0.6 MHz for
$V_{+0}^{ll}$, as can be seen in Fig. 3 of the main text. The most localized
“optical” vibrational mode (11% localization on the three carbon atoms) with
doubly degenerate “E” symmetry is at 160.7 meV and exhibits spin-phonon
coupling coefficients of 0.34 MHz for $V_{+-}^{ll}$ and 0.07 MHz for
$V_{+0}^{ll}$. We note that it is not necessarily the case for other systems
that the energies involved in the reduction of relaxation rate integral (Eq. 3
of the main text) to a limited sum (Eqs. 4 and 5 of the main text) will be
associated with especially strongly coupled phonon modes. It may also be the
case that the energies involved in the reduction are associated with a
particularly large phonon density of states, and that there are no especially
strongly coupled modes at those energies.

Figure S10:  Comparison of the relaxation rates predicted by ab initio
calculations (light, solid lines) and fits of the proposed model to the ab
initio rates (darker, dotted lines).

###  VIII.5 Fit of the proposed model to the ab initio predicted relaxation
rates

In order to compare the salient characteristics of the experimental and ab
initio relaxation rates, it is useful to fit both data sets to the same model
so that the fit parameters may be compared. Here, we fit a version of the
proposed model without sample-dependent constants (i.e., a sum of two Orbach-
like terms) to the relaxation rates predicted by the ab initio simulations up
to 5000 K. The results are shown in Fig. S10, and the fit parameters are
tabulated to 3 significant digits in Table S3. The figure shows the ab initio
predicted rates in lightly colored solid lines (matching the presentation from
Fig. 2 of the main text) and the fit rates in darker dotted lines.
Interestingly, we see that the activation energies extracted from this
procedure are around 5-10% lower than would be expected based on the spectral
function and the fit of the proposed model to the experimental data. We
attribute this effect to the larger occupation numbers, and thus larger
contributions to the calculated relaxation rates, associated with lower-energy
features of the spectral function. Using a toy model spectral function, we
observe that the magnitude of the activation energy underestimation increases
with the width of the Gaussian used for broadening the spectral function. We
therefore anticipate that larger supercell calculations which allow for
narrower Gaussian broadening may demonstrate less bias in the activation
energies extracted from the fit.

$A_{1}$ (s-1) |  $B_{1}$ (s-1) |  $\Delta_{1}$ (meV) |  $A_{2}$ (s-1) |  $B_{2}$ (s-1) |  $\Delta_{2}$ (meV)
---|---|---|---|---|---
$0.0700\text{\times}{10}^{3}$ |  $0.910\text{\times}{10}^{3}$ |  62.5 |  $0.169\text{\times}{10}^{3}$ |  $2.94\text{\times}{10}^{3}$ |  139
Table S3: Version of Table S1 for the relaxation rates predicted by the ab
initio calculations. The ab initio rates are fit according to the model
described by Eqs. (4) and (5) of the main text without sample-dependent
constants.

##  IX Generality of the theoretical arguments and methods

The theory presented in this work can be generalized to other spin systems. In
particular, second-order interactions will dominate Raman scattering in
systems that have no low-lying excited states and feature phonons whose
energies are much higher than the strengths of their respective couplings to
the spin. We expect systems in which both the defect and lattice consist of
light elements with strong covalent bonds to commonly meet these criteria.
Other likely candidates are defects with ground state orbital singlets, due to
the absence of the spin-orbit interaction. While our computational approach is
not limited to spin-1 or $C_{3v}$ symmetric systems, we highlight that our
arguments and methods could be applied with little modification to several
other spin-1 crystal defects with similar symmetries to the NV center, such as
divacancy centers in silicon carbide [21, 22] or the boron vacancy center in
diamond [23]. In contrast, the neutrally charged silicon vacancy center in
diamond is likely not a good candidate due to the presence of a low-lying
excited state [24], and the negatively charged silicon vacancy center in
diamond and analogous defects involving elements from the same group are also
likely poor candidates due to their large spin-orbit couplings [25].

##  X Complete set of experimental data

Fig. S11 and table S4 display the complete set of experimentally measured
rates collected for this work. In the table, NVs are indexed first by sample
and then numerically to distinguish single NVs in sample B.

Figure S11:  Semi-log and log-log versions (panels a and b respectively) of Fig. 2 from the main text showing the full temperature over which data was collected for this work. The rates predicted by the ab initio calculation are not shown.  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---
NV | $T\text{ (K)}$ |  $\Omega$ (s-1) |  $\gamma$ (s-1) | NV | $T\text{ (K)}$ |  $\Omega$ (s-1) |  $\gamma$ (s-1)
NVA |  8.9 |  $1.7(9)\text{\times}{10}^{-2}$ |  $5(3)\text{\times}{10}^{-2}$ |  NVA |  390.3 |  $1.65(7)\text{\times}{10}^{2}$ |  $2.6(2)\text{\times}{10}^{2}$
NVA |  50.0 |  $2(2)\text{\times}{10}^{-2}$ |  $9(6)\text{\times}{10}^{-2}$ |  NVA |  401.6 |  $1.79(7)\text{\times}{10}^{2}$ |  $3.5(2)\text{\times}{10}^{2}$
NVA |  84.6 |  $4(2)\text{\times}{10}^{-2}$ |  $5(2)\text{\times}{10}^{-1}$ |  NVA |  415.6 |  $2.13(10)\text{\times}{10}^{2}$ |  $3.7(3)\text{\times}{10}^{2}$
NVA |  124.1 |  $1.00(9)$ |  $3.4(4)$ |  NVA |  427.1 |  $2.43(10)\text{\times}{10}^{2}$ |  $3.6(2)\text{\times}{10}^{2}$
NVA |  148.4 |  $2.64(16)$ |  $7.0(6)$ |  NVA |  440.1 |  $2.54(11)\text{\times}{10}^{2}$ |  $4.1(3)\text{\times}{10}^{2}$
NVA |  158.6 |  $4.0(3)$ |  $1.03(8)\text{\times}{10}^{1}$ |  NVA |  454.0 |  $2.62(11)\text{\times}{10}^{2}$ |  $4.3(3)\text{\times}{10}^{2}$
NVA |  160.9 |  $4.0(3)$ |  $1.21(9)\text{\times}{10}^{1}$ |  NVA |  465.5 |  $3.04(13)\text{\times}{10}^{2}$ |  $5.2(3)\text{\times}{10}^{2}$
NVA |  172.9 |  $6.2(2)$ |  $1.66(8)\text{\times}{10}^{1}$ |  NVA |  471.4 |  $3.37(15)\text{\times}{10}^{2}$ |  $5.2(3)\text{\times}{10}^{2}$
NVA |  184.9 |  $8.3(4)$ |  $2.09(13)\text{\times}{10}^{1}$ |  NVB3 |  8.9 |  $1.6(6)\text{\times}{10}^{-2}$ |  $3.4(8)\text{\times}{10}^{-1}$
NVA |  196.9 |  $1.14(5)\text{\times}{10}^{1}$ |  $2.32(16)\text{\times}{10}^{1}$ |  NVB3 |  50.0 |  $0.0(8)\text{\times}{10}^{-2}$ |  $2.3(9)\text{\times}{10}^{-1}$
NVA |  208.5 |  $1.45(7)\text{\times}{10}^{1}$ |  $3.9(2)\text{\times}{10}^{1}$ |  NVB2 |  99.3 |  $2.3(7)\text{\times}{10}^{-1}$ |  $1.1(3)$
NVA |  221.4 |  $1.91(10)\text{\times}{10}^{1}$ |  $5.2(3)\text{\times}{10}^{1}$ |  NVB2 |  148.4 |  $2.7(2)$ |  $9.1(10)$
NVA |  233.5 |  $2.23(10)\text{\times}{10}^{1}$ |  $5.7(4)\text{\times}{10}^{1}$ |  NVB2 |  196.9 |  $1.33(13)\text{\times}{10}^{1}$ |  $3.0(4)\text{\times}{10}^{1}$
NVA |  244.7 |  $2.85(13)\text{\times}{10}^{1}$ |  $7.2(4)\text{\times}{10}^{1}$ |  NVB2 |  244.7 |  $2.7(2)\text{\times}{10}^{1}$ |  $7.3(7)\text{\times}{10}^{1}$
NVA |  256.4 |  $3.54(17)\text{\times}{10}^{1}$ |  $7.9(5)\text{\times}{10}^{1}$ |  NVB1 |  295.0 |  $5.0(11)\text{\times}{10}^{1}$ |  $1.2(4)\text{\times}{10}^{2}$
NVA |  268.8 |  $4.5(2)\text{\times}{10}^{1}$ |  $9.1(7)\text{\times}{10}^{1}$ |  NVB5 |  295.0 |  $5.1(5)\text{\times}{10}^{1}$ |  $1.31(15)\text{\times}{10}^{2}$
NVA |  281.4 |  $5.4(3)\text{\times}{10}^{1}$ |  $1.05(8)\text{\times}{10}^{2}$ |  NVB4 |  295.0 |  $6.5(5)\text{\times}{10}^{1}$ |  $1.23(13)\text{\times}{10}^{2}$
NVA |  295.0 |  $6.0(3)\text{\times}{10}^{1}$ |  $1.28(7)\text{\times}{10}^{2}$ |  NVB5 |  344.3 |  $1.06(9)\text{\times}{10}^{2}$ |  $2.2(2)\text{\times}{10}^{2}$
NVA |  301.5 |  $6.0(2)\text{\times}{10}^{1}$ |  $1.32(7)\text{\times}{10}^{2}$ |  NVB4 |  344.9 |  $7.6(12)\text{\times}{10}^{1}$ |  $1.9(4)\text{\times}{10}^{2}$
NVA |  308.1 |  $7.3(3)\text{\times}{10}^{1}$ |  $1.50(10)\text{\times}{10}^{2}$ |  NVB5 |  393.0 |  $1.60(13)\text{\times}{10}^{2}$ |  $3.4(4)\text{\times}{10}^{2}$
NVA |  323.8 |  $8.7(5)\text{\times}{10}^{1}$ |  $1.30(12)\text{\times}{10}^{2}$ |  NVB4 |  393.6 |  $1.77(16)\text{\times}{10}^{2}$ |  $2.9(4)\text{\times}{10}^{2}$
NVA |  326.1 |  $7.7(3)\text{\times}{10}^{1}$ |  $1.57(10)\text{\times}{10}^{2}$ |  NVB5 |  440.1 |  $2.23(19)\text{\times}{10}^{2}$ |  $4.9(5)\text{\times}{10}^{2}$
NVA |  328.3 |  $9.6(4)\text{\times}{10}^{1}$ |  $1.59(11)\text{\times}{10}^{2}$ |  NVB4 |  441.1 |  $2.7(3)\text{\times}{10}^{2}$ |  $5.2(8)\text{\times}{10}^{2}$
NVA |  337.6 |  $1.01(5)\text{\times}{10}^{2}$ |  $1.95(13)\text{\times}{10}^{2}$ |  NVB5 |  473.3 |  $3.2(3)\text{\times}{10}^{2}$ |  $5.8(7)\text{\times}{10}^{2}$
NVA |  353.7 |  $1.14(5)\text{\times}{10}^{2}$ |  $2.27(15)\text{\times}{10}^{2}$ |  NVB5 |  473.3 |  $3.5(3)\text{\times}{10}^{2}$ |  $5.0(7)\text{\times}{10}^{2}$
NVA |  368.7 |  $1.32(6)\text{\times}{10}^{2}$ |  $2.51(17)\text{\times}{10}^{2}$ |  NVB4 |  473.5 |  $3.9(3)\text{\times}{10}^{2}$ |  $5.9(6)\text{\times}{10}^{2}$
NVA |  380.4 |  $1.38(6)\text{\times}{10}^{2}$ |  $2.8(2)\text{\times}{10}^{2}$ |  |  |  |
Table S4: Complete set of experimentally measured rates collected for this
work.

## References

  * Toyli  _et al._ [2012] D. Toyli, D. Christle, A. Alkauskas, B. Buckley, C. Van de Walle, and D. Awschalom, Measurement and control of single nitrogen-vacancy center spins above 600 k, Physical Review X 2, 031001 (2012).
  * Myers  _et al._ [2017] B. A. Myers, A. Ariyaratne, and A. B. Jayich, Double-quantum spin-relaxation limits to coherence of near-surface nitrogen-vacancy centers, Physical Review Letters 118, 197201 (2017).
  * Gardill  _et al._ [2020] A. Gardill, M. C. Cambria, and S. Kolkowitz, Fast relaxation on qutrit transitions of nitrogen-vacancy centers in nanodiamonds, Physical Review Applied 13, 034010 (2020).
  * Cambria  _et al._ [2021] M. Cambria, A. Gardill, Y. Li, A. Norambuena, J. Maze, and S. Kolkowitz, State-dependent phonon-limited spin relaxation of nitrogen-vacancy centers, Physical Review Research 3, 013123 (2021).
  * Ryan  _et al._ [2010] C. A. Ryan, J. S. Hodges, and D. G. Cory, Robust decoupling techniques to extend quantum coherence in diamond, Physical Review Letters 105, 200402 (2010).
  * Souza  _et al._ [2011] A. M. Souza, G. A. Alvarez, and D. Suter, Robust dynamical decoupling for quantum computing and quantum memory, Physical Review Letters 106, 240501 (2011).
  * Shields  _et al._ [2015] B. J. Shields, Q. P. Unterreithmeier, N. P. de Leon, H. Park, and M. D. Lukin, Efficient readout of a single spin state in diamond via spin-to-charge conversion, Physical Review Letters 114, 136402 (2015).
  * Redman  _et al._ [1991] D. Redman, S. Brown, R. Sands, and S. Rand, Spin dynamics and electronic states of N-V centers in diamond by EPR and four-wave-mixing spectroscopy, Physical Review Letters 67, 3420 (1991).
  * Takahashi  _et al._ [2008] S. Takahashi, R. Hanson, J. Van Tol, M. S. Sherwin, and D. D. Awschalom, Quenching spin decoherence in diamond through spin bath polarization, Physical Review Letters 101, 047601 (2008).
  * Jarmola  _et al._ [2012] A. Jarmola, V. Acosta, K. Jensen, S. Chemerisov, and D. Budker, Temperature- and magnetic-field-dependent longitudinal spin relaxation in nitrogen-vacancy ensembles in diamond, Physical Review Letters 108, 197601 (2012).
  * Taminiau  _et al._ [2014] T. H. Taminiau, J. Cramer, T. van der Sar, V. V. Dobrovitski, and R. Hanson, Universal control and error correction in multi-qubit spin registers in diamond, Nature Nanotechnology 9, 171 (2014).
  * Lin  _et al._ [2021] S. Lin, C. Weng, Y. Yang, J. Zhao, Y. Guo, J. Zhang, L. Lou, W. Zhu, and G. Wang, Temperature-dependent coherence properties of NV ensemble in diamond up to 600 K, Physical Review B 104, 155430 (2021).
  * Bar-Gill  _et al._ [2013] N. Bar-Gill, L. M. Pham, A. Jarmola, D. Budker, and R. L. Walsworth, Solid-state electronic spin coherence time approaching one second, Nature Communications 4, 1743 (2013).
  * Herbschleb  _et al._ [2019] E. Herbschleb, H. Kato, Y. Maruyama, T. Danjo, T. Makino, S. Yamasaki, I. Ohki, K. Hayashi, H. Morishita, M. Fujiwara, _et al._ , Ultra-long coherence times amongst room-temperature solid-state spins, Nature Communications 10, 3766 (2019).
  * Abobeih  _et al._ [2018] M. H. Abobeih, J. Cramer, M. A. Bakker, N. Kalb, M. Markham, D. J. Twitchen, and T. H. Taminiau, One-second coherence for a single electron spin coupled to a multi-qubit nuclear-spin environment, Nature Communications 9, 2552 (2018).
  * Norambuena  _et al._ [2018] A. Norambuena, E. Muñoz, H. Dinani, A. Jarmola, P. Maletinsky, D. Budker, and J. Maze, Spin-lattice relaxation of individual solid-state spins, Physical Review B 97, 094304 (2018).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Physical Review B 47, 558 (1993).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B 54, 11169 (1996).
  * Bodrog and Gali [2013] Z. Bodrog and A. Gali, The spin–spin zero-field splitting tensor in the projector-augmented-wave method, Journal of Physics: Condensed Matter 26, 015305 (2013).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Son  _et al._ [2006] N. T. Son, P. Carlsson, J. ul Hassan, E. Janzén, T. Umeda, J. Isoya, A. Gali, M. Bockstedte, N. Morishita, T. Ohshima, and H. Itoh, Divacancy in 4H-SiC, Physical Review Letters 96, 055501 (2006).
  * Koehl  _et al._ [2011] W. F. Koehl, B. B. Buckley, F. J. Heremans, G. Calusine, and D. D. Awschalom, Room temperature coherent control of defect spin qubits in silicon carbide, Nature 479, 84 (2011).
  * Umeda  _et al._ [2022] T. Umeda, K. Watanabe, H. Hara, H. Sumiya, S. Onoda, A. Uedono, I. Chuprina, P. Siyushev, F. Jelezko, J. Wrachtrup, _et al._ , Negatively charged boron vacancy center in diamond, Physical Review B 105, 165201 (2022).
  * Rose  _et al._ [2018] B. C. Rose, G. Thiering, A. M. Tyryshkin, A. M. Edmonds, M. L. Markham, A. Gali, S. A. Lyon, and N. P. de Leon, Strongly anisotropic spin relaxation in the neutral silicon vacancy center in diamond, Physical Review B 98, 235140 (2018).
  * Hepp  _et al._ [2014] C. Hepp, T. Müller, V. Waselowski, J. N. Becker, B. Pingault, H. Sternschulte, D. Steinmüller-Nethl, A. Gali, J. R. Maze, M. Atatüre, _et al._ , Electronic structure of the silicon vacancy color center in diamond, Physical Review Letters 112, 036405 (2014).