# Quantifying the magnetic noise power spectrum for ensembles of P1 and NV
centers in diamond

Ethan Q. Williams<EMAIL_ADDRESS>Chandrasekhar Ramanathan
<EMAIL_ADDRESS>Department of Physics and Astronomy,
Dartmouth College, Hanover, New Hampshire 03755, USA

###### Abstract

We use Carr-Purcell-Meiboom-Gill (CPMG) dynamical decoupling to measure the
magnetic noise power spectra for ensembles of P1 and NV centers in diamond
using pulsed electron paramagnetic resonance (pEPR) at 2.5 GHz. The
stroboscopically detected pEPR experiments on NV centers were performed on an
HPHT (high pressure, high temperature) diamond sample at 13 mT and 190 mT,
while the experiments on P1 centers were performed on a CVD (chemical vapor
deposition) diamond sample at 89 mT. All power spectra show two distinct
features, a broad component that is observed to scale as approximately
$1/\omega$, and a prominent peak at the 13C Larmor precession frequency. The
broad $1/\omega$ behavior is consistent with an inhomogeneous distribution of
Lorentzian spectra due to clustering of P1 centers, which has recently been
shown to be prevalent in HPHT diamond. However, it is unknown if such
clustering occurs in CVD diamond. The maximum rate at which we can apply $\pi$
pulses is higher than the 13C frequency at 13 mT, but is lower than the 13C
frequency at 89 mT and 190 mT. We develop techniques that utilize the higher
harmonics of the CPMG filter function to improve our estimate of the 13C
contribution to the power spectrum at the higher fields. Surprisingly, the 13C
peak, when measured with higher harmonics of the CPMG filter, appears larger
than expected based on measurements with the lower harmonics. We assess the
robustness of our methods in the presence of finite pulse widths and flip
angle errors. These techniques could be used in a variety of ac magnetometry
and noise spectroscopy measurements such as chemical sensing and nanoscale
nuclear magnetic resonance.

diamond, nitrogen vacancy, NV, P1, electron spin paramagnetic resonance, ESR,
EPR, dynamical decoupling, pulse sequence, CPMG, noise spectroscopy, 1/f
noise, quantum sensing

††preprint: APS/123-QED

##  I Introduction

Defect electronic spins in diamond are proving to be a powerful platform for
quantum information and sensing technologies. The widely-studied nitrogen-
vacancy (NV) center has enabled new measurement modalities for studying
biological systems [1, 2, 3, 4], microelectronics [5], condensed matter
physics and chemistry [6, 7, 8, 9], nanoscale nuclear magnetic resonance (NMR)
[10, 11, 12, 13, 14], quantum information [15], and fundamental physics [16,
17]. The spin states of the NV center exhibit long coherence times at room
temperature and can be both initialized and read out optically, for either
ensemble or single spin measurements. The substitutional nitrogen defect, or
P1 center, also exhibits long coherence times but cannot be initialized or
read out optically except via the NV centers under special conditions [18, 19,
20]. While the P1 centers can be a dominant source of decoherence for the NV
centers [21], they are also a promising platform for studying the many-body
physics of interacting spin systems [20, 22] and may be used as qubits for
quantum sensing and computation [19].

Bulk single-crystal diamond is a disordered system in that each defect
electron spin has a unique configuration of other electron spins and 13C
nuclear spins in its vicinity. The coherence properties of the electron spin
are determined by its interactions with these neighboring spins. Figure 1(a)
shows a cartoon of the local environment of an NV center surrounded by both P1
centers and 13C spins. The NV center coherence times measured in optically
detected Hahn echo experiments [23] are primarily limited by the concentration
of nitrogen impurities for concentrations of $[\text{N}]\gtrsim 0.5\text{
ppm}$ and are limited by interactions with 13C nuclei and $T_{1}$-related
processes at lower nitrogen concentrations [24, 25, 21].

The $T_{2}$ of the NV center has previously been used to estimate the local P1
concentration in the high-concentration regime via ensemble experiments on
isotopically engineered CVD (chemical vapor deposition) diamond [25]. However,
such scaling relationships need to be used carefully as the homogeneity of the
P1 distribution is not always uniform. For example, it has recently been shown
that high pressure, high temperature (HPHT) diamond with a reported nitrogen
concentration of [N] $\approx$ 200 ppm that was irradiated and annealed for
the creation of NV centers had local P1 concentrations varying as widely as 10
to 300 ppm [26]. Clustering of P1 spins in HPHT diamond has also been observed
via pulsed electron paramagnetic resonance (pEPR) [27, 28] and gives rise to
multiple mechanisms for dynamic nuclear polarization [29]. Characterizing the
properties of bulk samples could better inform our understanding of the
variety of spin environments found in diamond.

Figure 1:  (a) Cartoon of the central spin model, depicting the NV center as
the central spin of a disordered system. (b) Illustration of the assembly of
the diamond sample, loop-gap resonator (cutaway view), optical fiber (emitting
green 532 nm light), microwave antenna, and sapphire dielectric used for
tuning the resonator. (c) Absolute values of resonance frequencies of the
electron spin transitions targeted in this work and the 2.5 GHz operating
frequency of the spectrometer. (d)–(g) show examples of data from an echo-
detected field sweep (“EDFS”) of the NV resonance near 190 mT. (d) Example of
the time trace of a spin echo, taken at the middle of the resonance line $\sim
191.2$ mT. The two lines show the in-phase and quadrature component of the
echo, which are demodulated from their carrier frequency and phase-corrected
to maximize the real signal. (e) Image plot of stacked echo time traces
obtained at different magnetic fields. The brightness is given by the log of
the signal magnitude, and the color represents the phase of the signal. (f)
The lower pair of lines shows the Fourier transform of the echo time trace in
(d), while the upper pair of lines shows echo amplitudes (taken as the
amplitude of the middle frequency of the echo Fourier transform) as a function
of magnetic field converted to relative frequency units. (g) Image plot of
stacked Fourier transforms of echo time traces seen in (e). More echo-detected
field sweep data for the other spin species addressed in this work are shown
in Appendix B.  Figure 2:  (a) Schematic of the experiment protocol used for
NV experiments including the initializing laser pulse, CPMG microwave pulse
sequence, and the echoes that are formed. The $N^{\text{th}}$ echo occurs at
detection time $t_{d}=N\tau$. In a single shot, we acquire each echo that
forms until the signal has diminished below the noise floor. The decay of the
echo amplitudes is represented by the coherence curve. (b) Examples of
coherence decays with fixed $\tau$. (c) Examples of coherence decays with
fixed $N$. The fits for the $N=1$ (Hahn echo) curves are described in Appendix
A. The data in (b) and (c) are extracted from the same two data sets (one
acquired at 13 mT and one at 190 mT). Each coherence data point represents the
magnitude of the echo formed at that time and is normalized for each dataset
such that the exponential fits to CPMG-1 have an amplitude of 1, i.e., $L=1$
at $t_{d}=0$. In (b) the vertical green and purple lines correspond to curves
in (c), indicating where the $N=1$ and $N=16$ coherence curves intersect with
the fixed $\tau$ curves respectively. Likewise, in (c) the red and blue
vertical lines show where the $\tau=9.83\ \mu$s and $\tau=13.8\ \mu$s curves
intersect with the $N=1$ and $N=16$ coherence curves.

Dynamical decoupling (DD), first developed and implemented in nuclear magnetic
resonance experiments [30, 31], has become an important tool in characterizing
the noise environments of qubits and when using qubits for sensing ac signals
in general [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. DD-
type pulse sequences have seen considerable use with the NV center in diamond
for characterizing decoherence [47, 48, 49, 50, 12, 51, 52, 53, 54, 11, 55,
13, 56, 14] and reading and controlling the state of 13C nuclear spins [57,
58, 59].

Here we characterize the average magnetic noise of these defect centers in
bulk diamond samples using DD-based noise spectroscopy measurements with pEPR.
Stroboscopic detection during pEPR allows us to acquire the spin signal after
each microwave pulse in a single shot of a DD experiment, as in the original
NMR DD experiments [30, 31]. We characterize ensembles of defect spins in two
diamond samples, both from Element 6. Sample A is a nitrogen-rich, HPHT
diamond ([N] $\approx$ 90 ppm), which has undergone electron irradiation for
the creation of NV centers with a concentration of [NV] $\approx$ 2 ppm.
Sample B is a CVD diamond plate with low nitrogen concentration ([N] $\approx$
0.4 ppm), with which we probe P1 electron spins. Both samples have a natural
isotopic abundance of carbon. Additional sample details are provided in
Appendix A.

The DD sequence used is the Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence
[30, 31]: $(\pi/2)_{x}-(\tau/2-\pi_{y}-\tau/2)^{N}$, where $\tau$ is the time
between the midpoints of two consecutive cycle pulses, and $N$ indicates the
number of $\pi_{y}$ pulses applied. CPMG is favorable for its simple structure
and short cycle length. Although DD sequences in the XY family make pulse
errors isotropic and better protect arbitrary initial states, CPMG typically
is comparable or slightly outperforms (at extending coherence times) the XY
sequences for known initial states [48, 60]. The XY sequences can also give
rise to spurious harmonics, which, although they can be used as a resource for
enhanced spectral resolution in certain conditions [54, 11], can obscure the
spectrum with subharmonics.

We show that interactions with 13C nuclei dominate and obscure the broad
features of the NV power spectrum at 13 mT. However, by increasing the
magnetic field to 190 mT, the effect of these 13C interactions is reduced, and
the background is found to exhibit a $1/\omega$ spectrum. The $1/\omega$
behavior was also observed with P1 centers in the CVD sample. We compare this
with the $1/\omega^{2}$ high-frequency tail from a Lorentzian power spectrum
seen in previous measurements on single spins and small ensembles [47, 50, 61,
55, 56].

We also characterize the resonant 13C Larmor frequency feature for the P1
centers at 89 mT and the NV centers at 190 mT by introducing a filter
function-based protocol for measuring sharp features in the power spectrum at
frequencies that are greater than the maximum CPMG modulation frequency
accessible with the system hardware.

The manuscript is organized as follows. In Section II we describe the basic
framework for our experiments. In Section III we summarize the theoretical
model describing coherence decay of an ensemble of central spins under
application of a DD pulse sequence. Section IV is the key results section,
with Section IV.1 providing an overview of the raw experimental data, Section
IV.2 discussing the analysis of the broad background noise, and Section IV.3
discussing the characterization of the resonant 13C peak. Section V summarizes
our findings and provides an outlook to future experiments.

##  II Experiment Overview

The system Hamiltonian of an isolated spin-1/2 electron in a magnetic field
$B_{0}$ is the Zeeman Hamiltonian $H_{S}=-\omega_{0}S_{z}$. For the isolated
P1 center, $\omega_{0}=\gamma_{e}B_{0}$, where $\gamma_{e}=-28.0\times 2\pi$
rad GHz/T. For the negatively charged NV center, which is spin-1, we treat the
two lowest energy levels as an effective spin-1/2 system with the mapping
$\ket{0}\rightarrow\ket{1/2}$ and $\ket{-1}\rightarrow\ket{-1/2}$, for which
$\omega_{0}=D+\gamma_{e}B_{0}$, where $D=2.87\times 2\pi$ rad GHz is the
energy of the zero-field splitting. We ignore the small difference in
$g$-factor between the two defects. To detect the electron spins at low and
intermediate fields, we use a microwave frequency of $|\omega_{0}|=2.5\times
2\pi$ rad GHz, which provides three resonance conditions as shown in Figure
1(c). For the NV center, the resonance frequency crosses 2.5 GHz at $B_{0}=13$
mT and 190 mT. For the P1 center the resonant field occurs at 89 mT. Echo-
detected field sweep spectra of the NV and P1 centers can be found in Figure 1
and in Appendix B.

Hyperfine interactions with the host 14N nucleus (spin-1), result in three
nuclear spin manifolds, and we set our magnetic field such that our pulses are
on resonance with the central ($m_{I}=0$) manifold. For the NV center, the
three hyperfine lines are within the bandwidth of our $\sim$50 ns microwave
pulse as evidenced by the echoes and spectra in Figure 1(d)–(g). For the P1
center, this hyperfine interaction is stronger, and the $m_{I}=+1$ or $-1$
resonance lines are outside the bandwidth of the pulse. For the NV
experiments, Sample A is oriented so that the [111] axis is aligned with
$\vec{B}_{0}$. For the P1 experiments, Sample B is oriented so that the [100]
axis is aligned with $\vec{B}_{0}$. A Jahn-Teller distortion defines a
hyperfine quantization axis along one of the N-C bond directions for the P1
center, generating up to 4 spectral lines for each state $m_{I}$ of the host
nitrogen [62, 18]. In our CPMG experiments with the P1 center,
$\vec{B}_{0}\parallel[100]$, resulting in an overlap of the 4 hyperfine lines
for each nuclear spin state.

A schematic of the CPMG sequence is shown in Figure 2(a). The NV centers in
the diamond sample are initialized to the $\ket{m_{s}}=\ket{0}$ ground state
with 532 nm laser light transmitted via an optical fiber. The initialization
laser pulse is 10 ms in duration. The microwave signal from the spin ensemble
is detected inductively using a loop-gap resonator [63, 64, 65]. We use a
capacitive antenna coupling to the resonator which houses the diamond sample
as depicted in Figure 1(b). Details of the experimental setup are provided in
Appendix C. No optical excitation is used in the P1 center experiments as they
cannot be optically hyperpolarized directly. In a single shot of a pulse
sequence of fixed $\tau$, we stroboscopically acquire each echo formed. As
shown in the sequence schematic in Figure 2(a), we acquire the echo formed at
$\tau/2$ after every $\pi_{y}$ pulse. The P1 center data was averaged over
65,536 shot repetitions, due to the low thermal equilibrium spin polarization,
while only 128 shot repetitions were needed for the optically hyperpolarized
NV centers. An example acquisition for the P1 data is shown in Figure 10 in
Appendix C.

The data in Figure 2 exemplify how the rate of coherence decay can be
significantly different for two sequences with only slightly different $\tau$
values. Figure 2(b) shows the stroboscopically acquired decays of echo
amplitudes for two CPMG sequences with different $\tau$ values at both 13 mT
and 190 mT. The time at which an echo is detected is $t_{d}=N\tau$. Since we
collect these decays for a wide range of pulse spacings, we can rearrange the
data to present a coherence curve for CPMG-$N$, where $N$ refers to
acquisition of the $N^{\text{th}}$ echo, which occurs after the
$N^{\text{th}}$ $\pi$ pulse. Figure 2(c) shows the results obtained at the two
fields for a CPMG-1 (Hahn echo) and a CPMG-16 sequence. The data at 190 mT
show only subtle variations with small changes in $\tau$, while the 13 mT data
exhibit a strong periodic dependence. The decoherence rate has a periodic
dependence on $\tau$ and the Larmor precession frequency of the 13C nuclear
spins in the form of electron spin echo envelope modulations (ESEEM), which
are explained in a noise spectroscopy framework with Figure 3. Without the
context of this periodic behavior, the 13 mT fixed-$\tau$ representation of
the data shown in Figure 2(b) may seem counterintuitive since the sequence
with the shorter $\tau$ of 9.83 $\mu$s decays more quickly than the sequence
with $\tau=13.8\ \mu$s.

##  III Dynamical decoupling-based spectroscopy in diamond

###  III.1 System and decoherence model in diamond

Since the concentrations of the measured spins – NV centers in Sample A and P1
centers in Sample B – are low, the dipolar interactions between identical
measured spins can be ignored. The underlying physics for both systems can
thus be described in terms of the central spin model [66]. Here the
environment consists of a spin bath as we can ignore interactions with lattice
phonons on timescales shorter than the spin-lattice relaxation time $T_{1}$.

For an electronic central spin in diamond, the total Hamiltonian can be
expressed as

$H_{\text{Tot}}(t)=H_{S}+H_{C}(t)+H_{E}+H_{SE}$ (1)

where $H_{S}$ is the aforementioned spin-1/2 system Hamiltonian. The control
Hamiltonian is $H_{C}(t)=-\gamma_{e}\vec{S}\cdot\vec{B}_{1}(t)$, where
$\vec{S}$ is the central spin operator vector $(S_{x},S_{y},S_{z})$ and
$\vec{B}_{1}(t)$ is the magnetic microwave control field. The environment
Hamiltonian $H_{E}$ consists of the Zeeman and spin-spin interactions for bath
spins in the vicinity of the central spin. The system-environment Hamiltonian
$H_{SE}$ includes the interactions between the central spin and the bath
spins. The bath consists of 13C nuclear spins, and unlike electron spins. When
the NV is the central spin, the unlike electron spins are the P1 spins. When
the P1 is the central spin the unlike electron spins are P1 spins of
spectrally-separated 14N hyperfine manifolds. Since the experiments on the P1
centers in Sample B are performed on the $m_{I}=0$ hyperfine manifold, the
unlike electrons that constitute the bath are those in the $m_{I}=\pm 1$
states.

Figure 3:  (a) Time domain illustration of dynamical decoupling for two
sequences (red and blue) with different interpulse delay times
$\tau_{\text{r}}$ and $\tau_{\text{b}}$. The microwave pulses, spin echoes,
and coherences $L$ are shown for both sequences in the first and third rows.
The associated modulation functions $f_{z}$ are shown in the second and fourth
rows. The fifth row shows an example of what the effective field noise
$b_{z}(t)$ may look like for a particular spin. The 13C oscillation is
exaggerated to highlight the resonance condition with the modulation function.
Below that are the rectified noise traces $b_{z}(t)f_{z}(t,\tau)$ for the red
and blue sequences applied to the given noise example. Note that the time
averaged value of the rectified noise is less than zero in the red example and
approximately zero in the blue example. The bottom line shows the accumulated
phase,
$\phi_{i}(t)=-\gamma_{e}\int_{0}^{t}b^{i}_{z}(t^{\prime})f_{z}(t^{\prime},\tau)dt^{\prime}$
for 100 spins with unique realizations of the phase and amplitude of the noise
$b_{z}^{i}(t)$ under application of either pulse sequence. The coherence for
the red sequence decays faster than the coherence for the blue sequence
because the red sequence satisfies the resonant condition
$\tau=m\pi/\omega_{\text{13C}}$, where $m$ is odd ($m=3$ in this case), and
$\omega_{\text{13C}}$ is the dominant frequency in $b_{z}(t)$, leading to a
large non-zero average value of the rectified noise. For the blue sequence,
$m=4$ and the rectified noise is time-averaged to nearly zero. (b) The overlap
integral, Equation 13, shown as a cumulative integral, i.e.,
$\int_{0}^{\omega}\chi^{\prime}(\omega^{\prime})d\omega^{\prime}$, where
$\chi^{\prime}=t_{d}^{2}S(\omega)F(\omega,t_{d},\tau)$. The cumulative
integral shows that the contribution of the fundamental peak is larger for the
blue sequence, but the contribution of the $m=3$ peak, which overlaps the
prominent Gaussian feature in the power spectrum for the red sequence, gives
the largest contribution to $\chi$, explaining the more rapid decay for the
red sequence. (c) The power spectrum $S(\omega)$ (black line, corresponding to
the right side vertical axis) and filter functions $F(\omega)$ (corresponding
to the left side vertical axis) for the red and blue sequences with a
detection time of $t_{d}=24\tau_{\textrm{r}}=18\tau_{\textrm{b}}$. The
prominent peaks in the filter function occur at $\omega=m\pi/\tau$, with
$m=1,3,5...$.

For purely dephasing interactions between the central spin and the environment
spins, the system-environment Hamiltonian is

$H_{SE}=S_{z}\left[\sum_{i}\left(A_{\perp}^{i}I_{x}^{i}+A_{\parallel}^{i}I_{z}^{i}\right)+\sum_{i}J_{\parallel}^{i}S_{z}^{i}\right].$
(2)

The first summation in the brackets is over the nuclear spins, where
$A_{\perp}^{i}$ and $A_{\parallel}^{i}$ are the transverse and longitudinal
hyperfine coupling components between the central spin and the $i^{\text{th}}$
13C nuclear spin (with spin operators $I_{x}^{i}$ and $I_{z}^{i}$). The second
summation is over the unlike electron spins in the central spin’s environment.
$J_{\parallel}^{i}$ is the coupling strength between the central spin and the
$S_{z}^{i}$ component for the $i^{\text{th}}$ unlike electron spin. The non-
secular transverse components of the electron-electron coupling are omitted
for reasons explained below.

In the interaction frame defined by the system and environment Hamiltonians
($H_{S}+H_{E}$), the system-environment Hamiltonian $H_{SE}$ becomes time-
dependent due to non-commuting terms. Notably, the transverse hyperfine
component transforms as $A_{\perp}S_{z}I_{x}\rightarrow
A_{\perp}S_{z}e^{-i\omega_{\text{13C}}I_{z}t}I_{x}e^{i\omega_{\text{13C}}I_{z}t}=A_{\perp}S_{z}(I_{x}\cos(\omega_{\text{13C}}t)+I_{y}\sin(\omega_{\text{13C}}t))$.
That is, the precession of the 13C nuclei gives rise to a modulation of the
hyperfine interaction with the central spin at the 13C Larmor frequency
$\omega_{\text{13C}}$. Similar terms with oscillating transverse spin
components also arise for the coupling to the unlike electron spins. However,
these terms can be ignored – and the transverse components of the electron-
electron coupling are omitted from Equation 2 – because the unlike electrons’
Larmor frequencies are very high compared to the coupling strengths, so their
effects are averaged out on short timescales [67]. Other terms in $H_{E}$ such
as dipolar couplings among environment spins do not commute with $H_{SE}$ and
give rise to fluctuations in the full system that are effectively random when
we trace over the environment. Performing such a trace over the bath spins
results in an effective spin Hamiltonian for the central spin [66]:

$H^{(E)}_{SE}(t)=-\gamma_{e}b_{z}(t)S_{z}\ ,$ (3)

where $b_{z}(t)$ is the $z$-component of the time-varying magnetic field noise
at the site of the central spin. It is the sum of the magnetic fields produced
by the fluctuating environment spins and contains periodic components at the
13C Larmor frequency. Figure 3(a) shows a cartoon example time trace of the
noise $b_{z}(t)$ that is dominated by an exaggerated oscillatory component.

###  III.2 The CPMG experiment and noise spectroscopy

The application of a $\pi$-pulse flips the sign of $H^{(E)}_{SE}$ of Equation
3 in a toggling frame representation. For DD sequences, which consist of a
train of $\pi$ pulses, the repeated flipping of the sign can be expressed by a
square-wave frame modulation function $f_{z}(t,\tau)$. The effective toggling
frame Hamiltonian for the dephasing noise acting on the central spin under
application of DD is expressed as

$\tilde{H}_{SE}(t,\tau)=-\gamma_{e}f_{z}(t,\tau)b_{z}(t)S_{z}.$ (4)

Figure 3(a) shows the frame modulation functions $f_{z}(t,\tau)$ for two
values of $\tau$ and illustrates how the rectified noise
$f_{z}(t,\tau)b_{z}(t)$ depends on $\tau$.

Prior to the first microwave pulse, the state of the system is given by the
density operator, which can be approximated at room temperature as

$\rho_{\text{pre}}=\frac{\mathbbm{1}}{2}-\epsilon S_{z},$ (5)

where $\epsilon$ is a polarization factor that is proportional to the signal
strength. If the system is in thermal equilibrium prior to the application of
microwave pulses, then $\epsilon=\epsilon_{\text{th}}$, which depends on the
temperature and resonance frequency as

$\epsilon_{\text{th}}=\tanh\left(\frac{-\omega_{0}\hbar}{2k_{B}T}\right),$ (6)

where $k_{B}T/\hbar$ is the Boltzmann energy converted to frequency units. For
$\omega_{0}=-2.5\times 2\pi$ rad GHz, and $T=300\text{ K}$ the polarization is
$\epsilon_{\text{th}}\approx 2.0\times 10^{-4}.$ For the NV center
experiments, hyperpolarization with green light achieves
$\epsilon=\epsilon_{\text{hyp}}\approx 120\epsilon_{\text{th}}$.

Upon application of the initial $(\pi/2)_{x}$ microwave pulse, the density
operator becomes

$\displaystyle\rho_{0}$ $\displaystyle=e^{-i\pi
S_{x}/2}\rho_{\text{pre}}e^{i\pi S_{x}/2}$ (7) $\displaystyle\rho_{0}$
$\displaystyle=\frac{\mathbbm{1}}{2}-\epsilon S_{y}.$ (8)

In the CPMG sequence, a series of $\pi_{y}$ pulses is applied such that the
$N^{\text{th}}$ $\pi_{y}$ pulse occurs at time $t_{N}=\tau(N-1/2)$ and the
$N^{\text{th}}$ echo is detected at time $t_{d}=N\tau$. The amplitude of the
$N^{\text{th}}$ echo is given by the expectation value for $S_{y}$:

$\displaystyle\langle S_{y}\rangle(t,\tau)$ $\displaystyle=\langle
e^{-i\phi(t,\tau)S_{z}}\rho_{0}e^{i\phi(t,\tau)S_{z}}S_{y}\rangle$ (9)
$\displaystyle=-\frac{\epsilon}{2}\cos(\phi(t,\tau))$ (10)

where
$\phi(t,\tau)S_{z}=\int_{0}^{t}dt^{\prime}\tilde{H}_{SE}(t^{\prime},\tau)$.

The coherence $L(t,\tau)$ that we analyze is the ensemble expectation value of
$S_{y}$ normalized by its value at time $t=0$, at which point
$\phi(0,\tau)=0$.

$L(t,\tau)=\frac{\langle S_{y}\rangle_{\text{ens}}(t,\tau)}{\langle
S_{y}\rangle_{\text{ens}}(0,\tau)}=\langle\cos(\phi(t,\tau))\rangle_{\text{ens}}$
(11)

The bottom row of Figure 3(a) shows the accumulation of many $\phi(t,\tau)$
for an example ensemble under application of two CPMG sequences with different
$\tau$ values.

For a disordered system composed of many weak couplings to the central spin,
it is useful to make the approximation that $b_{z}(t)$ is drawn from a
stationary Gaussian random distribution with zero mean. Then $\phi(t,\tau)$
for the ensemble is a Gaussian random variable, i.e., the ensemble
distribution of phases at any given time $t$ is
$P_{\phi}(t,\tau)=(e^{-\phi^{2}/(2\sigma(t,\tau)^{2})})/(\sqrt{2\pi}\sigma(t,\tau)^{2})$
[66]. In this case, the coherence can be shown to be
$L(t,\tau)=e^{-\langle\phi(t,\tau)^{2}\rangle/2}$, where
$\langle\phi(t,\tau)^{2}\rangle$ is the variance $\sigma(t,\tau)^{2}$ of the
ensemble phase distribution. We set
$\langle\phi(t,\tau)^{2}\rangle/2=\chi(t,\tau)$ such that the coherence can
simply be expressed as $L(t,\tau)=e^{-\chi(t,\tau)}$. We can compute $\chi$ at
the peak of the $N^{\text{th}}$ echo, when $t=N\tau\equiv t_{d}$, as

$\begin{split}\chi(t_{d},\tau)=\frac{\gamma_{e}^{2}}{2}\bigg{\langle}\int_{0}^{t_{d}}dt_{1}\int_{0}^{t_{d}}dt_{2}\qquad\qquad\\\
\qquad\qquad
b_{z}(t_{1})f_{z}(t_{1},\tau)b_{z}(t_{2})f_{z}(t_{2},\tau)\bigg{\rangle}_{\text{ens}},\end{split}$
(12)

which can be expressed in terms of a dimensionless filter function $F$ and the
power spectral density $S(\omega)$, which is the Fourier transform of the
autocorrelation function of $b_{z}(t)$ [66]

$\chi(t_{d},\tau)=\frac{t_{d}^{2}}{2}\int_{-\infty}^{\infty}d\omega\,S(\omega)F(\omega,t_{d},\tau).$
(13)

The dimensionless filter function $F$ is

$F(\omega,t_{d},\tau)=\frac{1}{t_{d}^{2}}\left|\frac{1}{\sqrt{2\pi}}\int_{0}^{t_{d}}f_{z}(t,\tau)e^{-i\omega
t}dt\right|^{2},$ (14)

and the power spectral density $S(\omega)$ is

$S(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}g(t)e^{-i\omega t}\,dt$
(15)

where $g(t)=\gamma_{e}^{2}\langle
b_{z}(t^{\prime})b_{z}(t^{\prime}+t)\rangle$.

Figure 3(c) shows the filter functions corresponding to the two example frame
modulation functions and an example noise power spectrum $S(\omega)$. Figure
3(b) shows the integration of Equation 13 in the form of a cumulative integral
for the two $\tau$ examples.

The challenge of noise spectroscopy is to calculate $S(\omega)$, given the
measured $\chi$ values and the form of the filter $F$ that describes the
control sequence. For DD sequences consisting of $N$ instantaneous $\pi$
pulses equally spaced by the interval $\tau$, the filter function is [66]

$F(\omega,N,\tau)=\frac{8\sin^{4}(\omega\tau/4)\sin^{2}(N(\pi-\omega\tau)/2)}{\pi\omega^{2}N^{2}\tau^{2}\cos^{2}(\omega\tau/2)}.$
(16)

In the large $N$ limit, the filter function can be approximated as the sum of
$\delta$-functions

$F(\omega,N,\tau)\approx\frac{4}{N\tau^{3}}\sum_{m=\pm\textrm{odd}}^{\infty}\frac{\delta(\omega-m\pi/\tau)}{\omega^{2}}.$
(17)

which results in a discrete sum expression for $\chi$:

$\chi(N,\tau)=\frac{4N\tau}{\pi^{2}}\sum_{m=1,3,5,...}^{\infty}\frac{S(\omega_{m})}{m^{2}},$
(18)

where $\omega_{m}=m\pi/\tau$. If $S(\omega_{1})\gg S(\omega_{m})$ for all
$m>1$, it is appropriate to approximate $F$ as a single peak at the
fundamental ($m=1$) frequency $\omega=\pi/\tau$ and ignore the contributions
of higher harmonics. We refer to this approximation as the “fundamental”
approximation of $S$:

$S(\pi/\tau)=\frac{\pi^{2}\ \chi(N,\tau)}{4N\tau}.$ (19)

For a smooth spectrum that monotonically decays with $\omega$, a sequence with
a shorter $\tau$ value decouples the system from a greater range of low-
frequency noise and therefore extends coherence more than a sequence with a
longer $\tau$ value does. However, if a peak in the filter function $F$
overlaps with a sharp peak in the power spectrum $S$, this will cause a sharp
dip in the coherence. Equivalently, if $\tau=mT/2\,\,(m=1,3,5...)$, where $T$
is the oscillatory period of a prominent noise source, such as the 13C Larmor
precession, then the coherence decays much more quickly than it otherwise
would. The 190 mT decay curves shown in Figure 2 are well-described by the
smooth spectrum scenario (since the periodic modulations are very subtle). The
13 mT decay curves exemplify the prominent peak scenario, as the periodic
modulations dominate the decay. The $\tau=9.83\ \mu$s and $\tau=13.8\ \mu$s
curves for the 13 mT dataset in Figure 2(b) correspond approximately to the
red and blue curves in Figure 3. In Figure 3, the red sequence decays more
quickly than the blue sequence, despite having $\tau_{r}<\tau_{b}$. This is
because we chose, for the sake of illustration, $\tau_{r}=3T/2$ and
$\tau_{b}=2T$, where $T$ is the period of the dominant oscillation in the
noise in Figure 3(a), corresponding to the central frequency of the prominent
peak in the power spectrum in Figure 3(c).

Figure 4:  (a), (b) and (c) are image plots of $\chi$ data from CPMG
experiments. (a) and (b) show data obtained with Sample A on NV centers at 13
mT and 190 mT respectively. (c) shows data obtained with Sample B on the P1
centers at 89 mT. $\chi=0$ indicates the maximum echo amplitude for the
dataset. These data are $T_{1}$-corrected as described by Equation 20. The
broad character of the power spectrum is captured by the transparent black and
white dashed lines, which are fits of contours of constant $\chi$ to a power
law form given by Equation 25. The corresponding $\chi$ values are labeled at
the lower right ends of the lines. The red vertical lines at the bottoms of
the plots indicate $\tau$ values where coherence revivals are expected:
$\tau=m\pi/\omega_{\text{13C}}$, where $m$ is even and $\omega_{\text{13C}}$
is the 13C Larmor frequency. For the 13 mT in (a), the fits are applied only
to the local $t_{d}$ maxima of the $\chi$ contours, which occur at these
revival $\tau$ values. (d) and (e) show the parameter values from the $\chi$
contour fits for the power law model of the noise spectrum
$S(\omega)=C\omega^{-\alpha}$. The error bars indicate the 95% confidence
interval of the fit. In the 13 mT dataset, the scarcity of local maxima lead
to large error bars for $\chi<1$ and an absence of fits for $\chi\leq 0.6$.
See Section IV.2 for further discussion.

##  IV Results and Discussion

###  IV.1 Overview of data

The observed transverse decay is a product of the decay due to spin bath
interactions and, to a lesser extent, $T_{1}$ relaxation,
$L_{\text{obs}}(t_{d})=\exp[-(t_{d}/(2T_{1})-\chi(t_{d})]$. In the absence of
a spin bath, we expect that the $T_{2}$ relaxation would be dominated by spin-
lattice interactions and $T_{2}=2T_{1}$. Accordingly, to correct for the
effects of $T_{1}$ relaxation, such that $\chi$ only reflects decoherence due
to the bath spins, we calculate $\chi$ as

$\chi(t_{d})=-\log(L_{\text{obs}}(t_{d}))-\frac{t_{d}}{2T_{1}}.$ (20)

The measured $T_{1}$ times are on the order of milliseconds and have little
effect on the NV $\chi$ calculation, for which the coherence times are on the
order of 100 $\mu$s. The $T_{1}$ correction is slightly more significant for
the P1 center data in Sample B, which achieves coherence times around
$T_{1}/2$ in the short $\tau$ limit. Details of $T_{1}$ and $T_{2}$
experiments and the resulting values are provided in Appendix A.

Figures 4 (a), (b), and (c) show image plots of $\chi$ calculated from
Equation 20. The vertical axes of the image plots indicate the time at which
an echo occurs. Therefore, a single-shot experiment with fixed $\tau$ is
represented by a vertical line ascending from bottom to top of the image plot.
As the echo amplitudes decay, $L$ decreases from 1 to the noise floor, as
limited by the signal-to-noise of the experiment, and $\chi$ increases from 0
to some $\chi_{\text{max}}$, which is about 5.5 for the NV center datasets and
about 3.5 for the P1 center dataset. In the lower right areas of Figures 4
(a), (b), and (c), one can see sections of data that are uniform in the
vertical direction, and the boundaries of the sections are slanted from lower
left to upper right. This is because each section corresponds to a particular
echo number $N$. The lower rightmost section appears almost as a triangle,
which is all $N=1$ (Hahn echo), with the $N=2$ section directly above it, and
so on. The image plot is constructed by defining the $\tau$ and $t_{d}$ axes
first and then populating the 2D array with $\chi$ values from the data set
that most nearly match for each point $(\tau,t_{d})$.

Much of our analysis focuses on the structure of fixed $\chi$ contours
throughout the data sets. A commonly used approach, as in References [35, 40,
68], is to focus on the contour $\chi=1$, which is encountered when
$t_{d}=T_{2}$, resulting in Equation 19 taking the form
$S(\pi/\tau)=\pi^{2}/(4T_{2})$. We extend the contour tracing approach to the
analysis of contours of densely sampled $\chi$ over the available range of
values for two reasons. First, we find that this provides a simple way to
quantitatively capture the broad character of the coherence data. Second,
since larger $\chi$ values generally occur at later detection times, the
$\delta$-function approximation of the filter function has greater validity,
enabling use of Equations 18 and 19. This notion is consistent with the
observation that longer CPMG trains lead to greater depth enhancement of ESEEM
[43]. Conversely, in the low $N$ limit, the peak in $F$ is wider and is
shifted from $\pi/\tau$, as can be seen in Figure 12(d).

The $\chi$ contour traces are obtained from the data as follows. For each
$\tau$, we convolve the coherence decay $\chi(t_{d})$ with a Gaussian
function: $\chi_{\text{conv}}(t_{d})=\int_{-\infty}^{\infty}dt_{d}^{\prime}\
\chi(t_{d})A\exp[-(t_{d}-t_{d}^{\prime})^{2}/(2\sigma^{2})]$, with $\sigma=4\
\mu$s and $A$ being a normalization constant to preserve the magnitude of
$\chi$. This smooths out fluctuations for short $\tau$ decays and leaves
longer $\tau$ decays unaffected. The $\chi$ contour consists of a list of
pairs of $N$ and $\tau$ values, where $N$ for a given $\tau$ is the first echo
for which $\chi_{\text{conv}}$ is greater than the $\chi$ for the contour
being traced.

The $\chi$ data exhibit a broad background and sharper modulations due to the
13C transverse hyperfine interactions. Accordingly, our analysis involves
characterization of the background spectrum $S_{B}$ and the sharp peak $S_{P}$
due to the 13C transverse hyperfine interaction. We model the total spectrum
as the sum of these contributions:

$S(\omega)=S_{B}(\omega)+S_{P}(\omega).$ (21)

As can be seen in Figure 4, in the 13 mT NV data, the modulations due to the
13C interaction dominate the data and prevent a clear measurement of $S_{B}$.
However, this also means that at 13 mT, $S_{P}$ can be measured with the
fundamental approximation of the filter (Equation 19). To obtain the
background at 13 mT, we analyze the maximal $t_{d}$ values of the $\chi$
contours, which, for larger $\chi$, occur near the revival $\tau$ values
($\tau=m\pi/\omega_{\text{13C}}$, where $m=2,4,6...$). At 89 mT and 190 mT,
the 13C modulations are more subtle, allowing a clearer characterization of
$S_{B}$, which is then used in a background subtraction step to characterize
$S_{P}$ with higher harmonics of the CPMG filter function.

Figure 5:  (a) and (b) both show the fundamental (“Fund.”) spectra, i.e., that
obtained by using Equation 19 with the $\chi$ contours 1 and 3.5 as indicated.
The solid straight lines show the power law fits (“Pow. Law”) obtained in the
manner described in Section IV.2. The reconstructed spectra (“Recon.”) are
discussed in Section IV.2 and the methods of reconstruction are explained in
Appendix F.2. Due to the contour trace consisting of experimentally realized
values of $N\tau$, which are discrete, the low-frequency (large $\tau$) limit
of the fundamental spectrum exhibits a jagged appearance, while the high-
frequency (short $\tau$) limit appears as more of a continuum. Each of the
datasets exhibits subharmonics that arise due to overlap of the CPMG filter
harmonics with the 13C peak. Only at 13 mT is the 13C peak at low enough
frequency to be present in the fundamental spectrum. The cyan dashed line in
(b) shows the Gaussian fit to that peak.

###  IV.2 Broad spectrum noise characterization

The electronic spin bath surrounding an NV center has been described well by
an Ornstein-Uhlenbeck (O-U) process [47, 50, 61, 46, 56]. The form of the
power spectral density of such a process is a Lorentzian centered at
$\omega=0$,

$S(\omega)=\frac{2\Delta^{2}\tau_{c}}{1+\omega^{2}\tau_{c}^{2}},$ (22)

which has a high-frequency power law tail of the form $1/\omega^{2}$. We
therefore model the broad spectrum as a power law
$S_{B}(\omega)=C\omega^{-\alpha}$, which enables the inclusion of harmonics
that provide a small correction to the calculation of $S_{B}$ [34, 40]. The
discrete approximation of Equation 18 becomes

$\chi(N,\tau)=C^{\prime}(\alpha)\,N\,\tau^{1+\alpha},$ (23)

where

$C^{\prime}(\alpha)=\frac{4\,\zeta(2+\alpha)}{\pi^{2+\alpha}}\left(1-\frac{1}{2^{2+\alpha}}\right)C,$
(24)

with $\zeta$ being the Riemann zeta function. We can find the exponent
$\alpha$ by first rearranging Equation 23 to the form

$t_{d}=\frac{\chi(\tau,N)}{C^{\prime}(\alpha)\tau^{\alpha}}.$ (25)

According to Equation 25, a $\chi$ contour should follow a straight line on
the $\log(t_{d})$ vs. $\log(\tau)$ axes if the noise spectrum obeys a power
law of the form $S_{B}(\omega)=C\omega^{-\alpha}$. The straight black and
white dashed lines in Figures 4 (a), (b), and (c) are obtained from fitting
the contours to Equation 25. The fit parameters $C$ and $\alpha$ from the
$\chi$ contours are plotted in Figures 4 (d) and (e), and the average values
over the plotted ranges are given in Table 1.

Sample |  | A | A | B
---|---|---|---|---
Species |  | NV | NV | P1
$B_{0}$ (mT) |  | 190 | 13 | 89
$\log_{10}(C)$ |  | 10(1) | 8(2) | 8(1)
$\alpha$ |  | 0.9(2) | 0.7(2) | 0.7(2)
Table 1: Average values of power law ($C\omega^{-\alpha}$) spectrum fit
parameters shown in Figure 4(d) and (e).

The origin of the $1/\omega$ scaling for the noise is not definitely
understood, but we suspect that the high concentration of nitrogen and the
heterogeneity of the nitrogen concentration throughout the sample play
significant roles. Lorentzian noise has been observed in single NV spins in
samples with low nitrogen concentration [47, 50, 61, 46]. Within an HPHT
diamond sample with reported nitrogen concentration of [N] $\approx$ 200 ppm
that has been irradiated and annealed for the creation of NV centers, the
localized concentrations of P1 centers can vary over ranges as wide as 10 to
300 ppm [26]. One possibility is that the high P1 concentration, which will
have much stronger P1-P1 dipolar couplings and exchange couplings [28], leads
to a fundamentally different type of dynamics and the O-U model no longer
applies. Another possibility, motivated by References [69, 70], is that the
high concentration P1 baths still undergo O-U dynamics, but their correlation
times $\tau_{c}$ vary with concentration, giving the appearance of $1/f$ noise
in our ensemble measurement. This would happen if the varied P1 concentrations
throughout the sample give rise to an effective distribution of correlation
times that can be approximated as $P(\tau_{c})=P_{0}/\tau_{c}$ between some
range of cutoff correlation times $\tau_{1}<\tau_{c}<\tau_{2}$. Then the
ensemble average power spectral density would be

$\displaystyle\langle S(\omega)\rangle_{\text{ens}}$
$\displaystyle=2\Delta^{2}P_{0}\int_{\tau_{1}}^{\tau_{2}}\frac{d\tau_{c}}{1+\omega^{2}\tau_{c}^{2}}$
(26)
$\displaystyle=\frac{2\Delta^{2}P_{0}}{\omega}(\arctan(\tau_{2}\omega)-\arctan(\tau_{1}\omega)),$
(27)

which takes the form of $1/\omega$ for $1/\tau_{2}\ll\omega\ll 1/\tau_{1}$.
Typically, $T_{2}\propto[\text{N}]^{-1}$, as has been shown with simulations
[50] and experiments [24, 25]. Considering also that $T_{2}\propto\tau_{c}$,
this suggests that the $1/\omega$ noise could arise for a heterogeneous
population of P1 concentrations ([N]) with an approximate distribution
throughout the sample of $P([\text{N}])\propto[\text{N}]$ over some range of
[N] corresponding to $\tau_{1}$ and $\tau_{2}$. Further research would be
needed to confirm or rule out either of these explanations.

Similar heterogeneity effects might explain the roughly $1/\omega$ spectrum
measured with the P1 centers in Sample B, however P1 clustering has not been
investigated to the same degree in CVD diamonds as it has in HPHT diamonds. In
Sample B, CPMG was effective enough to increase $T_{2}$ of the P1 centers
(seen in Figure 4(c) as detection times $t_{d}$ where $\chi=1$) from the
nominal Hahn echo value of 0.18(4) ms to values over 1 ms. For diamonds with
$0.1\text{ ppm}<\text{[N]}<1\text{ ppm}$, the Hahn echo $T_{2}$ of P1 centers
is limited by a combination of P1-P1 spin interactions and by the $T_{2}$ of
the 13C spins, which has been reported to be 190(10) $\mu$s [71]. It is
therefore possible that the $1/\omega$ spectrum results from a combination of
P1 clustering and 13C interactions.

Figures 5(a) and (b) display the “fundamental” calculations of the power
spectral density. In the fundamental spectrum, each $\tau$ value corresponds
to a point at the frequency $\omega=\pi/\tau$. We calculate the fundamental
spectrum using Equation 19 with $\chi$ fixed as the $\chi$ value of the
contour and with $N$ as that found for each $\tau$ value along the contour.
The solid straight lines in the fundamental spectra of Figures 5(a) and (b)
are the power law spectra generated from fitting the $\chi$ contours. Since
the power law fits for the 13 mT dataset were applied to the local maxima of
$t_{d}$ in the $\tau$ domain, the power laws that appear in the frequency
domain follow the minima between the subharmonics of the 13 mT fundamental
spectrum.

The power laws are seen in Figures 5(a) and (b) to match the respective
fundamental spectra in the mid- to high-frequency region. At lower
frequencies, particularly for the 190 mT dataset, the power law is
significantly higher than the fundamental spectrum data. Although the power
spectrum must ultimately approach a finite value in the limit
$\omega\rightarrow 0$, we find that we can reliably reconstruct the
fundamental spectrum from the power law model and a numerically computed CPMG
filter function. That is to say that the power law model is indistinguishable
from an alternative spectrum model (such as a stretched Lorentzian) that has a
finite value of $S$ at $\omega\rightarrow 0$. The details of the
reconstruction are given in Appendix F.2.

Nevertheless, we can estimate the finite magnitude of the power spectrum in
the limit $\omega\rightarrow 0$ by measurements of $T_{2}$ obtained from Hahn
echo experiments [35]. In the limit $\tau\rightarrow\infty$, the peaks of the
filter function approach $\omega\rightarrow 0$. In this limit, Equation 18
leads to $S(0)=2\chi/(N\tau)$, which differs only slightly from the
fundamental approximation by a factor of $\pi^{2}/8\approx 1.2$. The Hahn echo
$T_{2}$ values are provided in Appendix A. The resulting $S(0)$ are 41(11) kHz
for the NV at 13 mT, 82(1) kHz for the NV at 190 mT, and 11(2) kHz for the P1
at 89 mT. These values are in good agreement with the low-frequency levels of
the respective spectra in Figures 5(a) and (b).

The overall lower level of the spectrum for the P1 centers in Sample B is
primarily due to that sample having a much lower concentration of nitrogen
[24, 25]. Comparing the 13 mT and 190 mT data, we see that the 13 mT spectrum
dips lower than the 190 mT spectrum, indicating elevated low-frequency noise
at higher field. The cause of this may be increased axial alignment of the P1
bath spins at the higher field. The 13C may also play a significant role in
this, since at high field, the Zeeman interaction dominates the 13C
Hamiltonian, but at low field the NV hyperfine interaction is comparable to or
dominates the 13C Hamiltonian.

We also have investigated experimentally and numerically how finite pulses and
over/under rotation errors affect the power law characterization. These
investigations are presented in Appendix F. We find that the finite pulses in
our experiments do not significantly affect the results of the power law
characterization. This is because the finite pulses most significantly affect
the higher harmonics of the filter function, and in the context of a power law
spectrum, these harmonics contribute only a very small amount to the $\chi$
integral in Equation 13.

Figure 6:  (a) and (b) show the high-frequency spectra for the 89 mT P1 and
190 mT NV experiments respectively. These spectra are obtained via the
harmonic analysis method presented in Appendix D. The insets show detail of
the 13C peak fit to a Gaussian curve (dotted lines) and the expected 13C
Larmor frequency (vertical black lines). (c), (d), and (e) show the fit
parameters from applying Gaussian fits
($A\exp[-(\omega-\omega_{\text{13C}})^{2}]/(2\sigma^{2})]$) to the 13C peak.
The fits for the 89 mT and 190 mT data are those shown in the insets of (a)
and (b). For the 13 mT data, the fit was applied to the “fundamental”
$\chi=3.5$ spectrum indicated by the cyan dashed line in Figure 5 (b).

###  IV.3 Characterizing the 13C Larmor frequency

The transverse component of the 13C hyperfine interaction contributes a peak
$S_{P}$ to the power spectrum at the 13C Larmor precession frequency. We can
use the fundamental approximation of Equation 19 to observe this peak at 13
mT. However, at 89 mT and 190 mT, the 13C Larmor precession period is
significantly shorter than the shortest $\tau$ we can use in our CPMG
experiments, and we only see the periodic modulation of the coherence curves
with increasing $\tau$ as higher harmonics of the filter function overlap with
the 13C peak. We have developed an analysis procedure that utilizes all
applicable harmonics of the filter function to estimate the properties of
$S_{P}$ at high frequencies. The basic approach is outlined below and details
are provided in Appendix D. Note that $S_{P}$ could also be generalized to a
sum of sharp peaks $\sum_{i}S_{P}^{i}$ for contexts wherein the goal is to
obtain the spectrum of a more complex system.

Inserting the composite spectrum expression of Equation 21 into the discrete
$\chi$ calculation of Equation 18, we obtain

$\chi(N,\tau)=\\\
\frac{4N\tau}{\pi^{2}}\sum_{m=1,3,5...}^{\infty}\frac{S_{B}(m\pi/\tau)+S_{P}(m\pi/\tau)}{m^{2}}.$
(28)

For a sharp peak, we can assume that $S_{P}(\omega)=0$ for all $m$ except near
the sharp feature which is peaked at a frequency $\omega_{P}$. In the vicinity
of the peak we define the scanning frequency $\omega_{s}$ and scanning
harmonic $m_{s}$ with $\omega_{s}\equiv m_{s}\pi/\tau$. Then the summation of
$S_{P}$ collapses to only the $m=m_{s}$ term:

$\chi(N,\tau)=\frac{4N\tau}{\pi^{2}}\ \ \times\\\
\left[\left(\sum_{m=1,3,5...}^{\infty}\frac{S_{B}(m\pi/\tau)}{m^{2}}\right)+\frac{S_{P}(\omega_{s})}{m_{s}^{2}}\right].$
(29)

If the broad background component $S_{B}$ is known, for example by using
Equation 25, we can solve for $S_{P}$.

$S_{P}(\omega_{s})=m_{s}^{2}\ \ \times\\\
\left[\frac{\pi^{2}\chi(N,\tau)}{4N\tau}-\left(\sum_{m=1,3,5...}^{\infty}\frac{S_{B}(m\pi/\tau)}{m^{2}}\right)\right].$
(30)

We can express $S_{P}$ as a difference between the measured value $\chi$ at a
particular $N$ and $\tau$ combination and the expected background $\chi_{B}$
in the vicinity of that point.

$S_{P}(\omega_{s})=\frac{\pi^{2}m_{s}^{2}}{4N\tau}\left[\chi(N,\tau)-\chi_{B}(N,\tau)\right]$
(31)

When we perform this calculation in our analysis, we include an additional
term in the pre-factor to account for finite pulse widths as explained in
Appendix F.1. For a given combination of $N$ and $\tau$, we first compute
$\chi_{B}$ by determining which power law model (corresponding to the
different $\chi$ contour fits) most nearly approaches the given $N$ and $\tau$
point. This is the same as determining by an image plot in Figure 4, which
black and white dashed line most nearly crosses a particular ($\tau,t_{d}$)
point. That power law model is then used to compute the expected value of
$\chi_{B}$ using Equation 23.

A spectrum $S_{P}$ is thus obtained by performing the calculation of Equation
31 for all applicable combinations of $N$ and $\tau$ in the dataset in a
harmonic analysis procedure presented in Appendix D. In the harmonic analysis
procedure, we define a scanning window in the frequency domain. The maximum
spectral peak width that can be measured is the width of this scanning window.
This width is also the inverse of the largest $\tau$ that can be used in the
analysis. For one location of the window, we use Equation 31 to calculate a
segment of the spectrum. The window is stepped through a wide range of
frequencies, and the many spectral segments are combined. The resulting
spectra for the 89 mT and 190 mT data are shown in Figure 6(a) and (b). We
note that while negative values are not permitted in a power spectrum, they
appear here because obtaining the spectra $S_{P}$ via Equation 31 is
effectively a form of baseline subtraction.

Figures 6(c), (d) and (e) show the fit parameters obtained from fitting a
Gaussian to the reconstructed 13C spectral line. The error bars indicate the
95% confidence intervals. For the 13 mT experiment, we use the $\chi=3.5$
contour to obtain the fit. The data to which the fit is applied is the portion
of the fundamental spectrum between 90 kHz and 200 kHz (indicated by the
thicker blue line), and the fit curve is the cyan dashed line in Figure 5(b).
For the 89 mT and 190 mT analyses, the fit is shown in the insets of Figures
6(a) and (b) respectively.

The fits for the 13C peak widths $\sigma$ in Figure 6(c) show consistency
between the P1 measurement at 89 mT and the NV measurement at 190 mT. However,
the width of the NV peak at 13 mT is approximately 70% greater than the width
of the other two measurements. The reason for this is believed to be that at
13 mT, the 13C nuclear spins are in a regime where the NV hyperfine and Zeeman
terms are equal for a 13C nuclear spin [72]. Reinhard et al. [72] define a
critical radius around an NV center, where 13C spins outside this radius are
quantized along the direction of the external Zeeman field while 13C spins
inside this radius are quantized along the local hyperfine axis. We estimate
that the number of 13C nuclei within the critical radius is
$N_{c}(B_{0})=0.08\text{ mT}/B_{0}$, which is 6 at 13 mT and 0.4 at 190 mT.
The broader peak observed at 13 mT is due to the contributions of these
hyperfine-shifted spins.

Figure 6(d) shows the percentage error of the 13C frequency identified by the
Gaussian fitting of the spectral line. We calculate the expected
(“literature”) value as $\omega_{\text{lit.}}=B_{0}\gamma_{\text{13C}}$, and
we calculate $B_{0}$ as $|(D\pm\omega_{\text{spec}})/\gamma_{e}|$ with $\pm$
being $-$ for 13 mT and $+$ for 190 mT. The fractional uncertainty in the
frequency measurement is $\delta
D/D+\delta\omega_{\text{spec}}/\omega_{\text{spec}}$. Based on the
temperature-dependent studies of the NV resonance spectrum in Reference [73]
and the temperature fluctuations in our laboratory, we use a value of
$D=2.866(10)$ GHz for the 13 mT measurement and $D=2.865(10)$ GHz for the 190
mT measurement, which was slightly warmer due to heat from the magnet coils at
the higher field. $D=0$ for the P1 center. The uncertainty in the effective
spectrometer frequency depends on the centering and symmetry of the spectral
line which forms the echo signal. Example echoes obtained in echo-detected
field sweep spectra are shown in Figure 9. We used
$\omega_{\text{spec}}=2.490(3)\times 2\pi$ rad GHz in the NV experiments and
$\omega_{\text{spec}}=2.500(3)\times 2\pi$ rad GHz in the P1 experiment. For
the P1 center, the fractional uncertainty of 13C Larmor frequency is 0.0012,
and for the NV centers it is 0.005. These uncertainties are represented by the
gray bands in the background of Figure 6(d). At all fields, the measured value
is accurate within the uncertainty of the calculation. At 13 mT and 89 mT, the
calculated uncertainties are comparable to the 95% confidence intervals of the
fits. At 190 mT, the frequency measurement has notably higher precision,
demonstrating the utility of the harmonic analysis for frequency measurements
at intermediate magnetic fields.

Characterizing the amplitude of the 13C spectral peak is less straightforward
than characterizing its frequency position or width. Consider the
reconstructed $\chi=3.5$ fundamental spectrum of the 13 mT experiment in
Figure 5(b). To reconstruct the fundamental spectrum, we used the Gaussian fit
parameters from the primary peak at 140 kHz. The reconstructed subharmonic at
48 kHz appears to have only about half the amplitude of the subharmonic
obtained directly from the experimental data. The subharmonics of the 89 mT
and 190 mT data in Figures 5(a) and (b) show this trend as well. The size of
the peak as measured with different harmonics is shown in Figure 16. We
emphasize that the spectra obtained from the harmonic analysis shown in
Figures 6(a) and (b) represent the average effective 13C peak over the valid
range of inter-pulse delays $\tau$ and echo numbers $N$, as explained in
Appendix D. Some of the results on a single NV center in Reference [68] also
seem to suggest that the measurements of the 13C peak with higher harmonics of
the filter function return a larger value for the amplitude than measurements
with lower harmonics.

Measuring sharp peaks in the power spectrum is significantly affected by flip
angle (over/under rotation) errors of the CPMG $\pi$ pulse, which cause an
inversion in the middle of the spectral peak [58] and can shift the apparent
frequency. It is also significantly affected by finite pulse widths, which
effectively reduce the amplitudes of the peaks of the higher harmonics in the
filter function. We have taken these considerations into account and explored
their effects numerically and experimentally in Appendix F.

##  V Summary and outlook

Using CPMG pulse sequences, we have characterized the spectrum of dephasing
noise for bulk measurements of NV centers in HPHT diamond with [NV] $\approx$
2 ppm and [N] $\approx$ 90 ppm, and for P1 centers in CVD diamond with [N]
$\approx$ 0.4 ppm. In our analysis, we first obtained the broad spectral
character. We then used that background characterization to improve the
measurement of the narrow peak caused by the 13C transverse nuclear hyperfine
interaction.

The clearest picture of the broad character of the dephasing noise spectrum
was obtained for NV centers at 190 mT. Interestingly, it exhibits a frequency
dependence best approximated as $1/\omega^{\alpha}$ with $\alpha=0.9(2)$.
Measurements of the NV center at 13 mT and of the P1 center at 89 mT seem to
also agree with this model. For the HPHT sample, we suspect that the
$1/\omega$ noise is caused by the high concentration and heterogeneity of P1
centers in the sample, which is supported by recent studies on clustering [27,
29, 26, 28]. In the CVD diamond, it is likely that the approximately
$1/\omega$ spectrum emerges from a mix of P1 and 13C nuclear interactions and
that P1 clustering plays a significant role. However, P1 clustering has not
been studied to the same extent in CVD diamond. To test these hypotheses, one
could perform a more systematic survey of how the power spectrum varies from
site to site throughout a sample, as well as how the spectrum varies with
measured ensemble size ranging from the single spin limit to the bulk sample
limit.

By performing the harmonic analysis, we have demonstrated more precise
measurements of the 13C frequency in the 89 mT and 190 mT experiments than
could be obtained in the 13 mT experiment. The accuracy of these measurements
is limited not by the methodology, but by uncertainties in the zero-field
splitting parameter $D$ and in the effective resonant frequency of the central
spins. We emphasize that the harmonic analysis enables the measurement of
peaks at frequencies of $\omega>1/\tau_{\text{min}}$, where
$\tau_{\text{min}}$ is the shortest $\tau$ that can be applied with the
experimental hardware.

Finite pulse effects and flip angle errors do not greatly affect the broad
spectrum characterization, but can have non-negligible effects on the
characterization of sharp peaks in the power spectrum. Finite pulses make CPMG
less sensitive to peaks at higher frequencies. Flip angle errors will cause an
asymmetric splitting of a spectral peak.

DD is used to characterize noise in a wide variety of quantum systems. The
methods presented in this work can be applied generally when CPMG is used to
characterize noise. These methods applied to diamond could very directly
inform strategies for quantum sensing, e.g., for performing nanoscale NMR [10]
with electron spins as proxy sensors.

##  VI Acknowledgements

We thank Bhargava Thyagarajan, Linta Joseph, and James Logan for helpful
discussions. We thank Lihuang Zhu for helping design the pEPR spectrometer. We
are grateful to Kevin Villegas Rosales and Yoav Romach of Quantum Machines for
helping set up and program the OPX as well as helpful discussions. We thank
Johan van Tol for help with diamond sample characterization that was performed
at the National High Magnetic Field Laboratory in Tallahassee, FL. We thank
Gajadhar Joshi and Jonathan Friedman for helpful insight on the loop-gap
resonator design. We thank Dwayne Adams and Christopher Grant for help
building the EPR probe and resonator. This work was partially supported by
funding from the National Science Foundation under grant CHE-2203681 and
cooperative agreement OIA-1921199 and by the Gordon and Betty Moore Foundation
by grant GBMF12251. We acknowledge support of a QISE-NET Triplets Award with
EQW funded by NSF award DMR-1747426.

##  Appendix A Diamond samples

There are two diamond samples used in this work. Sample A is a high pressure,
high temperature (HPHT) Type Ib diamond from Element 6 with a stone-type cut
with approximate dimensions $4.2\times 5.7\times 1.3$ mm3. Sample B is a
chemical vapor deposition (CVD) diamond plate also from Element 6 with
approximate dimensions $3.25\times 3.17\times 0.28$ mm3. Here we discuss the
$T_{1}$ and $T_{2}$ experiments that we performed to characterize these
samples. Data are shown in Figure 7, and resulting fit values are listed in
Table 2.

Sample A has been processed for the generation of NV centers. The sample
received a 1 MeV electron irradiation dose of $2\times 10^{18}$/cm2 per side
followed by 2 hours of annealing at 800∘ C performed by US Diamond
Technologies. This sample is reported by Element 6 to have a nitrogen
concentration of [N] $\lesssim$ 200 ppm. van Wyk et al. [24] measured the Hahn
echo $T_{2}$ of P1 centers to depend on nitrogen concentration as
$T_{2}([\text{N}])=71\ \mu\text{s}\cdot\text{ppm}/[\text{N}]$. We have
performed Hahn echo experiments on the P1 center (central ($m_{I}=0$) 14N
hyperfine manifold) and fit the decay to a bi-exponential of the form.

$L_{\text{HE}}(t_{d})=M_{1}e^{-t_{d}/T_{2}^{1}}+M_{2}e^{-t_{d}/T_{2}^{2}}.$
(32)

A bi-exponential was found to give the best fit for early and late decay
times, which is indicative of the bulk sample containing a heterogeneous
mixture of nitrogen concentrations with some more dense clustering. Recent
studies suggest clustering of P1 centers is prevalent in HPHT Type Ib samples
[26, 29, 27, 28]. Measurement of the P1 Hahn echo decay curve is shown in
Figure 7(d). The dominant relaxation term in the fit, $T_{2}^{1}$, would imply
a concentration of $[\text{N}]=87(8)$ ppm, with the reported uncertainty based
on the 95% confidence bound of the $T_{2}^{1}$ fit. To estimate the
concentration of NV centers in Sample A, we multiply the nitrogen
concentration by the conversion efficiency of the NV synthesis process. The
conversion efficiency, estimated at 2%, was determined by integrating the
respective P1 and NV spectral lines of a continuous wave EPR spectrum obtained
on another HPHT diamond sample that was in the same NV synthesis batch with
Sample A. This spectrum was performed at the National High Magnetic Field
Laboratory in Tallahassee, FL. We estimate the NV concentration in Sample A to
be 1.7(3) ppm.

Sample | A | A | A | B
---|---|---|---|---
Species | NV | NV | P1 | P1
$B_{0}$ (mT) | 190 | 13 | 89 | 89
Conc. (ppm) | 1.7(3) | 1.7(3) | 87(8) | 0.39(9)
$T_{1}$ (ms) | 10(1) | 6.3(4) | 2.1(2) | 1.9(2)
$M_{1}$ (arb.) | 1.00(5) | 1.00(9) | 4.4(4) | 1.00(5)
$T_{2}^{1}$ ($\mu$s) | 24.3(3) | 49(13) | 0.82(5) | 181(34)
$M_{2}$ (arb.) | - | - | 0.16(2) | -
$T_{2}^{2}$ ($\mu$s) | - | - | 12(1) | -
Table 2:  Fit parameters for $T_{1}$ and $T_{2}$ experiments performed on the
samples in this work. $T_{2}$ data for the NV centers are shown in Figure
2(c). All other data are shown in Figure 7. The concentration
characterizations for the P1 centers are based on $T_{2}$ measurements and the
van Wyk [24] characterization. The concentration characterization for the NV
centers is based on the P1 concentration in Sample A and continuous wave EPR
spectra obtained at the National High Magnetic Field Laboratory. Uncertainties
are based on the uncertainty of Hahn echo $T_{2}$ fits. We emphasize that, due
to the heterogeneity of nitrogen [26, 29, 27, 28], these concentrations are
only rough averages over the whole sample. $T_{1}$ is in reference to
Equations 33 and 34. For the NV centers in Sample A and the P1 centers in
Sample B, the parameters $M_{1}$ and $T_{2}^{1}$ refer to fitting the Hahn
echo decay to an exponential $M_{1}e^{-t_{d}/T_{2}^{1}}$. For the P1 centers
in Sample A, the fit parameters $M_{1}$, $T_{2}^{1}$, $M_{2}$, and $T_{2}^{2}$
refer to the fit with Equation 32. Figure 7:  $T_{1}$ and $T_{2}$ data.
(a) Inversion depolarization, Sample A, NV center, 13 mT.
(b) Inversion depolarization, Sample A, NV center, 190 mT.
(c) Inversion recovery, Sample A, P1 center, 89 mT.
(d) Hahn echo, Sample A, P1 center, 89 mT.
(e) Inversion recovery, Sample B, P1 center, 89 mT.
(f) Hahn echo, Sample B, P1 center, 89 mT.

We have also fit the Hahn echo decays of the NV center in the HPHT sample.
These are the CPMG-1 decays shown in green in Figure 2(c) and are extracted
from the same data that are analyzed in the main text. These are adequately
represented by a simple exponential decay $L_{\text{HE}}=Me^{-t_{d}/T_{2}}$.
For the NV center at 13 mT, we apply the fit only to the extracted peaks of
the hyperfine modulations. According to Bauch et al. [25] who performed a
survey of ODMR experiments on NV center ensembles in diamonds of widely
varying nitrogen concentration, the NV $T_{2}$ is given by $T_{2}=(160\pm
12)\,\mu\textrm{s ppm}/[\textrm{N}]\textrm{ ppm }$. Based on this
characterization, we would expect coherence times of only 0.8(2) $\mu$s for
the NV center. However, as noted earlier, this relationship may not work well
in the presence of heterogeneous P1 distributions. Our longer $T_{2}$ values
may be measured from centers that have lower local P1 concentrations. Results
from a recent theoretical investigation [74] suggest a $T_{2}$ that is closer
to our measured value but is still lower by about a factor of 10.

Sample B is reported to have $\text{[N]}<1$ ppm. Based on the van Wyk
characterization [24], the 180(30) $\mu$s $T_{2}$ value we measured indicates
a concentration of [N] = 0.39(9) ppm.

To characterize $T_{1}$ in our samples we performed inversion recovery (IR)
experiments: $\pi-\tau_{1}-\pi/2-\tau_{2}-\pi-\tau_{2}-\text{acq.}$, where
$\tau_{1}$ is varied over a wide range of delay times, $\tau_{2}$ is fixed at
a small value for efficient readout, and “acq.” indicates the acquisition of
the echo signal. For the P1 characterization, we fit the echo amplitude data
to

$S_{\text{IR}}(t_{d})=M(1-2\sin^{2}(\theta/2))\exp\left[-t_{d}/T_{1}\right],$
(33)

where $t_{d}=\tau_{1}+2\tau_{2},$ $M$ is a normalization term, and $\theta$ is
the inversion flip angle whose deviation from $\pi$ indicates imperfect
inversion of the spin ensemble. For the NV centers, the same pulse sequence is
used, however it is preceded by the laser hyperpolarization pulse. The NV
centers’ thermal equilibrium signal is less than 1% of the magnitude of the
hyperpolarized signal. Accordingly, we refer to the NV $T_{1}$ experiment as
inversion depolarization (ID), and fit the decay to

$S_{\text{ID}}(t_{d})=M\exp\left[-t_{d}/T_{1}\right].$ (34)

For the experiments on Sample A, the diamond stone was oriented with
$[111]\parallel\vec{B}_{0}$, whereas with the P1 experiments on Sample B, the
diamond plate was oriented with $[100]\parallel\vec{B}_{0}$. The $T_{1}$ and
NV concentration values are in fair agreement with those reported in [75],
albeit at higher magnetic fields.

##  Appendix B Echo-detected field-swept spectra

Figure 8:  Echo-detected field sweeps. Each column shows a different way of
visualizing the data for the experiment (sample and spin species) specified by
the label at left. The samples are described in Appendix A. Column (i) shows
the time trace of the echo signal with phase indicated by the color and log of
the magnitude indicated by the brightness. Column (ii) shows the negative log
of the signal magnitude, i.e. $\chi$, indicated by the color. Columns (iii)
and (iv) are likewise visualizations of the phase and magnitude, but they show
the Fourier transforms of each echo time trace. The vertical axes indicate the
magnetic field. Figure 9:  (a) Example time traces of echoes extracted from
the middles of the echo-detected field sweep spectra shown in Figure 8. The
real and imaginary component of each echo is shown, phase-corrected so that
the real component has maximum amplitude. (b) Fourier transforms of the echo
time traces in (a). (c) Echo-detected field sweep spectra. These are plots of
the real and imaginary components along the line Frequency = 0 of the image
plots shown in Column (iii) of Figure 8. The Equivalent frequency offset is
calculated by multiplying the magnetic field displacement by the gyromagnetic
ratio for the electron.

Examples of echo-detected field-swept spectra are shown in Figures 8 and 9.
Time traces of echoes obtained at a range of different magnetic fields are
shown in Figure 8 Columns (i) and (ii) which show the phase and log magnitude
of the signal respectively. In the time traces, Time = 0 is the middle of the
echo. A Blackman window is applied to the time traces to suppress noise at the
edges of the acquisition buffer. Columns (iii) and (iv) show the phase and log
magnitude of the Fourier transforms of each time trace.

Figure 8 Row (a) shows a wide field sweep over the NV center around 190 mT (a
cropped version is shown in Figure 1(e)). The CPMG experiments, referred to
throughout this work as the 190 mT experiments, were performed by setting the
magnetic field to the middle of the main spectral line, which occurs around
191.2 mT. The outer peaks near 189 mT and 194 mT are due to hyperfine coupling
between NV centers and nearby 13C nuclei. The primary NV spectral line can be
seen in detail in the narrower field swept range of 13 mT data of Row (b). The
NV echo time traces show a ripple, which, by examining the Fourier transforms,
is understood to be an interference pattern caused by the excitation of the
three spectrally separated lines arising due to the hyperfine interaction with
the host 14N nuclei, clearly visible in the frequency trace shown in Figure
9(b).

Figure 8 Row (c) shows the Sample A P1 spectrum. The four outermost spectral
lines are those arising due to the hyperfine interaction with the host 14N
nucleus. Due to the Jahn-Teller distortion [62, 25] and Sample A being in the
$\vec{B}_{0}\parallel[111]$ orientation, the ratio of the magnitudes of the
inner to outer hyperfine-shifted lines is 3 to 1. Figure 8 Row (d) shows the
Sample B P1 spectrum over a narrower field range than the Sample A P1
spectrum.

In Figure 9(a) and (b), the example time traces and corresponding Fourier
transforms are presented for signals acquired in the middle of the primary
peaks of the respective spectra. Figure 9(c) shows the spectrum obtained by
taking the dc value of the Fourier transformed time-traces at each magnetic
field. The comparison of Figures 9(b) and (c) reveals how the spectral line is
embedded in the frequency content of the echo signal.

##  Appendix C Spectrometer and data collection

In order to explore the coherence properties of bulk electron spin systems we
have constructed a pulse electron paramagnetic resonance (pEPR) spectrometer
with an operating frequency of 2.5 GHz. It is worth emphasizing that whereas
most other NV diamond experiments use optically detected magnetic resonance
(ODMR), we use microwave detection, which is conventional in pEPR. We measure
a microwave signal that gives the magnitude and phase of the transverse
component of the net magnetization of the sample [76]. This is achieved with
the aid of a loop-gap resonator [63, 64], which houses the diamond sample as
depicted in Figure 1(b). The loop-gap resonator is tuned by inserting a shard
of dielectric (sapphire) into the gap. The coupling is adjusted by positioning
an antenna a few mm above the top of the resonator [65]. We are able to adjust
the $Q$-factor of the resonator over a wide range this way. For the P1
experiments with Sample B, $Q=2000$, and for the NV experiments with Sample A,
$Q=190$. The characteristic ring-down time is $2Q/\omega_{0}$, which is
approximately 250 ns and 24 ns for the respective experiments. The large $Q$
was used for the Sample B experiments to enhance sensitivity. The finite ring-
down effects for these different $Q$ values are discussed in Appendix F.

Figure 10:  Example recordings of a CPMG experiment with $\tau=5.4\ \mu$s
recorded in the 89 mT data on the P1 centers of Sample B. The time of the echo
$t_{d}$ is indicated in the bottom plot and is common to all 4 plots. The top
two plots show the time traces of the echo, with 0 on the vertical axis
indicating the center of the echo. The third plot shows the Fourier transform
of each echo time trace. For the first and third plot, the color indicates
signal phase, and the brightness indicates signal magnitude. The signal
magnitude for the time traces is also indicated in the second plot. The echo
magnitude, shown in the fourth plot, is the integral of the Fourier
transformed time trace.

The schematic of our pEPR spectrometer is shown in Figure 11. This is the
instrument that was used to perform the experiments on Sample A, the HPHT
diamond. This includes all of the NV experiments, as well as the P1
characterization shown in Figures 8 and 9. Here we provide details of the
spectrometer’s construction and performance. The experiments on the P1 centers
in Sample B, the CVD diamond, were performed on an earlier version of the
instrument with a few differences mentioned later.

Figure 11:  Microwave circuit schematic of the pEPR spectrometer. See text for
details.

The microwave source is a Windfreak SynthUSB3, which outputs a continuous
local oscillator (LO) signal at 2.8 GHz and +8.00 dBm. The LO signal is
amplified (Amp #1, Mini-Circuits ZX60-6013E-S+) and then split (Splitter,
Mini-Circuits ZAPD-30-S+) between the Transmitter and Receiver arms of the
spectrometer. In the Transmitter, the LO is mixed (IQ Mixer, Marki IQ1545LMP)
with the analog outputs of the OPX (Quantum Machines) in a single-sideband
mixing configuration. The OPX outputs pulses defined on a 300 MHz carrier to
the I and Q ports of the mixer. The lower-frequency sideband coming out of the
mixer contains the pulse on a 2.5 GHz carrier. The pulse is amplified (Amp #2
same as #1) then sent through Low-Pass Filter #1 (Mini-Circuits VLF-2250+) and
the Band-Pass Filter (Anatech Electronics, AB2593B765). Although the single-
sideband mixing configuration with corrective biasing should eliminate
unwanted sidebands, LO leakage, and harmonics, the extra filtering ensures the
attenuation of any signal other than the desired pulse. The pulse signal is
then amplified by the power amplifier (Amp #3, Mini-Circuits, ZHL-2425-250X+),
which we restrict to a maximum power of +49 dBm. Switch #1 (DBwave
Technologies, PASR0200500600B) transmits the pulse when the spectrometer is in
transmit mode, and routes amplifier noise into a 50 $\Omega$ termination at
all other times. The switching is controlled by an OPX TTL line. The pulse
signal enters port 1 of the 4-port circulator (Ditom Microwave, D4C2327) and
is routed out port 2 to the loop-gap resonator. The echo signal returns to
port 2 and is routed to the receiver. Port 4 of the circulator is capped with
a 50 $\Omega$ termination to isolate the transmitter and damp out stray
reflections. Limiter #1 (Fairview Microwave, FMMT1021) is inserted at the
start of the receiver chain to protect Switch #2 (Mini-Circuits, ZFSWA2-63DR+)
from the reflections of the high-power pulses. When receiving an echo signal,
Switch #2 routes the echo to the amplifier chain, and at other times it routes
signals to a 40 dB, 50 W attenuator (Mini-Circuits BW-N40W50+) followed by a
50 $\Omega$ termination. After Switch #2, the echo signal is picked up by two
low-noise amplifiers in series: Amp #4 (RF-Lambda, RLNA02G02G) followed by Amp
#5 (Mini-Circuits, ZQL-2700MLNW+). Limiter #2 (Mini-Circuits, ZFLM-252-1WL-S+)
prevents any damage from pulse reflections that leak through Switch #2. The
echo signal is amplified once more (Amp #6, same as Amps #1 and #2) before
being downconverted from 2.5 GHz to 300 MHz by the receiver IQ mixer (same
kind as transmitter IQ mixer). The I and Q outputs of the receiver IQ mixer
are sent through Low-Pass Filters #2 (Mini-Circuits, SLP-1650+) to remove the
5.3 GHz sideband. The filtered in-phase and quadrature components of the echo
signal are passed through equal-length transmission cables to bias tees (Mini-
Circuits, ZFTB-4R2G+) for AC coupled analog signal acquisition by the OPX. The
OPX also controls the triggering of the laser controller (CivilLaser,
LSR532H-1W-FC). The 532 nm laser light is transmitted to the sample via an
optical fiber (Thorlabs, FG365UEC), the ferrule of which is mounted in the
3D-printed assembly as depicted in Figure 1(b).

The experiments on P1 centers in the CVD diamond sample were performed on an
earlier version of this spectrometer. In the earlier version, the power
amplifier in place of Amp #3 (Spectrian, 20 W 2.4 GHz) output 30 dBm, and in
place of Switch #1 were two Mini-Circuits ZFSWA2R-63DR+ switches. It did not
include any laser optics. The earlier version did not have an OPX. Instead, a
Tektronix AWG7052 was used for pulse signal generation on a 100 MHz carrier
which was up-converted with a 2.4 GHz LO to 2.5 GHz. The AWG7052 was also used
for triggering the switches and the digitizer (SP Devices ADQ214). A SpinCore
PulseBlaster ESR-PRO-500-PCI was used to trigger the AWG7052.

A schematic of the CPMG sequence is shown in Figure 2(a). The NV centers in
the diamond sample are initialized with 532 nm laser light. The initialization
laser pulse is 10 ms in duration. 10 ms was found to be the laser pulse
duration at which the NV signal magnitude saturates. The laser is controlled
simply by triggering the laser controller. The controller has a maximum
modulation rate of 30 kHz. Accordingly, we include a 100 $\mu$s delay after
the end of the laser pulse before starting the microwave pulses. It is worth
noting that this initialization process is significantly different from ODMR
NV experiments, which typically have laser initialization pulses and delays on
the order of a few $\mu$s, achieved with acousto-optic modulators. The most
significant reason for the difference is that we are sending unfocused light
into an entire diamond sample, instead of focused light into a microscopic
region, so the intensity of light in our sensitive volume is accordingly much
lower than that in an ODMR experiment.

pEPR measurements are weak measurements in that the coupling between an
electron spin and the resonator is weak enough that there is no back-action
that could disrupt the state of the spin. This enables us to perform
stroboscopic acquisitions of the CPMG echoes. For a single $\tau$ value, in a
single shot, we capture every echo out to a maximum detection time of about 3
ms, at which point the echoes have largely faded into the noise. For each
acquisition, the signal is demodulated and integrated in order to obtain the
amplitude of the echo. Examples of the echo time traces and processing steps
for a CPMG experiment on the P1 centers in Sample B are shown in Figure 10. In
this example, the recording buffer is 3 $\mu$s in length. We apply a Blackman
window to the time trace before integration. For the CPMG experiments on NV
centers in Sample A, the acquisition buffer is 400 ns. The demodulation is
performed onboard the OPX, and we record just a pair of values (I and Q), such
that the echo magnitude is $\sqrt{\text{I}^{2}+\text{Q}^{2}}$. After the
microwave sequence is finished, the NV centers are re-initialized with the
laser, and the sequence is repeated for averaging. For the NV center
experiments we perform 128 averages.

For the P1 measurements, there is no direct hyperpolarization mechanism, so no
laser light is shown on the sample. The P1 signal magnitude is proportional to
the polarization of the P1 center in the applied $B_{0}$ field at thermal
equilibrium (see Equation 6). We wait 5 ms, $\approx 2.5\times T_{1}$, in
between shots. For the P1 center data presented here, 65,536 averages were
performed.

The minimum values of the interpulse delay $\tau$ that we can achieve are
limited by a variety of factors. For the P1 CPMG experiments, the minimum
$\tau$ presented in this work is 4 $\mu$s. At lower $\tau$ values, the 250 ns
ring-down of the resonator (discussed in the first paragraph of this Appendix)
begins to interfere with the echo signal. In the experiments on the HPHT
diamond, the minimum $\tau$ value presented is 1.4 $\mu$s. The rise/fall times
of the switches plus the finite pulse and recording durations put the minimum
possible $\tau$ around 1 $\mu$s. Additionally, we find that Limiter #1
saturates for $\tau$ of 1.4 $\mu$s after 3 ms of applying $\pi$ pulses. We
restrict the $\tau$ and sequence duration to these values to avoid saturating
the limiter.

We performed the NV CPMG experiments at $\tau$ values ranging from 1.4 $\mu$s
to 560 $\mu$s. As explained in Appendix D, the harmonic analysis for the 13C
peak characterization is restricted to $\tau\leq 10\ \mu$s. For the 190 mT
experiment, $\tau$ was sampled linearly up to $\tau=10\ \mu$s, with a step of
8 ns. Above $\tau=10\ \mu$s, $\tau$ was sampled logarithmically for
efficiency. For the 13 mT experiment, the linear to logarithmic sampling
transition was at 6.8 $\mu$s. In executing an experiment, the $\tau$ values
were selected in a random order. All 128 averages are acquired for the CPMG
decay curve with one $\tau$ value before continuing to the next $\tau$ value.
One acquisition loop iteration consists of the laser initialization (10.1 ms),
the microwave pulse sequence (3 ms), and a hardware recovery delay (60 ms),
such that acquiring all 128 averages for a given $\tau$ takes about 9.3 s.
Between runs for different $\tau$ values, the data is saved, and the magnetic
field control loop is executed, which typically takes 4.5 s. So, on average,
it takes 13.8 s to acquire the data for a given $\tau$ and proceed to the
next. For the 13 mT experiment, a total of 4356 $\tau$ values were sampled,
taking about 16.7 hours. For the 190 mT experiment, a total of 6447 $\tau$
values were sampled, taking about 24.7 hours.

For the P1 CPMG experiments, $\tau$ was sampled linearly from 4 to 41.3 $\mu$s
with a step of 100 ns. No laser was used. The microwave pulse sequences were
applied for 5 ms, followed by a spin system recovery delay of 5 ms (about
$2.6\times T_{1}$). This sequence was repeated until the digitizer’s memory
limit was reached, at which point each individual time trace was transferred
to the lab computer, and field correction was performed. This was repeated
until the 16,384 averages of the CPMG decay were obtained for a given $\tau$
value and a given global phase of the microwave pulses. This was repeated for
global phases of the control pulses of 0, $\pi/2$, $\pi$, and $3\pi/2$, and
the signals were added with appropriate phase corrections for a total of
65,536 averages. This phase-cycled averaging eliminated phase-independent
systematic noise. Due to the much higher number of averages needed and the
overhead time of data transfer, this experiment took 31 days to complete with
an average of 2 hours to acquire data for each of the 373 $\tau$ values.

##  Appendix D Harmonic Analysis

Noise spectroscopy strategies such as approximating the filter function as its
fundamental peak (Equation 19) to compute $S_{P}(\omega)$ or using Álvarez and
Suter’s matrix inversion method [40] enable one to measure the power spectral
density up to the frequency $\omega=\pi/\tau_{\text{min}}$, where
$\tau_{\text{min}}$ is the shortest $\tau$ one can implement with their
experiment hardware. In applications of ac magnetic field sensing, it may
occur that there are peaks in the power spectrum at frequencies a few times
greater than $\pi/\tau_{\text{min}}$, whose effects can still be seen in
experiments in which $\tau$ is swept with a small step size. It is desirable
to obtain a wide spectrum that includes these peaks, so that the peaks can be
measured with sufficient contrast to the background. Here we present an
analysis method for obtaining a wide spectrum that reveals prominent features
at frequencies greater than $\pi/\tau_{\text{min}}$.

Figure 12:  Illustrations of filter considerations in the harmonic analysis.
See text for details.

In essence, the harmonic analysis breaks the range of high frequencies into a
series of blocks of width $W$, which are instances of a scanning window. For
each instance of the window, Equation 31 is computed for each valid
combination of $\tau$ and echo number $N$ for which a harmonic $m_{s}$ of the
filter function is present within the window. Each calculation gives a point
in the spectrum at frequency $\omega=m_{s}\pi/\tau$. Interpolation and
averaging are used to convert the set of calculated spectral points into one
continuous spectrum. Results from this method are presented in Figures 6(a)
and (b) and 20(g) and (h).

In our case, the shortest interpulse delay $\tau$ that we can apply is 1.4
$\mu$s, as explained in Appendix C. The fundamental ($m=1$) frequency of a
CPMG filter function is $\omega=\pi/\tau$. So the highest frequency that we
can fundamentally measure is $\omega\approx 360\times 2\pi$ rad kHz, which can
be seen as the maximum frequency of the spectra in Figures 5(a) and (b) and
20(e). In those spectra however, the subharmonics of the 13C peak are clearly
visible. These subharmonics, which correspond to CPMG sequences with faster
decays, occur at frequencies for which a higher harmonic of the filter
function overlaps with the 13C peak, i.e., when
$\omega=\omega_{\text{13C}}/m$, where $\omega_{13\text{C}}$ is the 13C nuclear
Larmor frequency and $m$ is an odd integer, as illustrated in Figure 3(a).

The valid subharmonics are ones which meet certain criteria for avoiding
sensitivity limitations, avoiding distortions caused by finite widths of the
peaks in the filter function, and avoiding distortions caused by multiple
filter peaks overlapping with a single spectral line. These three criteria
define a maximum and minimum detection time $t_{d}$ and a maximum $\tau$
respectively. These limits can be thought of as defining a sensing region. For
the 89 mT and 190 mT data analyzed in this work, the sensing regions are those
indicated by the red dashed rectangles in Figures 14(a) and 15(a). The
criteria for avoiding sensitivity and distortion limitations inform the
following 6-step analysis protocol. Steps 1 through 5 explain the frequency
scan procedure over one instance of the window, focusing on an example where
the window contains a prominent peak, as depicted by the cartoon examples in
Figure 12. Step 6 explains how to iterate this procedure to obtain a wide
spectrum.

Figure 13:  The top plot shows a zoomed in section of the image coherence plot
shown in Figure 4(b). The bottom plot shows the value of $\chi$ along the
$\chi=2$ power law fit line. This figure serves as a visual example of the 13C
peak $S_{P}$ calculation described in Step 5 of the harmonic analysis
(Appendix D).

1\. Determine the width $W$ of the frequency-scanning window. $W$ is
determined by the requirement that only one harmonic of the filter function
overlap with the peak. This is the assumption that allows the isolation of
$S_{P}$ from the summation in Equation 29. $W$ is the maximum peak width that
can be resolved by the analysis process. $W$ sets the maximum $\tau$ value
that will contribute to the calculation as $\tau_{\text{max}}=2\pi/W$. The
center of the scan window is denoted $\omega_{c}$, and the upper $(+)$ and
lower $(-)$ bounds of the window are $\omega_{\pm}=\omega_{c}\pm W/2$. With
$\tau\leq\tau_{\text{max}}$, the minimum frequency difference between two
consecutive peaks of the filter function is greater than or equal to $W$, as
depicted in Figure 12(a). We note that in some systems $\tau_{\text{max}}$ may
be limited by a hardware or coherence limit, in which case $W$ should be
defined from $\tau_{\text{max}}$.

2\. Determine the range of harmonics, i.e., peaks of the filter $F$, that can
be used as the scanning harmonic $m_{s}$ (where $m_{s}$ can be odd integers in
reference to $m$ of Equation 18). For $m_{\text{min}}$ we want the lowest
filter peak available. With the shortest $\tau$ that can be applied with the
hardware, $m_{\text{min}}$ is the lowest peak that can still reach
$\omega_{+}$. For the example in Figure 12(b), which shows the filter for CPMG
with $\tau_{\text{min}}$, we get $m_{\text{min}}=7$ because the $m=5$ peak
falls short of $\omega_{+}$. Formally, $m_{\text{min}}$ is the lowest odd
integer greater than $\omega_{+}\tau_{\text{min}}/\pi$. As $\tau$ increases,
the peaks in $F$ contract toward 0. The maximum harmonic $m_{\text{max}}$ is
then determined by requiring that only one harmonic can be in the range
$[\omega_{-},\omega_{+}]$. Specifically, when the $m_{\text{max}}$ peak is at
$\omega_{-}=\pi m_{\text{max}}/\tau_{\text{max}}$, the $m_{\text{max}}+2$ peak
should be at a frequency greater than $\omega_{+}$:

$\frac{(m_{\text{max}}+2)\pi}{\tau_{\text{max}}}\geq\omega_{-}+W.$ (35)

This inequality is more succinctly expressed in terms of $W$. We say that
$m_{\text{max}}$ is the highest odd integer subject to the constraint

$m_{\text{max}}\leq\frac{2\omega_{-}}{W}.$ (36)

Figure 12(c) shows the filters for $\tau_{\text{min}}$ (light blue) and
$\tau_{\text{max}}$ (darker blue).

3\. For a given scan harmonic $m_{s}$ determine the starting and finishing
$\tau$ values to iterate over. Assuming we are sweeping $m_{s}$ in the
direction of increasing $\omega$ (as inferred by the increasing transparency
of the blue peaks in Figure 12(e)), then the $\tau$ values are swept in
decreasing order, $\tau_{\text{start}}>\tau_{\text{fin}}$. The corresponding
frequencies $\omega_{\text{start}}=m_{s}\pi/\tau_{\text{start}}$ and
$\omega_{\text{fin}}=m_{s}\pi/\tau_{\text{fin}}$ will not necessarily coincide
with $\omega_{-}$ and $\omega_{+}$. Rather, $\omega_{\text{start}}$ and
$\omega_{\text{fin}}$ should be within the bounds of $\omega_{-}$ and
$\omega_{+}$.

$\displaystyle\omega_{\text{start}}$ $\displaystyle\geq\omega_{-}$ (37)
$\displaystyle\frac{m_{s}\pi}{\tau_{\text{start}}}$
$\displaystyle\geq\omega_{-}$ (38) $\displaystyle\tau_{\text{start}}$
$\displaystyle\leq\frac{m_{s}\pi}{\omega_{-}}$ (39)

$\tau_{\text{start}}$ is the highest available $\tau$ subject to this
constraint. Similarly,

$\displaystyle\omega_{\text{fin}}$ $\displaystyle\leq\omega_{+}$ (40)
$\displaystyle\frac{m_{s}\pi}{\tau_{\text{fin}}}$
$\displaystyle\leq\omega_{+}$ (41) $\displaystyle\tau_{\text{fin}}$
$\displaystyle\geq\frac{\pi m_{s}}{\omega_{+}}$ (42)

$\tau_{\text{fin}}$ is the lowest available $\tau$ subject to this constraint.

Figure 14:  Visualizations of parts of the harmonic analysis procedure for the
89 mT P1 data. (a) Image plot of the coherence data. The red dashed rectangle
indicates the “sensing region” over which it is valid to use the $\delta$
approximation of the filter function to characterize the 13C peak. This is the
region that plots (b), (c) and (d) show. (b) Number of times each point of the
data is used in the harmonic analysis. (c) Log10 of the average harmonic
weighting prefactor in the calculation of Equation 31. $A(\omega_{s})$ is the
finite filter correction explained in Appendix F.1. (d) Average contribution
of each data point in the calculation of $S_{P}$ in the harmonic analysis. (e)
Spectral segments described in Steps 5 and 6 of the harmonic analysis. A
segment is obtained using one harmonic $m_{s}$ and averaged over all available
$N$. The vertical gray dotted line represents the 13C Larmor frequency. (f)
Detail of the area outlined by the gray dashed box in (e) with the legend
indicating which filter harmonic each segment was obtained with. (g) Maxima of
the peaks shown in (f). The black line is a quadratic fit. See Appendix E for
discussion. (h) Overlapping segments of width $W$ obtained when averaging the
contributions from all available harmonics as discussed in Steps 5 and 6 of
the harmonic analysis. These are then averaged to produce the spectrum in
Figure 6(a). Figure 15:  Visualizations of parts of the harmonic analysis
procedure for the 190 mT NV data. See the caption to Figure 14 for
descriptions. (f) does not include a legend as in Figure 14 because the
spectral segments are too dense to separately distinguish. The averaged
segments in (h) are averaged to produce the spectrum in Figure 6(b).

4\. For a given $m_{s}$ and associated range
$[\tau_{\text{fin}},\tau_{\text{start}}]$, determine the range of echo numbers
$[N_{\text{min}},N_{\text{max}}]$. To determine $N_{\text{min}}$, it must be
taken into account that the peaks in the filter function $F$ have finite
width, which becomes narrower with increasing $N$. Therefore, $N_{\text{min}}$
should be set by requiring that the harmonic peak is narrower than the sharp
peak needing to be resolved. Without knowing the sharp peak’s width
beforehand, one can instead require that the harmonic peak width be much
narrower than the scan frequency range $W$. With $\sigma=\sqrt{2\pi}/(N\tau)$
being the width of the scan harmonic peak (obtained from the Gaussian
approximation of the peaks in the filter $F$ [66]), we require
$\sigma<\epsilon W$, with $\epsilon\ll 1$. This leads to

$N_{\text{min}}>\frac{\sqrt{2\pi}}{\epsilon W\tau_{\text{fin}}}.$ (43)

$N_{\text{min}}$ is the lowest $N$ available subject to this constraint. In
Figure 12(d), filters are shown for a fixed $\tau$ and $N=1$ through 16. The
horizontal green bar has width $2\epsilon W$. The red filters are $1\leq N\leq
9$ and fail the width criterion. The blue filters are $10\leq N\leq 16$ and
satisfy the width criterion. $N_{\text{max}}$ should be determined by some SNR
threshold or $T_{1}$ limit. We use $N_{\text{max}}\tau_{\text{start}}<1$ ms,
which is sufficiently shorter than the 1.9 ms $T_{1}$ of the P1 center. Note,
in the definitions of $N_{\text{min}}$ and $N_{\text{max}}$, we have inserted
$\tau_{\text{fin}}$ and $\tau_{\text{start}}$ respectively, to place the more
restrictive constraints on the $N$ range. This ensures that the $N$ range
meets the requirements for the entire $\tau$ range and the given $m_{s}$.
Since the detection time is $t_{d}=N\tau,$ these width limits set a maximum
and minimum $t_{d}$, which are the tops and bottoms of the red dashed
rectangles in Figures 14(a) and 15(a).

5\. Compute $S_{P}(m_{s}\pi/\tau)$ using Equation 31. For each combination of
$m_{s}$ and $N$, a segment of the spectrum $S_{P}(\omega)$ is obtained by
sweeping $\tau$ from $\tau_{\text{start}}$ to $\tau_{\text{fin}}$. These
segments are averaged over the available $N$, then each of those segments are
averaged over the available $m_{s}$. Figure 13 provides example data from the
190 mT experiment to illustrate how Equation 31 is calculated. The effective
bounds of the scanning window represented in the $\tau$-domain by
$\tau_{\text{start}}$ and $\tau_{\text{fin}}$ are indicated by the vertical
dotted lines. In the middle of this window, we see more rapid decay due to a
filter harmonic overlapping with the 13C peak. $\chi(N,\tau)$ is the exact
value of $\chi$ at the point in question. To determine $\chi_{B}$ for a given
combination of $\tau$ and $N$, we determine which $\chi$ contour fit most
closely crosses over that point, and use the parameters from that fit to
calculate the precise value of $\chi_{B}$. We solve Equation 23 for the
various $\chi$ contour characterizations, and find which $\chi$ contour fit
predicts the closest $N$ to the point in question with the given $\tau$. For
example, the point at $\tau=4.15\ \mu$s and $t_{d}=N\tau=230\ \mu$s has the
$\chi=2$ contour fit (dashed line) crossing over it, so $\chi_{B}=2$ for that
point. The lower plot in Figure 13 shows the value of the measured $\chi$
along the power law fit of the $\chi=2$ contour. So the subtraction term in
brackets of Equation 31 is $\chi-\chi_{B}=2.5-2=0.5$. In experiments, it may
be the case that $\tau$ is sampled linearly with hardware-limited precision.
Therefore, the spectrum of $S_{P}(\omega)$ has $1/\omega$ sampling. The array
of $S_{P}$ values can be resampled, using interpolation, onto a linearly
spaced array of $\omega$ values with a small step size to create a segment of
the spectrum. The linearly resampled spectrum segments can then be added
together for each combination of $m_{s}$ and $\tau$ for averaging.

6\. Perform Steps 1-5 for a wide sweep of the window position $\omega_{c}$,
while maintaining the fixed window width $W$. The spectral segments of
adjacent ranges $[\omega_{-},\omega_{+}]$ can then be stitched together in
order to generate the spectrum over a wide frequency range. Furthermore, to
average out interpolation artifacts that may arise at the edges of the
scanning window, one can overlap the ranges such that any range starts at
$W/2$ greater than the $\omega_{-}$ of the previous range as in a bricklaying
pattern.

Figures 14 and 15 show intermediate steps of the harmonic analysis for the 89
mT and 190 mT data respectively. Figures 14(a) and 15(a) indicate the sensing
area with the red dashed rectangles. We set the maximum detection time $t_{d}$
to be 1 ms for both analyses to stay well within decoherence limitations. We
set $\tau_{\text{max}}=10\ \mu$s, which gives $W=2\pi\times 100$ rad kHz,
approximately 5 times wider than the width ($2\sigma$) of the 13C peak. We set
$W\epsilon=2\pi\times 3$ kHz, which is sufficiently narrower than the 13C peak
and corresponds to a minimum detection time of $t_{d}\approx 133\ \mu$s.
Figures 14(b), (c) and (d) and 15(b), (c) and (d) reveal how each data point
in the sensing region is incorporated into the harmonic analysis. The (b)
figures simply show how many times each point was used in the calculation of
Equation 31. In order to avoid the sensitivity limits of decoherence, points
were omitted if the $\chi$ value was greater than a noise threshold: 2 for the
89 mT data and 4.5 for the 190 mT data.

##  Appendix E Amplitude of the 13C peak

The challenge with measuring the 13C peak amplitude is that in the 13 mT data,
the 13C interaction is so strong that it causes complete decoherence at low
$N$, whereas in the 89 mT and 190 mT data, the peak seems to be greater in
magnitude when measured with larger $\tau$ (i.e., when a larger filter
harmonic $m_{s}$ overlaps with the peak).

The same $\tau$ and $t_{d}$ limits that constrain the harmonic analysis (red
dashed rectangles in Figures 14 and 15) apply to measuring the 13C peak in the
13 mT data as well. The width of the filter peak should be significantly less
than the width of the 13C peak, which is about 10 to 20 kHz. Approximating the
filter peak as a Gaussian gives $\sigma=\sqrt{2\pi}/t_{d}$. Accordingly, to
make use of the $\delta$-function approximation of the filter for the 13C peak
measurement at 13 mT, detection times of $t_{d}\gg\sqrt{2\pi}/(10\text{
kHz})\approx 25\ \mu$s ought to be used. It is clear however from Figure 4(a)
that the $\chi$ is nearing the noise threshold ($\chi\gtrsim 4$) at such
detection times in the middle of the first coherence dip ($\tau\approx 3.5\
\mu$s). In Figures 5(b) and 6(c), (d) and (e) we have presented and analyzed
the $\chi=3.5$ contour for the 13 mT data since it is mostly observed at
detection times later than 25 $\mu$s.

Figure 16:  (a) The 13C peak amplitude at 13 mT measured by fitting
fundamental spectra for different $\chi$ contours to a Gaussian of the form
$Ae^{-(\omega-\omega_{P})^{2}/(2\sigma^{2})}$ with $A,\omega_{P},$ and
$\sigma$ as fit parameters. (b) Peak amplitudes from all three main data sets
versus the filter harmonic $m_{s}$ that was used to scan the peak.

Figure 16(a) shows the amplitudes of the 13C peak obtained by applying
Gaussian fits to the different $\chi$ contours of the 13 mT data. For fits of
$\chi<3.5$, the detection times are early, and there is considerable
uncertainty due to the roughness of the $\chi$ contour at low $N$ (as in
Figure 5(a) with $\chi=1$). For fits of $\chi>3.5$ the signal-to-noise
(roughly $\log(5)/\log(\chi)$) becomes too small. Accordingly, the 13C peak
amplitude around $\chi=3.5$ is likely the most accurate measurement from the
13 mT data.

The 13C peak amplitudes presented in Figure 6(e) for the 89 mT and 190 mT data
represent the average over the sensing areas (red dashed rectangles of Figures
14(a) and 15(a)). These average 13C amplitudes are more heavily weighted
toward the measurements with lower harmonics $m_{s}$, which are obtained with
the shorter $\tau$ values. This is because there is a higher density of
acquisitions at shorter $\tau$ values. Figures 14(d) and 15(d) show the net
contributions of the data to forming the spectrum from the harmonic analysis.
Figures 14(g) and 15(g) show the amplitude of the 13C peak as obtained with
each available scan harmonic $m_{s}$. There is an apparent increase of the 13C
peak amplitude as measured with increasing harmonic number $m_{s}$. However,
this is likely largely due to the compounded uncertainty in the harmonic
analysis procedure. Based on Equation 31, and approximating
$\delta\chi\approx\delta\chi_{B}$, the uncertainty for the sharp peak,
$S_{P}$, measurement is

$\delta S_{P}\approx\frac{m_{s}^{2}\pi^{2}}{4N\tau}\sqrt{2}\,\delta\chi_{B}.$
(44)

Due to the $m_{s}^{2}$ term, $S_{P}$ tends to grow quadratically with $m_{s}$.
We have shown the peak amplitude fits as a function of the harmonic $m_{s}$
used to measure the peak for all three main data sets in Figure 16(b). To show
that the peaks in $S_{P}$ grow quadratically with $m_{s}$ due to uncertainty,
we have fit the 89 mT and 190 mT peaks as
$S_{P}^{\text{(peak)}}(m_{s})=p_{2}m_{s}^{2}+p_{0}$, where $p_{n}$ is the
coefficient of the $n^{\text{th}}$ order polynomial term. In the limit that
$\delta S_{P}=0$, our model predicts that $p_{2}=0$, and
$S_{P}^{\text{(peak)}}=p_{0}$. That is, in a classical Gaussian noise
spectrum, the measured amplitude of the peak should not depend on which filter
harmonic $m_{s}$ is used to measure it (taking into account the appropriate
weighting factors). These fits are depicted in Figures 14(g), 15(g), and 16(b)
and show approximate agreement with the measurements. For the 89 mT data we
obtain $p_{2}=5.4\pm 1.4$ kHz and $p_{0}=-160\pm 240$ kHz. For the 190 mT data
we obtain $p_{2}=1.9\pm 1.1$ kHz, and $p_{0}=400\pm 800$ kHz.

The difficulties of measuring the 13C peak amplitude are apparent. We can
place a lower bound on the peak amplitude of about 500 kHz. Likely, the most
accurate Gaussian approximation of the peak is slightly higher than this,
given by the lowest $m_{s}$ measurements from the 89 mT and 190 mT data. In
this sense, the values reported in Figure 6(e) are the most accurate
estimates. It is possible that the magnitude of the 13C noise varies with
magnetic field, but our measurements do not show that it is substantial.

##  Appendix F Numerical methods

###  F.1 Theory

The harmonic analysis method explained in this work relies on the high-
frequency content of the CPMG filter function to perform spectroscopy. In the
time domain, this high-frequency content is associated with sharp,
instantaneous transitions between $+1$ and $-1$ in the modulation function. In
real experiments, finite pulses smooth these transitions, creating a gradual
change between $+1$ and $-1$. This smoothing in the time domain corresponds to
a reduction of the high-frequency content of the filter. The modulation
function also takes on a perpendicular component during the transition. This
transverse component of the modulation function is also added to by flip angle
errors. Here we explain how we numerically obtain the leading order correction
to the decoherence calculation. We examine the case for CPMG with finite
pulses and flip angle errors. For a generalized treatment of noise
spectroscopy with multi-axis control, see Reference [77].

Figure 17:  Example calculations from the numerical investigations into finite
pulse effects. (a), (c), and (e) show the $x$ components and (b), (d), and (f)
show the $z$ components of the CPMG modulation function with $\tau=4\ \mu$s,
comparing the instantaneous pulses with no flip angle error (“Ideal”) to
finite pulses without errors (“Finite”) and with (“Error”) flip angle errors.
See Appendices F.2 and F.3 for details of the different sequence conditions.
(a) and (b) show a zoomed in view of the first pulse. (c) and (d) show the
entire modulation functions $f_{x}$ and $f_{z}$ for a 32-pulse sequence with
detection at 128 $\mu$s. The black lines in the background highlight the
approximately sinusoidal and cosinusoidal modulations to $f_{x}(t)$ and
$f_{z}(t)$ caused by the flip angle errors. (e) and (f) show a zoomed in view
of $f_{x}(t)$ and $f_{z}(t)$ for the first 6 pulses to highlight the
cumulating error effect of the flip angle error. (g) shows the overlap
integral of Equation 13 as a cumulative integral for the ideal and various
non-ideal pulses. Here we used a model power spectrum of $S(\omega)=5\times
10^{10}(\text{Hz})^{2}\omega^{-1}$, and a control sequence with $\tau=4\
\mu$s. The “$\delta$ Peak” data points are the cumulative sum from computing
the $\chi$ in the discrete form with Equation 18. (h) shows filter functions
calculated numerically from the ideal and non-ideal modulation functions. The
$2/(\pi\omega^{2}\tau^{2})$ dots are the filter peaks in the ideal limit. The
inset shows the fundamental (first) peak in detail, revealing the asymmetric
peak splitting caused by the flip angle errors.  Figure 18:  Results from
numerically computing $\chi$ arrays from a model power law spectrum
$S(\omega)=5\times 10^{10}(\text{Hz})^{2}\omega^{-1}$ and ideal and non-ideal
filter functions. $\alpha$ in (a) and $C$ in (b) are the recovered fit values
from analyzing the arrays with the same procedure that was performed for the
real experimental data.

In CPMG the initial density operator, following the $(\pi/2)_{x}$ pulse, is
$\rho_{0}=\frac{\mathbbm{1}}{2}-\epsilon S_{y}$. The control Hamiltonian for
CPMG is $H_{c}(t)=\omega_{1}(t)S_{y},$ where $\omega_{1}(t)$ is the amplitude
of the control. The control unitary is

$\displaystyle U_{c}(t)$
$\displaystyle=\exp\left[-iS_{y}\int_{0}^{t}dt^{\prime}\omega_{1}(t^{\prime})\right]$
(45) $\displaystyle=\exp\left[-iS_{y}\Phi(t)\right],$ (46)

where $\Phi(t)$ is the cumulative phase induced by the control pulses. In the
interaction picture defined by the Zeeman and environment Hamiltonians, the
effective system-environment Hamiltonian
$H^{(E)}_{SE}=-\gamma_{e}b_{z}(t)S_{z}$. In the toggling frame, which is the
interaction picture defined by the control pulses, the effective system-
environment Hamiltonian becomes
$\tilde{H}_{SE}(t)=U_{c}^{\dagger}(t)H^{(E)}_{SE}U_{c}(t).$ The multi-axis
filter function is
$f_{\mu}(t)=2\text{Tr}\\{S_{\mu}U_{c}^{\dagger}(t)S_{z}U_{c}(t)\\},$ which, in
the case of exclusively $y$-pulses, simplifies to

$\vec{f}(t)=(-\sin(\Phi(t)),0,\cos(\Phi(t))).$ (47)

We can then express the system-environment Hamiltonian as
$\tilde{H}_{SE}(t)=-\gamma_{e}b_{z}(t)(-\sin(\Phi(t))S_{x}+\cos(\Phi(t))S_{z}).$
To estimate the decoherence due to this interaction between initialization and
detection time ($t_{d}$), we approximate the Hamiltonian as a time-independent
Hamiltonian using the Magnus expansion. The leading order term of the Magnus
expansion is

$\displaystyle\tilde{H}_{SE}^{(0)}(t_{d})$
$\displaystyle=\frac{1}{t_{d}}\int_{0}^{t_{d}}\tilde{H}_{SE}(t)\,dt$ (48)
$\displaystyle\begin{split}&=-\frac{1}{t_{d}}\bigg{(}\int_{0}^{t_{d}}\gamma_{e}b_{z}(t)f_{x}(t)S_{x}\,dt\\\
&\ \ \ \ \ \ \ \ \ \
+\int_{0}^{t_{d}}\gamma_{e}b_{z}(t)f_{z}(t)S_{z}\,dt\bigg{)}\end{split}$ (49)
$\displaystyle=\frac{1}{t_{d}}(\phi_{x}S_{x}+\phi_{z}S_{z}),$ (50)

where $\phi_{\mu}$ is the net phase accumulated about the axis $\mu$ due to
the noise. Higher-order terms of the Magnus expansion contain nested
commutators of the form $[\tilde{H}_{SE}(t_{1}),\tilde{H}_{SE}(t_{2})]$. In
the case of perfect, instantaneous pulses $f_{x}(t)=0,$ so these commutators
are 0 and the leading order term is exact. With real pulses, the higher-order
terms will contribute to additional decoherence; however, in the weak-noise
limit, their effect is small enough to be ignored [77].

With this, we calculate the signal as the zeroth order approximation of the
normalized ensemble average value of $S_{y}$ at time $t_{d}$:

$\displaystyle\begin{split}&\langle s^{(0)}(t_{d})\rangle_{\text{ens}}\\\
&\qquad=\left\langle\frac{\text{Tr}\\{e^{-i\tilde{H}_{SE}^{(0)}(t_{d})t_{d}}\rho_{0}e^{i\tilde{H}_{SE}^{(0)}(t_{d})t_{d}}S_{y}\\}}{\text{Tr}\\{\rho_{0}S_{y}\\}}\right\rangle_{\text{ens}}\end{split}$
(51)
$\displaystyle\qquad=\langle\cos(\phi_{x})\cos(\phi_{z})\rangle_{\text{ens}}$
(52)
$\displaystyle\qquad=\int_{-\infty}^{\infty}d\phi_{z}\int_{-\infty}^{\infty}d\phi_{x}\cos(\phi_{x})\cos(\phi_{z})P(\phi_{x},\phi_{z}),$
(53)

where $P(\phi_{x},\phi_{z})$ is the distribution of the random variables
$\phi_{x}$ and $\phi_{z}$. Assuming $b_{z}(t)$ is a Gaussian random variable,
then so too are $\phi_{x}(t_{d})$ and $\phi_{z}(t_{d})$. We then obtain

$\langle
s^{(0)}(t_{d})\rangle_{\text{ens}}=e^{-\langle\phi_{x}^{2}\rangle/2-\langle\phi_{z}^{2}\rangle/2}$
(54)

leading, via the convolution theorem, to the overlap integral

$\displaystyle\chi$
$\displaystyle\equiv(\langle\phi^{2}_{x}\rangle+\langle\phi^{2}_{z}\rangle)/2$
(55)
$\displaystyle=\frac{t_{d}^{2}}{2}\left[\int_{-\infty}^{\infty}d\omega^{\prime}S(\omega^{\prime})(F_{x}(\omega^{\prime})+F_{z}(\omega^{\prime}))\right]$
(56)

where $S(\omega)$ is the power spectrum of $b_{z}(t)$, and the $\mu$-axis
filter function is

$F_{\mu}(\omega)=\frac{1}{t_{d}^{2}}\left|\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f_{\mu}(t)e^{-i\omega
t}\,dt\right|^{2}.$ (57)

We numerically compute $F$ by constructing $\Phi(t)$ as a finely sampled array
for a real pulse sequence based on our physical system, and then take the
square of the magnitude of the discrete Fourier transform.

Finite pulses cause the peaks of the filter function in the high-frequency
limit to be smaller than that given by Equation 17. Our strategy for
determining the correction involves numerically computing the approximate
filter based on the profile of our pulse. We begin with the realization that
the realistic modulation function can be represented as the ideal modulation
function convolved with some function $\vec{K}(t)$,

$f_{\mu}(t)=\int_{-\infty}^{\infty}d\tau f_{\text{Ideal}}(t)K_{\mu}(t-\tau)\
,$ (58)

where $f_{\text{Ideal}}(t)$ is the ideal modulation function and $K_{\mu}(t)$
is the convolving function for the $\mu$ axis. We do not need to know the
exact form of $K_{\mu}(t)$, since we will be computing $f_{\mu}(t)$ directly
from $\Phi(t)$. However, $K_{z}(t)$ is roughly Gaussian shaped, and $K_{x}(t)$
is roughly the derivative of a Gaussian, both with width comparable to the
pulse duration. We can use the convolution theorem to write

$\mathcal{F}\\{f_{\mu}\\}(\omega)=\mathcal{F}\\{f_{\text{Ideal}}\\}(\omega)\mathcal{F}\\{K_{\mu}\\}(\omega),$
(59)

where $\mathcal{F}\\{f\\}$ denotes the Fourier transform of $f$. Taking the
magnitude squared and summing over the components $\mu$ gives us the realistic
filter function.

$\displaystyle F(\omega)$ $\displaystyle=\sum_{\mu}F_{\mu}(\omega)$ (60)
$\displaystyle=\sum_{\mu}\frac{1}{t_{d}^{2}}\left|\mathcal{F}\\{f_{\mu}\\}(\omega)\right|^{2}$
(61) $\displaystyle=\sum_{\mu}F_{\text{Ideal}}(\omega)A_{\mu}(\omega)$ (62)
$\displaystyle=F_{\text{Ideal}}(\omega)A_{\text{Finite}}(\omega)$ (63)

where we have defined
$A_{\mu}(\omega)\equiv|\mathcal{F}\\{K_{\mu}\\}(\omega)|^{2}$ and
$A_{\text{Finite}}(\omega)\equiv\sum_{\mu}A_{\mu}(\omega)$. Figure 17 shows
examples of modulation functions and the resulting filter functions with
finite pulses.

To obtain $A_{\text{Finite}}(\omega)$, we compute $F(\omega)$ starting from
defining $\Phi(t)$ and using Equations 47 and 57, then divide by
$F_{\text{Ideal}}(\omega).$ Importantly, $A_{\text{Finite}}(\omega)$ is
independent of the interpulse delay $\tau$, so once we have found it for one
sequence (CPMG-1 for simplicity), the finite pulse filter $A(\omega)$ can be
applied in the analysis of all sequences without the need to reconstruct the
modulation function each time.

In the harmonic analysis, we include $A_{\text{Finite}}(\omega)$ as a
correction to Equation 31 as follows.

$\displaystyle\chi(N,\tau)$
$\displaystyle=\frac{t_{d}^{2}}{2}\int_{-\infty}^{\infty}d\omega
F_{\text{Ideal}}(\omega)A_{\text{Finite}}(\omega)$ (64)
$\displaystyle\begin{split}&=\frac{4N\tau}{\pi^{2}}\times\\\
&\quad\quad\Biggl{[}\left(\sum_{m=1,3,5...}^{\infty}\frac{S_{B}(m\pi/\tau)A_{\text{Finite}}(m\pi/\tau)}{m^{2}}\right)\\\
&\quad\quad+\frac{S_{P}(\omega_{s})A_{\text{Finite}}(\omega_{s})}{m_{s}^{2}}\Biggr{]}\end{split}$
(65)
$\displaystyle=\chi_{B}+\frac{4N\tau}{\pi^{2}}\frac{S_{P}(\omega_{s})A_{\text{Finite}}(\omega_{s})}{m_{s}^{2}}.$
(66)

The background term $\chi_{B}$ is acquired by the curve fitting described in
Section IV.2 and therefore needs no explicit correction. We thereby obtain the
corrected version of Equation 31.

$S_{P}(\omega_{s})=\frac{\pi^{2}m_{s}^{2}}{4N\tau
A_{\text{Finite}}(\omega_{s})}\left[\chi(N,\tau)-\chi_{B}(N,\tau)\right]$ (67)

###  F.2 Numerical reconstruction of spectra

We obtain the realistic filter function by construction of $\Phi(t)$ as
follows. First, we convolve the programmed pulse waveform profile with the
resonator impulse response to approximate the $B_{1}(t)$ pulse envelope. The
impulse response we approximate as $e^{-t|\omega_{0}|/(2Q)}$, where
$|\omega_{0}|=2.5\times 2\pi$ rad GHz is the resonant frequency and $Q$ is the
resonator quality factor. In the NV experiments on Sample A, the programmed
pulse waveform was a 64 ns Gaussian with $\sigma=13$ ns, and the resonator had
a $Q=190$; this pulse corresponds to the modulation and filter functions
labeled “Finite A” in Figure 17. In the P1 experiment on Sample B, the
programmed pulse waveform was a square pulse with length 210 ns, and the
resonator had a $Q$ of 2000 (for more sensitive detection of the weak echo);
this pulse corresponds to the modulation and filter functions labeled “Finite
B” in Figure 17. Since $\omega_{1}(t)\propto B_{1}(t)$, and we know the
calibrated delays and pulse amplitudes for implementing a $\pi$ pulse, we
normalize and center the pulse such that $\Phi(0)=0$, $\Phi(\tau/2)=\pi/2$ and
$\Phi(t_{f})=\pi$, where time $\tau/2$ is the middle of the pulse, and $t_{f}$
is some time well after the ring-down of the pulse. Applying Equation 47 then
produces the modulation function.

The $x$ and $z$-components of modulation functions with finite pulses are
shown in Figure 17(a) through (f). The lines labeled Finite A and Finite B are
for exact $\pi$ pulses with conditions corresponding to the NV and the P1
experiments respectively. Figure 17(h) compares the realistic filters against
the ideal (“Ideal”) filter. At low frequencies $\omega<1/t_{p}$, the finite
filters and the ideal filter match. At high frequencies $\omega>1/t_{p}$, the
finite filters decay rapidly. The smooth pulse transition – analogous to
sending a square wave signal through a low-pass filter – makes the sequence
insensitive to high-frequency noise. However, at frequencies near $1/t_{p}$,
the finite filter is slightly greater than the ideal filter by as much as 10%,
peaking in Figure 17(h) for the Finite A filter around 5 MHz and for the
Finite B filter around 0.7 MHz.

In our numerical reconstruction of the fundamental spectra of Figures 5(a) and
(b), we generate an array of $\chi$ values for different $\tau$ and detection
time $t_{d}$ values, similar to the experimental data presented in Figure
4(a), (b), or (c). The values of $\chi$ are calculated with Equation 56 using
the filter construction method described above and using the power spectrum
properties obtained from the power law and 13C peak fitting procedures
described in Sections IV.2 and IV.3. We then perform the contour tracing and
power law fit analysis on the numerically generated $\chi$ array in the same
manner as performed on an experimentally obtained $\chi$ array.

###  F.3 Modeling

It is convenient to use the $\delta$-function peak approximation of the filter
(Equation 17) in CPMG analyses, but how strong of an approximation is this in
the presence of finite pulses or flip angle errors? Figure 17 shows the
modulation functions, filters, and $\chi$ integrals of sequences with finite
pulses and flip angle errors. The flip angle error analysis was achieved by
altering $\Phi(t)$ such that each pulse increases $\Phi(t)$ by an increment of
$\beta\neq\pi$. The modulation and filter functions labeled “Error A” are
produced from a pulse that has the same duration as “Finite A”, but has a
greater amplitude such that there is an over-rotation of $\beta=1.05\pi$. The
functions labeled “Error B” likewise correspond to “Finite B” but have an
under-rotation of $\beta=0.95\pi$. Figure 17(f) shows how $f_{z}(t)$ for both
Error A and B reach the same level during the interpulse delay, and that this
level is incremented with each pulse. On the other hand, Figure 17(e) shows
that the functions $f_{x}(t)$ during the interpulse delays have opposite sign
for the over- and under-rotations. Figures 17(c) and (d) show $f_{x}(t)$ and
$f_{z}(t)$ over a longer time interval. With a small error, the modulation
functions for Error A and B can be seen producing an additional modulation to
the ideal $f_{z}(t)$ modulation function that is roughly sinusoidal for
$f_{x}(t)$ and cosinusoidal for $f_{z}(t)$.

$\displaystyle f_{x,\text{Err}}(t)$ $\displaystyle\approx
f_{z}(t)\sin\left(\omega_{\text{mod}}t\right)$ (68) $\displaystyle
f_{z,\text{Err}}(t)$ $\displaystyle\approx
f_{z}(t)\cos\left(\omega_{\text{mod}}t\right),$ (69)

where $\omega_{\text{mod}}=(\beta-\pi)/\tau$. In the frequency domain this is
equivalent to the ideal filter being convolved with a pair of
$\delta$-functions.

$\displaystyle\begin{split}F_{x,\text{Err}}(t)\approx&\frac{1}{t_{d}^{2}}\bigg{|}\int_{-\infty}^{\infty}\mathcal{F}\\{f_{z}\\}(\omega-\omega^{\prime})\
\times\\\
&\frac{1}{2i}\left[\delta\left(\omega-\omega_{\text{mod}}\right)-\delta\left(\omega+\omega_{\text{mod}}\right)\right]d\omega^{\prime}\bigg{|}^{2}\end{split}$
(70)
$\displaystyle\begin{split}F_{z,\text{Err}}(t)\approx&\frac{1}{t_{d}^{2}}\bigg{|}\int_{-\infty}^{\infty}\mathcal{F}\\{f_{z}\\}(\omega-\omega^{\prime})\
\times\\\
&\frac{1}{2}\left[\delta\left(\omega-\omega_{\text{mod}}\right)+\delta\left(\omega+\omega_{\text{mod}}\right)\right]d\omega^{\prime}\bigg{|}^{2}\end{split}$
(71)

Further distortions arise due to the asymmetry of the pulse, and due to the
modulation of the $x$-component of the modulation function being a sine
modulation at the modulation frequency $\omega_{\text{mod}}$. The inset of
Figure 17(h) shows how the flip angle error results in a splitting of the
peaks in the filter function. This feature appears as a general result of the
control, without any knowledge of the system-environment Hamiltonian. It has
been shown recently with Floquet analysis that over/under-rotations can be
used for polarization and readout of a single 13C via a single NV center [58,
78].

Figure 19:  (a) The fundamental ($m=1$) spectra of the 13C peak at 13 mT
obtained from NV experiments where the cycle pulse amplitudes were varied. The
numbers in the legend are $\beta/\pi$, where $\beta$ is the rotation angle of
the CPMG cycle pulse. (b) The mean (plotted value) and standard deviation
(error bars) of each peak – sampled in the shaded region of (a) – are shown as
a function of $\beta/\pi$. Figure 20:  Analyses similar to those presented in
Figures 4, 5, and 6 of the main text are applied to CPMG experiments on NV
centers at 190 mT with 10% over- and under-rotations. Image plots of $\chi$
are shown for the 10% under-rotation (a) and the 10% over-rotation (b). Black
and white dashed lines represent the power law fits, and the red vertical
stripes indicate revival $\tau$ values as in Figure 4. (c) and (d) compare the
$\alpha$ and $C$ power law ($S(\omega)=C\omega^{-\alpha}$) fit values for the
under- and over-rotation experiments ($\beta/\pi=0.9$ and $\beta/\pi=1.1$
respectively) to the experiment with no flip angle error ($\beta/\pi=1$). (e)
shows the fundamental spectra (“Fund.”) for the $\chi=3.5$ contour, the
corresponding fit power law spectra (“Pow. Law”), and the numerically
reconstructed spectra (“Recon.”). (f) shows the boxed portion of (e) in
greater detail, revealing the asymmetric peak-splitting effects caused by
under- and over-rotations. (g) and (h) show the spectra of the 13C peak
obtained with the harmonic analysis procedure for the under- and over-rotation
experiments respectively. (i), (j), and (k) compare the harmonic analysis
characterization of the 13C peak width, frequency percentage error, and
amplitude.

We also performed a numerical calculation to assess whether the flip angle
errors affect the power law spectrum measurement. Using a model spectrum of
$S(\omega)=C\omega^{-\alpha}$ with $C=5\times
10^{10}\,(\text{Hz})^{(1+\alpha)}$ and $\alpha=1$, and filter functions
calculated in the method described above, we generated $\chi$ arrays and
performed the power law fit analysis (Equation 25). The results are shown in
Figure 18. The four lines correspond to different conditions for generating
the filter function. For the line labeled “Ideal”, the modulation function
consisted of sharp transitions from instantaneous pulses. For the line labeled
“Finite”, the pulse duration was the same as conditions for Finite A, listed
above. For the line labeled “Error”, the modulation function was generated
with the same finite width as for “Finite”, but with an over-rotation of 10%.
Lastly, for the line labeled “$\delta$ Peak”, the filter was generated purely
from using the $\delta$-function approximation of Equation 17. The same
fitting procedure was applied to the four simulated data sets. The power law
fitting equation (Equation 25) is based on the $\delta$-approximation of the
filter and does not take into account the finite pulse or flip angle error
effects. The fact that the $C$ and $\alpha$ fits all significantly overlap
suggests that the finite pulse widths and flip angle errors do not drastically
affect the power law fitting results. However, the numerical construction of
the filters in this investigation was carried out only to the leading order
term of the Magnus expansion, and it is uncertain if incorporating the higher-
order terms would lead to significant deviations.

###  F.4 Experiments

We have examined the peak-splitting effect due to flip angle errors
experimentally by performing CPMG at 13 mT and varying the amplitude of the
cycle pulse. We varied the value of $\tau$ to sweep the fundamental modulation
frequency over the 13C peak. Results of the fundamental spectrum with $\chi=3$
are shown in Figure 19. We characterize the 13C peak amplitude as the mean of
the points within $\pm 5\%$ of the 13C Larmor frequency. The results suggest
that for flip angle errors of $\lesssim 5\%$, the peak is only slightly
diminished. However, this might be partially obscured by the fact that at the
13C peak, the data points in the spectrum are obtained at low echo number $N$.
Nevertheless, we see that for flip angle errors of $\gtrsim 10\%$, the peak is
strongly split, with the inversion in the middle being reduced to as little as
20% of the full peak amplitude.

To explore how the power law characterization and the 13C peak harmonic
analysis are affected by flip angle errors, we have performed two experiments
at 190 mT with over and under rotations of 10% achieved by adjusting the
amplitude of the pulse. Results are shown in Figure 20. We have performed the
same analyses as we have performed on the experiments presented in the main
text with no flip angle error. Recall that this analysis assumes the
$\delta$-function approximation of the filter in fitting the power law, and it
includes the finite pulse correction factor $A_{\text{Finite}}(\omega)$ in the
harmonic analysis but no flip angle correction. Figures 20(c) and (d) show the
power law fit parameters. At low $\chi\lesssim 1$, the sequences with flip
angle errors appear to measure a spectrum with larger $C$ and $\alpha$.
However, for the majority of $\chi$ values, there is mostly overlap between
the fits, suggesting no significant deviation of the $C$ and $\alpha$
measurement in the presence of small flip angle errors. For the 13C peak
characterization, the differences in fit parameters are slightly more
significant. Consistent with the realistic filter simulation – which is seen
in the inset of Figure 17(h) to widen the filter peak – the measurement with
the $\beta\neq\pi$ sequences indicates greater width and lower amplitude of
the 13C peak than does the measurement with no flip angle error. The
asymmetric splitting of the peak appears to lead to a shift in frequency
estimation of the 13C peak. However, the deviation in the frequency estimate
obtained from the harmonic analysis is within the uncertainty of the frequency
value based on experimental conditions as discussed in Section IV.3.

## References

  * Zhang  _et al._ [2021] T. Zhang, G. Pramanik, K. Zhang, M. Gulka, L. Wang, J. Jing, F. Xu, Z. Li, Q. Wei, P. Cigler, and Z. Chu, ACS Sensors 6, 2077 (2021).
  * Kazi  _et al._ [2021] Z. Kazi, I. M. Shelby, H. Watanabe, K. M. Itoh, V. Shutthanandan, P. A. Wiggins, and K.-M. C. Fu, Physical Review Applied 15, 054032 (2021).
  * Aslam  _et al._ [2023] N. Aslam, H. Zhou, E. K. Urbach, M. J. Turner, R. L. Walsworth, M. D. Lukin, and H. Park, Nature Reviews Physics 5, 157 (2023).
  * Barry  _et al._ [2016] J. F. Barry, M. J. Turner, J. M. Schloss, D. R. Glenn, Y. Song, M. D. Lukin, H. Park, and R. L. Walsworth, Proceedings of the National Academy of Sciences 113, 14133 (2016).
  * Nowodzinski  _et al._ [2015] A. Nowodzinski, M. Chipaux, L. Toraille, V. Jacques, J. F. Roch, and T. Debuisschert, Microelectronics Reliability Proceedings of the 26th European Symposium on Reliability of Electron Devices, Failure Physics and Analysis, 55, 1549 (2015).
  * Monge  _et al._ [2023] R. Monge, T. Delord, N. V. Proscia, Z. Shotan, H. Jayakumar, J. Henshaw, P. R. Zangara, A. Lozovoi, D. Pagliero, P. D. Esquinazi, T. An, I. Sodemann, V. M. Menon, and C. A. Meriles, Nano Letters 23, 422 (2023).
  * Ho  _et al._ [2022] K. O. Ho, Y. Shen, Y. Y. Pang, W. K. Leung, N. Zhao, and S. Yang, Functional Diamond 1, 160 (2022).
  * Lee  _et al._ [2023] J. Lee, K. Arai, H. Zhang, M. J. H. Ku, and R. L. Walsworth, npj Quantum Information 9, 1 (2023).
  * Lamichhane  _et al._ [2023] S. Lamichhane, K. A. McElveen, A. Erickson, I. Fescenko, S. Sun, R. Timalsina, Y. Guo, S.-H. Liou, R. Y. Lai, and A. Laraoui, ACS Nano 17, 8694 (2023).
  * Schwartz  _et al._ [2019] I. Schwartz, J. Rosskopf, S. Schmitt, B. Tratzmiller, Q. Chen, L. P. McGuinness, F. Jelezko, and M. B. Plenio, Scientific Reports 9, 6938 (2019).
  * Lang  _et al._ [2017] J. E. Lang, J. Casanova, Z.-Y. Wang, M. B. Plenio, and T. S. Monteiro, Physical Review Applied 7, 054009 (2017).
  * Staudacher  _et al._ [2013] T. Staudacher, F. Shi, S. Pezzagna, J. Meijer, J. Du, C. A. Meriles, F. Reinhard, and J. Wrachtrup, Science 339, 561 (2013).
  * Smits  _et al._ [2019] J. Smits, J. T. Damron, P. Kehayias, A. F. McDowell, N. Mosavian, I. Fescenko, N. Ristoff, A. Laraoui, A. Jarmola, and V. M. Acosta, Science Advances 5, eaaw7895 (2019).
  * Silani  _et al._ [2023] Y. Silani, J. Smits, I. Fescenko, M. W. Malone, A. F. McDowell, A. Jarmola, P. Kehayias, B. A. Richards, N. Mosavian, N. Ristoff, and V. M. Acosta, Science Advances 9, eadh3189 (2023).
  * Bradley  _et al._ [2019] C. E. Bradley, J. Randall, M. H. Abobeih, R. C. Berrevoets, M. J. Degen, M. A. Bakker, M. Markham, D. J. Twitchen, and T. H. Taminiau, Physical Review X 9, 031045 (2019).
  * Hensen  _et al._ [2015] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, Nature 526, 682 (2015).