\subsection*{\Large Supplemental Material}
\normalsize

\setcounter{equation}{0}
\setcounter{figure}{0}
\setcounter{table}{0}
\makeatletter
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\thetable}{S\arabic{table}}

\section{supplemental note 1: Sample preparation and details}

The diamond substrate used in this work is a 150-$\mu$m-thick plate with optically resolvable near-surface NV centers (depth $\sim 7$ nm). This same substrate is also used and described in Refs.~\cite{Myers2017b,Ariyaratne2018b,Bluvstein2019}. Optical access in the confocal setup is through the 150-$\mu$m thick diamond plate, and a metal waveguide patterned on the diamond is used to transmit microwaves.

The diamond begins as a polished, commercial Element Six electronic grade (100) diamond substrate of lateral dimensions 2 mm $\times$ 2 mm and a thickness of 0.5 mm. To reduce the thickness, the diamond is sliced and then polished from the sliced side to a thickness of 150 $\mu$m. Polishing damage is then mitigated by etching 1 $\mu$m with ArCl$_2$ plasma, and the substrate is then cleaned for 40 minutes in boiling acid H$_2$NO$_3$:H$_2$SO$_4$ 2:3. A 50-nm-thick layer of 99.99\% $^{12}$C is then grown on the diamond face by plasma-enhanced chemical vapor deposition. To form NV centers, $^{14}$N ions are then implanted with a dosage of $5.2\times 10^{10}$ ions/cm$^2$ at 4 keV and a 7$\degree$ tilt, yielding an expected N depth of 7 nm (calculated by Stopping and Range of Ions in Matter (SRIM)). Subsequently, the sample is annealed in vacuum ($< 10^{-6}$ Torr at max temperature) at 850$\degree$ C for 2.5 hours with a 40-minute temperature ramp. After annealing, the sample is cleaned in HClO$_4$:H$_2$NO$_3$:H$_2$SO$_4$ 1:1:1 for 1 hour at 230-240 $\degree$C. Materials processing techniques similar to those used here are detailed and characterized in Reference~\cite{Sangtawesin2018}.

Approximately three years span between the initial sample preparation and the measurements performed in this work, during which time several rounds of surface cleaning are performed. The measurements performed in this work are performed several months after a standard surface preparation protocol: the diamond is acid cleaned with Nanostrip (a Piranha analog) at 80$^{\circ}$C for 12 minutes, and then the diamond is oxygen annealed at 400$^{\circ}$C for 4 hours. Nanostrip is chosen because it preserves the metal waveguide on the diamond.

The NV centers' depths are measured by proton NMR and range between $\sim$4 and 17 nm; further details are given in the supplement of Ref.~\cite{Bluvstein2019}. To increase the photon collection efficiency from the NV centers, tapered nanopillars with a diameter of 400-nm are patterned by e-beam lithography and etched by O$_2$ plasma to a height of 500 nm. Under conventional spin-dependent photoluminescence readout, we measure spin contrasts of $\approx$ 35\% and a saturation fluorescence of $\approx$ 500 k counts/s. Pillars containing single NV centers are identified by second-order correlation measurements and then confirmed by measuring photon statistics from the negative and neutral NV charge states (see supplement of Ref.~\cite{Bluvstein2019}). By measuring the fraction of pillars hosting a single NV center and assuming the number of NVs per pillar follows a Poisson distribution, we estimate a conversion efficiency of implanted N to NV of $\approx 5\%$, consistent with previous reports \cite{Toyli2010}.

On day-to-month timescales, on some NVs we observe order-unity variations in the NV coherence time as well as the NV coupling strength to the surface spins. Accordingly, for Figs.~2 and 3(c) of the main text we interleave drive and no-drive sequences so that we are insensitive to sub-kHz changes in the environment. We note that for the DEER data in Figs.~1(c) and 1(d) of the main text, the NV $T_2$ was slightly lower during the measurement in Fig.~1(c), and so the quantitative values of NV coherence in Fig.~1 should not be directly compared.

\section{Supplemental Note 2: Surface spin density estimate}
Here we describe an estimate of the surface spin density by using the measurements presented in this work. One may use the DEER coupling strength to estimate the surface spin density; however, depending on the models and assumptions applied (\textit{e.g.} about correlation times of the noise), similar DEER coupling strengths can result in estimated densities that vary by orders of magnitude. Accordingly, in this section we present statistical arguments to set a lower bound on the surface spin density of $\approx 0.01 / \text{nm}^2$, and we use knowledge of the surface spin linewidth to set an upper bound of $\approx 0.1 / \text{nm}^2$. The bounded range of $0.01 - 0.1 /\text{nm}^2$ is consistent with the densities extracted by various different methods in Refs.~\cite{Sushkov2014,Myers2014,DeWit2018}.

\subsection*{2.1: Lower bound on surface spin density: statistical argument}

We can set a lower bound on the surface spin density by using statistical arguments, based on the sample of NVs we measure. A low spin density will lead to large variations in the surface spin signal that the individual NV centers feel: if we assume the extreme scenario that all observed variation between different NVs is due to a variation in the number of surface spins in the NV sensing volume, then we can estimate a density of surface spins which would be a lower bound. 

Figure~3(c) of the main text plots $\Gamma_\text{decoupled} = 1/T_{2} - 1/T_{2,\text{drive}}$ for eight individual centers of different depths. Six of the centers are between depths of approximately 7.5 to 12.5 nm; for these six centers, the mean decoupling rate mean($\Gamma_\text{decoupled}$) = 7.75 ms$^{-1}$ and the standard deviation std($\Gamma_\text{decoupled}$) = 2.9 ms$^{-1}$. From the plot, it appears a majority of the variation for these six centers originates from their varied depth, but for this estimate we assume all variation in $\Gamma_\text{decoupled}$ is due to the variation in the number of surface spins in the NV sensing volume.

For this estimate we assume that $\Gamma_\text{decoupled}$ is proportional to the mean-squared magnetic field $B^2_\text{rms}$ from the surface spins. $B^2_\text{rms}$ is proportional to $N$, where $N$ is the number of surface spins in the sensing volume, which here we assume follows a Poisson distribution. Then std($\Gamma_\text{decoupled}$)/mean($\Gamma_\text{decoupled}$) = std($B^2_\text{rms}$)/mean($B^2_\text{rms}$) = $\sqrt{N} /N = 1 /\sqrt{N}$. We observe std($\Gamma_\text{decoupled}$)/mean($\Gamma_\text{decoupled}$) $\approx 0.37$, and accordingly $N \approx 1/(0.37)^2 \approx 7$. Approximately, the relevant sensing area on the surface is a circle with area $\sim \pi r^2$, where $r \approx 10$ nm is the average depth of the six NV centers we consider here. These estimates correspond to a surface spin density of $7 / \pi (10 ~\text{nm})^2 \approx 0.02 / \text{nm}^2$, and to be conservative we estimate a lower bound of $\approx 0.01 / \text{nm}^2$.

This argument becomes invalid if the surface spins are not stationary over the course of the entire measurement of $\Gamma_\text{decoupled}$ on an NV, which lasts minutes to hours. If the spins are not stationary, then an arbitrarily low spin density would not necessarily result in any variation between different NV centers. However, the spins are likely stationary on the surface, as Sushkov \textit{et al}. \cite{Sushkov2014} image the locations of these surface spins by varying the magnetic field direction.

\subsection*{2.2: Upper bound on surface spin density: dipolar-limited surface spin linewidth}

We can set an upper bound on the surface spin density by noting the typical surface spin linewidths we observe, because dipolar coupling between the surface spins will broaden the line. In this work, we observe a surface spin $1/T_{2,\text{Rabi}} \approx 5$ MHz, roughly corresponding to the half width at half maximum (HWHM) of the resonance linewidth. Note that HWHMs of approximately 10 MHz are observed in Figs.~1(c) and 3(b) of the main text, but both of these HWHMs are Fourier broadened by approximately 5 MHz. If we assume that the measured 5 MHz decay constant is exclusively due to dipolar interactions between the surface spins, then we can set an upper bound on the surface spin density.

Roughly, the linewidth will be of the order of the dipolar coupling strength $A$. Consider, for example, the case of two spins with dipolar interaction strength $A$ and Zeeman interaction energy $\gamma B$, where $\gamma B \gg A$. The transition frequency of $|\downarrow \downarrow \rangle \rightarrow|\uparrow \downarrow \rangle$ will be $\gamma B - A$, and the transition frequency of $|\downarrow \uparrow \rangle \rightarrow|\uparrow \uparrow \rangle$ will be $\gamma B + A$. Having large numbers of dipolar-coupled spins thus creates a broadening with a HWHM of order $A$. Slichter calculates the second moment (\textit{i.e.} the variance) of the linewidth given $N$ spins limited by dipolar broadening \cite{Slichter1990}, and this variance is given by

\begin{equation}
\langle\Delta f^2\rangle = \left(\frac{1}{2\pi}\right)^2 \left(\frac{\mu_0}{4\pi}\right)^2\frac{3}{4} \gamma^4 \hbar^2 S(S+1) \frac{1}{N} \sum_{j,k} \frac{(1-3\cos^2{\theta_{jk}})^2}{r^6_{jk}},
\label{slichter}
\end{equation}

where $\gamma$ is the electron gyromagnetic ratio, $S$ is the spin (1/2 for the spin-1/2 surface spins), $N$ is the number of spins, and the sum adds the contribution from all spins, which are at a relative separation $r_{jk}$ and point at a relative angle $\theta_{jk}$.

As in Slichter \cite{Slichter1990}, we now consider the case where all spins are located in equivalent positions, such that the sum is independent of $j$. This then gives $N$ equivalent sums, canceling with the $1/N$. We arrive at

\begin{equation}
\Delta f_\text{rms} \equiv \sqrt{\langle\Delta f^2\rangle} \approx \text{39 MHz nm$^3$} \sqrt{\sum_{k} \frac{(1-3\cos^2{\theta_{k}})^2}{r^6_{k}}}.
\label{slichter2}
\end{equation}

We approximate the sum for our surface spins by considering a two-dimensional square array of spins. We use the nearest neighbor approximation (\textit{i.e.}, we sum the contribution from the four nearest neighbors to the central spin), and we consider a spin z-axis which is tilted 54.7$^{\circ}$ from the axis normal to the array (the direction of the applied magnetic field for our (100) diamond). These considerations then yield an upper bound for the surface spin density given by

\begin{equation}
\sigma \approx \frac{1}{\text{nm}^2} \left(\frac{\Delta f_\text{rms}}{125~ \text{MHz}}\right)^{2/3}.
\label{sigmahigh}
\end{equation}

If we take $\Delta f_\text{rms}$ to be our experimentally measured HWHM of 5 MHz, then Eq.~\ref{sigmahigh} gives an upper bound to the surface spin density of $\approx 0.1 / \text{nm}^2$.We emphasize that, by our calculations here, a surface spin density of $1/\text{nm}^2$ would give a FWHM of 250 MHz, far inconsistent with the observed linewidths in this work and in previous reports \cite{Mamin2012,Grinolds2014,Sushkov2014}.

This argument assumes that the spins lie perfectly in a 2D plane, without any height variation: however, rapid nm changes in height over nm lengthscales would lead to lower average dipolar coupling strengths than if the spins are all in the same plane. We also note that one should be careful in the interpretation of linewidth. A Lorentzian lineshape, for example, has infinite standard devation but a well-defined HWHM. In the case of the two-dimensional square array of spins, the HWHM is approximately the same value as the standard deviation and so we can use the method of moments safely.

\section{supplemental note 3: Radio-frequency control}

\subsection*{3.1: Radio-frequency equipment}
Three separate radio-frequency (RF) signal generators are used for controlling the electronic spin states of the NV spin qutrit, formed by $|m_s=0\rangle_{\text{NV}}$ and $|m_s=\pm1\rangle_{\text{NV}}$, as well as the surface spin qubits, formed by $|\uparrow \rangle_{\text{SS}}$ and $|\downarrow\rangle_{\text{SS}}$. A single microwave waveguide patterned on the diamond is used to deliver all RF signals. All three RF generators (SRS SG384) are gated by their own switch (Mini Circuits ZASWA-2-50DR+) and are then combined before an RF amplifier (Mini Circuits ZHL-16W-43+) which feeds into the patterned waveguide. Independent IQ modulation of each SG384 generator is used as an extra layer of isolation. Measurements are done in an applied magnetic field of 100s of Gauss, and so the surface spin RF frequency is $\sim$ GHz and passes through the amplifier.

\subsection*{3.2: Preparation, control, and readout of double-quantum basis}

The $|m_s=-1\rangle_{\text{NV}}$ and $|m_s=+1\rangle_{\text{NV}}$ states are not coupled by microwaves, so to manipulate superpositions in the DQ basis we employ microwave control of the SQ bases $\{0,-1\}_{\text{NV}}$ and $\{0,+1\}_{\text{NV}}$. We use subsequent single-frequency pulses to prepare, control, and readout DQ superpositions. For example, to prepare a DQ superposition, we apply an $|m_s=0\rightarrow-1\rangle_{\text{NV}} ~\pi/2$ pulse to prepare a SQ superposition, and then apply an $|m_s=0\rightarrow+1\rangle_{\text{NV}} ~\pi$ pulse. The DQ Hahn echo sequence we use in Fig.~2 of the main text is depicted in Fig.~\ref{DQecho}, which includes use of a differential measurement (see below). Dual-frequency pulses have been shown to also be an effective way to control the DQ basis \cite{Mamin2014}, at the cost of experimental complexity; however, dual-frequency pulses would likely be superior for use in measurements with long trains of RF pulses, such as CPMG-type measurements.

\subsection*{3.3: Differential measurement}
For all measurements plotted in the main text and supplement of this work, we use a common-mode rejection technique known as a differential measurement, used throughout the NV literature and described in Refs.~\cite{Myers2017b,Ariyaratne2018b,Bluvstein2019}. The differential measurement helps remove technical noise, can alleviate some effects of NV charge state conversion in the dark \cite{Bluvstein2019}, and also simplifies analysis by normalizing the signal so that the signal goes to $0$ as $\tau \rightarrow \infty$ in an echo sequence.

In the differential measurement scheme, at the end of a pulse sequence we readout the NV photoluminescence, denoted PL in Fig.~\ref{DQecho}. We then repeat the pulse sequence, but immediately before photoluminescence readout we apply a resonant microwave $\pi$ pulse to swap the $|m_s=0\rangle_{\text{NV}}$ and $|m_s=-1\rangle_{\text{NV}}$ populations and measure PL$_\text{swap}$ (Fig.~\ref{DQecho}). Subtracting PL and PL$_\text{swap}$ thus removes background PL signals and yields the population difference between $|m_s=0\rangle_{\text{NV}}$ and $|m_s=-1\rangle_{\text{NV}}$.
Explicitly, the overall signal we measure is given by 

\begin{equation}
\text{PL} - \text{PL}_\text{swap} = C \left(\rho_0 - \rho_{-1}\right),
\label{Diff}
\end{equation}

where $C$ is a constant describing the spin-dependent PL contrast, and $\rho_0$ and $\rho_{-1}$ are the $|m_s=0\rangle_{\text{NV}}$ and $|m_s=-1\rangle_{\text{NV}}$ populations after the projective $\pi/2$ pulse. Throughout this work, we define NV coherence as $\left(\text{PL} - \text{PL}_\text{swap}\right)/C$.

\begin{figure}
\includegraphics[width=152mm]{SuppFig1_DQsequence.pdf}
\caption{Pulse sequence used for performing Hahn echo measurements in the double-quantum basis, which includes use of a differential measurement technique by the additional $\pi$ pulse at the end of the $\text{PL}_\text{swap}$ measurement. The NV coherence is proportional to the differential measurement signal $\text{PL} - \text{PL}_\text{swap}$.
}
\label{DQecho}
\end{figure}

\section{supplemental note 4: The microwave analog of the AC Stark effect}

In this section we show that the drive field used for addressing the surface spins can cause an appreciable shift in the NV spin transition frequencies, even if far detuned from the NV spin transition. We identify this effect as a microwave analog of the AC stark effect, common to optical transitions in atomic physics. We observe this effect in our DEER measurement sequences as a slow oscillation in the NV phase as a function of the surface spin pulse duration, and we show in this section that this effect can be readily removed with the appropriate fitting function. We note that we also observe this effect as a frequency shift in Ramsey experiments (not plotted in this work). We emphasize that these results indicate that amplitude fluctuations in the surface spin drive field, even though far detuned from the NV transition, would lead to NV dephasing. Similarly, spatial inhomogeneity in the drive field would lead to dephasing in NV ensemble measurements.

\subsection*{4.1: Theoretical description of AC Stark effect on NV spin triplet}

The AC Stark effect is common to optical transitions in Atomic physics, and is described in detail in \textit{Theoretical Atomic Physics} by Friedrich \cite{Friedrich2017}. Here we adapt the calculation to the magnetic transitions of the NV ground state spin triplet. To be explicit, the NV interaction with the Rabi field can be understood by expanding the Zeeman interaction:

\begin{equation}
H^{\prime} = \gamma\boldsymbol{S} \cdot \boldsymbol{B} \cos(\omega t) = \left[\gamma B_z S_z + \left(\frac{1}{2} \gamma \left(B_x - i B_y\right) S_+ + \frac{1}{2} \gamma \left(B_x + i B_y\right) S_-\right)\right]\cos(\omega t),
\label{Stark1}
\end{equation}

where $\gamma$ is the electron gyromagnetic ratio, $\boldsymbol{B}$ is the magnetic field amplitude vector, $\cos(\omega t)$ describes the oscillation of the magnetic field, $\boldsymbol{S}$ is the spin-1 operator, and $S_{\pm} = S_x \pm i S_y$ are the spin raising and lowering operators. Analagous to the calculation in Friedrich \cite{Friedrich2017}, the energy shifts (up to second-order) to the NV spin triplet states due to an oscillating magnetic field are thus given by

\begin{equation}
\begin{split}
\Delta E_n &= \frac{1}{4} \left[ \sum_{E_m + \hbar \omega \neq E_n} \frac{|\langle\psi_n|\gamma\boldsymbol{S} \cdot \boldsymbol{B}|\psi_m\rangle|^2}{E_n - E_m - \hbar \omega} \right. \\ &+ \left. \sum_{E_m - \hbar \omega \neq E_n} \frac{|\langle\psi_n|\gamma\boldsymbol{S} \cdot \boldsymbol{B}|\psi_m\rangle|^2}{E_n - E_m + \hbar \omega} \vphantom{\sum_{E_m + \hbar \omega \neq E_n}} \right],
\label{Stark2}
\end{split}
\end{equation}

where $\psi_n$ are the NV $m_s = 0,+1,-1$ spin states, and $E_n$ is the energy of the $n^{\text{th}}$ level. Expanding the Zeeman interaction we arrive at

\begin{equation}
\begin{split}
\Delta E_n &= \frac{1}{4} \left[ \sum_{E_m + \hbar \omega \neq E_n} \frac{|\langle\psi_n|\left(\frac{1}{2} \gamma \left(B_x - i B_y\right) S_+ + \frac{1}{2} \gamma \left(B_x + i B_y\right) S_-\right)|\psi_m\rangle|^2}{E_n - E_m - \hbar \omega} \right. \\ &+ \left. \sum_{E_m - \hbar \omega \neq E_n} \frac{|\langle\psi_n|\left(\frac{1}{2} \gamma \left(B_x - i B_y\right) S_+ + \frac{1}{2} \gamma \left(B_x + i B_y\right) S_-\right)|\psi_m\rangle|^2}{E_n - E_m + \hbar \omega} \vphantom{\sum_{E_m + \hbar \omega \neq E_n}} \right],
\label{Stark3}
\end{split}
\end{equation}

where the $S_z$ term is cancelled out because the inner product will be $0$ unless $n=m$ for $S_z$, but if $n=m$ then the two terms cancel. In our experiment, we are interested in the relative energy shift induced between two spin levels. Here we will focus in particular on the relative energy shift between $m_s = 0$ and $m_s = -1$, which we use as our single-quantum qubit. Calculation with Eq.~\ref{Stark3} results in

\begin{equation}
\hbar\omega_{\text{stark}} \equiv
\Delta E_{-1} - \Delta E_{0} = \hbar \Omega^2_\text{NV} \left(\frac{1}{2\Delta_{-1}} + \frac{1}{2\Delta_{-1} + 4 \omega} + \frac{1}{4\Delta_{+1}} + \frac{1}{4\Delta_{+1} + 8 \omega}
\right)
\label{Eshift}
\end{equation}

where $\Omega_\text{NV} = (1/\sqrt{2})\gamma \sqrt{B^2_x + B^2_y}$ describes the NV coupling to the Rabi field, and $\Delta_{\pm1} = \omega_{\pm 1} - \omega$ are the detunings of the drive frequency $\omega$ from the $|m_s = 0 \rightarrow \pm 1 \rangle_{\text{NV}}$ transition frequencies $\omega_{\pm 1}$. Note that for $\Delta_{-1} \ll \omega, \Delta_{+1}$, Eq.~\ref{Eshift} becomes approximately $\hbar \Omega^2_\text{NV} / 2\Delta_{-1}$, which is commonly quoted as the AC Stark shift.

In this work we measure $\Omega_\text{SS}$ as a proxy for $\Omega_\text{NV}$, where $\Omega_\text{SS}$ is the surface spin coupling to the Rabi field. We emphasize that $\Omega_\text{NV} = \sqrt{2} ~\Omega_\text{SS}$, assuming the magnetic field is aligned to the NV axis. This $\sqrt{2}$ factor is due to the nature of the NV being a spin-1 particle, compared to the surface spins which are spin-1/2. Explicitly, the extra factor of $\sqrt{2}$ arises because $S_+ |S=\frac{1}{2}, m_s = -\frac{1}{2} \rangle = \hbar |S=\frac{1}{2}, m_s = +\frac{1}{2} \rangle$, whereas $S_+ |S= 1, m_s = -1 \rangle = \sqrt{2}~ \hbar |S=1, m_s = 0 \rangle$.

In summary, the theoretical result derived in this section shows that a relative energy shift is induced between the NV spin levels due to an off-resonant RF drive. Under the experimental parameters associated with the drive field used to address the surface electron spins in this work, this energy shift can be appreciable. For example, even at a large detuning of $\Delta_{-1}/2\pi =$ 1000 MHz, a drive field of strength $\Omega_\text{NV}/2\pi = 20$ MHz yields a Stark shift of approximately $(1/2\pi)~ \Omega^2_\text{NV} / 2\Delta_{-1} = 0.2$ MHz, which is easily observable with an NV $T_2$ of 10-100 microseconds. Further, these results indicate that the stability of the surface spin drive field is important for maintaining NV coherence. For example, for a 0.2 MHz shift, RF intensity variations of $\sim$ 10\% would lead to $\sim$ 20 kHz variations in the transition frequency and limit the NV coherence to $\sim$ 50 $\mu$s.

\subsection*{4.2: Experimental observation of AC Stark effect and its effect on NV measurements}

In this section, we focus on the AC Stark effects on DEER in particular, and we show that the AC Stark effect can be removed by an appropriate fitting function. Figure~\ref{starkfig}(a) shows the DEER measurement pulse sequence. We start the surface spin pulse 100 ns after the NV $\pi$ pulse ends due to technical reasons, \textit{e.g.} we have observed that sending both pulses at the same time through the RF circuit can cause artefacts which reduce pulse fidelity.

\begin{figure}
\includegraphics[width=159mm]{SuppFig2_ACStark.pdf}
\caption{Microwave analog of the AC Stark effect and its effect on NV measurements. (a) DEER pulse sequence used in this work. The 100-ns offset is used for technical reasons and does not affect measurements of surface spin properties. (b) Rabi oscillations of surface spins measured by the DEER sequence in (a). The data has a  downward trend caused by the AC Stark effect, which we fit by including an additional factor of $(1 - (\omega_\text{osc} t)^2/2)$, where $\omega_\text{osc}$ is a free fit parameter corresponding to the $\cos(\omega_\text{osc} t)$ signal observed in (d). (c) Data from (b) after dividing by the fitted $(1 - (\omega_\text{osc} t)^2/2)$. In both (b) and (c), the fitted $\Omega_\text{SS} / 2\pi = 13.7(1)$ MHz. (d) Same measurement as (b), with the same experimental parameters but for longer pulse lengths $t$. A coherent phase oscillation that persists far beyond the surface spin coherence time is observed and is fit to $A \cos(\omega_\text{osc} t)$. The fitted oscillation rate $\omega_\text{osc}$ agrees with the predicted oscillation frequency $\omega_\text{stark}$, which is calculated independent of the data in (d). To calculate $\omega_{\text{stark}}$ with Eq.~\ref{Eshift}, we use the experimental values: the applied drive field is $\omega/2\pi$ = 790 MHz, the state detunings are $\Delta_{-1} / 2\pi = 1292$ MHz and $\Delta_{+1} / 2\pi = 2866$ MHz, and the NV coupling to the Rabi field is $\Omega_\text{NV} / 2\pi = \sqrt{2}~ \Omega_\text{SS} /2\pi= \sqrt{2} ~13.7(1)$ MHz. (e) Measured $\omega_\text{osc}/2\pi$ for varied detunings $\Delta_{-1}/2\pi$, taken at $\Omega_\text{NV}/2\pi = 5.5$ MHz. The data are consistent with the plotted curve, which is calculated independent of the plotted data (no fit parameters).
}
\label{starkfig}
\end{figure}

In Fig.~\ref{starkfig}(b) we plot the DEER signal $\text{PL} - \text{PL}_\text{swap}$ (which is proportional to NV coherence, see Supplemental Material Note 3.3) when applying resonant microwaves to the surface spins at frequency $f_\text{SS}$. Rabi oscillations of the surface spins are observed, but there is also a downward trend of the data. At longer pulse lengths $t$, we observe in Fig.~\ref{starkfig}(d) that the downward trend observed in Fig.~\ref{starkfig}(b) is part of a $\cos(\omega_\text{osc} t)$ oscillation in the NV signal due to a coherent phase accumulation at rate $\omega_\text{osc}$. Thus, the downward trend in Fig.~\ref{starkfig}(b) can be removed by dividing by $\cos(\omega_\text{osc} t)$. In practice, we typically fit the DEER signal with the Taylor expanded $(1 - (\omega_\text{osc} t)^2/2)$; dividing this expression from the data yields the plot in Fig.~\ref{starkfig}(c).

The coherent phase accumulation observed in Fig.~\ref{starkfig}(d) is due to the AC Stark effect. We fit the data in Fig.~\ref{starkfig}(d) to $A\cos(\omega_\text{osc} t)$, where $A$ is in k counts/s and the fitted $\omega_\text{osc}/2\pi = 0.260(1)$ MHz is the oscillation frequency. Using Eq.~\ref{Eshift}, we calculate that the AC stark effect should yield an oscillation frequency $\omega_{\text{stark}}/2\pi = 0.264(4)$ MHz, consistent with the fitted $\omega_\text{osc}/2\pi = 0.260(1)$ MHz. To calculate $\omega_{\text{stark}}$ with Eq.~\ref{Eshift}, we use the experimental values: the applied drive field is at frequency $\omega/2\pi$ = 790 MHz, the state detunings are $\Delta_{-1} / 2\pi = 1292$ MHz and $\Delta_{+1} / 2\pi = 2866$ MHz, and the NV coupling to the Rabi field is $\Omega_\text{NV} / 2\pi = \sqrt{2}~ \Omega_\text{SS} /2\pi= \sqrt{2} ~13.7(1)$ MHz, where $\Omega_\text{SS}$ is fit in Figs.~\ref{starkfig}(b) and \ref{starkfig}(c). Here we use the surface spin Rabi oscillations in order to measure the NV coupling to the Rabi field, possible because the magnetic field is aligned to within a few degrees of the NV axis. Again, we emphasize the $\sqrt{2}$ factor that arises from careful consideration of the difference between the spin-1 NV and the spin-1/2 surface spins. The agreement here between theory and experiment is thus strong experimental evidence that these surface spins are spin-1/2 particles.

In Fig.~\ref{starkfig}(e) we plot measured values of $\omega_\text{osc}/2\pi$ as a function of drive detuning from the $|m_s=0 \rightarrow -1 \rangle_{\text{NV}}$ transition (here the microwaves are not resonant with the surface spins). At these comparatively small detunings, we find good agreement of $\omega_\text{osc}$ with $\Omega^2_\text{NV} / 2\Delta_{-1}$, where here we measure $\Omega_\text{NV}$ by measuring NV Rabi oscillations on resonance. Figure~\ref{starkfig}(e) plots the average of the measured oscillation rates at positive detuning $\Delta_{-1}$ and negative detuning $\Delta_{-1}$, which renders the measurement first-order insensitive to contributions from the $|m_s=0 \rightarrow +1 \rangle_{\text{NV}}$ transition and changes in $\Omega_\text{NV}$ as a function of frequency (due to nonlinearities in the RF circuit).

\section{supplemental note 5: contributions to NV dephasing}

\subsection*{5.1: Quantitative comparison of SQ and DQ coherence times}

We show here that the ratio of the single-quantum (SQ) and double-quantum (DQ) coherence times observed in Fig.~2(c) of the main text indicates there is significant common-mode noise decoupled by use of the DQ basis \cite{Kim2015}. As discussed in the main text, under a large applied magnetic field $\boldsymbol B = B_z \hat{z}$, the SQ transition frequencies are given by

\begin{equation}
    f_{0\rightarrow\pm1} \approx D + d_{\|} \Pi_{\|} / h \pm \left(\frac{\gamma}{2\pi} B_z + \frac{1}{2} \frac{\left(d_{\perp} \Pi_{\perp}/h\right)^2}{(\gamma/2\pi) B_z}\right),
    \label{SQFrequenciesSupp}
\end{equation}

and the DQ transition frequency is given by

\begin{equation}
    f_{-1\rightarrow+1} \approx 2 \left(\frac{\gamma}{2\pi} B_z + \frac{1}{2} \frac{\left(d_{\perp} \Pi_{\perp}/h\right)^2}{(\gamma/2\pi) B_z}\right).
    \label{DQFrequenciesSupp}
\end{equation}

where $h$ is Planck's constant, $D = 2.87$ GHz is the crystal-field splitting, $B_z$ is the magnetic field magnitude along the NV axis, $\gamma/2\pi = 2.8$ MHz/G is the NV gyromagnetic ratio, $d_{\|}/h = 0.35$ Hz$\cdot$cm/V and $d_{\perp}/h = 17$ Hz$\cdot$cm/V are the components of the NV's electric dipole parallel and perpendicular to the $z$-axis, and $\Pi_{\|}$ and $\Pi_{\perp}$ are the parallel and perpendicular components of the effective electric field, where $\boldsymbol{\Pi} = (\boldsymbol{E} + \boldsymbol{\sigma})$ has both electric field $\boldsymbol{E}$ and appropriately scaled strain $\boldsymbol{\sigma}$ terms. For simplicity, we define $\Tilde{E}$ and $\Tilde{B}$ such that

\begin{equation}
\begin{aligned}
    f_{0\rightarrow\pm1} = \Tilde{E} \pm \Tilde{B}\\ f_{-1\rightarrow+1} = 2 \Tilde{B}.
    \end{aligned}
\label{Freqs}
\end{equation}

Fluctuations in $\Tilde{E}$ and $\Tilde{B}$ lead to fluctuations in accumulated phase between measurements, and hence cause dephasing. As calculated in Ref.~\cite{Kim2015}, for the spin echo this results in coherences $C(T)$ as a function of the phase accumulation time $T$ given by

\begin{equation}
    \begin{aligned}    
    C_{\text{SQ}}(T) &= \exp\left[-\frac{\langle(\delta \phi_{\Tilde{E}})^2\rangle + \langle(\delta \phi_{\Tilde{B}})^2\rangle]}{2}\right] \approx \exp\left[-\left(\frac{T}{T_{2,\text{SQ}}}\right)^{n_\text{SQ}}\right] \\ 
    &C_{\text{DQ}}(T) = \exp\left[-2 \langle( \delta \phi_{\Tilde{B}})^2\rangle \right] \approx \exp\left[-\left(\frac{T}{T_{2,\text{DQ}}}\right)^{n_\text{DQ}}\right],
    \end{aligned}
\label{Coherences}
\end{equation}

where $\langle(\delta \phi_{\Tilde{E}})^2\rangle$ and $\langle(\delta \phi_{\Tilde{B}})^2\rangle$ are the variances in accumulated phase induced by $\Tilde{E}$ and $\Tilde{B}$, $T_{2,\text{SQ}}$ and $T_{2,\text{DQ}}$ are the coherence times of the single- and double-quantum bases, and $n_\text{SQ}$ and $n_\text{DQ}$ are stretch factors. Note that this calculation assumes that $\delta \phi_{\Tilde{E}}$ and $\delta \phi_{\Tilde{B}}$ are uncorrelated Gaussian random variables \cite{Kim2015}. If we assume that the NV is only dephased by magnetic field fluctuations, then $\langle(\delta \phi_{\Tilde{E}})^2\rangle = 0$, and if we assume that $n_\text{SQ} \approx n_\text{DQ} \equiv n$, then we can write

\begin{equation}
\left(\frac{T_{2,\text{DQ}}}{T_{2,\text{SQ}}}\right)^n \approx 0.25, \hspace{10pt}\text{for}~ \langle(\delta\phi_{\Tilde{E}})^2\rangle = 0.
\label{ratio}
\end{equation}

And if $\langle(\delta\phi_{\Tilde{E}})^2\rangle \neq 0$, then $\left(T_{2,\text{DQ}}/T_{2,\text{SQ}}\right)^n > 0.25$. For the data in Fig.~2(c) of the main text, we observe $(T_{2,\text{DQ}} / T_{2,\text{SQ}})^n = 0.48(4)$, where the fitted $n \approx 1.6$. This measured value indicates there is substantial common-mode noise eliminated by use of the double-quantum basis.

\subsection*{5.2: CPMG $T_2$ and Ramsey $T^*_2$}

In addition to Hahn echo measurements, we also perform CPMG-8 and Ramsey sequences with and without surface spin driving. By driving the surface spins we observe a CPMG-8 $T_{2,\text{SQ}}$ increase for three out of three NV centers measured, where the fractional increase on each center is similar to the Hahn echo $T_{2,\text{SQ}}$ increase. Future work can carefully map out the decoupled noise spectrum to evaluate the efficacy of surface spin driving as a function of the number of $\pi$ pulses, in both the SQ and DQ bases.

We do not observe a change in the Ramsey $T^*_{2,\text{SQ}}$ with surface spin driving, which is consistent with the fact that the $1/T^*_{2,\text{SQ}}$ values we measure (100s of kHz) are much greater than the DEER coupling rates we measure (10s of kHz). We attribute the other contributions to $1/T^*_{2,\text{SQ}}$ to low-frequency noise that is independent of the surface spins, such as electric field noise, magnetic noise from P1 centers, and drift of the applied magnetic field ($B_0 \sim 300$ G).

\end{document}

