Optimal frequency measurements with quantum probes

Optimal frequency discrimination using a quantum probe

As a diagnostic tool, NMR can be used to answer “yes–no” questions such as whether a certain toxin or metabolite is present in the sample. As sketched in Fig. 1a, the task is then to discriminate between two known spectra based on their frequency components. We define discrimination error as the error to decide on the wrong spectrum and our goal is to obtain a minimal discrimination error or equivalently a minimal discrimination time. A typical method compromises sampling the signal with consecutive, synchronized measurements and correlating the individual outcomes, e.g., by applying Fourier analysis, to obtain a spectrum. For sufficient recording time, the resolution can be high enough for an almost error-free discrimination. To illustrate, consider a simplified problem in which one wishes to discriminate between two spectra, each containing only a single frequency (ω1 or ω2) with the same amplitude B (in units of angular frequency). Naively, the method described above is Fourier limited, i.e., the time required for discrimination is \(T = \frac{1}{{|\omega _{\Delta}|}}\), where ωΔ = ω2 − ω1. Using a more sophisticated data analysis, such as Bayesian interference or machine learning, which can be applied for known B, the discrimination time lowers to \(T \sim \frac{1}{{B^{2/3}\omega _{\mathrm{{\Delta}}}^{2/3}}}\) (see Supplementary Note 1 and ref. 20) However, we show in the following that, given a sufficient coherence time, the discrimination time can be further reduced to:

$$T = \frac{\pi }{{2\sqrt B \sqrt {|\omega _{\mathrm{{\Delta}}}|} }}$$

(1)

The key idea is to drive the sensor such that the angle between the states, \(|\psi \left( {\omega _1} \right)\rangle\) and \(|\psi \left( {\omega _2} \right)\rangle\), is maximal, to guarantee a minimal error probability. Once orthogonality is achieved, quantum projection noise can theoretically be eliminated by measuring in the appropriate basis and it is possible to determine, with a certainty limited by readout fidelity, which frequency is present. Assuming a perfect readout, we show that this method is in general optimal, even if orthogonality cannot be achieved.

a The sensor state is tailored to depend on each frequency so that the sensor is driven in orthogonal directions. Readout of the state provides a “Yes/No” answer. The task is to correctly identify the frequency with highest fidelity in the shortest time. b Phase accumulation for frequency discrimination. For each Hamiltonian \({\cal{H}}_i\left( t \right) = B\;{\mathrm{sin}}\left( {\omega _i{\mathrm{t}} + \theta } \right)\sigma _{\mathrm{Z}} = {\mathrm{H}}_i\left( t \right)\sigma _{\mathrm{Z}}\), the sensor accumulates a different phase. Left: without control, the phase difference \(\varphi _{\Delta}\) oscillates and increases only slowly. Right: optimal control in this scenario implies applying π-pulses whenever the sign of H1–H2 changes, leading to a monotone phase increase with t2. c Geometrical picture of the optimal protocol. The states of the sensor can be thought of as two runners, where the aim is to maximize the gap between them, which is equivalent to the angle between the states. The speed of each runner is proportional to \({\mathrm{H}}_i\left( t \right)\). As soon as the initially slower runner becomes faster, we change their direction of motion to ensure an increasing gap. Note that the sketch just demonstrates the idea and the runners actually move on circular orbits. d Error probability as a function of time for three different strategies: optimal control and a single measurement (solid blue), this strategy is optimal and sets the fundamental error limit (see Eqs. 3 and 4). Correlated measurements (dashed red), in this illustration measurement is applied every \(\frac{{2\pi }}{{\omega _1 + \omega _2}}\). No control (dotted green). More details on this comparison are found in Supplementary Note 1.

Full size image

Given free evolution of the sensor, the discrimination time is \(T \sim \frac{\omega }{{B\omega _{\mathrm{{\Delta}}}}}\) (see Supplementary Note 1), which is basically the same scaling as the Fourier limit. As sketched in Fig. 1b, the phases \(\varphi (\omega _1)\) and \(\varphi (\omega _2)\), accumulated by the sensor under ω1 and ω2, move apart and get closer to each other such that their difference \(\varphi _{\mathrm{{\Delta}}} = \varphi (\omega _1) - \varphi (\omega _2)\) oscillates and only slowly increases. However, T can be significantly reduced by applying a suitable control, which is also shown in Fig. 1b: whenever the states start to get closer to each other (φΔ reduces), a control π-pulse can be used to change the direction of motion such that their distance increases instead. Formally, the distance between two states can be described by the angle \(\varphi _{\mathrm{{\Delta}}}/2\,{\underline{\underline {{\mathrm{def}}}}}\, \alpha \left( {\omega _1,\omega _2} \right)\), which evolves as:

$$\frac{{d\alpha \left( {\omega _1,\;\omega _2} \right)}}{{dt}} = \left( {{\mathrm{H}}_1\left( t \right) - {\mathrm{H}}_2\left( t \right)} \right)$$

(2)

where \({\cal{H}}_i\left( t \right) = B\;{\mathrm{sin}}\left( {\omega _i{\mathrm{t}} + \theta } \right)\sigma _{\mathrm{Z}} = {\mathrm{H}}_i\left( t \right)\sigma _{\mathrm{Z}}\) are the corresponding signal Hamiltonians, θ is the initial signal phase, and \(\sigma _Z\) is the Pauli spin-z operator. As a consequence, maximizing \(\alpha \left( {\omega _1,\;\omega _2} \right)\) implies applying a π-pulse whenever \({\mathrm{H}}_1\left( t \right) - {\mathrm{H}}_2\left( t \right)\) changes sign and it follows (see Supplementary Note 1):

$$\alpha \left( {\omega _1,\;\omega _2} \right) = {\int \nolimits_0^t} \frac{{\mu _{{\mathrm{max}}} - \mu _{{\mathrm{min}}}}}{2}dt^{\prime} \approx B\frac{2}{\pi }\omega _{\mathrm{{\Delta}}}t^2$$

(3)

where \(\mu _{{\mathrm{max}}}\;{\mathrm{and}}\;\mu _{{\mathrm{min}}}\) are the maximal and minimal eigenvalues of \({\cal{H}}_1 - {\cal{H}}_2\). Application of a series of π-pulses with a spacing \(\frac{{2\pi }}{{\omega _1 + \omega _2}}\) constitutes optimal control in this scenario. As orthogonality requires \(\alpha \left( {\omega _1,\;\omega _2} \right) = \pi /2\), the minimal discrimination time is \(T_{{\mathrm{opt}}} = \frac{\pi }{{2\sqrt B \sqrt {|\omega _{\mathrm{{\Delta}}}|} }}\).

For an illustrative explanation, one can think of \(\left| {\psi \left( {\omega _1} \right)} \right\rangle ,\;\left| {\psi \left( {\omega _2} \right)} \right\rangle\) as two runners, where the goal is to maximize the gap between them (see Fig. 1c). As long as the same runner (suppose \(\left| {\psi \left( {\omega _2} \right)} \right\rangle\)) is faster (speed proportional to the current amplitude), their separation gets larger as desired. However, once \(\left| {\psi \left( {\omega _1} \right)} \right\rangle\) starts to be faster than \(\left| {\psi \left( {\omega _2} \right)} \right\rangle\), we prevent a reduction of the gap by flipping the direction they run (equivalent to a π-pulse).

The phase acceleration in Eq. (4) coincides with the fundamental speed limit derived in ref. 21 and this implies a minimal discrimination error for every t. Hereafter, for all error analysis, we assume that the prior probability for both frequencies is 1/2, namely symmetric hypothesis testing. Given any two Hamiltonians \({\cal{H}}_1\left( t \right)\) and \({\cal{H}}_2\left( t \right)\), we prove in Supplementary Note 1 that the error probability when distinguishing between these Hamiltonians, optimized over all possible strategies, is lower bounded by:

$${\it{\epsilon }} \ge \left\{ {\begin{array}{*{20}{l}} {0.5\left( {1 - \sin \left( {\alpha _{{\mathrm{max}}}} \right)} \right)} & {\alpha _{{\mathrm{max}}} \le \frac{\pi }{2}} \\ 0 & {\alpha _{{\mathrm{max}}} > \frac{\pi }{2}} \end{array}} \right.$$

(4)

where \(\alpha _{{\mathrm{max}}} = {\int\nolimits_0^t} \frac{{\mu _{{\mathrm{max}}} - \mu _{{\mathrm{min}}}}}{2}dt^{\prime}\). This lower bound can always be saturated with a suitable control. Hence, even if \(\alpha _{{\mathrm{max}}} \le \frac{\pi }{2}\) at t, the strategy that minimizes the error probability is to apply the above discussed control and perform a measurement at t, which implies that even multiple correlated measurements at times shorter than t19,22 cannot beat this fundamental limit, as shown in Fig. 1d (see Supplementary Note 1).

Experiments

The experiments are performed using a single NV center in ultrapure diamond (Fig. 2a). Here we use a 400 G magnetic field, aligned along the NV symmetry axis, to lift the degeneracy of the three ground spin states and we use two of these states (here denoted as \(\left| 0 \right\rangle ,\;\left| 1 \right\rangle\)) as a qubit. We initialize our qubit into a coherent superposition state and map the sensor phase φ = 2α (compare Eq. (3)), acquired during interaction with the control sequence and the signals, into a population difference (such that the probability for the sensor to be projected to \(\left| 0 \right\rangle\) reads \(P = 0.5 \cdot \left( {1 - \sin \left( \varphi \right)} \right)\), which is subsequently read out optically.

a Experimental setup: the spin state of a single nitrogen-vacancy (NV) center in diamond responds to magnetic fields by a Zeeman shift of the spin levels, which can be optically initialized and readout with a confocal microscope. b Experimentally measured quantum phase difference \(\varphi _{\Delta}\) as a function of coherent interaction time and as log–log plot (inset). Solid lines are a fit to \(\varphi _{\Delta} = \frac{2}{\pi }B\omega _{\Delta}t^2\), where \(\omega _{\Delta} = |\omega _1 - \omega _2| = \left( {2\pi } \right) \cdot 2\) kHz. c Measured spin population as a function of coherent interaction time. The probability to be in state \(\left| 0 \right\rangle\) is: \({\mathrm{P}}\left( {\left| 0 \right\rangle \left| {\omega _1} \right.} \right) = \frac{1}{2} - \frac{1}{2}\sin \left( {\frac{{B\omega _{\Delta}t^2}}{\pi }} \right)\), \({\mathrm{P}}\left( {\left| 0 \right\rangle \left| {\omega _2} \right.} \right) = \frac{1}{2} + \frac{1}{2}\sin \left( {\frac{{B\omega _{\Delta}t^2}}{\pi }} \right)\). For \(\varphi _{\Delta} = \pi\), the sensor is driven to one of two orthogonal eigenstates. All error bars correspond to SD of several independent repetitions.

Full size image

In Fig. 2b, φΔ is plotted, when a magnetic field of frequency \(\omega _1 = (2\pi ) \cdot 0.999\) MHz or ω2 \(= (2\pi ) \cdot 1.001\) MHz, in the range of frequencies typical for nano-NMR experiments with NV centers and was measured with an XY8-N sequence with an interpulse spacing of 500 ns. The interaction time was extended by increasing the pulse number. A logarithmic plot of \(\varphi _{\Delta}\) shows the expected t2 increase until 44 µs when \(\varphi _{\mathrm{{\Delta}}} = \pi\), at which time the quantum sensor has been evolved into one of two orthogonal states, as \(|P\left( {\omega _1} \right) - P\left( {\omega _2} \right)| = 1\). By using a 90° phase-shifted \(\frac{\pi }{2}\) – pulse (compared to the initialization basis), the sensor phase can be mapped into a population difference and the resulting NV population as a function of interaction time is plotted (Fig. 2c). In Supplementary Note 4, we show that this balanced readout is also optimal. Experimental data for the dependency of the discrimination time on \(\omega _{\mathrm{{\Delta}}}\) and B are provided in the Supplementary Note 5. If the coherence of the sensor is not sufficient to achieve orthogonality, the described protocol still remains optimal (see detailed discussion in Supplementary Note 1).

Complex frequency discrimination

This method can be extended for complex spectra consisting of more than one frequency or amplitude component. Consider a Hamiltonian that takes the general form \({\cal{H}}_{1,2} = {\int} {f_{1,2}\left( \omega \right)\sin \left( {\omega t} \right)d\omega \;\sigma _Z = {\mathrm{H}}_{1,2}\sigma _Z}\) such that the two frequency distribution functions \(f_{1,2}\left( \omega \right)\) are known, but we do not know which of the Hamiltonians, \({\cal{H}}_1\) or \({\cal{H}}_2\), is present. As an illustrative example, we consider the 13C nuclear magnetic spectrum produced by either a sample of ethanol or propanol, both of which contain two chemically distinct carbon groups (Fig. 3a, b). A similar speed bound for discrimination applies to this case (see Supplementary Eq. 1) and is achievable with an analogous protocol. We can again define \({\mathrm{H}}_{\Delta} = {\mathrm{H}}_1 - {\mathrm{H}}_2\), where optimal control is achieved by applying π-pulses whenever \({\mathrm{H}}_{\Delta}\) switches sign (Fig. 3c). Importantly, for NMR detection of a polarized sample, the initial signal phase is known, as it is defined by a π/2-pulse applied to the sample. As a result, a sensor phase difference of π can be tailored to perform optimal discrimination. Of note, the sensor phase difference also increases according to t2, until the signals become completely out of phase with each other.

a, b Calculated magnetic field experienced by the NV center exposed to a solution of ethanol (red) or propanol (blue) in frequency (a) and time (b) domains, where a π/2 pulse is applied to the thermally polarized sample at zero time. c Calculated optimal control sequence of π-pulses applied to the NV sensor (purple) overlayed on the signal difference (green). d Sensor phase accumulation in presence of ethanol (red), propanol (blue), and the phase difference (green), showing quadratic scaling at short times.

Full size image

Implications of imperfect readout

Above, we showed that in principle, it is possible to discriminate two frequency distributions within a single measurement by eliminating projection noise. As perfect measurements are not possible in an experimental realization, additional readout noise has to be taken into account. For optical readout of NV centers, e.g., photon shot noise has to be considered. As a result, noise analysis is important to obtain a meaningful discrimination error. The readout of the NV center is performed by detecting the spin-dependent fluorescence emitted during a laser pulse (see “Methods”). The recorded photons are well described by a Poisson distribution with an average photon number of \(\lambda _0 = 0.084\) \(\left( {\lambda _1 = 0.07} \right)\) for the \(\left| 0 \right\rangle\) \(\left( {\left| 1 \right\rangle } \right)\) spin state (compare Supplementary Note 2). With this small contrast, the error of a single measurement, assuming again a symmetric hypothesis testing, is ~0.49, even for orthogonal states (see Supplementary Note 4), i.e., we have a nearly 50–50 chance to assign the frequency incorrectly.

The probability to make an incorrect decision can be reduced in two ways: first, by increasing the number of measurements, which we refer to as ensemble averaging. Then, the difference in the number of emitted photons increases, which leads to better discriminability. The second approach is based on improving the readout process itself by increasing the contrast between the states. This can be achieved by introducing an ancilla qubit, which acts as a quantum memory and stores the state of the NV center23,24. This method is usually referred to as single-shot readout (SSR). We experimentally implement these two methods and compare their performance using the error probability as figure of merit. We show that by benchmarking against the number of detected photons, ensemble averaging always performs at least as good as SSR, but when compared in terms of measurement time or number of coherent interaction periods, SSR has particular advantages.

Temporal ensemble averaging was performed using a 1.35 NA oil objective, to collect the NV fluorescence and repeating the sensing and readout Nens times (compare Fig. 4a). Hence, given the orthogonal states and Nens repetitions of the measurement, the discrimination between the states reduces to discrimination between two Poissonian distributions: \({\mathrm{Poi}}\left( {N_{{\mathrm{ens}}}\lambda _0} \right)\) and \({\mathrm{Poi}}\left( {N_{{\mathrm{ens}}}\lambda _1} \right)\). Assuming small contrast, namely \(\frac{{\lambda _0}}{{\lambda _1}}\) is close to 1, we observe that the error probability scales as \(\exp ( { - N_{{\mathrm{ens}}}0.5\left( {\sqrt {\lambda _1} - \sqrt {\lambda _0} } \right)^2})\) (Supplementary Note 4). For non-orthogonal states, the error probability also decays exponentially with Nens to 0 and the error exponent is \(\sim \frac{{\left( {\lambda _0 - \lambda _1} \right)^2\sin \left( \alpha \right)^2}}{{4\left( {\lambda _0 + \lambda _1} \right)}}\) (Supplementary Note 4), where α is the angle between the states. Hence, given that the duration of a single ensemble average is τ (including the overheads for readout and initialization), the discrimination time with this strategy goes as \(\sim \tau \frac{{4\left( {\lambda _0 + \lambda _1} \right)}}{{\left( {\lambda _0 - \lambda _1} \right)^2\sin \left( \alpha \right)^2}}\)

a Ensemble averaging. The complete sequence of sensor initialization, sensing period, and optical readout is performed Nens times. b Single-shot readout (SSR). After one sensing duration, the sensor state is mapped onto the ancilla readout qubit (13C nucleus). Optical readout of the ancilla spin is weakly perturbing, allowing its population to be probed NRR times with a quantum non-demolition measurement.

Full size image

SSR is implemented by using a weakly coupled 13C spin that forms a memory qubit and allows repetitive readout of the NV state with a quantum non-demolition (QND) experiment25 (compare Fig. 4b and Supplementary Note 2). The ancilla qubit allows more photons to be scattered before its state is destroyed and here we perform NRR repetitive readouts after one sensing time. The larger NRR, the closer we get to a perfect quantum measurement, limited ultimately by the lifetime of the ancilla. For NRR = 104, the photon statistics manifest digital step jumps, which allows high-fidelity readout of the nuclear spin state (see Supplementary Note 2). The discrimination between orthogonal states using this method is, in theory, a discrimination between \({\mathrm{Poi}}\left( {N_{{\mathrm{RR}}}\lambda _0} \right)\) and \({\mathrm{Poi}}\left( {N_{{\mathrm{RR}}}\lambda _1} \right)\), and thus given as \(N_{{\mathrm{RR}}} = N_{{\mathrm{ens}}}\), precisely the same as for ensemble averaging. Then, the only difference is the time required to perform readout of the nuclear spin and that SSR requires only one sensing period, which is discussed in more detail below. For non-orthogonal states however, this strategy is inferior as the error probability does not reduce to 0 in the limit of large NRR due to quantum projection noise. This additional uncertainty inherent to quantum measurements results in a scaling as \(0.5\left( {1 - |\sin \left( \alpha \right)|} \right)\) (see Supplementary Note 1).

We experimentally investigated the performance of both strategies by performing ensemble averaging and SSR for discrimination of the same signals described earlier. The error as a function of different parameters is plotted in Fig. 5a–e. In Fig. 5a, it is shown that for any angle α between the states (or equivalently, for any coherent interaction period), ensemble averaging performs better than SSR. However, in contrast to the theoretical expectation, we see that ensemble averaging yields a lower error than SSR also for orthogonal states (44 μs, \(\varphi _{\Delta} = \pi\)). The reason for this is that, although the nuclear state can be readout with high fidelity in our experiments, imperfect initialization of the NV center into the correct charge state and the finite ancilla T1 lifetime, limits the sensor readout fidelity to 0.8 (see Supplementary Note 2 and refs. 26,27). In Fig. 5a (red line), we plot the expected error achievable with SSR under perfect sensor initialization and control, which demonstrates the equivalence between the two readout methods when orthogonality is achieved.

a Experimental frequency discrimination error as a function of coherent interaction time. Data points are for sensor readout using a single SSR measurement (orange, NRR = 104) and ensemble averaging (lilac, Nens = 104). Fits are from Supplementary Eq. 89 and 91. The ideal SSR curve corresponds to perfect sensor initialization and mapping to ancilla. b Experimental frequency discrimination error as a function of detected photons for a single SSR measurement (orange) and ensemble averaging (lilac), for an interaction time of 44 µs \(\left( {\varphi _{\Delta} = \pi } \right)\). Imperfect sensor initialization limits the error achievable with a single SSR measurement. c Frequency discrimination error as a function of total measurement time for a single SSR measurement (orange, NRR = 104) and ensemble averaging (lilac, Nens = 104). Each data point corresponds to a different interaction time. d Frequency discrimination error as a function of total measurement time, with a fixed interaction time of 44 µs \(\left( {\varphi _{\Delta} = \pi } \right)\). Although for SSR (orange) the total time is incremented by performing additional measurements (each using NRR = 104), for ensemble averaging Nens is increased. Expected errors obtained from simulations for \(N_{{\mathrm{RR}}} = 10^2,\;10^3,\;10^4\) are included (solid blue, green, and yellow lines, see Supplementary Note 4). e Frequency discrimination error as a function of number of trials, i.e., number of times the signal is interrogated, with a fixed interaction time of 44 µs using many SSR (each using NRR = 104, orange) or ensemble averaging (lilac). Error bars are 1 SD.

Full size image

Figure 5b shows the discrimination error as a function of the number of photons detected with each readout technique, for a fixed interaction time of 44 μs (Δφ = π, dotted line in Fig. 5a). As the sensor is driven to an eigenstate, each photon should convey the same information about which frequency is present regardless of whether it was obtained from ensemble averaging or SSR. We observe that imperfect initialization reduces the information carried by each photon using SSR, resulting in the error converging to an offset of 0.2. By using charge-state detection to improve the sensor initialization, it is expected that a convergence in the two readout strategies would be observed.

However, SSR may still be useful if the total measurement time is taken into account. For ensemble averaging, every repetition cycle requires a duration of initialization \(\left( {t_{{\mathrm{init}}}} \right)\), interaction (t), and readout \(\left( {t_{{\mathrm{read}}}} \right)\); thus, the total measurement time is \(T_{{\mathrm{ens}}} = N_{{\mathrm{ens}}}\left( {t_{{\mathrm{init}}} + t + t_{{\mathrm{read}}}} \right)\). For our experiments, \(t_{{\mathrm{init}}} + t_{{\mathrm{read}}} \approx 1.5\;\mu {\mathrm{s}}\). For SSR, only one initialization and interaction period are required, while readout is performed NRR times. The initialization and readout time for SSR are different to ensemble averaging, however, as manipulation of the ancilla is required; thus, the total time is \(T_{{\mathrm{SSR}}} = t_{{\mathrm{init}}}^{{\mathrm{anc}}} + t + N_{{\mathrm{RR}}}\;t_{{\mathrm{read}}}^{{\mathrm{anc}}}\). For our experiments, \(t_{{\mathrm{init}}}^{{\mathrm{anc}}} \approx 100\;\mu {\mathrm{s}}\) and \(t_{{\mathrm{read}}}^{{\mathrm{anc}}} \approx 17\;\mu {\mathrm{s}}\), so for NRR > 10, most of the temporal overheads for SSR arise from readout of the ancilla spin via repeated mapping onto the NV electron spin (see “Methods”). For long interaction times (⪎17 µs), this overhead is less costly than repeating the interaction and SSR is faster than ensemble averaging, as shown in Fig. 5c.

In addition to performing only a single SSR measurement, we investigate the discrimination error when performing multiple SSR measurements. In Fig. 5d, the discrimination error is plotted as a function of the total time required by each readout strategy, including all measurement overheads. For SSR, the time is incremented by performing additional measurements, each using NRR = 104 repetitive readouts, which we in the following refer to as hybrid strategy, as many individual SSR are averaged, while for ensemble averaging just Nens is increased. Although about twice the number of photons is recorded for the hybrid strategy within the same time, ensemble average still performs better. The reason for this is that NRR is too large; we observe that by reducing NRR, the hybrid strategy can yield a smaller error as a function of time. For NRR = 103, we find a slightly better scaling for the hybrid strategy using simulations, taking our measurement parameters into account (see Supplementary Note 4), as projection noise is sufficiently reduced. By further reducing to NRR = 102, we find that the error exponent with the hybrid strategy is almost twice as large as the error exponent of ensemble averaging and there is a reduction of the overall error by a factor of 40 after measuring for 1 s. A more detailed analysis and optimization over NRR can be found in Supplementary Note 4. In addition to benchmarking against the total measurement time, we also compare in terms of the number of trials, (i.e. the number of times the signal is queried), which becomes critical if many interrogations are prohibited, e.g. due to sample contamination or degradation. As expected, the hybrid strategy performs significantly better, as a high confidence can be obtained from a single trial (Fig. 5e). In Fig. 5d, e the error decays exponentially with the total time (or equivalently number of queries) as expected, the error exponent for both methods is calculated explicitly in Supplementary Note 4.

It should be noted that these results apply to any problem with the aim to discriminate between two NV center states.

Qptimal frequency estimation using a quantum coherent probe

A related task is to estimate a single unknown frequency ω of a signal \({\cal{H}}\left( \omega \right) = B\;\sin \left( {\omega t + \theta } \right)\sigma _Z = H(\omega )\sigma _Z\) with minimum uncertainty \({\mathrm{{\Delta}}}\omega\). The variance \(\left( {{\mathrm{{\Delta}}}\omega } \right)^2\) is lower bounded by the inverse of the quantum Fisher information (QFI) \(I(\omega )\)28:

$$\left( {{\mathrm{{\Delta}}}\omega } \right)^2\, \ge \,\frac{1}{{I(\omega )}}$$

(5)

Notably, the QFI has a clear geometric meaning, as I(ω) can be written using the Bures distance dB between quantum states:

$$I(\omega ) = \mathop {{\lim }}\limits_{\omega _{\mathrm{{\Delta}}} \to 0} \frac{1}{{\omega _{\mathrm{{\Delta}}}^2}}d_{\mathrm{B}}\left( {|\psi \left( {\omega - \omega _{\mathrm{{\Delta}}}/2} \right)\rangle ,|\psi \left( {\omega + \omega _{\mathrm{{\Delta}}}/2} \right)\rangle } \right) \\ = \mathop {{{\mathrm{lim}}}}\limits_{\omega _{\mathrm{{\Delta}}} \to 0} \frac{{\alpha \left( {\omega - \frac{{\omega _{\mathrm{{\Delta}}}}}{2},\omega + \frac{{\omega _{\mathrm{{\Delta}}}}}{2}} \right)^2}}{{\omega _{\mathrm{{\Delta}}}^2}},$$

(6)

where \(\alpha \left( {\omega - \frac{{\omega _{\Delta}}}{2},\omega + \frac{{\omega _{\Delta}}}{2}} \right)\) is the angle between \(|\psi \left( {\omega \pm \omega _{\mathrm{{\Delta}}}/2} \right)\rangle\). Hence, maximizing the QFI is equivalent to optimizing the discrimination of \(\mathop {{{\mathrm{lim}}}}\nolimits_{\omega _{\mathrm{{\Delta}}} \to 0} \left| {\psi \left( {\omega \pm \frac{{\omega _{\mathrm{{\Delta}}}}}{2}} \right)} \right\rangle\). As shown above, this means maximizing the accumulation of \(\frac{{dH}}{{d\omega }} = \mathop {{{\mathrm{lim}}}}\nolimits_{\omega _{\mathrm{{\Delta}}} \to 0} \frac{{H\left( {\omega + \frac{{\omega _{\mathrm{{\Delta}}}}}{2}} \right) - H\left( {\omega - \frac{{\omega _{\mathrm{{\Delta}}}}}{2}} \right)}}{{\omega _{\mathrm{{\Delta}}}}}\), which is achieved by applying π-pulses whenever \(\frac{{dH}}{{d\omega }}\) changes sign. This optimal strategy is illustrated in Fig. 6a, where π-pulses are applied at the antinodes of the signal with a spacing close to \(\frac{\pi }{\omega }\). According to Eq. (3), it follows:

$$\alpha \left( {\omega - \frac{{\omega _{\mathrm{{\Delta}}}}}{2},\;\omega + \frac{{\omega _{\mathrm{{\Delta}}}}}{2}} \right) = B\frac{2}{\pi }\omega _{\mathrm{{\Delta}}}t^2 \Rightarrow I(\omega ) = \left( {\frac{2}{\pi }} \right)^2B^2t^4$$

(7)

The minimal uncertainty obtainable in a single experiment reads16:

$${\Delta}\omega = \frac{\pi }{{2Bt^2}}$$

(8)

and scales again as 1/t2 due to the phase acceleration of the sensor. This is a special case of the analysis performed in refs. 16,17,18,19 and partially realized in ref. 15 without taking readout noise into account.

a Maximal sensitivity to amplitude changes is obtained for π-pulses at the signal nodes \(\left( {\theta = \frac{\pi }{2}} \right)\), resulting in a linear accumulation of φ. Maximal sensitivity to frequency changes is obtained for π-pulses at the signal antinodes (θ = 0). b Experimental data of the phase difference for a 1 µT change in amplitude (red) and 10 kHz change in frequency (black). Frequency estimation results in a quadratic increase of φ, whereas amplitude estimation produces a linear increase. c Measured NV population as a function of π-pulse spacing \(\tau _\pi\), given \(\theta = \frac{\pi }{2}\) (blue) and θ = 0 (green). Fits correspond to Supplementary Eq. 39. d Measured derivative of the NV population as a function of π-pulse spacing \(\tau _\pi\), given \(\theta = \frac{\pi }{2}\) (blue) and θ = 0 (green). Fits correspond to Supplementary Eq. 41. e Measured frequency uncertainty as a function of interaction time, for X (dark purple) or Y (light purple) measurement basis. Red line corresponds to Eq. (8) and pink line to Eq. (9). f Measured frequency uncertainty as a function of interaction time, for X (dark orange) or Y (light orange) measurement basis using a single SSR measurement. Red line corresponds to Eq. (8). g Frequency estimation for short coherence signals (for ensemble averaging): measured frequency uncertainty given \(\theta = \frac{{\pi }}{2}\) (blue), θ = 0 (green), and unknown initial phase (gray) for a single measurement with interaction time shorter than the signal coherence time (200 µs). For longer times, the sensitivity is obtained from averaging multiple measurements, each with an interaction time of 97 µs.

Full size image

A further peculiarity is that the π-pulses should be applied close to the signal antinodes where \(\frac{{dH}}{{d\omega }}\) changes sign, which is in contrast to amplitude estimation. We demonstrate the underlying intuition by recording spectra for the two edge cases—control started close to a signal node \(\left( {\theta = \frac{\pi }{2}} \right)\) or antinode \(\left( {\theta = 0} \right)\). We again used XY8 sequences and chose the readout basis such that the final NV population reads \(P = 0.5 \ast \left( {1 + \sin \left( \varphi \right)} \right)\). When the π-pulses are placed at the signal nodes (Fig. 6b blue), the NV population P is maximal, but \(\frac{{\partial P}}{{\partial \omega }}\) vanishes; thus, although the sensor acquires maximal phase, it is insensitive to changes in the signal frequency. In contrast, for θ = 0 (Fig. 6b green), the sensor acquires minimal phase; however, small frequency changes lead to large population deviations (here, \(\left( {\frac{{\partial P}}{{\partial \omega }}} \right)_{{\mathrm{max}}}\, \approx 10^7\;{\mathrm{Hz}}^{ - 1}\), which corresponds to \(\left( {\frac{2}{{4\pi }}} \right)Bt^2\)).

Assuming perfect quantum measurements, any measurement basis in the X–Y plane would be optimal, as they all saturate the QFI (see Supplementary Note 7). The imperfect readout of the NV center not only worsens the estimate, but also provokes a dependency on the measurement basis29. When using conventional readout of single NV centers, the fluorescence contrast is very low and the measurement noise barely varies with the measurement basis (see Supplementary Note 8). As a consequence, it is optimal to measure (approximately) in the basis that yields the maximal slope. In our case, assuming an initialization in the X basis, a maximal slope is obtained by measuring in the Y-basis and it follows (in the limit of small contrast, see Supplementary Note 8):

$${\mathrm{{\Delta}}}\omega \approx \frac{\pi }{{\sqrt 8 Bt^2}}\frac{{\sqrt {\lambda _0 + \lambda _1} }}{{\left| {\lambda _0 - \lambda _1} \right|}}$$

(9)

Taking our experimental values of λ0, λ1 into account, a factor of ~40 in SD is lost.

We realize experimental frequency estimation by measuring an unknown frequency ω (=999.3 MHz) oscillating close to 1 MHz with known amplitude (1.7 µT). Applying XY8 sequences with a pulse spacing of 500 ns and locking θ = 0, we use ensemble averaging \(\left( {N_{{\mathrm{ens}}} \approx 1000} \right)\) to estimate the frequency several thousand times for a variety of interaction times t (controlled by tailoring the π-pulse number). By computing the SD and assuming shot noise \(\sqrt T\)-scaling, we are able to extract an estimate for the uncertainty \({\mathrm{{\Delta}}}\omega\) for a single measurement, which is plotted in Fig. 6c. As expected, there is a significant difference between X and Y readout (for initialization in X). Measurements in the Y basis saturate the bound in Eq. (9) and, in particular, \({\mathrm{{\Delta}}}\omega\) scales as t−2. In contrast, measuring in the X basis leads to much higher uncertainties, due to the very small slope \(\frac{{\partial P}}{{\partial \omega }}\). Interestingly, using SSR (here NRR = 104), the results almost coincide with the theoretical limits of perfect projective measurements (Fig. 6c lower), namely a close agreement in uncertainty using X/Y measurements and an SD of 1.6 µHz/Hz2 within a factor of 2 to the limit, 0.9 µHz/Hz2 set by the QFI (Eq. (8), red line in Fig. 6c).

Estimation: consequences of short coherence time

Finally, we address the relevance of this method to practical frequency estimation scenarios in NMR and communication. Although the described protocol is optimal if both, the signal and the sensor, are perfectly coherent, any deviations require further analysis of the estimation strategies.

For signals with short coherence time, which are especially relevant in nanoscale NMR where influences of molecular diffusion and nuclear couplings limit the coherence of the NMR signal6,30, the described protocol remains optimal, as it provides the most information per individual measurement (and correlations are not possible). As shown in Fig. 6d, the frequency can be estimated with an SD improving as t−2 during one sensing period, where a timing of θ = 0 (green) provides three orders of magnitude improvement in SD compared to θ = π/2 (blue). The estimate can be further improved by performing multiple measurements up to a total measurement time T, albeit with a reduced T−1/2 scaling, typical for incoherent averages. We compute sensitivities, which are here defined as \(\eta = {\mathrm{{\Delta}}}\omega \sqrt T\) of 58 Hz/Hz0.5 for θ = 0 and 23 kHz/Hz0.5 for θ = π/2. If the signal phase is not known a priori, we cannot apply pulses at the right timing. Nevertheless, using adaptive measurements, the t−2 uncertainty scaling can be preserved, albeit with a higher uncertainty \({\mathrm{{\Delta}}}\omega\). In Fig. 6d (gray), the frequency SD is plotted where the spacing \(\tau _\pi\) was updated for each interaction time and the measurement was averaged over all signal phases (see Supplementary Note 10). Using this strategy, the frequency can be estimated with an order-of-magnitude reduction in uncertainty and sensitivity (here η = 630 Hz/Hz0.5), but without the need to know the starting phase.

In the other limit, the coherence time of the sensor is much shorter than the coherence time of the signal. It has recently been shown that applying a chain of consecutive, correlated, measurements can enhance frequency estimation in this case19,22,31,32. Due to the stability of the phase of the signal, the amount of information gained from late measurements is much larger than the information gained from early ones. Hence, by employing the signal coherence, the T4 scaling of the QFI, which relies on the coherence of the sensor, degrades to T3. Without any control, the QFI in this regime reads (compare ref. 19):

$$I\left( \omega \right) = \mathop {\sum }\limits_\tau 4B^2t^2\tau ^2\;{\mathrm{sinc}}\left( {\omega t/2} \right)^2\cos \left( {\omega \tau + \frac{{\omega t}}{2} + \varphi } \right)^2 \\ \approx \frac{2}{3}B^2T^3t\;{\mathrm{sinc}}\left( {\omega t/2} \right)^2 \le 0.96\frac{{B^2T^3}}{\omega }$$

(10)

where the decay prefactor: \({\mathrm{sinc}}\left( {\omega t/2} \right)^2\cos \left( {\omega \tau + \frac{{\omega t}}{2} + \varphi } \right)^2\) stems from the oscillations of the signal slope. This term can be suppressed by applying π-pulses close to the antinodes and the ultimate limit of the QFI in this regime reads (compare Supplementary Note 11):

$$I(\omega ) = \left( {\frac{2}{\pi }} \right)^2B^2\mathop {\sum }\limits_\tau \left( {\left( {t + \tau } \right)^2 \,- \, \tau ^2} \right)^2 = \frac{4}{3}\left( {\frac{2}{\pi }} \right)^2B^2T^3t.$$

(11)

Hence, the scaling is the same as without control \(\left( {B^2T^3t} \right)\); however, there is a difference of \(\sim 0.56\omega t\), which can be significant for high frequencies (taking relevant experimental values, \(\omega t \sim 20\), it follows that a factor of ~10 is lost in the absence of control). Here again, the timing of the pulses is crucial. Applying resonant pulses generally leads to: \(I = \frac{4}{3}\left( {\frac{2}{\pi }} \right)^2B^2T^3t\cos \left( \theta \right)^2\) (see Supplementary Note 11). As a consequence, \(I(\omega )\) vanishes for θ = π/2 and is maximal for θ = 0. More generally, applying pulses with a general timing (θ) and a general detuning δ (defined as \(\omega - \pi /\tau _\pi\), where \(\tau _\pi\) is the spacing between the π-pulses), leads to:

$$\begin{array}{l}I(\omega ) \approx \mathop {\sum }\limits_\tau 4\left( \displaystyle {\frac{2}{\pi }} \right)^2B^2t^2\tau ^2\,{\mathrm{sinc}}\left( {\delta t/2} \right)^2\cos \left( {\delta \left( {\tau + \displaystyle\frac{t}{2}} \right) + \theta } \right)^2\\ \approx \displaystyle \frac{2}{3}\left( {\frac{2}{\pi }} \right)^2B^2T^3t\left\{ {\begin{array}{*{20}{c}} {2\cos \left( \theta \right)^2\delta T \ll 1} \\ {{\mathrm{sinc}}\left( {\delta t/2} \right)^2\delta T \gg 1} \end{array}} \right.\end{array}$$

(12)

It is noteworthy that in the limit of \(\delta t \ll 1,\;\delta T \gg 1\) (which was implemented in ref. 19), only a factor of 1/2 is lost, compared to the optimum, due to the detuning. Hence, in this regime, the maximal achievable advantage, compared to the results in ref. 19, is a factor of 2.