Introduction

The field of quantum sensing studies the precision in measuring various physical quantities using quantum protocols. It was shown that employing quantum schemes can greatly improve the precision in different problems, such as atomic clocks, magnetometry, and frequency estimation1,2,3,4,5. A typical settings in this field is when a measurable quantity; e.g., a frequency, is associated with a quantum register of which many copies are available. The registers can be individually addressed or be read by a strong measurement. When the quantum registers are qubits, and the measurable quantity is the Larmor frequency, the ultimate precision limit in a noisy environment is well established, and is achievable by performing Ramsey spectroscopy on each register6,7 (see left side of Fig. 1)

$$\Delta {\omega }_{{\rm{N}}}\ge \frac{1}{\sqrt{N}t},$$
(1)

where N is the number of copies and t is the measurement time. We henceforth refer to the equality in Eq. (1) as the standard quantum limit (SQL). Note, that usually the SQL is defined only by the scaling \(\Delta {\omega }_{{\rm{N}}}\propto \frac{1}{\sqrt{N}}\)1,4,5,8, while here it is defined as an equality. This distinction is key, since the prefactor can reduce the precision by orders of magnitude, and thefore it descriminates between non-optimal schemes to the fundamental percision limit.

Fig. 1: Possible measurement setups.
figure 1

Left — an ensemble of qubits (blue arrows); i.e., atoms, ions, superconducting circuits, with energy splitting ω that can be strongly measured (depicted above via optical access and an array of photo-detectors). The optimal precision of measuring ω, is achieved by performing Ramsey spectroscopy and is given by Eq. (1). Right — an ensemble of qubits, stores information about a desired measurable quantity; e.g., magnetic field, Larmor frequency. The desired information is accessible only through it’s weak interaction with an external quantum probe (red arrow). The interaction is illustrated by the color gradient, such that the amount of color is proportional to the interaction strength. The fundamental bound on the precision of a weakly interacting probe is unknown.

However, in some physical scenarios the registers cannot be measured directly, but can only be manipulated by global operations and measured through its weak interaction with an external probe (see right side of Fig. 1). It is unknown whether the limit (1) can be achieved in these settings, or more generally, if there is a tight bound on the precision of such measurements.

A prominent example of this scenario is nano-scale nuclear magnetic resonance (NMR), i.e., detecting the Larmor frequencies of nano samples of nuclear spins9,10,11,12,13,14,15. Since the nuclear spins cannot be measured directly, an external probe that interacts with the sample is manipulated and measured to retrieve information about the Larmor frequencies of the ensemble. The nitrogen-vacancy (NV) center is a natural candidate for an external probe, since in the past decade it has been shown to be an excellent magnetometer on the nano-scale16,17, and many successful experiments have been carried out in platforms similar to Fig. 29,10,11,12,13,14,15.

Fig. 2: Nano-NMR based on NV centers.
figure 2

The NV center (red) is situated at a depth d below the diamond surface. The nuclear spin ensemble (blue), which is located on top of the surface, is partially polarized due to the external magnetic field B. The local magnetic field at the NV’s position is only affected by the nuclear spins within a hemisphere of radius d centered above the NV’s position (dashed red line), since the dipolar interaction creates an effective cutoff. The sensor is initialized and then measured after a given time, in order to reproduce the probability distribution, which is analogous to the classic NMR signal (top right). However, as the sensor is brought into close proximity with the sample, strong back-action causes the signal to change (bottom right). This change dramatically affects the precision of the frequency estimation (the full expression of the probability distribution plotted on the right hand side is taken from Eq. (13)).

Here, we provide and analyze a protocol that in a general setting of an ensemble measurement, with a simplified interaction Hamiltonian between the ensemble and the quantum probe, achieves the SQL given by Eq. (1) (up to a small prefactor). The protocol is then extended to the physical settings of quantum NMR spectroscopy, where the interaction Hamiltonian is more involved. We show that our protocol still achieves the SQL up to a small prefactor. An open question in the field is whether spectroscopy of nuclear spins can be performed efficiently with shallow NVs, i.e., in the limit of strong back-action15,18. In our protocol the SQL is reached by utilizing the strong back-action caused by the entanglement formed between the sensor and the ensemble, and thus provides an affirmative (theoretical) answer.

Interestingly this convergence to the ultimate precision limit, despite of the weak interaction, stems from the superradiant nature of the interaction. In superradiance N atoms interact with the same probe, a single mode of radiation, and then even if the fluorescence of each atom is very weak the ultimate precision limit of a perfect projective measurement can still be reached given that N is large enough. In our case, N nuclear spins interact with the same quantum probe. The weak interaction corresponds to the weak fluorescence, and equivalently, given a large enough N, the SQL is achievable.

An alternative intuition to the behavior of the precision is given by the following reasoning. Whether the sensing is quantum or classic, the acquired signal is

$$S \sim NA\cos \left({\omega }_{{\rm{N}}}t\right),$$
(2)

where ωN is their Larmor frequency and A is the amplitude of the signal generated by a single nucleus. In the quantum case the signal is the population of the sensor, and the amplitude, A, is replaced by ϕ, the phase acquired from a single nucleus, which is a function of the interaction strength and the interrogation time. In the classical case the signal is just a classic magnetic field (i.e., measured by a current in coil). The uncertainty of the frequency measurement of the signal (2) scales inversely with the derivative and thus can be written as

$$\Delta {\omega }_{{\rm{N}}}=\frac{C}{Nt}=\frac{1}{Nt\phi },$$
(3)

where the first equality is correct for either classic or quantum sensing with some constant C, while the second is correct only for the quantum case. We henceforth refer to Eq. (3) as “Heisenberg scaling with N” or as the “weak limit” for reasons that will readily become apparent. Comparing Eqs. (3) and (1) implies that whenever \(\sqrt{N}\phi \ge 1\) the behavior of the signal must change and that the SQL might be achievable. Since the phase increases as the sensor is brought into close proximity with the sample, there is a critical distance dc, for which the behavior of the signal transforms from classical to quantum (see the right side of Fig. 2). State-of-the-art experimental setups can approach this critical distance (precise evaluations are provided in the following sections).

Results

Simplified model: constant coupling

Let us start with a simplified model that captures the essentials. The physical system consists of a quantum sensor, taken to be a two level system with energy spliting ω0, and an ensemble of N spin-1/2 nuclei with Larmor frequency ωN as described by the Hamiltonian

$${H}_{0}=\frac{{\omega }_{0}}{2}{\sigma }_{z}+\frac{{\omega }_{{\rm{N}}}}{2}\mathop{\sum }\limits_{i=1}^{N}{I}_{z}^{i},$$
(4)

where \({\sigma }_{i}/{I}_{i}^{j}\) is the Pauli matrix of the sensor/j-th nucleus in the i direction. The NV and the nuclei interact with a constant coupling

$${H}_{1}=g{\sigma }_{z}\mathop{\sum }\limits_{j = 1}^{N}{I}_{x}^{j}.$$
(5)

Our protocol for the estimation of ωN starts by initializing the system to the state \(\left|{\psi }_{i}\right\rangle ={\left|{\uparrow }_{X}\right\rangle }^{{\rm{S}}}\otimes \left|{\uparrow }_{X}\cdots {\uparrow }_{X}\right\rangle\) and the interaction is suppressed. Note, that we explicitly assume that the nuclear spins are completely polarized in the x direction. The assumptions of constant coupling, full polarizion, and interaction of the form σzIx are relaxed in the following sections.

The propagation in the interaction picture with respect to the sensor’s free Hamiltonian is

$$\left|{\psi }_{t}\right\rangle ={\left|{\uparrow }_{X}\right\rangle }^{{\rm{S}}}\otimes {\left(\cos \left(\frac{{\omega }_{{\rm{N}}}}{2}t\right)\left|{\uparrow }_{X}\right\rangle -{\rm{i}}\sin \left(\frac{{\omega }_{{\rm{N}}}}{2}t\right)\left|{\downarrow }_{X}\right\rangle \right)}^{\otimes N}.$$
(6)

Then H1 is turned on while H0 is turned off, so the propagation by an additional time τ will result in

$$\left|{\psi }_{\tau }\right\rangle =\frac{1}{\sqrt{2}}\left({\left|{\uparrow }_{Z}\right\rangle }^{{\rm{S}}}\left|{\psi }_{+}\right\rangle +{\left|{\downarrow }_{Z}\right\rangle }^{{\rm{S}}}\left|{\psi }_{-}\right\rangle \right),$$
(7)

where \(\left|{\psi }_{\pm }\right\rangle ={\left(\cos \left(\frac{{\omega }_{{\rm{N}}}}{2}t\right){{\rm{e}}}^{\mp {\rm{i}}g\tau }\left|{\uparrow }_{X}\right\rangle -{\rm{i}}\sin \left(\frac{{\omega }_{{\rm{N}}}}{2}t\right){{\rm{e}}}^{\pm {\rm{i}}g\tau }\left|{\downarrow }_{X}\right\rangle \right)}^{\otimes N}\). The reduced density matrix of the sensor after both steps is

$${\rho }_{{\rm{S}}}=\frac{1}{2}\left\{\begin{array}{*{20}{c}}1&\langle {\psi }_{-}| {\psi }_{+}\rangle \\ \langle {\psi }_{+}| {\psi }_{-}\rangle &1\end{array}\right\},$$
(8)

where the off diagonal element is given by

$$\langle {\psi }_{+}| {\psi }_{-}\rangle ={\left(\cos \left(2g\tau \right)+{\rm{i}}\sin \left(2g\tau \right)\cos \left({\omega }_{{\rm{N}}}t\right)\right)}^{N}.$$
(9)

To obtain the product it is convenient to write the complex number (9) in a polar representation

$$\Phi =N\arctan \left[\tan \left(2g\tau \right)\cos \left({\omega }_{{\rm{N}}}t\right)\right]\approx 2Ng\tau \cos \left({\omega }_{{\rm{N}}}t\right)$$
(10)
$$r={\left(1-{\sin }^{2}\left(2g\tau \right){\sin }^{2}\left({\omega }_{{\rm{N}}}t\right)\right)}^{N/2}\approx {{\rm{e}}}^{-2N{\left(g\tau \right)}^{2}{\sin }^{2}\left({\omega }_{{\rm{N}}}t\right)},$$
(11)

where the approximations are made using the assumption of weak coupling, gτ 1. We henceforth refer to the accumulated phase, Φ, as “the signal” since it corresponds to the “classical” signal (2) with the aforementioned accumulated phase ϕ = 2gτ, and we refer to r as the decay. The decay is caused by the entanglement between the NV and nuclear spin ensemble, as illustrated in Fig. 3. For times ωNt = nπ the dynamics is purely classical because the rotation around the x-axis is trivial, whereas for other times the rotation causes the sensor state to be entangled with the collective state of the ensemble. Therefore, the dimensionless quantity characterizing the decay, \(N{\left(g\tau \right)}^{2}\), is interpreted as the back-action. As the coupling constant g depends on the distance between the sensor and the sample, the transition from weak to strong back-action is dictated by the NV’s depth for a given τ. This transition is the one predicted in the introduction, as it occurs when \(\phi =2gt \sim \frac{1}{\sqrt{N}}\). Substituting Eqs. (10) and (11) into Eq. (9),

$$\langle {\psi }_{+}| {\psi }_{-}\rangle =r{{\rm{e}}}^{{\rm{i}}\Phi }\approx {{\rm{e}}}^{-2N{\left(g\tau \right)}^{2}{\sin }^{2}\left({\omega }_{{\rm{N}}}t\right)}{{\rm{e}}}^{{\rm{i}}2Ng\tau \cos \left({\omega }_{{\rm{N}}}t\right)}.$$
(12)

The measurable quantity associated with (12) is the probability distribution,

$${P}_{\left|{\uparrow }_{Y}\right\rangle }=\frac{1}{2}\left(1+{{\rm{e}}}^{-2{g}^{2}N{\tau }^{2}{\sin }^{2}({\omega }_{{\rm{N}}}t)}\sin (2gN\tau \cos ({\omega }_{{\rm{N}}}t))\right),$$
(13)

which is related to Eq. (12) by \({P}_{\left|{\uparrow }_{Y}\right\rangle }=\frac{1}{2}\left(1-{\rm{Im}}\left[\langle {\psi }_{+}| {\psi }_{-}\rangle \right]\right)\). The function (13) is the non-approximate form of Eq. (2). It is drawn on the right hand side of Fig. 2 with N = 3 105 and gτ = 10−6 (top) or gτ = 0.01 (bottom), which correspond to weak and strong back-action, respectively.

Fig. 3: The protocol.
figure 3

The NV and the nuclear spins are initialized at the x direction. The nuclei are then allowed, using an appropriate pulse sequence, to propagate according to the free Hamiltonian H0 for a time t, which results in a rotation by θ = ωNt around the z-axis. Then, by changing the external pulse sequence, the system propagates under the interaction Hamiltonian H1 for a duration of τ. This results in a rotation of the nuclei around the x-axis by an angle of ±ϕ = ± 2gτ depending on the sensor’s state. The sensor will experience an effective dephasing depending on the extent of its entanglement with the ensemble.

In order to determine the precision of the estimation of ωN we use the tools of quantum metrology. The Fisher information about a parameter g given a discrete distribution \({\left\{{p}_{i}\left(g\right)\right\}}_{i}\) is defined as \(I={\sum }_{i}\frac{{\left(\frac{d{p}_{i}}{dg}\right)}^{2}}{{p}_{i}}\)19. For a quantum system one can optimize over all possible measurement bases. This leads to the definition of quantum fisher information (QFI)1. For a density matrix \(\rho ={\sum }_{i}{p}_{i}\left(g\right)\left|{\Psi }_{i}\left(g\right)\right\rangle \left\langle {\Psi }_{i}\left(g\right)\right|\), the QFI about g is given by \({\mathcal{I}}={\sum }_{{p}_{i}+{p}_{j}\ne 0}\frac{2}{{p}_{i}+{p}_{j}}{\left|\left\langle {\Psi }_{j}\right|\frac{d\rho }{dg}\left|{\Psi }_{i}\right\rangle \right|}^{2}\)6. The precision of any measurement is bounded by the Cramer-Rao bound, \(\Delta g\ge \frac{1}{\sqrt{{\mathcal{I}}}}\). Since this is a tight bound, we use it henceforth to quantify precision. For the density matrix (8) the QFI can be expressed as the relevant Bures distance6,20:

$${\mathcal{I}}=\frac{{\left(\frac{dr}{d{\omega }_{{\rm{N}}}}\right)}^{2}}{1-{r}^{2}}\,+{r}^{2}{\left(\frac{d\Phi }{d{\omega }_{{\rm{N}}}}\right)}^{2},$$
(14)

thus, Eqs. (10), (11), and (14) lead to

$$\begin{array}{r}{\mathcal{I}}={\phi }^{2}{N}^{2}{t}^{2}{\sin }^{2}(\theta ){{\rm{e}}}^{-{\phi }^{2}N{\sin }^{2}(\theta )}\left[1+\frac{{\phi }^{2}{\cos }^{2}(\theta )}{1-{{\rm{e}}}^{-{\phi }^{2}N{\sin }^{2}(\theta )}}\right],\end{array}$$
(15)

where θ = ωNt. In the limit of weak back-action, \(N{\left(g\tau \right)}^{2}\ll 1\),

$${P}_{\left|{\uparrow }_{Y}\right\rangle }=\frac{1}{2}\left(1+\sin (N\phi \cos ({\omega }_{{\rm{N}}}t))\right),$$
(16)
$${\mathcal{I}}\approx {\phi }^{2}{N}^{2}{t}^{2}{\sin }^{2}\left(\theta \right)\left[1+\frac{{\cos }^{2}(\theta )}{N{\sin }^{2}(\theta )}\right].$$
(17)

The QFI in this limit is optimal when

$${\sin }^{2}\left({\theta }_{{\rm{opt}}}\right)=1,$$
(18)

since the derivative of the signal is then maximal and the decay is negligible. For the optimal time (18), the QFI (15) is

$${\mathcal{I}}={\phi }^{2}{N}^{2}{t}^{2},$$
(19)

which corresponds to the uncertainty \(\Delta \omega =\frac{1}{2Ng\tau t}\) as in Eq. (3) with ϕ = 2gτ. Thus, in the limit of weak back-action, we achieve an uncertainty that scales as N−1 in exchange for the large factor \({\left(g\tau \right)}^{-1}\).

Achieving better precision, however, requires Ng2τ2  1, where consequently, the decay starts to affect the quantum sensor. The optimal time (18), therefore, must change to account for the decay. The optimal time in this regime is approximately

$${\sin }^{2}\left({\theta }_{{\rm{opt}}}\right)=\frac{1}{N{\phi }^{2}},$$
(20)

for which the QFI (15) is

$${\mathcal{I}}\approx\, \frac{N{t}^{2}}{{\rm{e}}}$$
(21)

and corresponds to the uncertainty \(\Delta \omega =\frac{1}{t}\sqrt{\frac{{\rm{e}}}{N}}\). The decay, therefore, “corrects” the apparent Heisenberg scaling for high precision.

We note that the first step of the protocol, where the nuclei propagate according to H0, can be implemented by initializing the sensor in an eigenstate of SZ = 0, which eliminates H1 or by applying an external drive on the NV that suppresses the interaction. Therefore, the uncertainty is limited solely by the coherence time of the nuclei as in refs. 13,15,21,22,23,24 and it does not depend on the coupling constant g. Moreover, if we could manipulate and read-out each nucleus with a unit fidelity, the optimal QFI would be the SQL (1). Therefore, the presented protocol achieves the optimal precision, up to the numerical factor of e−1. Figure 4 compares the different QFI scalings to our protocol.

Fig. 4: Saturating the fundamental limit.
figure 4

The maximal value of the QFI (15) (red) for gτ = 0.01 and t = 1 as a function of N (log scale). The optimal QFI follows the weak back-action limit (green) given by Eq. (19) when N(gτ)2 1. As N increases the back-action corrects the scaling until N(gτ) > 1, where the QFI approaches the strong limit given by Eq. (21). In the strong back-action regime, the QFI of our scheme is only smaller by a factor of e−1 from the optimal QFI, Nt2, achieved by Ramsey spectroscopy (blue).

We show in Supplementary Note 2 that this behavior of the QFI is ubiquitous in supperradiant measurements. A similar behavior can be observed with N atoms that weakly interact with a single mode of radiation: given a large enough N the ultimate precision limit can be obtained despite the weak interaction. Hence this behavior is due to the superradiant nature of the interaction in Eq. (5).

To obtain this QFI we assume optimization of θ and the measurement basis, which requires knowledge of ωN. As shown in Fig. 5, the QFI increases with N(gτ)2, but it also becomes increasingly narrow, such that a more accurate estimation of ωN is required to achieve it. This, however, can be resolved by using an adaptive measurement protocol — a sequence of non-optimal measurements are performed to acquire an estimate of ωN, and the details of the next measurements are updated according to the outcomes of the previous measurements. This can be repeated until the estimate of ωN is good enough to attain optimal precision. An exmple of such a protocol is given in Supplementary Note 8 and further analysis of this model using multiple sensors can be found in Supplementary Note 7.

Fig. 5: Narrowing of the QFI for strong back-action.
figure 5

The QFI (15) for t = 1 and N = 105 for different values of gτ. As the back-action N(gτ)2 gets stronger the peak of the QFI increases until it reaches an optimum since it no longer depends on the interaction. The increase in QFI comes at the cost of it becoming increasingly narrow. Therefore, to achieve the optimal scaling (21) a good estimate of ωN is already required. This can be attained by using an adaptive protocol as explained in SI Note 8.

Spatially dependent coupling

In the previous section, we presented the simplified model with the interaction Hamiltonian (5). This is an approximate form of the dipole–dipole interaction Hamiltonian, which can be written as

$${H}_{{\mathrm{DD}}}={\sum \limits_{k = 1}^{N}}\,{\sum \limits_{i,j = \{x,y,z\}}}{\tilde{g}}_{i,j}^{k}{\sigma }_{i}{I}_{j}^{k}.$$
(22)

In the interaction picture with respect to sensor’s free Hamiltonian, after taking the rotating wave approximation, the remaining terms are

$${H}_{{\mathrm{DD}}}^{I}\approx {\sigma }_{z}{\sum \limits_{i = 1}^{N}}\left({g}_{0}^{i}{I}_{z}^{i}+{g}_{1}^{i}{I}_{x}^{i}+{g}_{2}^{i}{I}_{y}^{i}\right).$$
(23)

If we further assume that the nuclei can be driven in the x direction sufficiently fast compared to the interaction and the entanglement generation rate, we arrive at the Hamiltonian

$${\tilde{H}}_{1}={\sigma }_{z}{\sum \limits_{i = 1}^{N}}{g}^{i}{I}_{x}^{i}.$$
(24)

The difference between Eq. (24) and the toy model (5) is that the coupling differs from one nucleus to the other. It depends, in fact, on the position of the nucleus and, therefore, also on time. Further description of the interaction requires a specific realization of the sensor, henceforth we use the NV center. In a spherical coordinate system, where the NV is found at the origin, the NV’s magnetization axis coincides with the z-axis and the i-th nucleus position is denoted by \(\left\{{r}_{i},{\theta }_{i},{\varphi }_{i}\right\}\), the coefficients gi are given by

$${g}^{i}\left(t\right)\equiv g\left({{\bf{r}}}_{i}\left(t\right)\right)=-\frac{3}{2}J{\left[{r}_{i}\left(t\right)\right]}^{-3}\sin \left[{\theta }_{i}\left(t\right)\right]\cos \left[{\theta }_{i}\left(t\right)\right]\cos \left[{\varphi }_{i}\left(t\right)\right],$$
(25)

with the physical coupling constant \(J=\frac{{\mu }_{0}\hslash {\gamma }_{{\rm{e}}}{\gamma }_{{\rm{N}}}}{4\pi }=0.49\,{\rm{MHz}}\cdot {{\rm{nm}}}^{3}\), where μ0 is the vacuum permeability, is the reduced Plank constant and γe/N are the electronic/nuclear gyro-magnetic ratios.

Repeating the derivation of the previous section with the interaction (24), Eq. (9) becomes

$$\langle {\psi }_{+}| {\psi }_{-}\rangle ={\prod \limits_{j = 1}^{N}}\left(\cos \left(2{G}_{j}\right)+{\rm{i}}\sin \left(2{G}_{j}\tau \right)\cos \left({\omega }_{{\rm{N}}}t\right)\right),$$
(26)

where \({G}_{j}={\int \nolimits_{0}^{\tau }}dt{g}_{j}\left(t\right)\). In the limit of weak coupling, Gi 1, Eq. (26) can be approximated to

$$\langle {\psi }_{+}| {\psi }_{-}\rangle ={{\rm{e}}}^{-2\mathop{\sum }\nolimits_{i = 1}^{N}{G}_{i}^{2}{\sin }^{2}\left({\omega }_{{\rm{N}}}t\right)}{{\rm{e}}}^{2{\rm{i}}\mathop{\sum }\nolimits_{i = 1}^{N}{G}_{i}\cos \left({\omega }_{{\rm{N}}}t\right)}.$$
(27)

The QFI, therefore, depends only on \({s}_{1}={\sum }_{i}{G}_{i},\,{s}_{2}={\sum }_{i}{G}_{i}^{2}\) and it is approximately

$${\mathcal{I}}\approx 4{t}^{2}{s}_{1}^{2}\sin {\left(\theta \right)}^{2}\exp \left(-4{s}_{2}{\sin }^{2}\left(\theta \right)\right).$$
(28)

The quantity that determines the behavior of the optimal QFI is s2:

$${\mathcal{I}}\approx \left\{\begin{array}{ll}4{t}^{2}{s}_{1}^{2}&{s}_{2}\ll 1,\,{\sin }^{2}\left(\theta \right)=1\\ {t}^{2}\frac{1}{{\rm{e}}}\frac{{s}_{1}^{2}}{{s}_{2}}&{s}_{2}\gg 1,\,{\sin }^{2}\left(\theta \right)=\frac{1}{4{s}_{2}}\end{array}\right..$$
(29)

The regime of interest is s2 1, where, as we readily show, a modified SQL scaling can be achieved. Note that \({\mathcal{I}}\le {t}^{2}N/{\rm{e}},\) where equality is attained only for homogeneous couplings.

To obtain the QFI in a nano-NMR scenario, we need, therefore, to find the relevant s1, s2. The results will of course depend on the parameters of the problem. We denote by \(n=\frac{N}{V}\) the nuclei number density, D the diffusion coefficient of the sample, d the depth of the NV, and α the NV’s tilting angle, measured between the normal to the diamond surface and its magnetization axis. Let us first calculate s1 assuming N 1,

$${s}_{1}={\gamma }_{{\rm{e}}}\left\langle B\right\rangle \tau =n\tau \int\,{d}^{3}rg\left({\bf{r}}\right)=-\pi nJ\sin \left(2\alpha \right).$$
(30)

While s1 is proportional to the average magnetic field, it can be observed that s2 goes as the magnetic field auto-correlation,

$$\begin{array}{*{20}{l}}{s}_{2}={\gamma }_{{\rm{e}}}^{2}{\int \limits_{0}^{\tau }}dt {\int \limits_{0}^{\tau }}d{t}_{0}\int\,{d}^{3}r{d}^{3}{r}_{0}g\left({\bf{r}}\right)g\left({{\bf{r}}}_{0}\right)P\left({\bf{r}},{{\bf{r}}}_{0},t-{t}_{0}\right)\\ ={\gamma }_{{\rm{e}}}^{2}{\int \limits_{0}^{\tau}}dt {\int \limits_{0}^{\tau }}d{t}_{0}\left\langle B\left(t\right)B\left({t}_{0}\right)\right\rangle \end{array}$$
(31)

where \(\left\langle \cdot \right\rangle\) is an average over realizations and P(r, r0, t) is the stationary diffusion propagator from r, to r0 with time difference t. The correlation decays with a characteristic time of \({\tau }_{{\rm{D}}}\equiv \frac{{d}^{2}}{D}\), which dictates the behavior of s2. Although the full expression of s2 is involved25, the asymptotic behavior is rather simple,

$${s}_{2}=\left\{\begin{array}{ll}{\gamma }_{e}^{2}{B}_{{\rm{rms}}}^{2}{\tau }^{2}&\tau \ll {\tau }_{{\rm{D}}}\\ {\gamma }_{e}^{2}{B}_{{\rm{rms}}}^{2}\tau {\tau }_{{\rm{D}}}&\tau \gg {\tau }_{{\rm{D}}}\end{array}\right.,$$
(32)

where \({B}_{{\rm{rms}}}^{2}\) is the instantaneous fluctuation in the magnetic field:

$$\begin{array}{l}{\gamma }_{{\rm{e}}}^{2}{B}_{{\rm{rms}}}^{2}={\gamma }_{{\rm{e}}}^{2}\langle {B}^{2}\rangle =\int\,{d}^{3}r\ g{\left({\bf{r}}\right)}^{2}\\ =\pi n{J}^{2}\left(\frac{35-3\cos \left(4\alpha \right)}{256{d}^{3}}\right).\end{array}$$
(33)

The mean field, (30), and Brms, (33), can be estimated more explicitly,

$${\gamma }_{{\rm{e}}}\left|\left\langle B\right\rangle \right|\approx \frac{n}{{n}_{{\rm{water}}}}48{\rm{MHz}},$$
(34)
$${\gamma }_{{\rm{e}}}^{2}{B}_{{\rm{rms}}}^{2}\approx \frac{n}{{n}_{{\rm{water}}}}{\left(44.5{\rm{KHz}}\right)}^{2}{\left(\frac{10{\rm{nm}}}{d}\right)}^{3},$$
(35)

where we took α = 54. 7 and nwater = 33 nm−3 as the water number density. We use these quantities for the remainder of the article in our quantitative estimates.

The validity condition of the second order approximation in Gi can be put into physical terms with the definitions (30) and (33); e.g., Eqs. (27) and (28) are valid when \(\frac{{B}_{{\rm{rms}}}}{\left|\left\langle B\right\rangle \right|}\ll 1\), which means that the polarization should be larger than the statistic polarization. Hence, for increasingly shallow NVs, higher orders should be taken into account (see discussion in Supplementary Note 3).

We are now fully equipped to calculate the optimal QFI by substituting s1 and s2 (Eqs. (30) and (32), respectively) into Eq. (29). In the weak back-action regime, s2 1, the optimal QFI reads:

$${\mathcal{I}}=4{\gamma }_{{\rm{e}}}^{2}{\left\langle B\right\rangle }^{2}{\tau }^{2}{t}^{2}.$$
(36)

For strong back-action the optimal QFI depends explicitly on τD,

$${\mathcal{I}}\approx \left\{\begin{array}{ll}\frac{{\left\langle B\right\rangle }^{2}}{{\rm{e}}{B}_{{\rm{rms}}}^{2}}{t}^{2}\approx 17.5\frac{n{d}^{3}}{{\rm{e}}}{t}^{2}&\tau \ll {\tau }_{{\rm{D}}}\\ \frac{{\left\langle B\right\rangle }^{2}}{{\rm{e}}{B}_{{\rm{rms}}}^{2}}\frac{\tau }{{\tau }_{{\rm{D}}}}{t}^{2}\approx 17.5\frac{n{d}^{3}}{{\rm{e}}}\frac{\tau }{{\tau }_{{\rm{D}}}}{t}^{2}&\tau \gg {\tau }_{{\rm{D}}}\end{array}\right.,$$
(37)

where the optimal θ is given by

$${\sin }^{2}\left(\theta \right)\approx \left\{\begin{array}{ll}\frac{1}{4{\left({\gamma }_{{\rm{e}}}{B}_{{\rm{rms}}}\tau \right)}^{2}}&\tau \ll {\tau }_{{\rm{D}}}\\ \frac{1}{4{\left({\gamma }_{{\rm{e}}}{B}_{{\rm{rms}}}\right)}^{2}\tau {\tau }_{{\rm{D}}}}&\tau \gg {\tau }_{{\rm{D}}}\end{array}\right..$$
(38)

Note that the transition from Eq. (36) to Eq. (37) occurs when \(\frac{2}{\pi }{\gamma }_{{\rm{e}}}{B}_{{\rm{rms}}}{T}_{2}^{{\rm{NV}}}=1\), which defines the critical depth dc in terms of the physical parameters. The addional factor of π−1 derives from the dynamical decoupling sequence applied to the NV in order to produce H1 (see next section or a detailed derivation in Supplementary Note 4). For common parameters of liquids, an interrogation time of \({T}_{2}^{{\rm{NV}}}=1\,{\rm{ms}}\) and a fully polarized nuclear spin ensemble, we find that dc ~ 130 nm.

The result (36) and the short times limit of (37) are similar to those of the simplified model, with the difference that due to the dipolar interaction and the geometry, there is an additional dependence on the NV’s tilting angle. In Eq. (37) this is translated into a change in the total number of particles in Eq. (21) by an effective number N ≈ 17.5nd3, which is proportional to the number of particles in the effective interaction region.

The long times regime of Eq. (37) may look puzzling since for large enough τ we can get an arbitrarily large QFI, which seems paradoxical. This arbitrarily large QFI is due to the fact that we consider an infinite sample volume; i.e., an infinite amount of nuclear spins. Restricting ourselves to a finite volume, V, with N nuclei yields the restriction \(\frac{{\left\langle B\right\rangle }^{2}}{{\rm{e}}{B}_{{\rm{rms}}}^{2}}\frac{\tau }{{\tau }_{{\rm{D}}}}{t}^{2}\le N{t}^{2}\) by imposing the SQL. This limitation on τ can be taken to be ττV (see Supplementary Note 3), where τV = V2/3/D is the volumetric diffusion time; i.e., the characteristic time it takes for a particle to move from one of the volume’s boundaries to another. For sufficiently long times, ττV, the QFI scaling changes to

$${\mathcal{I}}\approx\, \frac{nV{t}^{2}}{{\rm{e}}}=\frac{N}{{\rm{e}}}{t}^{2}.$$
(39)

This is the equivalent of Eq. (21). In the simplified model the sensor is coupled equally to all the nuclei from the start; therefore, the QFI scales with the total number of nuclei, while in reality the dipolar interaction creates an effective cutoff, such that in short times the QFI scales as nd3 and only after a long time, ττV, when all the nuclei have passed through the interaction region, does the QFI scale with the total number of nuclei. As in Eq. (21), in this model the QFI is not limited by the coherence time of the NV, \({T}_{2}^{{\rm{NV}}}\), but only by the coherence time of the nuclei \({T}_{2}^{{\rm{N}}}\), which is usually longer by orders of magnitude.

Undriven nuclei

In the previous section, we assumed that we can drive the nuclei sufficiently fast to achieve the approximate dynamics given by Eq. (24). In some experimental settings this is unfeasible for practical reasons. Hence, in what follows, we drop this assumption. The approximate dipolar Hamiltonian (23) is rewritten as

$${H}_{1}={\sigma }_{z} {\sum \limits_{i = 1}^{N}}\left({g}_{0}^{i}{I}_{z}^{i}+{g}_{1}^{i}{I}_{+}^{i}+{g}_{2}^{i}{I}_{-}^{i}\right).$$
(40)

The NV is then driven with π−pulses every time τp so that the effective pulse frequency \({\omega }_{{\rm{p}}}=\frac{\pi }{{\tau }_{{\rm{p}}}}\) is close to ωN. This yields the effective Hamiltonian (see Supplementary Note 4)

$${H}_{1}\approx \frac{2{g}_{\pm }^{i}}{\pi }{\sigma }_{z}\left[{I}_{x}^{i}\sin \left({\varphi }_{i}\right)-{I}_{y}^{i}\cos \left({\varphi }_{i}\right)\right],$$
(41)

and the remaining free Hamiltonian

$${H}_{0}=\frac{\delta \omega }{2}{\sum \limits_{i = 1}^{N}}{I}_{z}^{i}$$
(42)

where δω = ωN − ωp and \({g}_{\pm }^{i}=-\frac{3}{2}J{\left[{r}_{i}\left(t\right)\right]}^{-3}\sin \left[{\theta }_{i}\left(t\right)\right]\cos \left[{\theta }_{i}\left(t\right)\right]\). Following the same protocol as before (see Supplementary Note 5) the signal is

$$\begin{array}{r}\Phi =\frac{2}{\pi }{\gamma }_{e}\left\langle B\right\rangle \tau \sin \left(\delta \omega t\right),\end{array}$$
(43)

where the average magnetic field is as in the driven case, (30). The decay, however, changes to

$$r=\exp \left\{-{\gamma }_{{\rm{e}}}^{2}\left[{\left({B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right)}^{2}{\cos }^{2}\left(\delta \omega t\right)+{\left({B}_{{\rm{rms}}}^{{\rm{UDC}}}\right)}^{2}\right]{\tau }^{2}\right\},$$
(44)

for times ττD, where

$${\left({B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right)}^{2}=\frac{4}{{\pi }^{2}}{B}_{{\rm{rms}}}^{2}\left(\frac{2{\sin }^{2}(\alpha )(-13\cos (2\alpha )+1)}{13\cos (4\alpha )+51}\right),$$
(45)
$${\left({B}_{{\rm{rms}}}^{{\rm{UDC}}}\right)}^{2}=\frac{4}{{\pi }^{2}}{B}_{{\rm{rms}}}^{2}\left(\frac{28\cos (2\alpha )+13\cos (4\alpha )+87}{52\cos (4\alpha )+204}\right)$$
(46)

are the quantum and classical contributions to the instantaneous fluctuation of the magnetic field and Brms is given by Eq. (33).

Comparing Eqs. (43) and (44) to the results of the previous section, Eqs. (30) and (32), note that the signal is the same upto a prefactor of order 1 and a constant known phase. The decay, however changes, so it is non-zero for any given time t. This is expected, since previously the interaction caused all the nuclei to rotate around the same axis, which resulted in an entanglement induced decay. Without the external drive, the interaction causes each nuclear spin to rotate around a different axis in the xy plane of the Bloch sphere, which induces additional classical dephasing (see Supplementary Figure 3). It is worth pointing out that the quantum contribution (45) can be negative, but the total Brms is always real and positive. When the decay is small, γeBrmsτ 1, we retrieve the result (36) up to a prefactor ~ 1, since it only depends on the signal.

If the decay is dominant, the classical and quantum decay compete when optimizing the QFI. On the one hand the strong back-action regime requires \({\gamma }_{{\rm{e}}}\left|{B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right|\tau >1\), which implies that τ has to be large enough. On the other hand, the classical dephasing causes an exponential decrease in the QFI, \({{\rm{e}}}^{-{\left({\gamma }_{{\rm{e}}}{B}_{{\rm{rms}}}^{{\rm{UDC}}}\tau \right)}^{2}}\), that depends only on τ. To limit this effect we require \({\gamma }_{{\rm{e}}}{B}_{{\rm{rms}}}^{{\rm{UDC}}}\tau <1\). Hence, when the nuclei are undriven fine tuning of τ is required in order to optimize these competing processes. The optimal protocol will no longer be universal and will depend on the physical parameters.

First, we assume strong back-action, in order to derive the analog protocol to the previous sections, and to emphasis our last statement regarding the decay. The optimal t satisfies

$${\cos }^{2}\left(\delta \omega t\right)=\frac{1}{{\gamma }_{{\rm{e}}}^{2}\left|{\left({B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right)}^{2}\right|{\tau }^{2}},$$
(47)

for which the QFI is

$${\mathcal{I}}\approx \frac{{\left\langle B\right\rangle }^{2}{t}^{2}}{{{\rm{e}}}^{{\rm{sgn}}\left({({B}_{{\rm{rms}}}^{{\rm{UDQ}}})}^{2}\right)}\left|{\left({B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right)}^{2}\right|}{{\rm{e}}}^{-{\gamma }_{{\rm{e}}}^{2}{\left({B}_{{\rm{rms}}}^{{\rm{UDC}}}\right)}^{2}{\tau }^{2}}.$$
(48)

As aforementioned, the QFI, (48), has a reduction by an exponential factor compared to the previous result (37), arising from the classical dephasing. As long as \(\left|{\left({B}_{{\rm{rms}}}^{{\rm{UDC}}}\right)}^{2}\right|\le \left|{\left({B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right)}^{2}\right|\) we can satisfy the strong back-action requirement, while limiting the classical decay by taking \({\tau }^{2}=\frac{1}{{\gamma }_{{\rm{e}}}^{2}\left|{\left({B}_{{\rm{rms}}}^{{\rm{UDC}}}\right)}^{2}\right|}\). This results in a reduction by a reasonable factor of e−1. This inequality is a function of the NV’s tilting angle alone, and is correct for \(\alpha \ge \frac{\pi }{3}\). Since the NV’s natural tilting angle is α ~ 0.32π this is approximately the case, where the slight deviation will lead to further reduction of the QFI by a small factor. Another issue with this strategy is that it is not guaranteed that \(\tau =\frac{1}{{\gamma }_{{\rm{e}}}{B}_{{\rm{rms}}}^{{\rm{UDC}}}}\le {T}_{2}^{{\rm{NV}}}\), and indeed for a typical number density and an NV depth of d = 140nm, the time required is \(\tau \approx 3\,{\rm{ms}}>{T}_{2}^{{\rm{NV}}}\).

Obtaining other scalings is possible by optimizing t and τ, depending on the parameters; for example, by choosing \({\cos }^{2}\theta =1,\,{\tau }^{2}=\frac{1}{{\gamma }_{{\rm{e}}}^{2}\left[{\left({B}_{{\rm{rms}}}^{{\rm{UDQ}}}\right)}^{2}+{\left({B}_{{\rm{rms}}}^{{\rm{UDC}}}\right)}^{2}\right]}\) we achieve

$${\mathcal{I}}=\frac{{\left\langle B\right\rangle }^{2}}{{\rm{e}}{B}_{{\rm{rms}}}^{2}}{t}^{2},$$
(49)

which is the same result as in Eq. (37). Note that this is the case of “critical back-action”, when γeBrmsτ = 1, therefore smaller τ will lead to the weak back-action limit and larger τ will cause an exponential decrease. Since the same typical parameters yield τ ≈ 2 ms, the lack of external drive reduces Eq. (37) by a factor of e−1. The results for ττD can be derived by using reasoning similar to the previous section. However, the long times limit of Eqs. (37) and (39) will typically no longer be achievable because the optimal τ will be much longer than \({T}_{2}^{{\rm{NV}}}\).

Partial polarization

The initial state of partially polarized nuclear spins is \({\left[\frac{1}{2}\left({\mathbb{1}}+p{I}_{X}\right)\right]}^{\otimes N},\) where p is the polarization (−1 ≤ p ≤ 1). Owing to this partial polarization, the ultimate precision limit is no longer \({\mathcal{I}}=N{t}^{2}\) (see Supplementary Methods 1), but

$${\mathcal{I}}=N{p}^{2}{t}^{2},$$
(50)

Hence the QFI is degraded by a factor of p2. We thus wish to inquire whether we can approach this limit using our protocol.

The effect of partial polarization is somewhat similar to the effect of the undriven case. The polarized nuclei behave as in the driven case and induce a quantum dephasing on the external probe, while the unpolarized dynamics, creates classical dephasing. The QFI for ττD is given by (see Supplementary Note 6)

$${\mathcal{I}}\approx 4{\gamma }_{{\rm{e}}}^{2}{\tau }^{2}{t}^{2}{\sin }^{2}(\theta ){{\rm{pol}}}^{2}\cdot {e}^{-4{\gamma }_{{\rm{e}}}^{2}{B}_{rms}^{2}{\tau }^{2}\left[\left(1-{{\rm{pol}}}^{2}\right)+{{\rm{pol}}}^{2}\cdot {\sin }^{2}\left(\theta \right)\right]}{\left\langle B\right\rangle }^{2},$$
(51)

where \({\rm{pol}}=\left|p\right|\) and \(\left\langle B\right\rangle\), Brms are given by Eqs. (30) and (33), respectively. Therefore, the results of the previous section also apply for a partially polarized driven ensemble with minor modificatios. The optimal strategy in the strong back-action regime is to take \({\sin }^{2}\left(\theta \right)=1\) and \({\tau }^{2}=\frac{1}{4{\gamma }_{{\rm{e}}}^{2}{B}_{{\rm{rms}}}^{2}}\), which yeilds a QFI of

$${\mathcal{I}}=\frac{{{\rm{pol}}}^{2}}{{\rm{e}}}{\left(\frac{\left\langle B\right\rangle }{{B}_{{\rm{rms}}}}\right)}^{2}{t}^{2}.$$
(52)

Therefore, even with finite polarization our protocol achieves the ultimate precision limit, as long as the polarization is greater than the statistical polarization. When the dephasing caused by the back-action is larger than the one inflicted by the unpolarized dynamics, other approaches can achieve the optimal QFI (see Supplementary Note 6). In that case it is required that pol ≥ 70% and the critical distance changes to \({d}_{c} \sim {{\rm{pol}}}^{2/3}{T}_{2}^{2/3}\), such that for pol = 70% it drops from 130nm to dc ~ 100 nm for T2 = 1 ms. This fine tuning might yield a certain advantage over the optimal time parameters provided above, since the QFI depends explicitly on the optimal time.

Discussion

We provide a protocol for nano-NMR that achieves the SQL up to a prefactor in the strong back-action regime. Moreover, the uncertainty does not depend on the coherence time of the sensor or the probe-nuclear coupling. Our anlysis implies that in the strong back-action regime, as long as the polarzation is greater than the statistical polarization, the optimal precision is achived by preforming a single measurment and not by using sequential measurement schemes. These features make it highly applicable for sensing very small samples, which is the ultimate goal of the field.