Nano-NMR based flow meter

D. Cohen,1 R. Nigmatullin,2 O. Kenneth,3 F. Jelezko,4 M. Khodas,1 and A. Retzker1 1Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Givat Ram, Israel 2Complex Systems Research Group, Faculty of Engineering and IT, The University of Sydney, Sydney, New South Wales 2006, Australia 3Dept. of Physics, Technion, Israel 4Institute for Quantum Optics, Ulm University, Albert-Einstein-Allee 11, Ulm 89081, Germany (Dated: May 28, 2019)

Introduction -In the past few decades there has been a surge in the use of microfludic devices. They are now a wellestablished platform, which has revolutionized bio-medical research, primarily as a result of the Lab on a Chip approach to chemical and biological analysis [1][2][3][4][5][6][7]. Microfludic channels present rich physical phenomena which depend on a large number of parameters (geometry, fluid type, fabricated material, etc.). Understanding the physics of these channels is key to the future development of the field, not only for the fabrication of such devices, but also because their research applications demand both high temporal and spatial measurement resolution, which require an in-depth understanding of the flow characteristics [7,8]. One of the main interests is the flow profile near the boundary of microfludic channels, since these boundary effects will play a major role when attempting to downsize these structures into the nanoscale regime. In recent years accumulating evidence has shown that the flow in such channels does not always obey the no-slip boundary condition [8][9][10], which is commonly used to solve the Navier-Stokes equations. Though there have been advances in the microscopic theory of this phenomenon [10], an experimental method to accurately measure the velocity near the surface has not yet been developed.
Current methods for velocity measurement in microfludic channels are based mostly on fluorescent beads or dyes. Typically, fluorescent molecules are injected into the channel and their propagation is tracked using a confocal microscope [11][12][13][14][15][16][17]. Alternatively, the velocity can be determined by light scattering from a periodic array [15]. Using fluorescent molecules can also be used to measure the diffusion coefficient of these molecules within a liquid [16]. These methods lack in performance and only achieve an accuracy of about 5% for the mean channel velocity due to technical issues which include laser beam focusing and the temporal resolution of the camera. Though progress in technology may help overcome these limitations in the future, these methods suffer from fundamental issues as well. First, the fluorescent molecules are often large, and in a small channel they can affect the flow profile. Second, the flow trajectories of these molecules generally pass through the middle of the tube (in a Poiseuille settings) because of the velocity gradient, which supplies little information about the flow profile near the surface. Finally, current methods can only be used to measure the diffusion coefficient of the fluorescent molecule within the liquid, whereas the interesting parameters are usually the self-diffusion and rotational diffusion coefficients of the liquid.
Here, we propose a non-intrusive measurement technique for the mean drift velocity and the self-diffusion coefficient, which is sensitive to surface effects, based on nano Nuclear Magnetic Resonance (nano-NMR). Nano-NMR is an emerging field, that attempts to adapt classic NMR techniques into the nano scale [18][19][20][21]. The essential difference is that the core assumption of a macroscopic number of molecules is no longer valid, and as a result the NMR signal is below the signal-to-noise ratio of classic NMR measurement devices. Recent experiments have demonstrated that Nitrogen-Vacancy (NV) centers are a promising platform for nano-NMR [18][19][20][21][22][23][24][25][26]. The main idea behind this method is the fact that the NV center is an excellent magnetometer at the nano scale that can read the magnetic field created by the nuclear spins effectively and thus replace the role of the coils in the regular NMR setting.
In this paper, we show that a nano-NMR approach for velocity estimation in microfludic channels outperforms current fluorescent based methods. The main idea is sketched in fig.1, where a flow through a microfluidic channel is sensed by an NV center ensemble via a similar setup to the one in [27]. The unpolarized spins motion induces a random magnetic field at the location of the NV centers. The power spectrum of the magnetic field noise can be estimated by optically probing the NV center. By analyzing the noise characteristics the flow properties can be deduced. We shall consider the NV center as a two level system (spin 1 2 ) with an energy gap of γ e B ext , where B ext is the externally applied magnetic field [28].
The power spectrum generated by the motion can be read by a decoherence spectroscopy method [29][30][31][32][33][34][35] or spin relaxation [36,37]. These methods exploit the fact that the transition rate or relaxation time in a two-level system is propor- Figure 1. The scheme-The nuclear spins (red arrows), which are randomly oriented, generate a random magnetic field on an ensemble of NV centers (blue arrows), which are located next to the surface of the diamond. The NVs sense the magnetic field, whose randomness is amplified due to the flow. By optically measuring the state of the NV centers, the power spectrum of the magnetic noise can be estimated, from which the velocity can be deduced.
tional to the spectral density evaluated at the transition frequency or the Rabi frequency of an external drive. In both cases, the relaxation rate is given by Γ = S(ω), where S(ω) is the power spectrum at frequency ω, which is either the external drive's Rabi frequency or the energy gap of the two level system. The relaxation rate and, therefore, the power spectrum can be efficiently estimated from the population decay.
Measurement scheme -Since all the effects of the velocity emerge from the magnetic noise induced on the NV, polarization is not essential here and it can be measured in unpolarized samples in a way that is similar to diffusion measurements [38]. ! < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > S(!) < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > ! = 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > Larmor frequency Dynamics information is in the width Figure 2. A diagram of the total spectrum that can be probed by the two terms. The T 0 term (blue) probes the zero frequency part of the spectrum, and the T ±1 terms (orange) probe the Larmor frequency part of the spectrum. The fluid drift velocity can be determined from the structure of the spectrum.
The different components of the magnetic field can be sensed by probing the different parts of the dipole-dipole interaction separately. We denote S i (I i ) as the spin operator of the NV (nuclear spin) at the i direction and S ± (I ± ) as the NV's (nuclear spin's) raising/lowering operators. The velocity can be sensed by the T 0 = 3 cos(θ ) 2 − 1 S z I z term or the T ±1 = 3 2 sin θ cos θ e ∓iφ (S z I ± + S ± I z ) terms (see for example [39]) as shown in fig. 2. The T 0 is the classical term, deriving from the thermal polarization fluctuations of nuclei within a volume of d 3 , which inflict a random field on the NV. In the case of a polarized fluid, the noise is decreased and this term vanishes. However, probing the dynamics is still possible if the nuclei are polarized in the x − y plane [40].
The second term, T ±1 , is a rotating term which manifests itself at high frequencies and is routinely used to probe the nano-NMR spectrum [18][19][20][21]. In order to probe this term, a rotation of the NV at the Larmor frequency must be introduced. The method of choice is to introduce driving by pulsed dynamical decoupling [20][21][22][23][24][25][26], which yields the following effective Hamiltonian: H 1 = 3 2 sin θ cos θ S z (I x cos(δt + φ ) + I y sin(δt + φ )) , where δ is the frequency difference between the Larmor frequency and the frequency of the pulses.
Sensitivity estimation -The velocity precision can be understood in a simplified model. In this nano-NMR setting there are two characteristic time scales: the diffusion time scale τ D = d 2 D and the drift time scale τ v = d v , where d is the NV's depth below the diamond surface, D is the self-diffusion coefficient and v is the mean drift velocity. These time scales control the correlation time of the NMR signal if only one of the processes is dominant. Respectively, they are the inverse characteristic widths of the NMR spectra. In the case when both processes are significant, the natural expectation is that the characteristic width is In order to estimate the accuracy we follow the commonly accepted assumption that the nano-NMR spectrum is a Lorentzian function [18,38], S(ω) = σ 2 S(ω=0) ω 2 +σ 2 . The peak is given by S(ω = 0) = γ 2 e B 2 RMS σ , and B 2 RMS = n µ 0h γ N 4π 2 1 d 3 , where γ e = 2π × 28 GHz T −1 is the electron gyromagnetic ratio, γ N is the nuclei gyro-magnetic ratio and n is the nuclei number density [38], which results in a B RMS of 350kHz/γ e for a 5 nm deep NV. The sensitivity is best described by the dimensionless quantities: The first denominator can be interpreted as a unit power spectrum taken as S(ω = 0, v = 0), whereas the second denominator is a unit velocity for which ω v = ω D , i.e.,ṽ = ω v ω D . We use these units henceforth in our estimations.
The decay rate Γ of the NV can be estimated with an accuracy of ∼ Γ for a single shot measurement within one experiment which lasts for t ∼ Γ −1 . Repeating this measurement for a total duration of T results in an accuracy of ∆Γ = Γ T . Since most NV measurements are not single shot, the accuracy is reduced by an order of magnitude. Since the standard regime is the one in which the velocity effect is dominated by the diffusion effect; i.e., ω v ω D the decay rate equals tõ S ≈ 1 −ṽ 2 2 , which is quadratic in the velocity; hence the sen-

Diffusion length equals drift length
Time hB(T + t)B(T )i < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > (c) Figure 3. (a) The structure of the power spectrum for diffusion only -The power spectrum is divided into two regimes which are defined by ω D . The regime ω ω D is characterized by the expected Lorenzian type decay. The ω ω D regime is described by a √ ω decay. (b) The structure of the power spectrum for diffusion and drift in normalized units -The power spectrum has three different regimes which are divided by ω v and ω D . The high frequency part ω > ω D is roughly equal to the power spectrum at zero velocity and damped by a factor of (1 −ṽ), while for the intermediate region ω v < ω < ω D this damping factor increases to approximatelyṽ. The low frequency part is nearly a constant, which equals S v=0 (ω = ω v ) (c) The structure of the correlation function -The decay of correlations changes at the point where the diffusion distance equals the drift distance, t = D v 2 . In the case that the correlation did not decay substantially within an NV coherence time, T 2 , the correlations in the NV measurements are proportional to the magnetic correlations. sitivity with respect to the velocity is substantially reduced in this regime. This leads to the natural conclusion that whenever the effect of diffusion is larger than the velocity effect, it is very challenging to estimate the velocity. This can be intuitively understood in the following way. The velocity is estimated from changes in the magnetic field due to the drift of the fluid, as in fig. 1. Diffusion, however, causes the magnetic field to change as well. When the molecules move away from the vicinity of the NV due to diffusion rather than drift, the effect of the drift is suppressed.
Before we proceed, let us quantify this effect. The sensitivity of the velocity is optimal at ω = 0. Therefore, the uncertainty in the velocity, √ This can be also rewritten as ∆ṽ = 1 vφ √ M , where φ is the phase which is collected during a time of τ D and M is the number of measurements. For the experimental parameters of water: D ∼ 2.3 · 10 3 nm 2 µs −1 , n ∼ 33 nm −3 [41] and d ∼ 10 nm, v = 1 mm s −1 , the sensitivity is approximately: T . Limited contrast and collection efficiency will reduce the sensitivity by a factor of 10, whereas collection from many NVs can increase it by a factor of 50 A µm 2 , where A is the area covered by the NV centers [27]. Thus in total A , (A is the area in units of (µm) 2 ) which requires a large NV area (150 µm) 2 to get reasonable accuracy. Differences in the distance from the surface in the NV ensemble result in negligible changes [40].
In the following we show that the Lorentzian approximation of the power spectrum is very crude since in fact the power spectrum is linear in the velocity or proportional to the square root of the velocity, depending on the frequency regime. It is possible to utilize these corrections to enhance the velocity sensitivity.
A heuristic sketch of the power spectrum without drift, based on Sec. VIII of [40], is shown in fig. 3a. While a Lorentzian type behavior is valid for large frequencies, the behavior for low frequencies deviates considerably from a Lorentzian and decays as 1 − ω/ω D . This non-analytic behavior derives from the fact that the correlation function at long times decays as 1 t 3/2 due to diffusion dynamics which inflicts a diverging derivative of the power spectrum at zero frequency. In the following we show that this non-analytic behavior has important implications for parameter estimation.
The velocity can be estimated either from the correlation function or from the power spectrum. Although these two quantities are related via the Fourier transform, the sensitivity obtained from each method is different, because of the modifications in the experimental procedure for probing them. The optimal strategy corresponds to different regimes of the coherence time of the NV center. In the case where the coherence time of the NV center is shorter than the correlation time of the signal, it is advantageous to probe the correlation function, since in this regime each measurement probes approximately a constant magnetic field. Therefore, the correlation between measurement results is proportional to the magnetic correlation function. In the opposite regime it is best to probe the power spectrum directly. Because an increased signal-tonoise ratio can be achieved by an ensemble measurement, a reasonable NV coherence time is 100 µs [18] after dynamical decoupling and should be compared to τ D when deciding on a correlation function or power spectrum probing.
Estimation via the correlation function -The general structure of the correlation function, based on Secs. V and VI of [40], is shown in fig. 3c. In the case where τ D is longer than the coherence time of the NV, which is valid for viscous fluids such as oil and NVs deeper than 5 nm, each measurement result is a function of the instantaneous magnetic field, thus making it advantageous to probe the correlation function. The time separation that yields infor-mation about the velocity can be found in the regime in which the drift length is larger than the diffusion length (vt > √ Dt). In this regime, the phase acquired by the NV during one measurement is φ ≈ γ e B(t)T 2 , and thus the correlation function is

RMS
. As the velocity in oil can reach 10 −2 mm s −1 , yielding an optimal depth of 25 nm and taking into account a limited collection efficiency, a fractional uncertainty of 2 is achieved. Thus, the total uncertainty is of the order of ∆v v = 1 A . The limitation is that this method works as long as T 2 is shorter than τ D and τ v , which only applies to very viscous fluids. For example, for an oil and NV depth of 25 nm, the diffusion time τ D ≈ ms [42]. For water on the other hand and for a shallower NV, τ D = 0.02 µs, which reduces the sensitivity by four orders of magnitude.
Estimation via the power spectrum -For lower viscosity fluids it is advantageuous to estimate the velocity directly from the power spectrum. The structure of the power spectrum, based on Secs. II and IX of [40], is shown in fig.3b and its behavior is divided into three regimes: the low frequency regime ω < ω v , the intermediate regime ω v < ω < ω D and the high frequency regime ω > ω D . In the high and intermediate frequency regimes the power spectrum is linear in the velocity and in the low frequency regime it has a square root dependence.
In the intermediate frequency regime, for ω ≈ ω D , the uncertainty is this constitutes an improvement by a large factor ofṽ −1 = ω D ω v in comparison with eq. (2).Substituting the parameters suiting water we estimate ∆v = 1 mm s s T for 15nm deep NV, which is the optimal depth as the power spectrum at this depth equals T −1 2 . For an ensemble of NVs we achieve the enhancement ∆v = 1 50 mm s s T (µm) 2 A . To appreciate this accuracy, consider that for water the dimensionless parameterṽ = 1 corresponds to a drift velocity of v = 100 mm sec . In the low frequency region the decay rate is equal tõ and the uncertainty is therefore, which is a substantial improvement by a factor of √ṽ compared to the intermediate range (eq. 3). This regime is applicable whenever T 2 1 ω v , which is attainable for low viscosity fluids as this time scale is around 10µs in that case. The accu- T v , is a decreasing function of v as expected. Using the same parameters of water, the accuracy in this is estimated by ∆v v = 10 −5 A . This result represents an improvement by more that four orders of magnitude in fractional uncertainty over state of the art methods even when demanding high temporal (second) and spatial resolution (micrometer).
Molecular dynamics - . . . In order to verify our analytical results we preformed a molecular dynamics simulation. Fig. 4 illustrates the simulation setup. The system consists of N dipolar particles confined to a simulation box of volume L x L y L z . Particles are interacting via the Lennard Jones (LJ) potential 4ε σ r 12 − σ r 6 with an interaction cut-off distance of r c = 2.5σ . The system is initialized into a thermal state at temperature T by applying a Langevin thermostat and each particle is assigned a random spin value I z ∈ {−1, 1}. After thermalization, the average particle velocity is subtracted from each particle, e.g. v j ←v j − v . The initialization concludes with the addition of the drift velocity v in the x-direction to each particle v x j ← v x j + v. The simulation dynamics is deterministic, following Newton's laws, and integrated using the Velocity-Verlet method with step size ∆t = 0.00125 ε/(mσ 2 ). During the simulation the z component of the magnetic field induced by the particles at the position of the NV center is computed using the equation where r i is the distance between particle i and the NV center, and θ i is the angle between r i and the NV's quantization axis.
The NV is placed a distance d below the simulation box in the center of the xy plane. Specular reflections are applied at the boundaries that are perpendicular to the z-axis and periodic boundary conditions are used in x and y directions. To emulate particles coming into and leaving the simulation box, we  flip the spin of particles that reach one of the periodic boundary walls with probability 1/2. √ v + cv, with b = −0.42 and c = 0.053 which yielded R 2 = 0.9926, thus corroborating eq. (4) and the following accuracy estimation (5). The difference b < 1 between Eq. (4) and the simulation fit should be accounted for by the missing prefactor of the √ṽ term, which cannot be calculated using our analytic tools.
Conclusion -We have proposed a novel setup to study flow properties in micro-and nano -fluidic structures which meets a crucial need in the fast growing field of microfluidics. The proposed setup is quantum inspired and is based on quantum sensing via color centers in solids.
The various parts of the dipole-dipole interaction could be used to sense the dynamics of the nuclei. A natural division of the different terms in the interaction is via the T 0,±1,±2 parts [1]. In the following, we mainly concentrate on the T 0 and T ±1 and analyze their efficiency to probe the dynamical quantities of polarized and unpolarized fluids. We will show that in complete contrast to all other NMR examples, polarization does not help in this case and all quantities could be estimated also in the unpolarized case with a similar signal-to-noise ratio.
A. The term T 0 = 3 cos 2 (θ ) − 1 S z I z This term could probe both the polarized and the unpolarized dynamics. The two cases will, however, result in the same efficiency.

Unpolarized nuclear spin ensemble
The Master equation in this case originates from the following stochastic Hamiltonian: where r i (t) is the vector connecting the NV center and the ith nucleus, θ i is the angle between between r i and the quantization axis, S z and I i z are the components of the spin operators of NV and ith nucleus respectively along the quantization axis of the NV. The total number of nuclei is denoted by N.
As the dynamics of the nuclei is probed via a noise measurement and is manifested in the dephasing of the NV we apply the Master equation starting with the equation of motion for the joint density matrix of the NV and the nuclei, It follows from Eqs. (1) and (2)ρ where · stands for average over realizations of nuclear spin polarization and their trajectories. Namely, · denotes tracing over the spin degrees of freedom, and averaging of the spatial ones. The density matrix χ of the nuclear spins and NV is assumed to be in the form χ = ρ NV ⊗ ρ Nuclei , where ρ Nuclei is approximately constant in time. Hereon, we denote the NV's density matrix as ρ for brevity. As all the spin operators commute, the nuclei spin operators, I i z can be assumed to be a classical uncorrelated random variables taking values x i = ±1 in the non-interacting limit. Equation (3) in the weak coupling regime giveṡ where The master equation (4) describes a dephasing process due to the random magnetic field induced by the nuclear spins at the NV's location. The dephasing is the strongest when the nuclei are immobile since both diffusion and drift, tend to smear the magnetic field thereby reducing the fluctuations. The correlation function, (5) could be further simplified as follows: where the last equality holds as different nuclear spins are uncorrelated in the non-interacting limit. We rewrite the expression (6) in terms of the nuclear density, n and the conditional probability P( r,t| r 0 ,t 0 )d 3 r for the nucleus located at r 0 at the time t 0 to be found in the volume element d 3 r centered at r at a later time instant t > t 0 [2], where since the correlation function (6) is symmetric with respect to the interchange of t and t , we have defined P ( r,t | r ,t) = P ( r,t| r ,t ) for t < t. In addition, we have assumed the steady state with the nuclear density being constant in space and time.
By introducing the dynamical decoupling we obtaiṅ where ω is the dynamical decoupling frequency. Equation (8) expresses the dephasing rate via the power spectrum, S BB (ω). And in this way the measurement of the dephasing rate allows one to study the effect of the flow on the power spectrum.

Statistically polarized nuclear spin ensemble
In the case of finite polarization, the polarization x i takes the value ±1 with the probability p and 1 − p, respectively. We have x i = 1 − 2p for the mean and (x i − x i ) 2 = α = 2 p(1 − p) for the STD. The finite mean average leads to an Overhauser field which is not important for us here, since it does not depend on the parameters v and D, characterizing the flow. Therefore, only the random fluctuations affect the noise. This will add a multiplicative prefactor of 4p(1 − p) in Eq. (6). It ensures that the dephasing is absent when the spins are fully polarized as expected.

Fully polarized nuclear spin ensemble
In the case of full polarization the calculation above produces a null result (p = 0 or p = 1) as in that case the NV will only feel the Overhauser field with no noise at all. This means that there is no classical noise. Noise, however, could still be induced quantum mechanically. This can be done by generating the nuclear spin polarization in the x − y plane and let the NV-nuclei interaction decohere the NV. This can be seen in the following way. Starting with the NV in the | ↑ x and the nuclei in the x direction (see fig. 2(a)), the state of the NV and the nuclei is NV center Nuclei initial state Nuclei final state In the polarized case all the nuclei start in a definite direction in the x − y plane, for example the x direction. (b) The interaction between the NV and the nuclei creates an entangled states of the GHZ type where the | ↑ z is entangled with the nuclei which rotate clockwise and | ↓ z is entangled with the nuclei which rotate counter-clockwise; although each nucleus rotates by a different angle which is dictated by the individual interaction strength. This entanglement is detected via the NV by a measuring its coherence (purity). The faster the flow the weaker the effective interaction between the NV the the nuclei, which can be estimated the by amount of dephasing. The same effect is detected for polarized and unpolarized nuclei.
Under the T 0 term this state evolves to (see fig. 2 where the θ i 's are the rotation angles from the x axis in the x − y plane due to the NV -nuclei interaction. The NV's density matrix is Measuring in the x basis we get: as The expectation value of the x-component of the NV spin polarization, which holds for immobile nuclei. The nuclei's motion, however, decreases the dephasing rate due to dynamical averaging, which means that instead of , where T is the total time. Assuming the Gaussian distribution of the θ i s we have cos θ i = exp − 1 2 θ 2 i . We then have from Eq. (14), Tr (ρ · σ x ) = exp (−S(ω = 0)T /2). The same can be of course rigorously derived from the master equation as in Eqs. (3) and (4).
We shall now show that this result also applies for an unpolarized ensemble of nuclear spins. First, lets examine a single nuclear spin polarized to | ↓ x . The joint state of the NV and the nuclease is | ↑ x NV | ↓ x , and it propagates in time to . It follows that the inner product (13) is not affected by the different initial polarization of the nuclear spin. For many nuclear spins, as |ψ 1,2 are tensor product states, (13) still remain unchanged if each spin's polarization is chosen at random to be | ↑ x or | ↓ x . This can be seen clearly by fig. 2, with simple modifications -initially each nuclear spin starts in | ↑ x / ↓ x , after some time, due to the interaction, the i th nuclease is found rotated by an angle θ i from it's initial position. Since the direction of the rotation is dependent on the state of the NV, and the value of θ i is dependent on the position of the nuclear spin we arrive at the same dephasing dynamics as before.
Finally, we note that since the unpolarized ensemble is invariant under spin rotations, the direction chosen to represent the individual polarization of each nuclease is irrelevant. This scenario is, therefore, equivalent to the unpolrized ensemble presented in the previous sections.
This term could also be used for sensing the parameters v and D. In order to turn this interaction resonant, a very specific magnetic field [3] or dynamical decoupling pulse sequence [4,5] has to be applied. This term allows for an exchange of one energy quanta between the NV and a nuclear spin. The optimal exchange is achieved for immobile spin ensemble, since the nuclei motion, due to diffusion or drift, shifts the interaction from resonance, and therefore decreases the efficiency of the process. By measuring the state of the NV center, the rate of the process can be estimated and the flow parameters v and D can be deduced.
The master equation for the NV's density matrix in this case iṡ This results in a dephasing dynamics with a dephasing rate which is equal to the power spectrum at the Larmor frequency. Namely, this result is the same as in the previous section, where the peak of the power spectrum is at the energy difference between the nuclei and the NV, which is approximately the NV's Larmor frequency. Eq. (15) holds both for polarized and unpolarized nuclear spin ensembles, as in both cases the NV's initial state can be chosen to allow this type of energy exchange.
C. The term T ±1 = 3 2 sin θ cos θ S z I x ± iI y Unlike the T 0 term that contributes to the low frequency part of the power spectrum, this term affect the power spectrum at the vicinity of the nuclear spins Larmor frequency as in [6][7][8][9]. Therefore, dynamical decoupling is always required for this interaction to be the dominant one. After the dynamical decoupling this term can effectively be written as T ±1 = 3 2 sin θ cos θ S z (I x cos (δt + φ ) + I y sin (δt + φ )), where δ is the detuning of the dynamical decoupling frequency from the nuclear Larmor frequency. At first sight it seems that this term is very much different from the previous ones, however, we will see that the results are very much similar.
First, we shall examine a polarized nuclear spin ensemble. The NV is initialized at the | ↑ x NV state and the nuclear spins are all initially oriented in the | ↓ z stat, see fig. 3. The T ±1 terms causes them to rotate around an axis in the x − y plane. The axis will be time dependent, but for a propagation time τ in which δ τ 1 the rotation axis is approximately a constant vector in the x-y plane, and therefore the total time propagation can be viewed as concatenated rotations around different axis in the x-y plane. The exact axis of rotation is not important, as entanglement between the NV state and the nuclei state will be generated for any rotation axis in the x-y plane. Thus, the same dephasing dynamics as in the previous sections is achieved. This can also be seen by comparing fig. 2 to fig. 3 -the two processes are the same as in both cases we get a GHZ type state in different bases. It can also be shown that quantitatively the result is similar.
By following the same reasoning as in part I A 3, the same dephasing dynamics can be shown to apply for an unpolarized spin ensemble. This can be seen by comparing fig. 3 and fig. 4 -the important parameter is the angle difference between the nuclei which are entangled to the | ↑ z and the ones which are entangled to the | ↓ z as in Eq. (13). This is left unaffected by the initial polarization of each nuclues.

NV center
Nuclei final state (a) + Nuclei initial state In the polarized case all the nuclei start in a | ↓ z state. (b) The interaction between the NV and the nuclei creates an entangled states of the GHZ type where the | ↑ z NV is entangled with the nuclei which rotates in the counter-clockwise direction and | ↓ z NV is entangled with the nuclei which rotates in clockwise direction with respect to an axis in the x − y plane (black arrow). This entanglement is detected via the NV by a decoherence measurement.

II. QUALITATIVE ESTIMATES OF THE POWER SPECTRUM BY DIMENSIONAL ANALYSIS
In the previous section we showed how the dephasing rate, which is proportional to the power spectrum could be measured, and we argued that the dynamical properties of the liquid can be estimated from this rate. Therefore, calculating the power spectrum and the temporal correlation function in terms of properties of the liquid is required for all our sensitivity estimates. In the following section we present crude estimation for the power spectrum based on dimensional analysis. We chose to present them first, since they are simple, independent of geometry, and most importantly they reproduce our main results while capturing their intuition.

A. Conditional probability as a Green function of the drift-diffusion equation
The noise caused by the drift/diffusion is characterized by the conditional probability P(r,t|r ,t ) of a particle being injected at time t at the locationr to be found at locationr at subsequent times, t > t . At times exceeding the inter-particle collision time the motion becomes diffusive and for these times P(r,t|r ,t ) can be found by solving the differential equation it satisfies. This equation at t > t follows from the Fick's law for the probability current, wherev(r) is the drift velocity and D is the diffusion coefficient. The drift velocityv(r) is assumed to be time independent and given. The probability conservation is given by Combining (16) and (17) we readily obtain, In the unpolarized case half of the nuclei start in a | ↓ z state and half in | ↑ z . (b) The interaction between the NV and the nuclei creates an entangled states of the GHZ type where the | ↑ z NV is entangled with the nuclei which rotate in the counter-clockwise direction, upper row. | ↓ z NV is entangled with the nuclei which rotates in clockwise direction with respect to an axis in the x − y plane (black arrow), lower row. This entanglement is detected via the NV by a decoherence measurement.
Assuming that the flow has no sources or drains in the region of interest, and that the fluid is incompressible we have, and eq. (18) simplifies to By definition, the conditional probability complies with the initial condition, (20) and (21) can be combined into the single equation, It can be shown that the power spectrum of the noise is a weighted Fourier transformation of the conditional probability [2], where the weighting functions f 1,2 (r) depend on the known form of the dipole-dipole interaction, presented in the previous section. The first equality is shown in the previous section (7). The second equality is due to the fact that the diffusion propagator is taken to be causal and coincides with [11]. The notation means integration over the volume occupied by fluid. The Fourier transform of the conditional probability defined in (23) satisfies according to (22) −iω − D∇ 2 +v ·∇ P ω (r,r ) = δ (r −r ) with the standard requirements on analyticity in the space of a complex frequency.

B. Universal features of the diffusion induced noise
In general, the geometrical constraints provide us with the typical distance d of the NV from the rest of the diffusing/drifting spins. The diffusion frequency scale is and we expect the power spectrum to depend on ω/ω D .
1. Static limit: ω = 0 In the static limit, ω = 0, the problem is equivalent to electrostatics, and (24) becomes the Poisson equation. We therefore have where we omitted the contribution of "image" charges, which are necessary in order to satisfy the boundary condition. This will be readily explained, while we first proceed with the approximation (26). Substituting (26) into (23) we obtain from dimensional analysis the following estimate, which means that where C 1 is a non-universal constant. The contribution of "image" charges may only affect the value of C 1 , because the typical distances fixing their location are again ∼ d. Therefore the "images" contribution can be eliminated. This is the same result obtained in [? ] with an appropriate choice of C 1 . This is used in the main text for S(ω = 0) of the Lorentzian model, and therefore as the normalization factor of the power spectrum in Eq. (1).

Low frequencies, ω ω D
Finite frequency yields a finite length scale, l ω = D/ω . Therefore, in the regime l ω d which is the same as ω ω D we may estimate roughly the suppression of the power spectrum by cutting off the distances that are larger than this characteristic length: The constants C 1,2 depend on geometry but other than that the result (29) is universal.

Large frequencies, ω ω D
Lets turn to (24), substituting v = 0, We estimate the spatial Fourier transform of (30), by concluding that the important values of q are of the order d. Therefore in the real space, which is the same as approximating (30) directly by considering the D∇ 2 ≈ −D/d 2 as a perturbation on top of the large frequency, ω D/d 2 . Returning to (32), since the typical distances are given by d, which is a typical damped oscillator behavior.
We could obtain this behavior starting from the Lorentzian curve, This crude dimensional analysis is validated later on in the spherical geometry. Equations (29) and (33) correspond to the behavior depicted in Fig. 3a in the main text.

C. Estimates of the effect of the drift on the power spectrum
We now make a crude but useful estimates of how the drift affects the power spectrum at different frequencies. We characterize the drift by the parameter similar to (25). To determine the change in the power spectrum we apply the idea of "Doppler shift" presented in [10]. Lets assume that the potential of the interaction between the NV and the spins in the fluid has a simple form, where the geometry is encoded in the scale d such that the typical q ∼ 1/d. For a harmonic wave in principle there are two contributions, cos(qx) = 1 2 e iqx + 1 2 e −iqx , which gives rise to the Doppler shifts ±∆ Doppler , that enter the propagator In the following, we expand (37) for "small" Doppler shift. We derive the exact requirement post-priori. There is an obvious concern about the cancellation of the linear terms in the expansion of (37). The cancellation appears because the two terms e ±iqx are present in (37) with equal weights. This cancellation will not occur if the two components were with different weights. This means that the magnetic fields f 1 (r) ∝ Y m l=2 and f 2 (r ) ∝ Y m l=2 enter (23) with different values of m, m = m . This is the case when the rotational symmetry is broken, i.e., when the NV's magnetization axis points in an arbitrary direction (exact definition is geometry dependent). This would bring about the breaking of the symmetry between the two senses of propagation. We note that even in the symmetric case, cancellation does not occur for ω = 0 since Hence, instead of (37) we proceed with We then Taylor expand (39) for sufficiently weak Doppler shift, In every given range of frequencies it implies different limitations on the parameters as we discuss below. By dimensional arguments (40) also implies the same relation for the power spectrum, Our claim is that all the results obtained in later sections are consistent with (41). We now consider different frequency regimes in turn. In this case the application of the Doppler shift equation (41) to the asymptotic expression (33) gives We can also deduce the relative change as which tells us that the expansion works in the range, Eq. (42) is used for the high frequency domain in Fig. 3b of the main text.

Low frequencies: ω ω D
In this case the application of the Doppler shift equation (41) to the asymptotic form (29) gives Lets understand when (45) holds. First, the correction has to be small relative to the first term. Therefore, we have to impose the inequality, Since we are already in the regime ω ω D , we must have ω v ω D . These conditions can be written in a compact manner together with (46), Yet, if we want to stop our series expansion, we must ensure that the next order contribution is negligible. Consider including one more term in the expansion (45) The absolute value of the ratio of the third and the second term is given by ω v 2ω , therefore, we have instead of (47) the following limitation, Eq. (45) is later verified for the planar geometry, and is used for the intermediate regime of Fig. 3b of the main text.

Low frequencies: ω ω v
To get an estimate in this case we assume that at low frequencies the only relevant scale is Doppler shift itself, ω v . Therefore, we can simply take the expression (45) and substitute in it ω = ω v , Although this is a very crude estimation, it agrees with the results of the molecular dynamics simulation presented in the main text, and it reproduces the results obtained in the spherical geometry shown in the next section. This result also seems to be universal for sufficiently slow varying flow profile. Eq. (51) coincides with Eq. 4 of the main text, and is used in Fig. 3b of the main text for the low frequency region

III. SPHERICAL GEOMETRY
After obtaining the generic behavior of the power spectrum we wish to consider a specific model. We start with a simple spherical geometry as a "toy model", which can be solved analytically and help develop intuition. Consider the fluid in the region outside a sphere of radius d with the NV placed in it's center. Moreover, we introduce the drift as a rotation of the fluid as a whole at the angular frequencyΩ = Ωẑ. We will first derive the propagator by solving Eq. (24) withv =Ω ×r. The drift term in this case,v ·∇P ω (r,r ) = (Ω ×r) · ∇P ω (r,r ) =Ω · (r ×∇P ω (r,r )).
If v is the velocity at the radius d in equatorial plane, we simply have henceforth we use the notation ω v instead of Ω.
Using the definition of the angular momentum,L =¯h ir ×∇ in spherical coordinates, we obtain for the drift term (52),v and (24) takes the form Recall that in spherical coordinates where the L-operator is divided byh. Substituting (57) in (56) and using the spherical harmonics completeness relation, We shall look for a solution in the form Substituting (67) into (58) results in Let us introduce the Doppler shifted frequency,ω (60) can then be rewritten as Defining the dimensionless distances we can write (60) as follows, where the sign of (±i) is the same as the sign ofω. For convenience we define, with the convention that the imaginary part is positive in both cases. Note that for any solution, F (z) of the spherical Bessel equation,    l ((±i) 1/2 x) form a pair of independent solutions of (64) for x = x . Because of the convention (65) the physically acceptable solutions will always be h (1) l ((±i) 1/2 x) at large x. The continuity at x = x and the convergence at x → ∞ tells us that we have a solution of the form, To fix the two constant A 1,2 we need two conditions, one is vanishing of the flux at r = d and the second is the jump of the derivative at r = r . Vanishing of the flux across the hard surface gives Lets introduce the notation, with a similar notation for h (2) l [(±i) 1/2d ]. Then by (68) we have the relation, To find A 2 we multiply (60) by (− r D ) and apply the operation lim ε→0 r +ε r −ε dr to both sides to get, which in terms of rescaled variables (63) translates into, Using (67) we write the last equation as Simplifying the last equation we are left with the following equation for A 2 , Notice that the Wronskian of the two independent spherical Hankel functions is which gives Then from (70) And from (67) Alternatively, It will be more convenient to express the result in terms of the dimensionless frequency, according to (63), With (80), (79) can be written as  Figure 5. The exact power spectrum, Γ(ω) as given by (89). The units of power spectrum are 1/(dD) and the power spectrum is plotted as a function of the dimensionless frequency, ω/ω D .

B. Asymptotic behavior in the limit of large and small frequencies
At lower frequencies the Taylor expansion of (89) yields in accordance with the general expectation (29). For finite fluid rotation frequency ω v = 0, (92) describes the Doppler shift of the diffusion smeared power spectrum. The Doppler shift holds of course for all frequencies as well. Note that at zero frequency we have a non-analytic square root behavior of the power spectrum, At large frequencies and at ω v = 0, the result (89) reduces to as again expected from (34). The question arises whether the singular suppression of the power spectrum as given by (93) is generic. The answer is unfortunately negative. In the next section we present a model with the angular frequency decaying to zero at large distances, for which the singular square root is eliminated by the combined action of the diffusion and frequency variations across the fluid volume, see Fig. 6. The model we previously presented allows for exact solution at finite drifts. It is, however, artificial inasmuch the drift gives rise to only a single Doppler shift, rather than a spectrum of many such frequencies shifts. Here we are try to overcome this limitation by looking at a different velocity profile in which the frequency decays with the distance, Repeating all the steps leading to (60), eq. (67) is still valid, but instead of (60) we now get We will focus on ω = 0 as otherwise the solution is not given in terms of known special functions. With x = r/d, and x = r /d, Eq. (96) implies that the function g lm (r, r ) in units of 1/(dD) satisfies where We can rewrite (97), To determine λ we may introduce the notation, (l + 1/2) 2 β = β , which for l = 2 gives β = β (25/4) l(l + 1) = (l + 1/2) 2 − 1/4, λ (λ + 1) = (λ + 1/2) 2 − 1/4.
We have the relation, where we assume β > 0 for definiteness, and that the square root has a positive imaginary part. The solution of (99) for x = x then reads, The boundary condition at the hard surface, x = 1, is Another boundary condition is which gives The previous equation can be solved for B, which finally gives, The coefficient A can then be calculated from (104) by substituting (108), The solution (103) can now be explicitly written using (108) and (109), The power spectrum can now be written as The integrals are easily evaluated, Substituting (112) into (111) we obtain This result is plotted in Fig. 6 as a function of the ratio ω v /ω D which is in terms of β for m = 1 reads ω v /ω D = β (25/4) using the definitions (103) with l = 2.

IV. SPECIAL NOTATIONS
In the next few sections we calculate integrals using series expansions of special functions. Here we clarify all the notations we use (Ω) -The spherical harmonics, Ω is the Solid angle. 14. p F q a 1 , · · · , a p ; b 1 , · · · , b q ; z -The generalized hypergeometric function.

V. TIME CORRELATION CALCULATION -PLANAR GEOMETRY WITHOUT DRIFT
In this section we derive a closed expression for the temporal auto-correlation function of the magnetic field, for a planar boundary condition, meaning, we assume that the NV is situated in a depth d within a planar diamond surface, and that the liquid fills the half space above the diamond surface.

A. Calculating the autocorrelation function in the time domain
Solution to the diffusion equation with a planar boundary condition -As in the previous sections, we are interested in correlation functions of the type (23), where 2 (r) are the spatial dependencies of the magnetic field, P(r 0 ) is the initial spatial distribution of the nuclear spin ensemble, assumed throughout this work to be uniform, and P (r,t|r 0 , 0) is the propagator from the pointr 0 tor given a time difference t. The first step in calculating these correlation functions is to find the appropriate propagator. Assuming, the ensemble is freely diffusing in the half space above the diamond surface, the propagator Ψ will be the solution to the diffusion equation in the half space with the initial condition Ψ (r, 0) = δ (r −r 0 ). We recall that the diffusion equation with the diffusion coefficient D is The equation can be solved for the whole space (by using the Fourier transform): We define a coordinate system in which the NV center is found at the origin and the diamond surface coincides with the z = d plane. The appropriate boundary condition is the vanishing of the flux through this plane, meaning: ∂ z Ψ| z=d = 0. We can enforce it by using a sum of two solutions to the diffusion equation (method of images): We can write the total correlation function of the magnetic field as a sum of the real spins contribution and that of the images, where we substituted f (m) (r) with the dipolar fields spatial dependence. We denoted by α (m) the squared dimensionless coefficient required to transform from the spherical harmonics representation of the dipolar interaction to the common angular notation, and by n the nuclear spin density. For simplicity, we assumed that the NV's magnetization axis is perpendicular to the diamond surface; namely, theẑ direction; Therefore, the only correlation functions that need to be considered are between f (m 1 ) and f (m 2 ) with m 1 = m 2 , due to the polar symmetry. We shall relax this assumption later on. Solution in cylindrical coordinates -In the following we calculate the correlation functions G (m) (t) analytically and without approximations. The main idea behind the derivation is to use the polar symmetry of the problem, therefore the important tools used are the 2D Fourier transform and the Anger-Jacobi expansion [13], where, J m (z) is the m'th Bessel function of the first kind. The main result of this section is √ π 2d 5 (122) Evidently, the correlation functions do not decay exponentially as usually assumed, but in a more complex fashion. The characteristic time τ D = d 2 D = (ω D ) −1 , dictates a change in the behavior of the function, rather then just the amount of decay. We now continue with the full derivation. Let us focus on the first part of the correlation function, (119), We solve this by using the Fourier transform with respect to (ρ −ρ 0 ) and then the Anger-Jacobi expansion (121), where θ /θ are the angles betweenk andρ/ρ 0 . Since we are in the plane we can write where ϕ is the free angle of thek vector. Substituting (127) back into (126) and integrating over the polar angles results in We shall now pursue the integration over ρ, ρ 0 for different values of m. Starting with m = 0, The integral over ρ is Substituting this result back into (132) Returning to (131) with m = 1, The integral over ρ will now be Using this result, the integral (135) reduces to Returning to (131) again with m = 2, The integral over ρ is Substituting back in (138) we arrive at Eqs. (134), (137) and (140) validate (123). Therefore, we turn our focus to m = 0. We now integrate over z 0 , z, noting that the integration limits are therefore [d, ∞), since the nuclei are confined to the upper half space, The integral over z is Substituting (145) back into (144) we can integrate over k, We now turn our attention to the other part of the correlation function, G re f lected , given by (125). The integration follows the same steps up to (134), The integration over z, z 0 yields Substituting (152) back into (147), We solve (153) by integrating each term individually, 5 .
Summing (154) and (155), we arrive at √ π 2d 5 (156) Finally, summing the two parts of the correlation function (146) and (156) leads the final result (122) Asymptotic behavior -We shall now examine the asymptotic behavior of the correlation function at short and long times. The result (159) presented in this section is used in Fig. 2c of the main text.
For t τ D we can expand eq. (122) as follows, This result is expected, since when taking the formal limit t → 0, the diffusion propagator approaches a delta function. Consequently, in this limitr →r 0 , and the integral (119), For t τ D we can expand eq. (122) to This also can be understood in terms of the integral (119). In this limit, the unnormalized propagator is approximately 1, therefore, These limits can be interpreted intuitively. In short times, the magnetic field squared goes as r −6 . Since the minimal distance from the NV is d, the main contribution comes from particles in about this distance from the NV. Their magnetic field decays as d −6 , and there are about nd 3 such particles in the effective interaction region, which is approximately a hemisphere of radius ∼ d on top of the diamond surface near the NV. Similarly, for long times the correlation ∼ r −3 r + √ Dt −3 . Taking the same proximity argument yields the long times limit.

VI. THE TIME CORRELATION OF DRIFTING AND DIFFUSING PARTICLES IN THE WHOLE SPACE
In the following section we derive the correlation function of diffusing nuclei subjected to constant driftv = vx. We shall further assume that the nuclei occupy the whole space, while the NV center is found at the origin. Though, this geometry is unrealistic, the asymptotic behavior (scaling) of the correlation function should hold for the planar geometry as well.

A. The drift-diffusion equation propagator
The drift-diffusion equation with constant drift is We solve it by first using the spatial Fourier transform, which implies where A is a constant. In order to retrieve the desired propagator we apply the inverse Fourier transformation, This is a Fourier transform of a Gaussian function with respect to the variablex +vt. Therefore, Adding the initial condition ψ (t = 0,x) = δ 3 (x −x 0 ) we arrive at as expected.

B. Explicit calculation of the correlation function
In the following we derive explicitly the correlation function G 0 . We show that for the specified geometry andv = vx, We start by substituting the propagator (167) into (114) together with the dipolar field dependency, (169) We define cylindrical coordinate, ρ = y 2 + z 2 , y = ρcosφ , z = ρsinφ , and rewrite eq. (169) as As in previous section, we use the 2D Fourier transform and the Anger-Jacobi expansion (121) to arrive at andφ is the free angle of the k vector. Substituting (172) into (171) and integrating overφ yields We now continue integrating over the angles and ρ, ρ 0 term by term. The first integral is The second integral is We note that therefore (180) can be written as The third integral is the same as the second with x 0 ↔ x, therefore, The froth integral is Summing the four results (176), (182), (183) and (189) leads to We now pursue the integration over x, x 0 . The integral over these coordinate is We solve the integration over x term by term: Again, we integrate term by term:

C. Asymptotic behavior of the correlation function
In the following we examine the asymptotic behavior of (168). First, we would like to check that the formal limit v → 0 recovers the result without drift (159), which is expected from dimensional analysis. These equations are used in Fig. 3c in the main text.

VII. THE DRIFT-DIFFUSION PROPAGATOR IN THE FREQUENCY DOMAIN
As explained in the main text, when T 2 > τ D , τ v , it is beneficial to estimate the drift velocity from the power spectrum. Since the calculation of the Fourier transform of the temporal correlation function is difficult, we chose a different approach -we calculate the correlation function in the frequency domain directly.
In the following section we take the first step towards this goal, when we derive the drift-diffusion propagator for a constant drift velocityv, We start our derivation with the drift-diffusion equation in the time domain, Fourier transform with respect to time yields We then use the transformation e −ā·r ∇eā ·r = ∇ +ā (209) in order to rewrite eq. (208), where in the last equality we used We chooseā =v 2D , so eq. (210) is simplified to We then defineψ = e −v 2D ·r ψ, and substitute into the previous equation, In order to find the propagator, we add a source term to eq. (213), The solution to this equation is given bỹ In order to solve the integral (217), we use the residue theorem. The poles of the integrand follow the equation where R and φ are defined by (206). The poles are, therefore, or alternatively, If ω > 0 the phase π 2 ≤ φ ≤ π so k 1 is in the first quarter and k 2 is in the third quarter. If ω < 0 the phase π ≤ φ ≤ 3π 2 so k 1 is in the second quarter and k 2 is in the fourth quarter. If ω = 0 the solutions are on the imaginary axis. The solution in all of these cases is the same, We now need to define the required integration paths for the residue theorem. For the first (left) integral to converge, we choose half a circle that goes through the positive imaginary axis, so the relevant singularity is k 1 . For the second integral, we choose half a circle that goes through the lower side of the complex plane, thus, the relevant singularity is k 2 . These are simple poles, so the integration is trivial, ψ = 1 4πD r −r k 1 2k 1 e ik 1 ·|r−r | + 1 4πD r −r k 1 2k 1 e ik 1 ·|r−r | = 1 4πD r −r e ik 1 ·|r−r | (222) Using the definition ofψ we arrive at (205).

VIII. ASYMPTOTIC BEHAVIOR OF THE POWER SPECTRUM OF DIFFUSING PARTICLES
In the following we derive the asymptotic behavior of the power spectrum of diffusing particles. We first show that for some constant C 1 , in accordance with (28). We then show that for ω ω D with ω D = D/d 2 , the power spectrum to the leading order in ω/ω D is for some constant C 2 , up to logarithmic accuracy in the small parameter. Eq. (224) agrees with the universal estimation (29).
A. Explicit expression for the propagator for v = 0 The first step in calculating the power spectrum is to derive an explicit expression for the propagator.
Substituting v = 0 at (220), Therefore, for ω ≥ 0, (205) can be written as In oder to enforce the boundary condition, ∂ ∂ z ψ z=d = 0, we shall use the method of images. The image propagator is Taking the derivatives of eqs. (226) and (227), it can be explicitly shown that the boundary condition is satisfied.

B. Approximate solution for low frequencies
In the following we explicitly derive the asymptotic behavior of the power spectrum, (224) at low frequencies ω ω D . We shall use the same notion of direct (indirect) correlation function for the "real" ("image") contribution as was previously introduced in (119). For brevity, we omit in the following the nuclear spin density n, which multiplies all the correlation functions.
Calculation of the direct correlation -The direct correlation function in the frequency domain is given by We start by estimating (228) at ω = 0. In terms of the dimensionless integration variables (r) i = (r) i d , (r 0 ) i = (r 0 ) i d , where in the last equality we used the spherical harmonics expansion for |r −r 0 | [12], and we omitted the prime notation for brevity. We continue by preforming the radial integration of (230), We solve the integrals of (231) for cos θ > cos θ 0 , The opposite case cos θ < cos θ 0 , from symmetry, will be the same with θ ↔ θ 0 . We can therefore define θ ≶ = θ 0 cos θ 0 ≶ cos θ θ else for and rewrite (230) as [I l (θ > , θ < |r > r 0 ) + I l (θ > , θ < |r < r 0 )] The summation and integration over θ , θ 0 is just a numerical factor which we can evaluate. It can be approximated by truncating sum , for example, at l = 4 since the terms go as ∼ 1 (2l+1)l and higher l terms contribution will be smaller by a large factor. We note that the scaling of (235) is in accordance with the long time limit of the correlation function (159), since it's truncated Fourier transform over times for which the approximation t τ D holds is which also scales like 1 Dd . We note that this coincides with the simple calculation (28). Now, we would like to go beyond the static limit ω = 0, and estimate the leading order of (228) in ω ω D . For this we use the series expansion [12] where cos Θ =r ·r 0 , in order to expand (228) into We define the dimensionless parameters In terms of the parameters (239) the integral (238) can be written as Lets examine the integration of (240) over the magnitude of x, while substituting (225), We define (243) So (242) can be written as where log is the natural logarithm and γ is the Euler constant. For x 1, it follows that [12] j l (x) ≈ x l (2l + 1)!!
Using the approximations (249),(250) and (251), the integral(248) evaluated by Integrating (253) with respect to √ ω yields, up to an integration constant, √ ω · (2l)! l!2 l · l cos θ 1 − cos θ cos θ 0 l (257) Substituting (255) and (257) into (244) we arrive at where we used Some of these approximations might seem hasty, since even if ωτ D 1 it does not necessarily mean that ωτ D cos θ 1. Our approximations are checked post-priori -if the angular integrals converge, we expect an error by at most a prefactor, but the scaling of physical constants to be correct. It can be shown that the integral does converge for l = 0 which is the most divergent contribution to the sum.
Calculation of the indirect correlation -We now evaluate in a similar fashion the indirect contribution to the correlation function, We define The limits of the integral (261) in terms of the new variables (262) are z ∈ [0, ∞) , z 0 ∈ (−∞, 0], and in spherical coordinates r , r 0 ∈ [0, ∞]. Omitting the prime notation from the radial coordinates for brevity, (261) can be written as We would like to expand (263) in a way which is similar to a multipole expansion. We can write, For l > 0 [12], Therefore, by substituting (265) into (264), We now use (266) in order to simplify the integral (263). For now, we will only concern ourselves with the integration over the magnitude of r, r 0 , Using the expansion (237) as in the previous part, the integral I (1) lk , defined in (268), can be written as We are interested in the limit ωd 2 D 1, for which we can use the expansions in order to simplify (269), (2n)! (2n − 1)!!n!2 n P n (cos Θ) We recall that the leading contribution of the direct correlation in the parameter ω ω D is ∼ ω ω D log ω ω D by eq. (258), while in (273) the leading contribution ∼ ω ω D , which is smaller. We therefore look for larger contributions in the remaining integrals. Using (237) in order to expand I We can estimate the leading order of (278) in which is also small compared to the leading contribution of the direct correlation. The last term I This is the radial integration of the direct correlation function (228) with cos θ = cos θ 0 = 1. So we know by (258), that it will diverges as log 2 with the same prefactors as (258). We note, that the angular dependency differs from (240) -apart from the distinction made above about the radial integration limits, we recall that the integrals are not defined over the same coordinate system. Since we have an explicit expression for the time correlation (122), and it does not depend (up to a numerical prefactor) on the value of m, we expect that for all m it will diverge as log 2 or all will diverge logarithmically.

IX. ASYMPTOTIC BEHAVIOR FOR LOW VELOCITY
In the following section we want to extend our previous analysis to include finite drift velocity. We will analyze the correlation functions in the regime v 2 Dω 1, which we name the "low drift limit". It signifies times where the motion of nuclei is "diffusion dominated", meaning, that the diffusion length scale is much larger than that of the drift -vt √ Dt → v 2 t D 1, which translates into (288) in the frequency domain. We first continue with assumption that the NV's magnetization axis is perpendicular to the diamond surface for simplification, and later on we show how it can be relaxed.

A. NV's magnetization axis coincides with boundary condition
In the following we show that in the low drift limit the correlation scales quadratically in the drift velocity; namely G . This corresponds to (37), where the linear dependency in the drift velocity is eliminated because of the symmetry.
Approximating the propagator -Since the velocity dependency enters the correlation function only through the propagator, the first step in taking the limit (288) is to approximate the propagator. We recall (205), ψ = 1 4πD |r −r 0 | e ik 1 ·|r−r 0 | e¯v where In the low drift limit (288) we can approximate (290) by ik 1 ≈ e i π We note that under the assumption that the magnetization axis is perpendicular to surface the first order will surely be eliminated due to azimuthal symmetry. Therefore, the leading term in the expansion will be of second order. We continue with the spherical Bessel expansion [12] of (296), (2l + 1) i l j l −i v 2D r −r P l (cos Θ) .
Only even values of l will survive the integration due to the symmetry, The direct correlation with the approximated propagator (299) is given by where we, again, omitted the multiplicative factor of the nuclear spin density n for brevity. We define the dimensionless integration variablesr We then write the integral (301) in terms of the variables (302) 1 + e i 5π The zeroth order in our small parameter, defined by (288), reproduces the correlation function without drift. The leading contribution of (304) will be given by l = 0, 1,