Parallel detection and spatial mapping of large nuclear spin clusters

Nuclear magnetic resonance imaging (MRI) at the atomic scale offers exciting prospects for determining the structure and function of individual molecules and proteins. Quantum defects in diamond have recently emerged as a promising platform towards reaching this goal, and allowed for the detection and localization of single nuclear spins under ambient conditions. Here, we present an efficient strategy for extending imaging to large nuclear spin clusters, fulfilling an important requirement towards a single-molecule MRI technique. Our method combines the concepts of weak quantum measurements, phase encoding and simulated annealing to detect three-dimensional positions from many nuclei in parallel. Detection is spatially selective, allowing us to probe nuclei at a chosen target radius while avoiding interference from strongly-coupled proximal nuclei. We demonstrate our strategy by imaging clusters containing more than 20 carbon-13 nuclear spins within a radius of 2.4 nm from single, near-surface nitrogen–vacancy centers at room temperature. The radius extrapolates to 5–6 nm for 1H. Beside taking an important step in nanoscale MRI, our experiment also provides an efficient tool for the characterization of large nuclear spin registers in the context of quantum simulators and quantum network nodes.


Single-spin free induction decay under periodic weak measurements
Weak measurements probe the transverse nuclear polarization (specically, the Î x component) linearly with the measurement strength β [1]. The measured free-induction decay (FID) signal for nuclear spin i as a function of time, that is, the measured transition probability for the spin sensor, has the form: x i (t) = 1 2 p 0 sin (β) cos (ω i t) e −Γ i t e −(t β /T 2,e) a (S1) where ω i = γn t t 0 B (t ) dt is the average precession frequency of nuclear spin i. γ n is the nuclear spin gyromagnetic ratio. For completeness, we have included the NV electronic spin coherence time T 2,e , where a is an exponent which describes the dephasing rate of the spin sensor (exponential or Gaussian). Because in our experiments t β T 2,e , sensor dephasing is not important for our study. p 0 describes the initial nuclear polarization.
Weak measurements simultaneously address all nuclear spins within the detector bandwidth of approximately t −1 β [1]. Because the measurements are weak, β ≈ a ⊥ t β 1, we can neglect nuclear-nuclear interactions mediated by the sensor spin. For n nuclear spins, the total FID signal is thus the superposition of all individual contributions, leading to Eq. (2) in the main text.

Amplitude and decay rate
The evolution of a quantum stateρ k describing a single nuclear spin subject to a series of sequential weak measurements of duration t β at a repetition period t s , is approximately given by [1]:ρ k (t) ≈Î e + Î x cos (ωt) +Î y sin (ωt) e −Γ β t + O(β 2 ) where t = kt s (k = 1 . . . K), ω is the nuclear Larmor frequency, β = a ⊥ t β /π is the measurement strength and Γ β = β 2 /(4t s ) is the measurement-induced dephasing due to quantum back-action.
Weak measurements map the expectation value Î x of the nuclear-spin state, proportionally to the measurement strength β onto the optically-readableŜ z sensor state [1]. The measured signal thus has an amplitude where p 0 is the initial polarization given by Î x (t = 0). This is Eq. (5a) in the main text.
At the same time, a series of weak measurements leads to an exponential decay of the nuclear coherencesÎ x andÎ y , ultimately leading to a fully mixed state. Detailed analysis of the spin-state projection shows [1] that the loss of coherence is, on average, β 2 /4 per weak measurement, leading to an exponential decay rate Γ β = β 2 /(4t s ). Further, owing to our o-resonant optical readout and re-initialization scheme, the NV center undergoes stochastic spin-state, electronic-state and potentially charge-state transitions during optical illumination [2,3]. This leads to stochastic jumps in the hyperne interaction experienced by the nuclear spin, leading to a readout-induced dephasing rate Γ optical = (a||t ) 2 2ts . Here, t is an experimentally-determined parameter that roughly equals the time that the NV center spin is in an undened m S state during laser readout. Finally, nuclear spins are subject to intrinsic dephasing Γ 0 = (T * 2,n ) −1 given by the intrinsic dephasing time T * 2,n . The total decay rate is then the sum of all contributions: This is Eq. (5b) in the main text.

Dynamic nuclear polarization (DNP)
DNP methods refer to the transfer of polarization from electron onto nuclear spins [? ]. For this to occur, the involved energy scales must be commensurate. Thus, the main idea is to match the diering energy scales of the involved spin species by engineering their interaction, such that polarization transfer can take place. If the resulting timescales are faster than spin diusion, repeated application can lead to hyperpolarization of the nuclear spin-bath [4,5].
We employ a ramped nuclear-spin orientation via electron spin-locking (NOVEL) [6] sequence ( Fig. S1), consisting on resonantly driving the electron-spin using a Rabi eld with (variable) amplitude 2B 1 . To rst order, the eective interaction has the form: for each nuclear spin.Ŝ ± ,Î ± are ladder operators for the electron (here in the x-basis) and nuclear spins, respectively. In our experiments, we work in the sub-space spanned by the electron-spin projections m S = {0, −1}. Eq. S5 has the form of a canonical Rabi-drive, meaning that polarization transfer is driven by a ⊥ at a rate Ω = ∆ 2 + a ⊥ /2 2 . ω n = γ n B 0 is the bare nuclear Larmor frequency and ω 1 = √ 2γ n B 1 is the Rabi frequency (the √ 2 factor accounts for the fact that we work with an S = 1 electron spin). ∆ = ω n + a || /2 − ω 1 is the detuning. We assume the hyperne couplings to be much faster than the timescales for spin-diusion.
Since ∆ is in principle dierent for each nuclear spin, the parallel hyperne couplings a || introduce z-disorder in the nuclear-spin bath. This is an advantage since so-called dark states [5,7], which prevent complete bath polarization, are expected in systems without z-disorder. However, even in the presence of nuclear spins with identical (a || ,a ⊥ ) couplings, we expect a high degree of polarization upon repeated application of the NOVEL sequence.

SENSITIVE SLICE
We dene the power signal-to-noise ratio (pSNR) as the ratio between the spectral peak amplitude and the standard deviation in a reference spectral region where there are no signals. For a single nuclear spin, the total signal after K weak measurements is given by [1]: where C represents the total photon counts, the optical contrast, and p 0 {0...1} is the nuclear polarization. Two processes contribute to the total noise in a single weak measurement, a Bernoulli process resulting from quantum projection noise and a Poisson process associated with photon shot noise [8]. In our experiments, the photon shot noise greatly exceeds the quantum projection noise, and the noise variance is given by: where the last approximation is valid for small 1, which is usually the case in our experiments. The power SNR is given by: To obtain the amplitude SNR per unit time, we divide by the duration t m = t pol + Kt s of one complete measurement cycle and take the square root: We now dene the sensitivity function S (β) by the amplitude SNR expressed as a function of β, dropping the count rate factor √ C and taking the approximation sin(β) ≈ β: (S10) This equation corresponds to Eq. (6) in the main manuscript. Eq. S10 is maximized in the regions of large a ⊥ and vanishing a || . In such case Γ optical = 0, and we can rewrite where we have set p 0 = 1 for simplicity. Eq. S11 reaches a maximum S max at ΓKt s ≈ 2π/5 when Γ 0 → 0. For nite Γ 0 , an optimum point is reached when the intrinsic and induced decay rates are commensurate [1], i.e., Γ 0 ≈ Γ β and thus Γ ≈ 2Γ β = β 2 / (2t s ) = 2Γ 0 , meaning that: Since β = a ⊥ /π t β and a ⊥ ∝ r −3 , assuming short sensor readout times t s ≈ t β the sensitivity radius at Γ 0 ≈ Γ β scales as r ∝ t 1/6 β . Specically, we obtain: where a 0 = µ 0 γ n γ e 3 cos (θ m ) sin (θ m ) / (4π) and θ m = 54.7 • is the magic-angle. is the reduced Planck constant and µ 0 the free-space permeability. To maximize signal acquisition (i.e. Kt s ) while minimizing decay, we set Γ 0 = Γ β in Eq. S11 and numerically nd the point at which the sensitivity has decayed by 1/e, i.e. solve S β (opt) = 1 − e −1 S max to nd Γ 0 Kt s ≈ 3.06 ≈ 0.5π (3 + π), and thus: (S14) To verify Eqs. S13 and S14 we have numerically maximized Eq. S10 as a function of β and evaluated S, r and Γ 0 /Γ β at β max for dierent values of K and t s ≈ t β . We nd good agreement with our analytical approximations.
To compute the combined sensitive slices plotted in Fig. 5(a,d) of the main manuscript we evaluate the sensitivity function S (β) for each t β value and take the maximum max S (β i ) at each (ρ, z) value.

Information criteria
The selection of appropriate approximate models is critical for statistical inference out of experimental data. A very general methodology for model selection and parameter estimation is the use of information criteria and likelihood concepts [9]. When applying information criteria for model selection one aims to measure the distance or information, which in turn can be linked to the concept of entropy maximization, between two models [9].
From a simplied perspective, the idea is to balance the goodness of t and the complexity of a model using a cost function of the form where G σ 2 accounts for the goodness of t and depends on an unknown variance estimator σ 2 of the t residues. To minimize notation overhead, in the following we will refer to estimators of any quantity (e.g. variance, standard deviation, mean) simply with a letter.
P (K, M ) can be regarded as a penalty which depends on the sample size K and the number of parameters M . In a likelihood framework, the goodness of t can be expressed in terms of a negative likelihood function [9], where L(θ) is the likelihood of the candidate model given the data, evaluated at the maximum likelihood estimate of the model parameters θ. Computing ∆G between two models is a measure known as the Kullback-Leibler information [9]. In isolation, Eqs. S15 and S16 are meaningless; only dierences in their values when calculated using dierent models are useful. Dierent choices of P (K, M ) have been proposed, being the Akaike Information Criteria (AIC) [10] and the Bayesian Information Criteria (BIC) [11] the most common ones: where M is the total number of estimated parameters in the model and K is the sample size.
The last term in Eq. S17 is a correction introduced to compensate the AIC tendency to overt for nite samples [12]. BIC has a strong penalty, and therefore can lead to undertting in large samples [13]. The Weighted Average Information Criteria (WIC) [14] provides a solution which performs well, independent of the sample size: (S19) We use P WIC as the the penalty term P (K, M ) in Eq. (8) of the main manuscript.
In the limit of large sample sizes K 1, and assuming normally-distributed residues (errors) with a constant variance, the negative log-likelihood in Eq. S16 becomes (up to a constant term) [10]: where σ 2 is the variance estimator for the true variance σ 2 0 of the residues. Note that the variance estimator is itself a function which depends on the data and the model parameters.
For least-squares estimation, the maximum-likelihood estimator for the variance σ 2 = Σ/K, where Σ is the residual sum of squares (RSS), is usually used. However, this is a biased estimator. Instead, using the sample variance: where k are the individual residues, provides an unbiased estimator which contributes an additional penalty term and further avoids overtting [13,15].
We now derive Eq. S20 for our specic FID models. Assume a random variable x whose probability distribution f 0 (x; θ) is dened on the parameters θ = {θ m }, where m = 1 . . . M and M is the number of t parameters. A set of observations x = {x k }, where k = 1 . . . K, can be regarded as K samples of the distribution f 0 (x; θ). In our experiment, θ are the three hyperne parameters for each spin plus the three global parameters p 0 , T * 2,n and t (see main manuscript), and x is the experimental FID trace. Assuming a model f (x; θ) to be a good approximation of the true distribution f 0 , we want to nd estimates for the model parameters θ which most likely produce the measured data record x. The likelihood function, which is to be maximized, is a function on the parameters θ and corresponds to the model f evaluated at the data points x, where θ are now parameters to be estimated. Thus, to dene a likelihood function we need to make some assumption about the distribution of the data.
We now assume that the observations x are independent and their noise normally distributed.
Thus, the approximating model becomes: where k are the residuals from the tted model and σ 2 is an estimator of their variance.
It is useful to work with logarithms because they are continuous and monotonous, thus the negative log-likelihood becomes: where the approximation holds for large K and we have dropped the term K 1 + ln (2π) because it is a constant oset [10,12]. In our case, the data record x consists of a signal and a noise contribution. The latter is dominated by photon shot noise, which in the limit of many collected photons converges to a normal distribution. Using Eq. S21 as variance estimator we thus write: where the functionx(θ) represents the chosen FID model for the set of parameters θ. Eq. (S26) represents the G(θ, x) term in Eqs. (8,9) in the main manuscript.
In our analysis, we t the complex Fourier spectra instead of the time-domain FID traces.
Denoting spectral quantities with a hat, we have for each complex spectral component The noise of the Fourier spectra is also normally distributed [8]. Since we assume the measurement points x to be independent, the noise in each the realâ j and imaginaryb j parts of the complex spectrum also follows a normal distribution, that is,â j ,b j ∼ N µ = 0, σ 2 0 for K 1.â j andb j are independent for j < K/2 [8].
In terms of the power spectrum, its components areŷ j = |x j | 2 =â 2 j +b 2 j , meaning that their noise follows a chi-squared distribution with two degrees of freedom (χ 2 2 ). Therefore, the residues from the tted modelˆ j =ŷ j −ŷ j (θ) would follow a displaced chi-squared where we use standardized units and the fact that the expectation value µ for a χ 2 2 distribution equals its standard deviation σ. The model becomes: where we have taken into account that K j=1ˆ j follows a normal distribution with zero mean (central-limit theorem for K 1). Note the normalization factor 1/σ multiplying the χ 2 2 distributions such that their integrals are still unity. Taking the negative logarithm we nd: where we express the standard deviation as the square root of the variance estimator σ = √ σ 2 and again dropped constant osets. Using again Eq. S21, the negative log-likelihood function for the power spectrum becomes: whereŷ (θ) = {ŷ j } j=1···K represents the power spectrum of the chosen FID modelx(θ).
For a number D of independent datasets (i.e. FID time traces), we then write the combined likelihood function for the spectral components and their square magnitude as: where the subscripts (·) re , (·) im and (·) psd label the real, imaginary and magnitude-squared parts of the power spectra. Using Eq. S21 for the variance estimators, we nally write: where Σ d,ic is the residual sum of squares for dataset d and ic ∈ {re, im, psd}. To nd maximum likelihood estimates using the WIC we thus add Eq. S30 and Eq. S19.

Generalized Simulated Annealing (GSA)
GSA is a stochastic approach [16] that combines the original method of Classical Simulated Annealing (CSA) [17] and Fast Simulated Annealing (FSA) [18]. It has proven very useful to nd global minima of complicated, multidimensional systems with large numbers of local minima such those found in quantum chemistry, genetics or non-linear time-series. Due to its statistical nature, local minima can be escaped much more easily than with steepest-descent or gradient methods [19]. The core idea is to combine importance sampling with an articial temperature which is gradually decreased to simulate thermal noise. The algorithm uses a distorted Cauchy-Lorentz visiting distribution, to sample the parameter space [16]: where D is the number of dimensions (i.e. t parameters). Eq. S31 is also known as a Tsallis distribution. Compared to the Boltzmann (CSA) or Cauchy (FSA) distributions, Eq. S31 has a longer tail and thereby combines a local search with frequent long jumps which facilitate detrapping from local minima. The shape of the distribution is controlled by the visiting parameter q v , which also controls the cooling rate [16]: At (computational) time t and articial visiting temperature T qv (t), GSA thus requires the generation of random jumps ∆x t distributed according to Eq. S31. Here, we use an exact D-dimensional Tsallis random number generator [20]. A new conguration x t = x t−1 +∆x (t) is accepted with a probability [16]: where q a is an acceptance parameter and β ∝ T a (t) −1 is inverse acceptance temperature. T a (t) is also controlled by q v according to the cooling rate given by Eq. S32. Eq. S33 is a generalized Metropolis algorithm. Congurations with a lower cost (∆E < 0) are thus always accepted. For q a < 0, zero acceptance probability is assigned to the cases where [16]: meaning that ∆E must be smaller than β/ (1 − q a ) for hill-climbing to occur. The choice of q a < 0 in the argument of the power-law acceptance function (Eq. S33) was originally suggested because it favors lower energies. However, alternative choices of p qa for q a > 1 have also been proposed [21]. Setting q v = q a = 1 recovers CSA, and q v = 2, q a = 1 recovers FSA [16]. As T → 0 (or t → ∞), the algorithm becomes a steepest-descent method [19].
For our calculations, we selected q v = 2.7, q a = −1 and a large initial acceptance temperature T a (1) ∼ 10 3 to facilitate the escape from local minima. A technique to accelerate convergence is to set: which is equivalent to a monotonously decreasing q a ∝ t −1 [19]. It is worth noting that the computational times associated with the visiting and acceptance temperatures need not to be equal. To improve convergence towards physically-relevant congurations, we modied the acceptance probability p qa → p qa using a step-function: where d ij is the estimated spatial distance between any two spins i and j at any computational time t, and r cc = 0.154 nm is the diamond bond-length. In other words, we reject any proposed nuclear-spin congurations where any two spins are closer to each other than the smallest 13 C -13 C distance in a diamond crystal.
We also assigned an individual visiting temperature T qv,m to each t parameter θ m in order to account for their dierent scales/units, and set their initial values large enough such that their entire search intervals were covered at t = 0 (see Sec. 3 3.4). Furthermore, we observed better results when normalizing the visiting and acceptance times to a few hundred Monte Carlo steps. Specically, we used t v = t/400 and t a = t/200. This re-scaling of the computational time eectively helps to more widely explore the parameter space at each temperature step. We also observed better convergence when updating one parameter at a time instead of the whole parameter set.
For a given number n of spins, we run the algorithm starting from ∼ 10 2 −10 3 initial random spin congurations, where each conguration is drawn by randomly selecting n points within a radius of 2.5 nm around the coordinate origin. We ran the algorithm over 5000 iterations and took the best result. As a guide, a single run (1 initial spin conguration) of the algorithm for n ≈ 20 took ∼ 3 h.
It is worth noting that the performance of GSA, and in general of annealing algorithms, is problem dependent and therefore some parameter tuning is usually required. Our selections for q v and q a lie within ranges where good performance has been reported [16].

Estimation of uncertainties
As mentioned in the main text, we simultaneously tted all spectra. To get an estimate for the uncertainties in the t parameters θ, we resampled the data using a bootstrapping protocol [22]. We calculated the mean and standard deviation of the t residues (green dots in Fig. 4b,c in the main text) for the Re[FFT] and Im[FFT] parts of each spectrum. We then used the resulting deviations to generate normally-distributed noise and thereby created a set of P = 100 new, articial datasets. Each newly generated dataset (set of 4 spectra) was minimized again (using a steepest-descent search) and the original-dataset solution as starting point. We thereby generated a set θ m,p ≡ {θ m } p of estimates for each t parameter θ m , where p = 1 . . . P . Finally, we calculated the t results as:

Selection of parameter intervals
To start our minimization routines, we chose random starting values from a pre-dened interval for each parameter.
For a ||,i , we dened the interval by looking at the spectral support across all power spectra and set it to 2π ×[2 (f 0 − f k ) . . . 2 (f 0 + f k )], where f 0 is the central-peak frequency (assumed to be around the bare Larmor frequency) and we selected f k such that all visible peaks in the power spectra where covered. Specically, we selected f k = 4.5 kHz and f k = 4 kHz for dataset 1 and 2, respectively.
where t is a phenomenological time constant that is roughly of the duration of the laser pulse. t is not correlated with the 13 C environment [2,3] and is a global free t parameter in our maximum likelihood estimation. The total decay rate is then the sum of all contributions The third parameter, the initial nuclear polarization p 0 , in principle diers for every 13 C nuclei, because the polarization rate is proportional to a ⊥ . However, because we repeat the polarization transfer for typically > 10 3 cycles, we expect that the polarization of all nuclei within a resonant slice become saturated. (This level is reduced for pairs or clusters of spins with near-identical (a || , a ⊥ ) coupling constants, however, not many such pairs are expected for our dilute 13 C concentration). This assumption is supported by the fact that spectra show little change in peak intensities once the number of cycles exceeds 10 3 (see Fig. 5) and is consistent with similar models using cross-relaxation induced polarization (CRIP) [23]. Therefore, we assume p 0 to become homogeneous within our sensitive slices and treat it as a global parameter.

INTERNUCLEAR COUPLINGS
Continuous weak measurements can be regarded as the equivalent of inductive detection in conventional NMR. In this way, they capture the full evolution of the nuclear environment, including potential nuclear spin-spin interactions. Such nuclear spin-spin interactions lead to additional peaks in the spectra, which can interfere with our minimization strategy.
To illustrate the eect of nuclear spin-spin interactions and analyze whether they aect our present experiment, we performed a density matrix simulation with n = 3 randomly positioned nuclear spins using a Master equation approach. To estimate whether such couplings are potentially present in our experimental spectra, we computed the average number of spin pairs yielding resolvable couplings within our spectral resolution given by Γ = max{1/ (Kt s ) , Γ 0 }. The FFT resolution is 156 Hz and 110 Hz for NV1 and NV2, and the estimated nuclear linewidth is Γ 0 = 1/T * 2,n = 86 ± 7 Hz for NV1 and 120 ± 11 Hz for NV2, limited by uctuations in the magnetic bias eld. We generate 10 4 random 13 C nuclear spin congurations within the sensitive volumes of NV1 and NV2 (Fig.   5a,d in the main text), respectively, and count the number of pairs with couplings larger than Γ. Considering the worst-case scenario (neglecting the inuence of measurement and readout-induced dephasing), we nd on average 1.14 ± 1.06 (69% probability that at least one spin-pair is present) and 2.14 ± 1.43 (89% probability) spin-pairs for NV1 and NV2 , respectively. Thus, although the presence of one or two spin pairs is likely, this will not aect the majority of peaks observed in the spectrum. Moreover, since the two spins forming a pair will have very similar hyperne coupling constants, their peaks will be spectrally close, and the peak splitting may not be resolved because of peak overlap.
As a control, we have inspected the reconstructed spatial positions of the mapped nuclei and computed their relative spatial distances. We nd an average distance between the closest pairs of 0.212 nm and 0.178 nm for dataset 1 and 2, respectively. This corresponds to homonuclear couplings of 544 Hz and 645 Hz. However, due to the inversion symmetry of the diamond crystal, each recovered nuclear position has a 50% probability to be mirrored along the coordinate origin.     Figure 4c. The parameters for the FID shown in Fig. 4a are those for t β = 4.944 µs. A detailed pulse diagram is given the Extended Data Fig. 1b  (upper bound) of 11.68 ± 1.03 ms. Errors for a || , a ⊥ and φ are calculated by bootstrapping. Errors for r, ϑ, φ and the global parameters are calculated using Monte Carlo error propagation. The estimated uncertainty volume δV (also depicted in Fig. 5c) is given in cubic Angstrom. * denote spins whose uncertainty volume is larger than the volume per carbon atom in the diamond lattice (5.69 Å 3 ).  Fig. 5f) is given in cubic Angstrom. * denote spins whose uncertainty volume is larger than the volume per carbon atom in the diamond lattice (5.69 Å 3 ).

Polarization Excitation Detection
Start of nuclear precession FIG. S1. Experimental sequence. Horizontal axis is time. Yellow and red blocks are microwave (MW) π and π/2 pulses, respectively, that act on the lower-energy (m S = 0 to −1) transition of the electronic spin. Black sinusoidal curves are radio-frequency (RF) pulses acting on the 15 N (frequency ∼ 0.86 MHz) and 13 C (frequency ∼ 2.1 MHz) nuclear spins. Green blocks are laser pulses with typical duration 1.5 µs. After each laser pulse we include a delay of ∼ 1 µs to allow the NV center to decay to the electronic ground state. Our protocol proceeds in three steps: Polarization, excitation and detection. The polarization step consists of M ∼ 10 3 repetitions of the NOVEL sequence [6,24]. The NOVEL sequence consists of a (π/2) X rotation followed by a spin-lock along Y of duration t novel . We use a linear amplitude ramp of the spin-lock pulse of ∼ 10% to improve robustness of the polarization transfer. Every four NOVEL cycles we add a c-NOT gate to polarize the 15 N nuclear spin to improve the polarization transfer eciency. The c-NOT gate consists of a selective MW π pulse (conditional on the 15 N spin state) followed by a selective RF π pulse (conditional to the electronic spin state) [25]. The excitation step consists of a single, non-selective RF pulse on the 13 C nuclear spins. The detection step consists of a train of K ∼ 10 3 weak measurement readouts [1,8,26,27]. Each weak measurement block consists of an XY8 sequence [28] of between 4-24 MW π pulses, sandwiched by two MW π/2 pulses that have a relative phase shift of 90 • [29]. The delay between π pulses (2τ ) is adjusted to the inverse of the expected 13 C Larmor frequency f , τ ≈ 1/(4f ). The XY8 sequence is followed by a delay time t d . During the delay, the electronic spin state is read out by applying a green laser pulse and turning on the photo-detector for ∼ 300 ns. An additional MW π pulse at t d /2 is used to average the hyperne interaction [29]. The sampling time t s = t β + t d is adjusted, via the delay time, to match an (j + 0.25) multiple of the expected 13 C Larmor precession period, where j is an integer and the 0.25 implements a time-proportional phase increment (TPPI) [30]. All measurement parameters are collected in Tables S1 and S2.  Complex Fourier spectra of the 13 C environment (NV2) for a series of interaction times t β . Shown are from top to bottom (vertically oset for clarity): power spectrum (PSD), real part of the complex spectrum (Re[FFT]), t residues for Re[FFT], imaginary part of the complex spectrum (Im[FFT]), and t residues for Im [FFT]. The horizontal axis shows the spectral shift relative to the 13 C Larmor frequency. To illustrate the agreement between our analytical FID model (Eq. 2 in the main text) and the numerical density-matrix simulations, we take the measurement parameters (Table S2) and the best-t parameters obtained for NV2 and simulate the resulting set of spectra (blue traces) using the density matrix formalism. The gray traces show the corresponding spectra obtained using our analytical model. The negligible residues (green traces), well smaller than those in Figs.4b,c in the main text, highlight the excellent match between the density-matrix simulations and our second-order analytical model. where τ is the spacing between π pulses (Fig. 2b in the main text). In the absence of 13 C -13 C couplings (top panel), we observe three peaks associated with the three nuclear spins. The oset from f c is proportional to the parallel hyperne coupling constant (grey arrows). In contrast, for non-vanishing internuclear couplings (bottom panel), splittings proportional to g ij appear but only those within the spectral resolution δf = 1/(Kt s ) = 156.25 Hz set by the Fourier limit of the time trace are resolved (orange arrows).