Power-law scaling of correlations in statistically polarised nano-NMR

Diffusion noise is a major source of spectral line broadening in liquid state nano-scale nuclear magnetic resonance with shallow nitrogen-vacancy centres, whose main consequence is a limited spectral resolution. This limitation arises by virtue of the widely accepted assumption that nuclear spin signal correlations decay exponentially in nano-NMR. However, a more accurate analysis of diffusion shows that correlations survive for a longer time due to a power-law scaling, yielding the possibility for improved resolution and altering our understanding of diffusion at the nano-scale. Nevertheless, such behaviour remains to be demonstrated in experiments. Using three different experimental setups and disparate measurement techniques, we present overwhelming evidence of power-law decay of correlations. These result in sharp-peaked spectral lines, for which diffusion broadening need not be a limitation to resolution.


I. INTRODUCTION
Nuclear magnetic resonance (NMR) spectroscopy is widely used in the life and material sciences. However, classical NMR techniques are not readily available on the single cell level.
Therefore, further scientific breakthrough in biosensing may be hindered by the invasive analysis techniques that most conventional approaches require in this regime. Existing methodologies entail tagging, e.g., with fluorescent nano-particles, cryogenic temperatures, or high magnetic fields [1][2][3][4][5]. Either causes substantial modifications to the sample, meaning that its natural properties cannot be determined accurately. Nuclear magnetic resonance at the nano-scale (nano-NMR) with spin sensors, such as nitrogen-vacancy (NV) centres, offers a label-free, room temperature approach capable of studying biologically relevant samples down to the molecular level without altering their properties [6][7][8][9]. NV centres have already been shown capable of detecting the magnetic field created by nano-sized distributions of nuclei in liquid samples at ambient temperature [4,[10][11][12][13][14]. However, a fundamental limitation to the ability of resolving spectral lines is thought to exist.
The nano-NMR approach with spin sensors relies on nano-scale sample sizes containing a sufficiently low number of nuclei, in which statistical fluctuations are significant enough to overcome the thermal averaging of nuclei orientation, leading to the emergence of time-correlated magnetic fields without resorting to sample polarisation. These fields can be detected at room temperature with quantum probes such as a shallow NV Ô nicolas.staudenmaier@uni-ulm.de; À These authors contributed equally centre, whose interaction region is typically considered a hemisphere with a radius of the order of the NV depth as shown in Fig. 1(b) [4,[15][16][17]. Sample sizes within the NV detection region are on the order of 10 −24 −10 −20 L for several nm deep NV centres. There, the statistical polarisation exceeds the thermal by more than three orders of magnitude in ambient conditions [10,[18][19][20][21][22]. Statistically polarised nano-NMR with spin sensors enables studying sample properties inaccessible with any other NMR protocol.
The main challenge for nano-NMR with statistically polarised liquid samples is posed by the diffusion of molecules out of the interaction region defined by the probe. Diffusion changes the spatial distribution of statistical polarisation, and leads to a decay of the correlations of the magnetic field signal in the probe, which is typically considered exponential (see Fig. 1(c)) [16,23]. Then, measurements performed beyond the characteristic time of the exponential decay yield no information about the spectral signatures of interest. When this time is short, not enough information can be gathered to allow for spectral reconstruction in post-processing [24], creating a resolution problem. The corresponding spectral line for correlations with exponential decay is Lorentzian, as shown in Fig. 1(d), where frequency resolution is typically defined by its full width at half maximum [25][26][27].
Recent theoretical analysis has claimed that diffusion induced decay of correlations follows a power-law at long times (see Fig. 1(c)), due to the specific dipole-dipole interaction of the shallow NV centre sensors and the nuclei in the sample [28]. Such correlations allow suitable measurement protocols to extract sufficient information beyond the characteristic decay time T D , so estimating frequencies smaller than 1/T D becomes feasible [24]. Moreover, they provide a more accurate route to measure the diffusion coefficient in liquid samples. Although previous work in liquid NV nano-NMR has shown an instance where a non-exponential function fits better the observed correlation, this was attributed to surface effects that led to reduced translational diffusion of the nuclei close to the surface [15]. Thus, experimental evidence of interaction dependent power-law behaviour has not been demonstrated yet, to the best of our knowledge. The main reason might be that the early time decay of correlations resembles an exponential function, and it is only the long-time decay beyond the diffusion time that is critical to show deviations from the exponential paradigm [29].
In this work, we provide compelling experimental evidence of a power-law decay of correlations at long times, in accordance with the theoretical prediction in Ref. [28]. We demonstrate our results using three distinct measurement settings and several independent statistical analysis tools. We include measurements of the magnetic field originated in the fluorine nuclear spins of the sample in one of the experiments. These nuclei are assumed to be separated from the diamond surface by an immobilised proton layer adsorbed onto the surface, and permit us to rule out surface effects as the source of a long correlation tail [15]. Our results pave the way for statistically polarised, room temperature nano-NMR experiments with biomedically relevant samples, with a scope that is not limited to NV centres but that, due to the power-law scaling stemming from dipole-dipole coupling, extends to any sensor that is based on such interactions.

A. Theory
The NV centre is a point defect in the diamond lattice composed of a substitutional nitrogen atom and an adjacent vacancy which, together, are described as a system whose ground state is a 3 A 2 spin triplet. Applying a bias magnetic field, the degeneracy of the |m s = ±1 spin levels is lifted, and the NV centre can be used as an effective two-level system, as shown in Fig. 1(a). Owing to its spin properties, an initial equal superposition state of two of the relevant states (e.g. {|0 , | − 1 }) accumulates a phase due to interaction with an external magnetic field. A rotation of the evolved state into the measurement basis reflects the accumulated phase as a population imbalance between the NV centre spin states, which can be measured thanks to the distinct fluorescence rates of each of the states. Such an initial state is most sensitive to the slowest frequency components of the (magnetic) noise. Therefore, dynamical decoupling (DD) sequences are routinely used to filter out the effect of slow noise, prolong coherence times, and enable sensing of specific frequency components of the signal [18,30,31]. The NV centre is located at a depth d below the diamond surface. In the standard diffusion model a hemisphere sensing volume (blue) of radius ∼ d is considered, in which the nuclei interact with the NV centre. Alternatively, we accurately take into account the 1/r 3 dependence of the dipole-dipole interaction. The colour intensity of the orange hemisphere indicates the interaction strength of the nuclear spins (red) with the NV centre. (c) Decay of the auto-correlation according to the exponential and the power-law model. With the accurate dipole-dipole model, decay is (nearly) exponential for short times (t TD) and a power-law for t TD with TD the characteristic diffusion time. (d) Power spectrum of the NMR signal. The power-law model shows a distinct, sharp peak whereas the exponential model resembles a Lorentzian lineshape.
In statistically polarised nano-NMR, the magnetic signal originates on an ensemble of nuclei within a small sensing volume oscillating at their Larmor frequency. Diffusion of molecules induces magnetic noise which we model as a stationary Gaussian process with zero mean. For all experiments, we estimate the possible spectral broadening effect from back-action on the NV centre or a gradient field due to the NV sensor [32] to be at least one order of magnitude smaller than the broadening caused by diffusion (see Appendix E for details). Therefore we consider it negligible for the ensuing analysis of the decay of correlations due to diffusion. Throughout the experiments described here, we will be probing the auto-correlation of the time evolution of the NV centre interacting with the nuclear spins, which is given by the noise covariance of the external signal, Here, Φ rms is the accumulated phase on the NV probe, which is proportional to the root-mean-square field B rms originating from the statistically polarised nuclei, and δ is the typically small undersampling frequency of the detected signal. The latter is equal to the difference between the Larmor frequency and the frequency defined by the sampling times in the measurement protocol.
The envelope C(t/T D ) describes the effect of noisy fluctuations due to diffusion, with characteristic diffusion time T D = d 2 D where D is the diffusion coefficient [16]. Its particular shape is determined by the specific interaction between the sample and the sensor. For NV centres, it is the magnetic dipole-dipole interaction with the nuclei. In the common paradigm, this interaction lasts for T D . Thus, nuclei that diffuse beyond a distance √ T D D ∼ d are too far to interact with the NV centre. This defines an interaction volume which is a hemisphere of radius ∼ d above the surface of the diamond, illustrated by the solid colour region in Fig. 1(b) [16]. Correlations of the magnetic field at the NV position reflect the diffusion of nuclei in or out of this hemisphere, and the assumption about the interaction duration means that correlations are lost when nuclei diffuse out of the interaction region. In this picture C(t/T D ) ∝ exp(−t/T D ). The exponential model is widely accepted in the field due to its success when the measurement times are short and the frequencies probed are easily resolved [15,16,33]. However, such approximation has profound implications over measurements about the diffusion coefficient D in microfluids [28,29]. Moreover, it results in a resolution problem which cannot be overcome [24].
A more careful analysis of the dipole-dipole interaction renders a significantly different behaviour [28]. Heuristically, it can be understood as follows: the magnetic field generated at the NV position as a result of a nucleus located at a point r is B(t) ∝ 1 r 3 . Hence, the correlation of the magnetic field is B(t)B(0) ∝ 1 where the coordinates r and r of a specific nucleus are related by its diffusive motion. This connection is usually observed in the second moment r 2 = x 2 + y 2 + z 2 = x 2 + y 2 + 4Dt + z 2 ≈ r 2 + 6Dt, where 2Dt is the variance in the nuclei positions in one dimension. Diffusion of molecules close to the surface is free along the surface direction x and y but limited in the orthogonal z. Hence, the second equality is always valid, while the last approximation is correct for t T D or t T D , where the exact geometry of the problem is less important. Substituting the second moment into the auto-correlation function yields T D only close nuclei contribute, the hemisphere region approximation is valid, and correlations replicate the exponential behaviour. At long times, the interaction region expands beyond the hemisphere paradigm, as shown in Fig. 2(b), and correlations decay as a power law B(t T D )B(0) ∝ 1 (6Dt) 3/2 . The full expression for C(t/T D ) in this scenario is given by G(t/T D ) in Eq. (12) (see Methods). A more detailed derivation of the asymptotic behaviour of the power spectrum due to these correlations shown in Fig. 1(d) can be found in Appendix C and in Ref. [28].

B. Experiments
We performed a series of experiments comprising three different setups: correlation spectroscopy measurements with single NV centres and with an ensemble of NV centres, and quantum heterodyne (Qdyne) with single NV centres. The experiments and results are reported in each subsection, altogether reinforcing the validity of the power-law scaling of the correlation's decay. Description of the measurement sequences can be found in Section IV (Methods). All experimental details and parameters are thoroughly described in Appendix D.

Correlation spectroscopy with a single NV centre
The first set of experiments is carried out on a shallow NV centre which is located at a depth ≈ 2.9 nm below the surface of an isotopically enriched (99.999% 12 C) diamond sample grown by chemical vapour deposition [34,35]. The shallow NV centre is addressed via a fluorescent confocal scan and is initialised and read out using a 532 nm laser pulse. We apply a bias field of ≈ 450 G along the NV axis using a permanent magnet to lift the m s = ±1 degeneracy (see Fig. 1(a)). Microwave pulses for coherent control of the NV centre are applied through a copper wire of 20 µm diameter strapped across the diamond, and we use state-dependent photoluminescence measurements to detect the population of the spin states. The interaction of the NV centre with hydrogen nuclear spins in the immersion oil (Fluka 10976, viscosity of 400 cSt) is measured.
The NV centre is initialised into its m s = 0 state, and read out optically. We use correlation spectroscopy to probe the diffusion spectrum (see Methods IV C for details). First, we prepare the system in state |Y = |0 +i|1 √ 2 by applying a π/2 pulse around x-axis and use two Knill dynamical decoupling (KDD4) sequences, separated by a waiting time difference T [36][37][38][39]. The centres of two subsequent π pulses in the sequence are separated by the time τ = 1/2f L , where f L is the Larmor frequency of the nuclear spins [36][37][38][39], which allows for maximum phase accumulation of the NV centre due to the interaction with the nuclear spins. In this case, the Larmor frequency is 1.9 MHz. The information about the accumulated phase during the first decoupling sequence is mapped onto the spin population by another π/2 pulse applied around the y-axis. This result is then correlated with the second interrogation starting at time t = N τ +T of the second KDD4 decoupling sequence, where N = 20 is the number of pulses in KDD4. The obtained signal is proportional to the correlation of the accumulated phases during the two sequences, and allows us to directly extract the auto-correlation of the signal.  Figure 2. Experimental data from correlation spectroscopy measurement with a single NV centre. (a) Signal of the correlation spectroscopy measurement as a function of the time difference between the two DD sequences (black dots) and fits of the exponential (blue line) and power-law (orange line) models. Labels show the goodness of fit. The power-law model, described by Eq. (12), offers a better fit to the data, which can be appreciated also in the C(t/TD) envelope (dashed lines). In both models, the initial amplitude for the fitting algorithm is fixed, obtained from an independent estimate in a power spectrum measurement. (b) Fast Fourier transform (FFT) of the experimental data together with FFT from the fitted data in (a) to the exponential and the power-law decay models, shown for illustrative purposes. The sharpness displayed by the power law model allows for a more precise frequency estimation than the exponential model. fitted to the exponential and the power-law decay models described by the correlation function Eq. (1) with envelopes Eqs. (11) or (12), respectively. Note that Eq. (12) comes from exact mathematical calculations such that the envelope is valid for any time and shows a power-law decay for longer correlation times. We restrict the parameter space for both models by estimating the initial value of the signal amplitude ∼ Φ 2 rms from an independent power spectrum measurement, which is typically used for finding the depth of the NV centre (see [16] and subsection D 2 in the Appendix). This is preferable, as the signal sampling times are chosen for better estimation of long-lived correlations. Thus, the estimate of the initial amplitude is not efficient and differs significantly between the two models, making comparison difficult. We then use the estimated value as a fixed input parameter for both models and apply non-linear least squares fitting of the signal. We start from 500 different initial conditions and take the best fit.
The obtained fits in Fig. 2(a) show that the power-law decay model performs better with an R 2 ≈ 0.93 in comparison to R 2 ≈ 0.90 for the exponential model. Its goodness of fit is evident especially at long times, e.g., between 50 − 100 µs, which correspond to more than three times the expected diffusion time T D ≈ 17 µs, highlighting the importance of long-lived correlations for an accurate estimation of the diffusion coefficient, and emphasising the need to estimate independently the initial contrast and the characteristic decay time in order to obtain accurate measurements of the diffusion coefficient.
In addition, Fig. 2(b) shows a comparison of the Fast Fourier transform (FFT) of the experimental data and the FFT from the fitted data to the exponential and the power-law decay models. The results show that the power-law decay model provides a better fit to the data than the exponential model. Its spectral linewidth is narrower than with the exponential fit, which is due to the long-lived correlations of the power-law model. Such a feature enables improved precision and resolution in sensing experiments due the sharply peaked spectrum of the power-law model (see also Fig. 1(d)). The slight broadening of the experimental FFT is mainly due to the fitting procedure, which uses zero padding, so we prolong the time domains of the fits of both theoretical models accordingly. This results in a sharper peak for the power-law model, as expected from theory. Such broadening can also be present due to extra noise sources which produce exponential decay at a much slower rate than that of diffusion noise, as described in Appendix E.

Correlation spectroscopy with an NV centres ensemble
The next set of experiments is carried out on an ensemble of NV centres. Similarly to the single NV experiment, we use correlation spectroscopy to probe the diffusion spectrum. Here, a perfluoropolyether oil (Fomblin Y, Sigma Aldrich 317926, viscosity of 60 cSt) is analysed. We detect the signal coming from nuclei in the Fomblin oil, which sits above a proton layer located on top of the diamond surface, preventing fluorine nuclei from sticking to the surface. Thus, we can rule out surface effects as the underlying cause of the power-law behaviour.
The correlation spectroscopy protocol is the same as in the single NV experiment, except that each of the two DD sequences is KDD4-4, i.e., KDD4, repeated four times, accounting for the larger average depth of the NV centres in the ensemble of about 9 nm and the higher Larmor frequency, for which the centres of the decoupling pulses are separated by τ = 1/2f L = 135.6 ns (f L = 3.687 MHz). Figure 3(a) displays the measured signal vs. the time t between the beginnings of the two DD sequences, together with its fit to the two decay models, exponential and power-law, as described.
We use the same data analysis approach as for the single NV centre experiments. Since the measurements focus on targeting the later times of the correlations, sampling times are chosen accordingly. Thus, instead of considering the signal amplitude as a free parameter, which is inefficient and leads to big differences between models (see subsection D 3 in the Appendix), we estimate it from independent power spectrum measurements and provide it as a fixed input parameter to each of the 500 non-linear least squares fittings of the signal, analogously to the single NV experiment, and consider the best fit. The obtained fits in Fig. 3(a) show that the power-law decay model provides a better fit to the experiment, in line with the results obtained with the single NV centre. The goodness of the fit of the power law model is especially better at long times, e.g., between 200 − 400 µs, similarly to the single NV centre experiment. Finally, Fig. 3(b) compares the FFT of experimental data and the FFTs of the exponential and power-law fits to the experimental signal. We observe that not only the power-law model provides a better fit to the power spectrum of the data but also a sharply peaked spectrum, which is in principle key for improved precision and resolution.

Qdyne with a single NV centre
Qdyne experiments rely on individual storage of each single readout performed at a constant rate [7,40], in contrast to conventional measurements (as correlation spectroscopy) where the acquired data is averaged. In post-processing, the auto-correlation of the obtained time trace reveals the accumulated signal. For that reason a single dynamical decoupling measurement sequence is repeated continuously. The experiments were performed with NV centres at a moderate depth between 8 and 15 nm. We performed six different Qdyne experiments spanning a total measurement time of ≈ 40 days. The Fluka immersion oil is used for five and another perfluoropolyether oil (Fomblin Y, Sigma Aldrich 317993, viscosity of 1508 cSt) for one measurement. XY8 dynamical decoupling is employed with different number of repetitions according to the depth of the NV centre where for all measurements the Larmor frequency is about f L ≈ 2 MHz. Full details of the experimental parameters are given in subsection D 4 in the Appendix. For each experiment, the resultant time-traces are divided into 15 minutes slices (∼ 10 7 measurements) for which the auto-correlation is calculated. Further noise reduction is achieved by averaging 20 of these auto-correlations.  . Experimental data from correlation spectroscopy measurement with an ensemble of NV centres. (a) Signal of the correlation spectroscopy measurement vs. the time difference between the two DD sequences (black dots) and fits of the exponential (blue line) and power-law models (orange line) with their goodness of fit and the corresponding envelopes shown as dashed lines. The power law model shows a better fit to the data, especially at long times. In both fits we fixed the initial amplitude of the fitted signal to an estimate from an independent power spectrum measurement. (b) Fast Fourier transform (FFT) of the experimental data together with the FFT from the fitted data to the exponential and the power-law decay models, shown here for illustrative purposes. Similarly to the single NV centre experiments, the power law model provides a better fit of the sharp peaked spectral line, thus allowing for better frequency estimation. a non-physical phase ϕ, included for reasons of stability in the numerical analysis. For each C(t/T D ) model considered throughout, we perform a local optimisation with fixed initial parameters, taking the initial frequency from FFT of the auto-correlation of the full time-trace, and a global optimisation repeated 500 times with random initial parameters drawn from uniform distributions, in which the best fitting is determined by the highest R 2 (see IV E in Methods for full details). The latter analysis mimics the procedure in the case when FFT yields no results, e.g when the signal decay is too fast.
Focusing on the first experiment (see Table I), the upper row of Fig. 4 displays the histogram of 18 frequency estimators for global optimisation of each model, with ϕ a free fitting parameter in Fig. 4(a) and kept fixed (ϕ = 0) in Fig. 4(b). Results are consistent with the detection of a frequency δ ≈ 900 Hz, on a sample with estimated T D = 400 µs, i.e. a frequency smaller than the noise bandwidth defined by the characteristic decay time of the auto-correlation. FFT analysis of the full signal's correlation (see Fig. 10) indicates a frequency centred at 970 Hz with a FWHM of 1205 Hz. In both instances of Fig. 4, the root-mean-square error (rmse) of the estimator is smaller in the case of power-law correlations fitting as compared to exponential; 292 Hz vs 468 Hz in Fig. 4(a) and 242 Hz vs 317 Hz in Fig. 4(b), with the standard deviation showing similar scaling.
Further evidence of power-law correlations is provided by generating artificial time-traces numerically. We simulate Qdyne experiments with parameters similar to those in the experiment shown in the upper row of Fig. 4, with two distinct noise models which produce two data sets with either exponential or power-law correlations as explained in IV F in Methods. Each of the data sets is simultaneously analysed with both models, and the resulting histograms for the frequency estimators displayed in the lower row of Fig. 4. In Fig. 4(c), we show the results for data generated with power-law correlations, where a peak at a frequency ≈ 900 Hz can be estimated. Power-law model fitting yields a rmse of 195 Hz, while exponential fitting of the power-law auto-correlations is slightly wider at 270 Hz. Fig. 4(d) displays the results corresponding to time-traces generated with exponential correlations. Here, no frequency is detected and the rmse is 368 Hz for exponential fitting and 385 for the power-law model. Notice that the search region limit is 404 Hz for a flat histogram rmse, indicating that both estimators are essentially featureless. Together with the detection of a frequency below the noise bandwidth, our results are consistent with power-law correlations.
Note that in Fig. 4(a), fitting to exponential correlations results in a histogram which goes to the extreme of the search region, signalling that the fitting generally fails. Only by reducing the number of parameters, restricting ϕ to 0, is the problem ameliorated. This does not happen for power-law model fitting. Results for local optimisation (displayed in Fig. 11) display similar behaviour and, for power-law correlations, there exists a trend in the rmse, which diminishes as the fitting model is refined, e.g. in local optimisation or with fixed parameters. Such trend is not found for exponential correlations fitting. That the fitting fails when an extra, meaningless parameter is added provides further evidence that it is the power-law model which best describes the experimental data.
We complete the analysis by considering the rmse ratio between power-law and exponential fittings. Here, most experiments have a frequency higher than the noise bandwidth (see Table I for details) and therefore frequency estimation poses no problem with either model [24]. Hence, we resort to compare which model provides the most accurate fitting according to a smaller rmse.  Histograms of frequency estimators. The data is obtained from non-linear least squares fitting of auto-correlations for experimental (upper row) and numerical (lower row) Qdyne time-traces. In blue, fitting function assumes exponential decay of correlations while in orange the power-law model is considered. Each estimator corresponds to the highest R 2 in 500 fittings of the same time-trace. In (a), the phase ϕ is left as a free fitting parameter while in (b) it is kept fixed to ϕ = 0. (c) displays the results for numerical data generated for auto-correlations with power-law decay while in (d) the data is generated for exponentially decaying auto-correlations. Root-mean-square error ratio between local optimisation for power-law and exponential model fittings. The ratio is plotted as a function of the frequency δ and the noise bandwidth defined by TD. Each dot is calculated with the rmse of histograms as those shown in Fig. 4. In black squares, experimental results. Orange circles display the ratio for data generated with power-law decay while blue stars correspond to data generated with exponential decay. Solid lines show the mean ratio for all dots while shaded areas correspond to the standard deviation for each mean. Note that experiment six in Table I is not shown but it is considered for the mean ratio. Figure 5 displays the rmse ratio for local optimisation of the experimental data, with the dashed line showing the average ratio standing at 0.62 ± 0.16, which means that power-law fitting of the experimental data has an average improvement of 40 % of the rmse with respect to exponential fitting. For global optimisation the improvement is 20 % (see Fig. 12). A numerical analysis of artificial time-traces generated with either correlation model and analysed with both models shows the striking difference existing when analysing the data with the incorrect model. For data generated with power-law correlations the root-mean-square error (rmse) ratio between local optimisation for power-law and exponential model fittings is 0.65 ± 0.27 while for data generated with exponential correlations the ratio is 1.42 ± 0.29, further confirming the experimental results of detection of power-law decay of correlations.

Diffusion coefficient estimation
Estimating the diffusion coefficient on the micro-and nano-scales is an important but challenging task, and current techniques are prone to errors [29]. Using NV centres to calculate D from the characteristic time T D in an exponential envelope of correlations leads to inaccuracies, as this model disregards correlations beyond the distance d, which nonetheless carry diffusion information, as shown by a more rigorous modelling [28]. The ability to measure power-law correlations, as we have done throughout this article, opens up the possibility to estimate D more accurately.
The diffusion time T D is defined as the characteristic time to diffuse at a distance d = √ DT D . For a hemispherical sensing volume, the number of nuclei is rms . Thus, the procedure to obtain D shall be as follows. The B rms can be independently estimated from, e.g., power spectrum measurements (as we do here) or rapid Ramsey measurements; hence the depth d of the NV centre can be obtained. T D is then estimated from a fitting procedure as described in the previous sections, with the B rms a fixed parameter. From these, the diffusion coefficient D can be obtained.
Following this procedure, we estimate D for the correlation spectroscopy measurements with single NV centres shown in Fig. 2. For the exponential model we The estimate of the diffusion coefficient of the oil in the experiment, based on its viscosity characteristics is of the order of magnitude [16] D ≈ 6 × 10 −13 m 2 s −1 , which is within the confidence intervals of the estimates of both models. While the reported confidence intervals are rather large, they could be substantially decreased in future experiments, specifically designed to improve the efficiency of the procedure, e.g., by improving the signal-to-noise ratio and obtaining a clear signal at even longer times, as expected theoretically from the power-law model. Investigation of surface effects on diffusion are also envisaged, e.g., by probing it with NV centres at different depths with surface effects expected to be more pronounced with shallower NV centres.

III. DISCUSSION
An exponential modelling for correlation's decay is the natural assumption for nano-NMR experiments. Yet this hypothesis has profound implications on the possible applications of any such experiments. In particular, exponential correlations lead to Lorentzian lineshapes for which spectral resolution is fundamentally limited. It is possible, however, to refine the original assumption by carefully examining the microscopic diffusion process which is the cause of correlations decay. Diffusion causes power-law decay of correlations at long times for the dipole-dipole interaction that is characteristic for the NV centre nano-NMR spectrometer [28]. This has substantial consequences, as spectral lineshapes become sharp-peaked, for which in theory there is no limit to resolution [24]. Yet, measuring deviations from the exponential paradigm and utilising them for improved spectral analysis is a challenging task due to experimental noise. Furthermore, by considering the full interaction model for correlations rather than the truncated exponential model, it shall be possible to measure more precisely the diffusion coefficient D in microfluids, a long sought goal, where current methods have only limited accuracy [2,41].
By probing the auto-correlation of the dipole-dipole interaction between NV centres and nuclear spins, and prolonging the measurement time to several times the characteristic exponential decay time, we show strong evidence supporting a power-law decay of correlations as they provide a better fit to the experimental results. Furthermore, with the Qdyne data frequency estimation yields a results, which is 40% more precise than with the exponential model. We can also exclude surface effects causing a long correlation tail by detecting the signal created by fluorine nuclei lying on top of a proton layer adhered to the diamond surface, in which such effects originate [15]. Additionally, we show that Qdyne measurements auto-correlations are compatible with the same power-law decay to the extent that in the regime where exponential correlations limit frequency resolution, a frequency can still be detected. Finally, we rule out the possibility of exponential correlations by showing with numerical analysis that such correlations in the same regime would lead to featureless histograms or altogether wide spectral lines in which frequency estimation is not possible.
Our results offer a way to significantly improve spectral resolution, and enable broad applications of nano-NMR with NV centres. The power-law scaling of correlations which originates in the dipole-dipole interaction between sensor and sample is key to this result, extending the applicability of our results to any quantum sensor based on this interaction, such as Rydberg atoms, squid based sensors, molecular quantum sensors, or alternative colour centres such as silicon carbide [42]. In addition, the ability to measure accurately correlations of nano-sized fluid samples, opens the door to study flow properties at these scales.
It has been recently shown that single trajectory auto-correlations differ from ensemble averaged (multiple trajectories) auto-correlations for anomalously diffusing fluids [43], which also show power-law scaling characterised by critical exponents [44,45]. These single auto-correlations contain crucial information about the sample fluid which might be missed on typical ensemble averages [46]. Our results show that it is possible to directly access the microscopic behaviour that underpins the properties of nano-fluids, with applications in a wide variety of fields.

A. Samples
All experiments are carried out with an isotopically enriched (99.999 % 12 C) diamond sample grown by chemical vapour deposition. The diamond substrate with natural abundance of 13 C was equipped with a 99.999 % 12 C enriched homoepitaxially grown diamond film with a thickness of about 150 nm in a home-built plasma enhanced chemical vapour deposition growth reactor [34,35,47] using 99.999 % enriched 12 CH 4 gas (Cambridge Isotope Laboratories) at a concentration of 0.2 % with respect to hydrogen.
Isolated shallow NV centres were then created by ion implantation with 15 N + at a dose of 5 × 10 8 N + cm −2 , using an acceleration energy of 2 keV (correlation spectroscopy), and 2.5 keV (Qdyne). Additionally, for Qdyne measurements utilising deeper NV centres, an ion dose of 1 × 10 11 N + cm −2 at an acceleration voltage of 2.5 keV was used, followed by processing the diamond with the indirect overgrowth method according to Ref. [47], burying the NV centres deeper and thus shielding them from noise sources at the surface of the diamond. To heal resulting radiation damage, mobilise vacancies, and eventually create the desired NV centres, the diamond samples are annealed in a home-built UHV furnace at 1000 • C for 3 hours, while ensuring extremely low process pressures < 1 × 10 −7 mbar [48].
The NV centre ensembles were created by implanting 15 N + ions with a dose of 1 × 10 12 N + /cm 2 and an energy of 2.5 keV, followed by annealing as described for the single NV centres. The NV ensemble creation yield was determined to be (1.0 ± 0.1) % which corresponds to an average NV concentration of 60 ppb taking into account the depth distribution of implanted nitrogen in reference [47]. After annealing, the diamond is boiled in a 1:1:1 mixture of sulphuric (97%), perchloric (70%) and nitric acid (65%) at 200 • C in a microwave reactor system (MWT AG, type ETHOS Lab) for 30 minutes to remove any (graphitic) residues from the surface.

B. Power spectrum measurement
We perform power spectrum measurements to determine the B rms of the nuclear spins at the NV centre position [16]. This is used to determine the initial contrast of the auto-correlation signal and fix the value as in Figure 2(b) and 3(b). We consider a two-state system, which is prepared initially in state |0 , e.g. by optical pumping. The power spectrum measurement consists of the microwave pulse sequence π/2(x) pulse -dynamical decoupling -π/2(∓x) pulse, where the coordinate in parentheses indicates the axis of rotation. For two alternating measurements, which differ by the axis of rotation of the last π/2-pulse, we take their signal difference. The advantage of using the alternating sequence is reduction of the effect of unwanted laser power fluctuations.
During the DD sequence a phase Φ is accumulated due to the Larmor precession of the nuclear spins. Φ follows a Gaussian distribution where the expectation of Φ is zero. Its variance depends on the strength of the interaction between the NV centre and the nuclear spins in the sensing volume, the filter function of the applied DD sequence and the interaction time. The expected value of the phase variance in case of statistical polarisation is given by Φ 2 = Φ 2 rms (see Appendix B) where γ e is the gyromagnetic ratio for the NV electron spin, B rms is the root-mean-square of the magnetic field of the statistically polarised sample, N is the number of the DD pulses, τ is their pulse separation and ω L is the Larmor frequency of the sensed spins (in angular frequency units). This allows us to estimate B rms and the corresponding depth of the NV centre [16].

C. Correlation spectroscopy
Next, we analyse the expected measurement results with correlation spectroscopy.
Again the two-state system is considered, which is prepared initially in state |0 . We perform the correlation spectroscopy measurement, which consists of two DD blocks, separated by a waiting time T . During the waiting time, the information about the accumulated phase is mapped onto a population difference, so it cannot exceed T 1 of the system. The pulse sequence is π/2(x) − DD sequence − π/2(y) − waiting time T − π/2(x) − DD sequence − π/2(±y), where the two versions of the last pulse π/2(±y) give the alternating projections on the |0 and |1 states, respectively. The signal is the difference between these two and is given by where Φ 1 and Φ 2 are the accumulated phases during the first and the second DD sequences, c max is the maximum readout contrast (see Appendix A for details), and we can neglect the effect of decoherence and relaxation because the duration of each dynamical decoupling sequence τ N T 2 and the waiting time T T 1 . Averaging multiple measurement readouts leads to where the last approximation is valid only for small Φ k , k = 1, 2. The phases Φ k follow the same Gaussian distribution as in the power spectrum measurement, i.e. centred at zero with a variance Φ 2 rms , and we obtain (see Appendix B) where t = N τ + T is the time separation between the beginnings of the two DD sequences and ω L = 2πf L with f L again the Larmor frequency of the sensed spins (or the respective undersampling frequency). In our analysis the correlation B(t )B(t + t) ≈ B 2 rms cos (ω L t)C(t), where the correlation function C(t) can be exponential [16] or have a power-law decay for long times [24,28].

D. Quantum heterodyne detection
Lastly, we focus on the Qdyne protocol [40]. As before, the initial state of the NV centre is the |0 state. In this case, the basic building block of the protocol is a DD sequence encased in between two π/2 pulses as π/2(x) − DD sequence − π/2(y).
The salient feature of the Qdyne protocol is that the projection of the NV state through the second π/2 pulse, is followed not only by state readout/re-initialisation but also by a fixed delay time that defines the sampling frequency with which measurements are repeated, thereby recording information about the phase of the target signal, a crucial point that permits correlating all measurement outcomes in post-processing.
No alternating measurement protocol is used here. The signal for a given measurement at time t in a Qdyne measurement sequence is In post-processing of the data the auto-correlation of the recorded time trace is calculated. The correlation between any two readout pairs that are separated by t n = nT Qd , where T Qd is the sequence length of a single measurement, is then With Φ t = Φ 1 and Φ t+tn = Φ 2 the same signal as in correlation spectroscopy (5) is measured.

E. Statistical analysis
Parameter estimation is done by numerically fitting the experimental auto-correlation functions to the theoretical mode where a 0 is a signal offset, a 1 ∼ Φ 2 rms is the signal amplitude, δ is the undersampling frequency, ϕ is the signal phase, and T D is the diffusion time. In the fits of Figs. 2 and 3 the parameter a 1 ∼ Φ 2 rms is fixed and estimated from an independent power spectrum measurement. Note that we have included an artificial phase ϕ. This is necessary for correlation spectroscopy experiments as the undersampling frequency δ is typically much smaller than the actual Larmor frequency ω L , which can result in a non-zero phase, e.g., when the first sampling time is smaller than the oscillation period at δ. The reason is purely numerical for Qdyne as the fitting algorithm is prone to crashing, especially when considering the power-law correlations with Eq. (12). Then, the phase aids in adding stability and thus saving computational time. In all of the instances where we include it as a free parameter, we check that the estimation result for ϕ is either 0 or 2π within numerical precision. For the correlation function decay described by the envelope C(t/T D ) we consider either that (exponential model) or (power-law model) as in [28] with z = t/T D and C(z) = G(z) to distinguish it from other possible power-law models. The latter is shown to fall off as C(z 1) ∝ z −3/2 , for which reason we call it the power-law model.
The fitting procedure utilises a non-linear least squares algorithm with finite-difference estimation of the gradient. In correlation spectroscopy data, due to measurement vectors being indivisible, we compare goodness of fitting for each of the two correlation models, using the R 2 as a figure of merit. For each signal and model, the fitting procedure is repeated 500 times and the best fitting according to the highest R 2 is selected. The procedure is akin to the global optimisation described below for the Qdyne data.
Qdyne data analysis deserves a more in-depth explanation.
Qdyne experiments are performed sequentially, therefore, each experiment yields a single time-trace of measurements spanning the total duration of the experiment. Here, we describe the processing of the data corresponding to the experiment shown in Fig. 4(a) and (b), as the rest are analogous. In this case, the total measurement duration is 90 hours, while each single measurement requires 49.740 µs. Contrary to the case of correlation spectroscopy, where the data recording procedure did not allow to slice the time-traces, in Qdyne we can partition the data in arbitrarily short vectors. For reasons of setup stability the slices are of 15 minutes worth of measurements. In between two of these 15-minute time-traces the NV centre is optically refocused. However, the resulting auto-correlations are too noisy for adequate fitting. Thus, we compromise between eliminating noise by averaging several auto-correlations and having sufficient statistics to build a meaningful histogram, i.e. to perform Bayesian analysis. For Fig. 4(a) and (b) that means averaging 20 auto-correlations from 15-minute time-traces, resulting in a total of 18 auto-correlations with which to perform statistical analysis. Each of the 18 resulting averaged auto-correlations is fitted to the model as explained above, yielding a total of 18 frequency estimators that compose the shown histograms.
In the case of local optimisation each of the 18 auto-correlations is fitted once to the model Eq. (1), with input parameters taken from the FFT of the total auto-correlation calculated over the full 90 hours time-trace for the frequency, and amplitude and decay time taken directly from the full auto-correlation. For global optimisation, the fitting procedure is repeated 500 times, each time with different initial parameters randomly drawn from uniform distributions of size equal to the allowed searching regions for the fitting procedure (e.g. from 200 to 1800 Hz for the frequency estimator). This procedure is repeated for each of the 18 time-traces and for each case the frequency of the highest R 2 fitting results in the estimator considered for statistics. In all cases the figure of merit is the root-mean-square error of the resulting histogram of estimators.

F. Numerical simulations
Qdyne signals are generated numerically employing data from molecular dynamics simulations [24,49]. We utilise molecular dynamics simulations data of N ≈ 46 × 10 3 dipolar particles within a simulation box of size L x,y,z = (50, 50, 24) nm with an NV centre located at depths ranging from 0.3 to 12 nm. The particles within the box diffuse according to a Lennard-Jones fluid with normalised parameters = σ = 1 (where is the depth of the potential while σ is the distance at which the potential is zero), starting from an initial thermal state at temperature T = κ B T (K) [28]. The magnetic signal is calculated as the magnetic field induced by the particles at the NV position. The data thus generated proves to have no trends and the standard deviation is scale invariant, i.e. depth independent, meaning that noise data generated for one depth can be scaled through the corresponding T D to suit a different depth [24]. Therefore, we are able to use all depths, thus, increasing the statistics.
The resultant time-traces contain the power-law noise model with correlations scaling as C(t/T D 1) ∼ (t/T D ) −3/2 and thus we use them directly to generate Qdyne time traces associated with this model. Exponential correlations are simulated by fitting the molecular dynamics data to Orstein-Uhlenbeck noise which process exponential decay, and using such fitting as the noise source. for experiment control and data processing, SoftwareX 6, 85 (2017).

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code used for obtaining the presented numerical results is available from the corresponding author upon reasonable request. We consider a two-state system, which is prepared initially in state |0 , e.g., by optical pumping. In all measurement schemes, we find the population of the NV centre via fluorescence collection during laser excitation, which relates the average number of detected photons to the population of the system because the fluorescence rate is higher for the NV centre being in the spin state |0 . To reduce laser power fluctuations, the photoluminescence (PL) of the NV centres at the first few hundred nanoseconds (signal window) of the laser pulse are normalised with the reference PL at the end of the optical pumping period (reference window; usually last 2 µs of a 5 µs long laser pulse). The signal window is adjusted such that highest SNR can be achieved. The fluorescence signal for readout of the state |i (i = 0, 1) is F i = ηi η ref where η i is the photon count rate during the signal window if the NV centre is in state |i and η ref is the photon count rate during the reference window that is independent of the spin state after long enough optical pumping. The maximum contrast is then In the experiments usually a maximum contrast c max ≈ 30 % is obtained.
We describe now the power spectrum measurement, which consists of the microwave pulse sequence π/2(x) pulse -dynamical decoupling -π/2(−x) pulse. This sequence keeps the system in state |0 in case of no phase accumulation during dynamical decoupling. The density matrix after the sequence is given by where Φ is the accumulated phase during the DD sequence and the population of state |0 and |1 is given by p 0 = cos Φ 2 2 and p 1 = sin Φ 2 2 , respectively.
Performing the alternating sequence π/2(x) pulsedynamical decoupling -π/2(x), which transforms our system to state |1 when Φ = 0, effectively exchanges the resulting populations p 0 and p 1 . The respective fluorescence signals of the two measurements are then where F 0 ( F 1 ) is the expected fluorescence after the first (alternative) sequence. Finally, we obtain the contrast for the power spectrum measurement where we used that cos (Φ) = cos Φ 2 2 − sin Φ 2 2 .
Multiple readouts are performed and we average the result, which leads to c ps = c max cos (Φ) . (A4) Assuming that the phase Φ follows a Gaussian distribution, centred at zero with a variance Φ 2 , we obtain where the last approximation is valid only for small Φ 2 . The variance of the phase Φ depends on the strength of the interaction between the NV centre and the nuclear spins in the sensing volume, the filter function of the applied DD sequence, and the interaction time. The expected value of the phase variance in case of statistical polarisation is given by where γ e is the gyromagnetic ratio for the NV electron spin, B rms is the root-mean-square of the magnetic field of the statistically polarised sample, N is the number of the DD pulses, τ is their pulse separation, and ω L is the Larmor frequency of the sensed spins (in angular frequency units). We note that we assumed that the coherence time of the sensed signal is much longer than the sensing time during the DD sequence, so B(t )B(t + t) ≈ B 2 rms cos (ω L t). If this assumption is not feasible the spectrum of the signal from the nuclear spins can be taken into account [16,24]. However, this effect is typically not large for our power spectrum experiments, so we will neglect in the present analysis. If the DD sequence is applied on resonance, i.e., τ = π/ω L , the variance of the accumulated phase takes the form The formulas for the phase Φ rms in Eqs. (A6) and (A7) allow us to estimate it by a power spectrum measurement, where we vary τ and keep N constant (see Fig. 6) or DD order scan measurement where we keep τ = π/ω L and vary the number of DD pulses N , respectively. The formula for the expected contrast of the power spectrum is given by where T 2 is the coherence time of the respective DD sequence. We have neglected the very small effect of T 1 decay during the usually very short duration of the DD sequence in comparison to the T 1 time. We note that DD is not efficient for compensation of high frequency noise and pulse errors might also lead to a loss of the contrast c max even without an external signal, i.e., for Φ = 0. )UHTXHQF\0+]
The experimental data points are fitted with a Lorentzian fit (blue line).

Correlation spectroscopy
We measure the population of the NV centre in the correlation spectroscopy experiments analogously to the power spectrum measurement. The pulse sequence is illustrated in Figure 7(a). We again apply two alternating measurements, where the spin state of the NV centre is projected on either |0 or |1 with a π/2-pulse around the y-or −y-axis (see Methods).
In the following we describe in more detail the state of the system at each step of the correlation spectroscopy experiment. We assume that the system is initially prepared in state |0 by optical pumping. Then, the density matrix after the first correlation spectroscopy sequence π/2(x) − DD sequence − π/2(y) is where Φ 1 is the accumulated phase during the first DD sequence and we assumed negligible relaxation during the sequence. The off-diagonal elements become zero due to decoherence during the waiting time T T * 2 and the density matrix afterwards is Finally, the density matrix after the second correlation spectroscopy sequence π/2(x) − DD sequence − π/2(y) is with the accumulated phase Φ 2 during the second DD sequence and we again assumed negligible relaxation during the sequence. The population of state |0 is then p 0 = 1 2 (1 + sin (Φ 1 ) sin (Φ 2 )) and the population of state |1 is p 1 = 1 2 (1 − sin (Φ 1 ) sin (Φ 2 )). The populations of the two states are exchanged with the alternative sequence π/2(x) − DD sequence − π/2(−y). Again, the measurement signal is the difference between them where the last approximation is valid for small accumulated phases during each DD sequence, which is typically the case. We perform multiple readouts and average the result, which leads to with the envelope C(t) of the correlations (see B for details).

Qdyne
In the Qdyne measurements the NV centre's PL is not normalised owing to the fact that the individual readouts are not averaged but stored individually. The photon count rates of our setup were of around η 0 = 0.04 and η 1 = 0.03, i.e. in 100 readouts in average four (three) photons are detected if the spin state is |0 (|1 ). Again, the readout window is optimised to achieve the highest SNR and all other photons during the laser pulse that are not within the readout window are discarded.
The measurement sequence (see Figure 7(b)) is similar to the first step of correlation spectroscopy, so the density matrix is the same as in (A9) with the accumulated phase Φ during dynamical decoupling. The population oscillates around 1/2 for which reason the signal in Eq. (8) of the main text is measured.

Appendix B: Modelling statistical polarisation
The signal we observe is due to statistical polarisation in the detection volume of the NV centre(s). We assume that the magnetic field, generated by the sample due to statistical polarisation is given by [50] where A, D ∼ N (0, B rms ) the normal distribution with zero mean and standard deviation B rms , so the variance of the magnetic field is B(t) 2 = B 2 rms . When we apply a DD sequence on resonance, i.e., by instantaneous π pulses at times π(1 + 2k)/ω, k = 0 . . . n we can sense only A cos (ωt). The accumulated phase takes the form [50] Figure 7. Pulse sequences for the used protocols. (a) Correlation spectroscopy pulse sequence. The NV centres are initialised into ms = 0 state using a green laser. The MW pulse sequence consists of two DD sequences with equal length N τ where N is the total number of π-pulses and τ = 1 2f L is the spacing between two subsequent π-pulses where fL is the Larmor frequency of the nuclei. The dynamical decoupling sequence is sandwiched between two π/2-pulses of orthogonal phases. The final spin state is optically measured using the green laser where the signal counts (first red window) are normalised with the reference counts (second red window). The measurement is averaged over n sweeps.(b) Qdyne pulse sequence. The measurement sequence with a single DD block is repeated continuously. The measurement result during each sweep is stored individually and the sequence synchronised before next sweep starts. The photon counts (red window) are not normalised.
where f (t) is a modulation function due to the DD pulses and N τ is the duration of the DD sequence with N the number of pulses and τ = π/ω the time between the centres of the pulses. Since A ∼ N (0, B rms ), we obtain Φ = 0, Next, we calculate the expected value of Φ 1 Φ 2 of two phases Φ 1 and Φ 2 , accumulated by two DD sequences, starting at times t 1 and t 2 = t 1 + t, where t is the difference between the starting times of the two DD sequences. First, it is useful to state that we assume now that A and D can be time dependent as the time separation t between the starting times of the two DD sequences can be much longer than the diffusion time of the sample, which results in a change in A and D between the two sequences. However, we also assume for simplicity that their values do not change during a DD sequence (N τ T D ), thus they depend only on the starting time of the respective sequence. Finally, we also assume The correlation function C(t 2 − t 1 ) = C(t) of the magnetic field envelope can have a different shape that we probe, e.g., exponential decay [16] or a power-law decay at long times [28]. We obtain B(t 1 )B(t 2 ) = A(t 1 )A(t 2 ) cos (ωt 1 ) cos (ω(t 1 + t)) + D(t 1 )D(t 2 ) sin (ωt 1 ) sin (ω(t 1 + t)) where we used A(t k )D(t l ) = 0 in the first equality. The accumulated phases during the first DD sequence then take the form We took t 1 = 0 for simplicity of presentation and without loss of generality. For the second DD sequence where N τ is again the duration of the DD sequence. Then, the expected value of Φ 1 Φ 2 takes the form where we used A(0)D(t) = 0 in the first equality.
Appendix C: Power spectrum shape analysis In the following, we derive the asymptotic behaviour of the power spectrum of the magnetic noise induced on the NV centre by diffusing particles in three spatial dimensions. This derivation is based on the supplementary material of [28].
We assume that the nuclei obey the diffusion equation. The diffusion is therefore described by the equation where P ( r, r 0 , t) is the conditional probability that a particle is found at position r at time t if it was initially at position r 0 at time t = 0, and δ is the Dirac delta function. The power spectrum is the Fourier transform of the correlation Since the time dependence enters only through the conditional probability P ( r, r 0 , t) taking the Fourier transform of (C2) yields At ω = 0 the solution to Eq. (C4) is equivalent to that of the electrostatic potential of a point charge Substituting the conditional probability (C5) into Eq. (C3) leads to where the lower bound for the integration is determined by the depth d of the NV centre. The only length scale in the problem is the NV's depth, therefore, by dimensional analysis of Eq. (C6), where C 1 is a non-universal constant that depends on the geometry of the problem. The first correction in frequencies to the spectrum Eq. (C7), in the limit ω ω D , where the diffusion time T D is the characteristic time to diffuse at a distance the characteristic angular frequency of diffusion noise, can be accounted for by considering the finite diffusion length scale l ω = 2πD/ω d. This creates an effective cut-off to the interaction, such that Eq. (C6) is corrected to where again C 2 is a non-universal constant. This reproduces the sharp-peaked behaviour around the Fourier peak of the power spectrum shown in Fig. 2(d) on the main text. The correction can be understood intuitively by considering that sensing an oscillating signal at a small angular frequency ω (or two signals with an angular frequency difference ω) typically requires some non-zero correlation C(t) at time t = t ω = 2π/ω, so we could detect at least one oscillation cycle of the signal. In the regime t ω t D , or equivalently ω ω D , this leads to a minimum time t ω where any oscillating component of the auto-correlation function at ω after this time contributes to a lower value of S(ω) in comparison to S(ω = 0) and thus a narrower linewidth. Spatially, this means that there is an effective cut-off length l ω = 2πD/ω d, where the spectral density at ω, i.e., S(ω) is smaller than S(ω = 0) due to an oscillating signal at ω coming from nuclei which diffuse more than the cut-off interaction length l ω .
In the other limit of large frequencies, e.g. ω ω D , which is equivalent to l ω d, we approximate D∇ 2 ≈ −D/d 2 ∼ −ω D as the small correction to the large frequency ω. Eq. (C4) then takes the form Hence, the spectrum, Eq. (C3), approaches (C10) Eq. (C10) recovers the well known Lorentzian behaviour of the spectrum.
Note that whether one considers the angular dependence and the geometry changes the prefactors but not the overall behaviour of the correlation function or the spectrum shape.

Appendix D: Additional experimental details
We measure the auto-correlation function of the detected signal by three different approaches: correlation spectroscopy with a single shallow NV centre, correlation spectroscopy with an ensemble of NV centres, and Qdyne with single NV centres. Next, we describe more details on the experimental setup and measurements results for each of these.

Experimental setup for single NV centre experiments
Experiments for the correlation spectroscopy (CS) and Qdyne measurements have been done on two different but conceptually equivalent setups. The shallow NV centres are addressed via a fluorescent confocal scan done  Figure 8. Experimental data from correlation spectroscopy measurement with a single NV centre. (a) Signal of the correlation spectroscopy measurement vs. time t between the beginnings of the two dynamical decoupling sequences in correlation spectroscopy (black dots) and fits of the exponential (blue line) and power-law models (orange line) with their goodness of fit. The power law model shows a better fit to the data, especially at long times. There is also a significant difference in the estimate of the initial amplitude of the oscillation a1. (b) FFT of the experimental data, FFT from the fitted data to the exponential and the power-law decay models (in (a)) is shown for illustrative purposes. The results show that the power law model allows for a more precise frequency estimation. (c) Optimum R 2 for different fixed values of the initial amplitude a1 of the correlation spectroscopy signal for the single NV experiment. The blue (orange) curves shows the best fits of the data to the function a0 + a1 cos (ωt + φ)C(t) for the exponential (power-law) decay model while the green lines show the estimated values of the initial amplitude a1 = 0.12, obtained from independent power spectrum measurements. The peak value of the R 2 curve is higher for the power-law model in comparison to the exponential decay model. In addition, the power-law model is substantially better when the initial amplitude is fixed to the estimate from the independent power spectrum measurement, i.e., comparing the R 2 of the two models for a1 on the green lines. The fits with fixed a1 and their FFT are shown in Figure 2 of the main text.
on a home-built confocal microscope. The NV centres can be initialised and read out using a 532 nm (CS) and 517 nm (Qdyne) laser pulse of about 5000 ns (CS) and 1000 ns (Qdyne) duration generated by a CW laser (Laser Quantum Gem532, Toptica iBeam smart 515) and chopped by an acousto-optic modulator (CS, Crystal Technology 3200-146) or by the pulsed-laser output itself (Qdyne

Correlation spectroscopy with a single NV centre
The NV centre for the correlation spectroscopy measurements is located ≈ 2.9 nm below the diamond surface. The depth was estimated by measuring the power spectrum of the hydrogen nuclei in the immersion oil with dynamical decoupling (see Fig. 6) [30]. We used the Knill dynamical decoupling (KDD4) pulse sequence for our experiments [36][37][38][39]. KDD has five π pulses with phases (χ + π/6, χ, χ + π/2, χ, χ + π/6) [36][37][38][39]. It is based on a composite pulse, first proposed by Tycko and Pines [52], which was recently shown to be one example from a set of universal composite pulses for population inversion [53]. We obtain KDD4 by nesting KDD in the XY4 sequence [37][38][39], i.e., the phase χ takes values (0, π/2, 0, π/2), resulting in a sequence of twenty pulses. We repeat it two times, which we label as KDD4-2, and vary the pulse time separation [36][37][38][39]. The power spectrum measurement allows estimation of B rms , so we can use it as an input parameter for the data analysis of the correlation spectroscopy measurements in Fig. 2. The T 2 time with the same DD sequence is measured 110 µs by varying the DD order where the spacing between the two adjacent π pulses was set off-resonant to 1/(f L ×1.3). The T * 2 is measured to 200 ns using Ramsey sequence. We note that the translation diffusion correlation time for the nuclear spins is given by T D = d 2 /D oil , where d is the depth of the NV centre and D oil = 5 × 10 −13 m 2 /s is the diffusion coefficient of the immersion oil. The spin lattice relaxation time, T 1 of the NV centre is 1.11 ms. Figure 8 includes additional data analysis of the experimental results from the single NV centre experiment. Figure 8(a) shows that the power-law decay model has a better fit to the data with an R 2 ≈ 0.96 in comparison to R 2 ≈ 0.93 for the exponential model. Its goodness of fit is better especially at long times, e.g., between 50 − 100 µs, which correspond to more than three times the expected diffusion time T D ≈ 17 µs. In contrast to Fig. 2 in the main text, the initial contrast is considered here a free parameter. We note that the  Figure 9. Experimental data from correlation spectroscopy measurement with an ensemble of NV centres. (a) Signal of the correlation spectroscopy measurement vs. time t between the beginnings of the two dynamical decoupling sequences in correlation spectroscopy (black dots) and fits of the exponential (blue line) and power-law models (orange line) with their goodness of fit. The power law model shows a similar overall fit to the data as the exponential model, nonetheless fitting the data better at long times. There is a significant difference in the estimate of the initial amplitude of the oscillation a1.
(b) FFT of the experimental data, FFT from the fitted data to the exponential and the power-law decay models (in (a)) is shown for illustrative purposes. The results show that the power law model allows for a more precise frequency estimation. (c) Optimum R 2 for different fixed values of the initial amplitude a1 of the correlation spectroscopy signal for the single NV experiment. The blue (orange) curves shows the best fits of the data to the function a0 + a1 cos (ωt + φ)C(t) for the exponential (power-law) decay model while the green lines show the estimated values of the initial amplitude a1 = 0.079, obtained from independent power spectrum measurements. The peak value of the R 2 curve is higher for the power-law model in comparison to the exponential decay model. In addition, the power-law model is substantially better when the initial amplitude is fixed to the estimate from the independent power spectrum measurement, i.e., comparing the R 2 of the two models for a1 on the green lines. The fits with maximum R 2 for either model are shown in (a) in this figure and the fits for fixed a1 = 0.079 is shown in the Fig. 3(a) in the main text.
two models give almost a factor of three difference in its estimate (see Fig. 8(c)). Figure 8(b) shows the FFT of the experimental data, and the FFT from the fitted data to the exponential and the power-law decay models (in (a)) for illustrative purposes. Figure 8(c) shows the estimation of the initial amplitude of the auto-correlation function from the correlation spectroscopy data of the single NV centre. The data analysis shows that the overall fit of the power-law model is better, i.e., the peak value of the R 2 curve is higher, in comparison to the exponential decay model. However, the optimal values of the initial correlation amplitude a 1 for the two models differs substantially. Because a 1 ∝ Φ rms ∝ B rms we can estimate it independently with a power spectrum measurement and use the estimate (a 1 = 0.12) as a fixed parameter in the fitting procedure. The result also shows that the power law decay model performs especially well at long times (see Fig. 2 in the main text).

Correlation spectroscopy with an NV ensemble
The NV ensemble is optically initialised and read out via a self-made widefield setup.
In order to increase the collection efficiency, a solid immersion lens with a diameter of 6 mm is placed between the diamond and the objective lens (NA: 0.63, effective focal length: 20 mm, Edmund Optics 85-301). The NV centres are excited by a diode-pumped solid-state laser at 532 nm. An acousto-optical modulator (AOM) (Crystal Technology 3200-146) is used to generate laser pulses.
The fluorescence is recorded by a silicon photomultiplier (Excelitas LynX-A-33-A50-T1-A) and measured by a digital oscilloscope (Adlink PXIe-9834). To eliminate back scattered light from the excitation beam, a bandpass filter is placed in front of the detector. Microwave (MW) pulses are first generated by an arbitrary waveform generator (Keysight M3202A) and then amplified. The MW pulses are delivered to the NV sensor through a copper coil placed on top of the diamond. The magnetic field of about 920 G is generated by an electromagnet and the magnetic field direction is aligned to one of the NV symmetry axis. Figure 9 includes additional data analysis of the experimental results from the NV centre ensemble experiment. Figure 9(a) shows that the power-law decay model has a similar fit to the data with an R 2 ≈ 0.96. Its goodness of fit is better especially at long times, e.g., between 200 − 400 µs, which correspond to more than three times the expected diffusion time. In contrast to Fig. 3 in the main text, the initial contrast is considered a free parameter, which is the reason behind the exponential model showing very similar fit at long times as the power-law model. Since the aim is to demonstrate the power-law model, our sampling is optimised towards improving the signal-to-noise ratio for long times, where the difference between models is more apparent. But this means that the sampled data does not allow us to efficiently estimate the signal amplitude a 1 at very short times (there are very few data points and the a 1 confidence interval is quite large). As a result, the fit of unconstrained exponential model tends to underestimate a 1 , allowing for an artificially increased signal at long times, similar to the power-law model. We note that the two models give almost a factor of two difference in its estimate (see Fig. 8(c)). Figure 9(b) shows the FFT of the experimental data and the FFT from the fitted data to the exponential and the power-law decay models (in (a)) for illustrative purposes. Figure 8(c) includes additional results for the estimation of the initial amplitude of the auto-correlation function from the correlation spectroscopy data of the ensemble measurements. Similarly to the single NV experiments, the data analysis shows that the overall fit of the power-law model is better, i.e., the peak value of the R 2 curve is higher, in comparison to the exponential decay model. The optimal values of the initial correlation amplitude a 1 for the two models differs substantially, similarly to the single NV experiment. Because a 1 ∝ Φ rms ∝ B rms we can estimate it independently with a power spectrum measurement and use the estimate (a 1 ≈ 0.008) as a fixed parameter in the fitting procedure. The result shows that the power law decay model performs especially well at long times (see Fig. 3 in the main text).

Qdyne with single NV centres a. Experimental description
For the Qdyne measurements different NV centres located a few nanometers between 8 and 15 nm below the surface of the diamond are used. The auto-correlation of the measured time trace is calculated using the statstool package in Python. It typically shows a decay on a time scale of a few ten milliseconds ( Fig. 10(a)), which probably results from fluorescing particles at low concentration that occasionally diffuse through the laser beam bath. To account for this decrease an exponential decay is added to the model function and the data is corrected by subtracting the exponentially decaying part. Figure 10(b) and (c) show the auto-correlation and the fitted model functions of the raw and the corrected auto-correlation data. The underlying frequency of the oscillation can be calculated as the undersampling frequency from the nuclear Larmor frequency f L and the sequence length T Qd as outlined in [40].
The measurement sequence, as illustrated in Figure 7(b), consists of optical initialisation, that at the same time acts as readout, and subsequent dynamical decoupling. For the dynamical decoupling a XY8 protocol is embedded in two π/2-pulses with orthogonal phases. The order of the XY8 block is set between 8 and 24 depending on the depth of the NV centre.
Details to the measurement parameters are shown in I. The results of measurement number 1 are presented and discussed in the main text and Figure 4. All other measurements are analysed in the same fashion where the results are shown in Figure 5. For the measurements number 1 to 5 fluorescence-free immersion oil from  Figure 10. Auto-correlation calculation with the Qdyne data. Correlations show exponential decay at the millisecond scale (a). That is removed by fitting an exponential together with the model function (b). In (c) the fit of the exponential and power-law models are shown for the corrected data. FFT of the data and the fitted models for exponential (power-law) models in blue (orange) in (d) with inset showing the peak region.
Fluka is used and for the last one the high-viscosity perfluoropolyether.
b. Additional results Figure 11 shows the complete histogram analysis for experiment one in I. Here we include both the global (upper row) and the local (lower row) optimisation results, both without (left column) and with a fixed phase ϕ (right column) for the non-linear least squares fitting. Subfigures (a) and (b) are also presented in Fig. 5 in the main text). As can be seen, only the power-law model does result in meaningful histograms in all instances, with root-mean-square errors ranging from 292 Hz in the case of global optimisation without fixing the phase Fig. 11(a) to 151 Hz for local optimisation with fixed phase in Fig. 11(d), with a clear trend of diminishing rmse as the fitting model becomes more accurate. In contrast, the exponential correlations model results in wider histograms. The smallest rmse is obtained in  Figure 11.
Histograms of frequency estimators from non-linear least squares fitting of 18 auto-correlations from experimental Qdyne time-traces. In orange, fitting to a correlation with power-law decay C(t/TD 1) ∼ (t/TD) −3/2 while in blue the fitting is to a correlation with exponential decay. (a) corresponds to a global optimization with 500 fittings per auto-correlation, with the best result decided by the highest R 2 and with ϕ a free fitting parameter. In (b) the phase ϕ = 0 is fixed. (c) displays the histogram for local optimization with fixed initial seeding frequency as taken from FFT of the original full time trace, while (d) does the same but with fixed ϕ = 0.
fixed phase does the exponential correlations model show meaningful histograms, with rmse, respectively 365 and 350 Hz for global and local optimisations. Note that the flat histogram limit lies at 404 Hz.
In Fig. 12 we consider the rmse ratio between global optimisation for power-law and exponential fittings, both for all Qdyne experiments as well as numerical simulations of time-traces generated with either of the correlation models and analysed simultaneously with both of them. We observe an experimental mean ratio of 0.81 ± 0.21, which again favours the interpretation of power-law decay of correlations, since it is this model that better fits the data. Numerical simulations are consistent with this analysis, yielding a ratio of 0.89 ± 0.16 for correlations generated according to the power-law model while this ratio is 1.12 ± 0.17 for correlations generated with the exponential model.  main text show the FFT of the experimental data together with the fits to the exponential and power-law decay models. We see that the latter provides a better fit to the data but nonetheless the experimental FFT has a broadening not directly accounted for by the model. In this section, we hypothesise on the origin of that extra broadening by studying a combined power-law and exponential decay model. • G(t/TD) Figure 12. Root-mean-square error ratio between global optimisation for power-law and exponential models fittings. Each dot is calculated with the rmse of histograms as in Fig. 11. In black squares, experimental results. Orange circles display the ratio for data generated with power-law decay while blue stars correspond to data generated with exponential decay. Solid lines show the mean ratio for all dots while shaded areas correspond to the standard deviation for each mean. Note that experiment six in I is not shown but it is considered for the mean ratio.
For a diffusion dominated noise, we expect the main contribution to the correlation function decay to be a power-law, as we have demonstrated in the main text. Other noise sources such as spin diffusion noise by nuclear dipolar interactions are already accounted for by the T 2 of the NV centre and, provided measurements do not get close to that time, as happens in our experiments, they do not contribute significantly. Yet there are possible contributions from a magnetic gradient of the NV centre (≈ 1 Gauss for correlation spectroscopy with 2.9 nm deep NV and ≈ 0.1 Gauss for Qdyne), magnetic instabilities, spatial inhomogeneities, or back-action (= γBN τ ≈ 0.1 < 1) in the NV centre, that can contribute when measurement times are sufficiently long. These various sources of noise cause exponential decay of the correlation function at a much slower rate than diffusion, and can therefore contribute to a slight Lorentzian broadening of the sharp-peak from the power-law decay.
To account for these extra sources, we consider here a model which includes the correlation function decay G(z) from diffusion noise (Eq. (12) in the main text), with an additional exponential decay at a much lower rate, such that C(z) = G(z) exp (−t/T E ), with T E T D . In Fig. 13, we show the FFT for such a mixed model, with T E = 50T D , together with the power-law decay model and the exponential model both with just with T D decay rate. We can see that such a mixed model reproduces well the observed extra broadening from the experiments. Note that considering each possible noise source would require studying the particular impact that each of them would have on the correlation function, which need not  Figure 13. FFT of the correlation function with exponential, power-law, and combined exponential and power-law decays. The first two have a decay rate TD while for the combined model the power-law decay is characteristic time is TD and the exponential decay characteristic time is TE = 50TD.
necessarily be exponential, and estimating their own characteristic time. Such an analysis is however beyond the scope of this work.
Appendix F: Frequency sensitivity for the power-law and exponential models Frequency resolution in a power spectrum is defined as the ability to differentiate two close spectral lines. When a physical process characterised by a frequency is affected by noise, the corresponding spectral line broadens, such that two lines closer than their characteristic linewidth resemble a single, broad spectral line [50]. In liquid state nano-NMR this broadening is mainly caused by diffusion. The origin of the spectral line lies in the time-correlated magnetic field generated by a statistically polarised nuclei sample. When such a signal is noisy, correlations decay, limiting the amount of information that can be extracted about the corresponding frequency. We quantify this information by means of the Fisher Information of classical parameter estimation, such that for each specific measurement protocol, and for each model considered for correlations, we have a particular Fisher Information which measures the ability to estimate the frequency, as shown in [24].
We can define frequency sensitivity for a specific protocol as the error on the frequency estimation divided by the square root of the total measurement time η = ∆δ √ Ttot . Theoretically, the frequency estimation error is bounded by the inverse square root of the total Fisher Information about the frequency in a given experiment ∆δ ≥ Thus, considering equal experimental parameters, we can compare the frequency sensitivity in light of each correlations model, i.e. exponential or power-law, through their respective Fisher Informations.
Considering the ratio of sensitivity in the exponential correlations η exp to power-law η pl yields, for correlation spectroscopy ηexp η pl ∼ 1 δT D , which diverges for small frequencies δ, meaning that having power-law correlations improves the frequency sensitivity, although globally, for correlation spectroscopy, a resolution problem still exists. For Qdyne, the ratio is ηexp η pl ∼ log (δTtot) δ 2 T 2 D , which not only diverges for small frequencies, but grows with the total measurement time. Thus, power-law correlations significantly improve frequency sensitivity and in principle allow to vanquish the resolution problem [24]. For different experimental protocols, e.g. between correlation spectroscopy and Qdyne, the photon shot noise has an impact on which protocol is optimal, as correlation spectroscopy is less affected. Here as well power-law correlations lessen the impact that such a noise has, as they provide more information for equal total measurement time.
We can as well estimate the impact on frequency sensitivity from additional noise sources which cause exponential decay at a much slower rate than diffusion. To do so, we calculate the Fisher Information for our combined model in which diffusion induces a power-law decay while other noise sources induce exponential decay of correlations with a characteristic time T E T D . In this case, for Qdyne measurements with Ei(x) the exponential integral function, which goes to 0 for small T E → 0 and yields log (δT tot ) for large T E → ∞. For frequencies such that δT D < 1 < δT E , provided the measurement time is sufficiently long, there is no significant decrease in sensitivity due to these additional noise sources. For frequencies δT E < 1 the resolution problem reappears and lim δ→0 = 0. However, note that these additional noise sources contributing to T E , such as, for example, magnetic field instabilities or magnetic defects, contrary to what happens with diffusion, can be mitigated through tweaking the setups, and therefore the resolution problem they lead to is not fundamental in the sense of the diffusion resolution problem which, in the case of exponentially decaying correlations, can only be avoided altering the sample itself, which would forfeit the purpose of nano-NMR. It is expected moreover, that future improvement on setups and experimental techniques leads to the ability to vanquish these noise sources.