Utilising NV based quantum sensing for velocimetry at the nanoscale

Nitrogen-Vacancy (NV) centers in diamonds have been shown in recent years to be excellent magnetometers on the nanoscale. One of the recent applications of the quantum sensor is retrieving the Nuclear Magnetic Resonance (NMR) spectrum of a minute sample, whose net polarization is well below the Signal-to-Noise Ratio (SNR) of classic devices. The information in the magnetic noise of diffusing particles has also been shown in decoherence spectroscopy approaches to provide a method for measuring different physical parameters. Similar noise is induced on the NV center by a flowing liquid. However, when the noise created by diffusion effects is more dominant than the noise of the drift, it is unclear whether the velocity can be efficiently estimated. Here we propose a non-intrusive setup for measuring the drift velocity near the surface of a flow channel based on magnetic field quantum sensing using NV centers. We provide a detailed analysis of the sensitivity for different measurement protocols, and we show that our nanoscale velocimetry scheme outperforms current fluorescence based approaches even when diffusion noise is dominant. Our scheme can be applied for the investigation of microfluidic channels, where the drift velocity is usually low and the flow properties are currently unclear. A better understanding of these properties is essential for the future development of microfluidic and nanofluidic infrastructures.

condition [33][34][35] , which is commonly used to solve the Navier-Stokes equations. Though there have been advances in the microscopic theory of this phenomenon 35 , an experimental method to accurately measure the velocity near the surface has not yet been developed. Our proposed setup is sensitive to surface effects and therefore meets this need.
The main idea is sketched in Fig. 1, where a flow through a microfluidic channel is sensed by an NV center ensemble, in a setup similar to the one that was used in prior works 23,24 . 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 NVs. By analyzing the noise characteristics, the flow properties can be deduced. It should be noted that velocimetry methods using classic NMR techniques exist [36][37][38] . In principle, these methods can also be applied in the nanoscale NMR settings. The implementation of this type of experiment is challenging, since it requires efficient polarization of the nuclear spins together with a strong, stable magnetic field gradient.
In this work, we show that a nanoscale NMR approach for velocity estimation in microfludic channels is possible without polarization and even when diffusion dominates the dynamics. The remainder of this paper is organized as follows. First, we discuss measurement protocols. Then, we estimate the sensitivity for the velocity measurement assuming a Lorentzian noise spectrum. We proceed by showing that in fact the spectrum deviates from the Lorentzian approximation due to the diffusion process and that the corrected spectrum leads to improved sensitivity even when the diffusion process is more dominant than the drift. Finally, we present molecular dynamics simulations that support our analytic results and compare our method to current velocimety techniques.

Results
Measurement scheme. We 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 1 . The power spectrum generated by the motion can be read by a decoherence spectroscopy method [39][40][41][42][43][44] or spin relaxation 45,46 . These methods exploit the fact that the transition rate or relaxation time in a two-level system is proportional 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 (see Fig. 2). The relaxation rate and, therefore, the power spectrum can be efficiently estimated from the population decay. Since all the effects of the velocity emerge from the magnetic noise induced on the NV, polarization is not essential here and the velocity can be measured in unpolarized samples in a way similar to diffusion measurements 22 .
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, S ± (I ± ) as the NV's (nuclear spin's) raising/lowering operators and θ i as the angle between the NV's magnetization axis to the vector connecting the NV to the i-th nucleus. The velocity can be sensed by the T S I (3cos( ) 1) z z (see for example 47 ) as shown in Fig. 2. The T 0 is the classical term derived 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. This is demonstrated in the Supplementary Information (SI).
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 16,18,20,21 . In order to probe this term, a rotation of the NV at the Larmor frequency must be introduced, for example, by pulsed dynamical decoupling 5,[15][16][17][18][19]48 , which yields the following effective Hamiltonian: H S I t I t sin cos ( cos( ) sin( )) , where δ is the frequency difference between the Larmor frequency and the frequency of the pulses. The velocity can be estimated also by measuring the correlation function of the magnetic field directly. The idea is presented in Fig. 2b and discussed further later on. 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 near the surface of the diamond. The NVs sense the magnetic field, whose randomness is amplified by 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. Created by Studio Mirage, Israel. , 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 where both processes are significant, the natural expectation is that the characteristic width is σ τ In order to estimate the sensitivity we follow the commonly accepted assumption that the nano-NMR spectrum is a Lorentzian function 20,22 where γ e∕N is the electron/nuclear gyromagnetic ratio, μ 0 is the vacuum permeability, ℏ is the reduced Plank constant and n is the nuclei number density. A 5 nm deep NV and water density of 33 nm −3 result in B RMS ≈ 250kHz∕γ e . The sensitivity is best described by the dimensionless quantities: Figure 2. (a) The Lorentzian power spectrum at different velocities -The power spectrum consists of multiple peaks associated with different parts of the dipole-dipole interaction between the NV and the ensemble. Shown here are the zero frequency term, T 0 ∝ S z I z term (dotted light blue), and the Larmor frequency terms T ±1 ∝ S z I ± (dotted purple). Each point on the spectrum is associated with a different decay rate of the NV's quantum state; for example the zero frequency peak is proportional to the NV's pure dephasing time T 2 . The velocity changes the width and the peak of the spectrum, and therefore it can be deduced by measuring them. (b) The correlation function at different velocities -the velocity can also be estimated by directly measuring the correlation function, the Fourier transform of the power spectrum. The correlation oscillates at the Larmor frequency, and decays in a rate that depends on the velocity and diffusion. The peak of the correlation is proportional to τ γ B c e RMS 2 2 , where τ c is the correlation time, which is also determined by the velocity. (c) While the spectrum (blue) at zero frequency can be measured in a Hahn-echo experiment, different parts of the spectrum require external driving of the NV. The Fourier transform of the pulse sequence acts as a window function on the the power spectrum, such that the dephasing time T 2 is given by a convolution between the window function and the power spectrum. Shown above are two examples of window functions at frequencies 1/2 1/2 ω = π τ (orange/ green) that sample different parts of the spectrum. The window functions peak sharply at the characteristic frequency of the external drive and they have smaller peaks at higher harmonics of that frequency. Both the window functions and the power spectrum are symmetric with respect to ω = 0.
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 . We use these units henceforth in our estimations. The decay rate Γ of the NV can be estimated with a sensitivity 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 a sensitivity of ∆Γ = Γ T . Since most NV measurements are not single shot, the sensitivity is reduced by an order of magnitude compared to this estimate. The standard regime is the one in which the velocity effect is small compared to the diffusion effect; i.e., ω v ≪ ω D , therefore the decay rate equals ≈ − S 1 , v 2 2 which is quadratic in the velocity; hence, the sensitivity 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 estimation is optimal at ω = 0. Therefore, the uncertainty in the velocity, S This can also be rewritten as where φ is the phase collected during a time τ D and M is the number of measurements. For the experimental parameters of water: . Limited contrast and collection efficiency will reduce the sensitivity by a factor of 10, therefore the sensitivity achieved by a single NV, henceforth denoted by a subscript s, is ≈ . This can be mitigated by collection from many NVs, which increases the sensitivity by a factor of 100 A which requires a large NV area ~(150 m) 2 μ to get reasonable sensitivity. Note that we assume an NV density of 10 12 cm −2 , which is about the current upper bound for shallow NV ensembles 50 . Therefore it should be interpreted as an upper limit. Differences in the distance from the surface in the NV ensemble result in negligible changes in sensitivity (see SI).
In the following we show that the Lorentzian approximation of the power spectrum is very crude since in fact the power spectrum is linear or proportional to the square root of the velocity, depending on the frequency regime (see Methods). This enhances the sensitivity (Eq. 2) by a factor of up to v 3/2  (see Eq. 5 in sec. 2.4). For water flowing in a microfluidic channel v 5 10 3  ≈ ⋅ − , which implies a decrease in the uncertainty by a factor of v 3 10 3/2 4  ≈ ⋅ − compared to eq. 2. A heuristic sketch of the power spectrum without drift is shown in Fig. 3. While a Lorentzian type behavior is valid for large frequencies, the behavior for low frequencies deviates considerably from a Lorentzian and decays as 1 . This behavior derives from the fact that the correlation function at long times decays as t −3∕2 due to diffusion dynamics which transforms to ω dependence in the power spectrum. This non-Lorentzian behavior has been observed experimentally for diffusing Rb vapor in N 2 gas 51,52 . 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, T 2 . 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.
The coherence time is dictated by external noise factors, such as magnetic impurities, that are not taken into account in our estimation. In shallow NVs the main source of noise is due to surface effects [53][54][55][56] , while the dynamics is controlled by magnetic impurities for bulk NVs [57][58][59] . Though the coherence time of single NVs has gradually increased in past years, due to surface engineering, annealing techniques and sophisticated pulse sequences, the task remains challenging for a dense ensemble 50 . This will not effect the estimation via the correlation function, as it assumes that T 2 is the smallest time scale of the system. It can, however, impact the power spectrum estimation, as it requires T 2 to be larger than τ D . Hence, when performing an experiment, an optimization will be required to achieve the best sensitivity.
We henceforth assume a reasonable NV coherence time is 100 μs 20 after dynamical decoupling. Since achieving such coherence times might suggest taking lower NV density, the numeric estimations via the power spectrum might deviate by an order of magnitude.

Scientific RepoRtS |
(2020) 10:5298 | https://doi.org/10.1038/s41598-020-61095-y www.nature.com/scientificreports www.nature.com/scientificreports/ estimation via the correlation function. The general structure of the correlation function, based on secs. V and VI of the SI and the Methods section, 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 information 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 As the velocity in oil can reach 10 −2 mms −1 , yielding an optimal depth of 25 nm and taking into account a limited collection efficiency, a fractional uncertainty of 5 10 is achieved. Thus, the total uncertainty is on the order of = ⋅ .
The limitation is that this method only 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 6 . 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 advantageous to estimate the velocity directly from the power spectrum. The structure of the power spectrum, based on secs. II and IX of the SI and the Methods section, 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 velocity and in the low frequency regime it has a square root dependence. Further exploration of the parameter space is shown if Fig. 4. In Figs. 4a,b the power spectrum as a function of ω and ω D is presented for d = 10nm and v = 0. The effects of finite drift are shown in Figs. 4c,d, where is plotted as function of the same parameters. Since the most likely experiment would be to measure first the spectrum without drift and then to perform a second measurement with finite drift and compare the two, the normalized contrast is the most relevant quantity.
In the intermediate frequency regime, for ω ≈ ω D , the uncertainty is This constitutes an improvement by a large factor of = In the low frequency region the decay rate is equal to and the uncertainty is therefore,  μ ∆ . This result represents an improvement by more then three orders of magnitude in fractional uncertainty over the Lorentzian model.

Molecular dynamics.
In order to verify our analytical results we performed a molecular dynamics simulation. The simulated system consisted of N interacting particles confined to a box and the NV was placed a distance d below the simulation box at the center of the x-y projection of the box.   with finite drift velocity as a function of the frequency and the diffusion frequency at d = 10 nm. The power spectrum is given by Eq. 23 for ω ≤ ω v and by Eq. 22 for ω v < ω with C = 1. The bounds of (c/d) match those of (a/b) and the drift velocity is taken to be 0 01/1 mm s .
. The contrast decays gradually as a function of the two frequencies. The strong contrast results in an improved sensitivity scaling. (2020) 10:5298 | https://doi.org/10.1038/s41598-020-61095-y www.nature.com/scientificreports www.nature.com/scientificreports/ Comparison to fluorescent based velocimetry. 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 [60][61][62][63][64][65] . Alternatively, the velocity can be determined by light scattering from a periodic array 66 . Fluorescent molecules can also be used to measure the diffusion coefficient of these molecules within a liquid 64 . These methods lack in performance and only achieve an accuracy of about 5% for the mean channel velocity due to various 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. The fluorescent molecules are often large, and in a small channel they can affect the flow profile. In addition, the flow trajectories of these molecules generally pass through the middle of the tube 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. Our proposed setup solves these problems since it is non-intrusive and sensitive to surface effects.
To verify our claim that our method outperforms the fluoresce based techniques we note that they all suffer from a limitation on the molecule size. The relative error that is caused by the Brownian motion can be estimated by 61 where T is the experiment time. Since the diffusion is an unbiased noise, this error can be reduced by averaging where N, is still limited by systematic errors. The distance between the molecules, provides these system a characteristic length, d c , which is bounded by the molecule's size. This enables us to write (7) in a similar way to our previous assessments, We note that the relative error, Eq. 8, has the same scaling in  v as the one in the Lorentzian model given by Eq. 2. Therefore, these methods will not be able to estimate velocities for which <  v 1. The lower limit will be given by v c D d c = . This is problematic as d c should be small enough, for the molecules not to change the flow profile in the channel, but it has to be large enough to eliminate the effect of diffusion. This poses a fundamental problem, which cannot be overcome. Our method provides an improved scaling in  v even when it is small. www.nature.com/scientificreports www.nature.com/scientificreports/ conclusion Here we proposed a new 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 proposal relies only on statistical polarization and therefore can be easily implemented. Our setup outperforms current velocimetry methods as it is sensitive to surface effects and it is effective even when diffusion dominates the dynamics.

Methods
The correlation function and power spectrum of diffusing particles. The dipole-dipole interaction can be decomposed into six parts associated with angular momentum exchange between the NV and a nucleus. Each of these terms rotates at a different frequency when the NV is driven by an external field and therefore contributes to a different part of the power spectrum.
This leads to the definition of the correlation functions, where where are the spatial dependencies of the magnetic field given by the Spherical Harmonics, 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 point r 0 to r given a time difference t. The dimensionless correlation functions G m m ( , ) 1 2 are equal to the correlation functions of the magnetic field up to a multiplicative factor of ( ) . 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 P will be the solution to the diffusion equation in the half space with the initial condition δ = = − P r t r r r ( , 0 , 0) ( ) 0 0 . Solving for the whole space and using the method of images we find ( 2) 4 If the NV is not tilted, meaning, it's magnetization axis coincides with the normal to the diamond surface, the only correlation functions that do not equal zero are ) Calculating Eq.11 by performing the integrals in Eq. 9 explicitly we find that  www.nature.com/scientificreports www.nature.com/scientificreports/ The asymptotic behavior of Eq. 13 is given by The instantaneous correlation defines the B RMS , but the long time limit, often assumed to decay exponentially, is actually given by a power law. The power spectrum around ω = 0 will be give by the Fourier transform of Eq. 13, www.nature.com/scientificreports www.nature.com/scientificreports/ as the special functions strongly diverge at large arguments. We deduce the this behavior by other arguments (see SI).
where the proportion constant cannot be determined by our estimation. Since Eq. 17 is consistent with the Lorentzian model, measuring the deviation from this model requires an accurate measurement of the power spectrum at zero frequency or at the Larmor frequency by Eq. 12. A heuristic sketch of the power spectrum is shown in Fig. 3. Detailed derivation of the results of this section and the following are found in the SI.
corrections to the power spectrum at low drift. The calculation of the Fourier transform of the temporal correlation function with finite drift is difficult. Hence, we chose to estimate the behavior of the power spectrum directly by taking the Fourier transform of Eq. 9, where P r r ( , , 0) 0 ω is the diffusion propagator in the Fourier domain. Assuming the nuclei flow in a constant velocity v , the propagator is given by The propagator, Eq. 19, can be expanded into a series in the low drift limit ω ωω  1 v D 2 . Taking advantage of the non-Lorentzian spectrum and requiring the series expansion to converge yields further restrictions - Assuming the direction of the velocity is restricted to the x − y plane, in this regime, the expansion of the spectrum of an untilted NV will be The linear contribution in the velocity vanishes, because of the polar symmetry of the Spherical Harmonics. Since the velocity scaling is similar to the Lorentzian model, no advantage is gained in the sensitivity when the NV is untilted. However, when the NV is tilted the polar symmetry is broken and the expansion reads where C is some constant that we were unable to evaluate. Extending Eq. 22 to ω ~ ω v by taking ω = ω v results in Analytic calculation of the power spectrum in another geometry and the molecular dynamics simulations support the scaling found in Eq. 23 (see SI). A heuristic sketch of the power spectrum in the low drift limit is shown in Fig. 3b. Figure 6 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 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 v v j j ← − < >. The initialization concludes with the addition of the drift velocity v in the x-direction to each particle v xj ← v xj + v. The simulation dynamics is deterministic, following Newton's laws, and integrated using the Velocity-Verlet method with step size t m 0 00125 /( ) 2 ε σ ∆ = .

Molecular dynamics.
. 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 Scientific RepoRtS | (2020) 10:5298 | https://doi.org/10.1038/s41598-020-61095-y www.nature.com/scientificreports www.nature.com/scientificreports/ 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. Figure 5 shows the power spectra computed using the LJ fluid model for d = σ and L x = 163.08σ. The diffusion time scale τ D is computed from the relation ω τ The integration was carried out for 3.2 × 10 6 time steps, which resulted in total integration time of 18.7τ D . The value of the power spectrum at ω = 0 for different velocities, shown in the inset of the figure, was fitted using least-squares algorithm.