Noise-robust quantum sensing via optimal multi-probe spectroscopy

The dynamics of quantum systems are unavoidably influenced by their environment, but in turn observing a quantum system (probe) can allow one to measure its environment: Measurements and controlled manipulation of the probe such as dynamical decoupling sequences as an extension of the Ramsey interference measurement allow to spectrally resolve a noise field coupled to the probe. Here, we introduce fast and robust estimation strategies for the characterization of the spectral properties of classical and quantum dephasing environments. These strategies are based on filter function orthogonalization, optimal control filters maximizing the relevant Fisher Information and multi-qubit entanglement. We investigate and quantify the robustness of the schemes under different types of noise such as finite-precision measurements, dephasing of the probe, spectral leakage and slow temporal fluctuations of the spectrum.


Results
Filter function approach and protocols for spectrally resolved sensing. Let us consider a two-level system, given by the two levels |0〉 and |1〉, that interacts with a fluctuating field E(t) via E(t)σ z , where σ z is the z-Pauli matrix -see also Fig. 1. Following 6,20,34 , we prepare the system in its lower state |0〉 and apply an initial π/2-pulse to bring it into the superposition state | 〉 + | 〉 ( 0 1 )/ 2. Then, we apply a sequence of π-pulses (i.e. population flips between |0〉 and |1〉) that we describe by a pulse modulation function y(t) ∈ {−1, 1} that switches sign at the position of the pulses.
The system then acquires a phase φ(t) leading to the state ( ) , and this phase can be measured via a population measurement after a final π/2-pulse that brings the system into the state , where T is the total operation time of the pulse sequence. The phase is determined by the fluctuating field and the modulation function, and can be written as If we average over many realizations of the sequence, then we can describe the dephasing by a decoherence function χ(t). The latter enters the probability p(t) to find the system in state |0〉 upon measurement as t ( ) and is given by is the autocorrelation function of the fluctuating field E(t). In Eq. (4), the brackets 〈·〉 mean averaging over the realizations of the stochastic field E(t).
The Fourier transform of the autocorrelation function denotes the power spectral density S(ω) of the field E(t): ωτ −∞ ∞ S g e d ( ) ( ) (5) i As shown in detail in the Methods, also for continuous pulse modulations, the decoherence function χ(T) can be written in terms of S(ω) and a so-called filter function F(ω), i.e.
T 0 In the following, we call T the filter operation time. The filter function is defined as T T i t 0 is the Fourier transform of the pulse modulation function. To simplify the notation, we will omit the time in the subscript and identify F ≡ F T and Y ≡ Y T . By engineering the pulse modulation function, we can design different filter functions F(ω) that select specific frequency ranges of the power spectral density S(ω). In order to estimate the functional behaviour of the power spectral density S(ω), we employ in the following a set of K filter functions F k (ω), k = 1, …, K, each of them generated by a different sequence of π-pulses and we measure the decoherence function after the application of each of these filters. We will then describe protocols that compose the functional behaviour of S(ω) from the single measurement outcomes, so that we can reconstruct S(ω) in a range ω ∈ [0, ω c ] for some cut-off frequency ω c . Let us first introduce a novel protocol, which is based on the orthogonalization of the filter functions F k (ω), k = 1, …, K. For the sake of brevity, from here on we will denote this new protocol with the acronym Filter Orthogonalization (FO) protocol. Then, we recall the Alvàrez-Suter (AS) protocol 34 , that will be used as a benchmark to test the performance of the new approach. Finally, we will also introduce a fidelity measure to evaluate the performance of the protocols.
Filter Orthogonalization protocol. The filter functions of the FO protocol in principle can be chosen arbitrarily.
Here we choose to base them on equidistant pulse sequences and design the latter by positioning the π-pulses at the zeros of Note that this includes also a sequence without any π-pulse (for k = 1) corresponding to a filter function sensitive to the frequency ω = 0. Instead, ω max is a constant that determines the bandwidth within which we can analyse the spectrum of the signal (the field E(t)), and the filter corresponding to the highest frequency is centered around (1 − 1/K)ω max . The population measurements at the end of the filter application will yield the coefficients k k 0 which is the value of the decoherence function χ(T) at the end of the application of the filter F k . We then calculate the K × K matrix A, whose matrix elements kl k l 0 c quantify the overlap in the frequency domain between the filter functions F k (ω) and F l (ω) (k, l ∈ 1, …, K). We truncate the integral at the cut-off frequency ω c , since we want to analyse S(ω) only in the interval [0, ω c ]. The matrix A is symmetric and can be orthogonalized by means of the following transformation: T where Λ ≡ diag(λ 1 , …, λ K ) are the eigenvalues of A and V is an orthogonal matrix. The K filters span a K-dimensional function space that has an orthonormal basis k l kl 0 c with the Kronecker-Delta δ kl . If we expand the spectral density S(ω) in this orthogonal basis, we obtain the coefficients Hence, we obtain an estimate ω  S( ) of the power spectral density S(ω) given by the expansion Alvàrez-Suter protocol. For the AS protocol we use as well K equidistant pulse sequences. However, due to the requirements of this protocol, this time we position the π-pulses at the zeros of As a consequence, the minimal frequency that can be resolved is ω max /K, while the maximal frequency that can be resolved is exactly ω max . Indeed, the protocol will yield the value of S(ω) at the K discrete points ω ω = k m ax k K . To reconstruct S(ω) pointwise, we then employ the AS protocol 34 , where long pulse sequences (ca. 30-100 pulses) produce a Dirac-delta-like shape of the filter functions. The FO protocol, instead, does not rely on this Dirac-delta-shape filters and thus can in principle work with a smaller number of pulses (smaller filter operation time) and reconstructs the power spectral density S(ω) as a continuous function. In the following, unless specified otherwise, for both protocols we choose K = 20 and ω c = 10. For the FO protocol, we choose ω max = 11.5 to avoid border effects near to the cut-off frequency. For the AS protocol, instead, we choose ω max = ω c = 10 so that the K values of S(ω) obtained by this protocol cover the desired bandwidth.
Estimation Fidelity. To compare the two estimation protocols and different sets of parameters, we introduce a fidelity measure. This is not straightforward, since the AS protocol reconstructs S(ω) only in a set of discrete points ω ω = k m ax k K , k = 1, … K, while the FO protocol reconstructs a continuous function. We choose to compare the values of S(ω) only on the mentioned set of discrete points, and define the fidelity as where ω  S( ) is the reconstructed power spectral density and the norm ω S( ) is given by for the power spectral density and analogously for its estimate ω  S( ).
Robust noise sensing under measurement noise and dephasing. The estimation of the power spectral density S(ω) relies on the measurement of the probability = − − p e (1 ) k c 1 2 k , which quantifies the overlap between S(ω) and each filter function F k (ω), k = 1, …, N, given by the coefficients . The values of the p k 's are measured via state population measurements, which are unavoidably affected by external sources of error, due to the presence of imperfections of the measurement device, statistical errors and detector noise contributions.
In this section, we will model two essential sources of error in the state population measurement. In particular, we assume that each measurement of p k is affected by an absolute error Δp k , which is independent from the value of p k (and thus its statistics is independent from k). This error Δp k is chosen randomly from a uniform distribution defined by the interval Δp k ∈ [−Δp max , Δp max ], where Δp max ≥ 0 is the maximum absolute value admissible by Δp k . Furthermore, we assume also that the relative phase between the two probe qubit states (|0〉 and |1〉, respectively) is additionally damped by a dephasing contribution (different from the influence of the field E(t)). This additional dephasing models the instability/natural imperfection of the probe system and is described by the rate Γ, such that for a filter operation time T the decoherence function increases by ΓT with respect to the case without dephasing. We have thus to take into account that under the presence of measurement noise and dephasing we will not effectively measure p k , but the probability p k defined as As a consequence, for a small error Δp k it holds that which has been obtained by performing the derivative of Eq. (17) with respect to c k . In Eq. (19), Δc k is the resulting error in the coefficient c k . Since we want to determine the functional shape of S(ω), the relevant error is the relative error Δc k /c k (the absolute value of the coefficients c k can usually be tuned by changing the distance of the probe to the source of the signal S(ω), which results in a change in the intensity of the field E(t) at the site of the probe). For this relative error we find k k c k T k k which takes its minimum for c k = 1, with a quite flat behaviour (about 10% increase) between about 0.5 and 1.5. Accordingly, if possible, the measurement has to be performed in this regime of maximum sensitivity and we will do this in the simulations. Let us observe, moreover, that the relative error scales exponentially with respect to the filter operation time T and the dephasing rate Γ.
In the remaining part of this paragraph, we have accordingly studied the performance of the FO protocol in comparison with the AS protocol by simulating the measurement of the power spectral density which is a spectrum that can arise for example from the spontaneous decay at two different frequencies. In the next figures, the operation time is a dimensionless number given in the unit of ωT, while for the error model we have assumed a measurement noise given by Δp max = 0.01. Then, we have chosen S 0 ∝ T in such a way that the probe operates in the maximum sensitivity range (c k ≈ 1) for all choices of T and most of the K filters (we consider the same S 0 and T for all K filters, but not all of them have the same overlap with S(ω)). Let us observe that for a given dephasing rate Γ the time-exponential increase of the measurement error Δc k /c k will favour smaller values of the filter operation time T. On the other side, if we choose T too small, then the filters lose their frequency selectivity as they flatten out in the frequency domain. Figure 2 shows how the fidelity scales with the filter operation time T for the FO protocol given a dephasing rate of Γ = 0.4. For each value of Γ and for both sensing protocols, we have prepared such a graph to determine the optimal value of T out of the tested values, e.g. T = 2 out of the values T = 1, 2, 3, 5, 7, 10 used in the example of Fig. 2. Then, in Fig. 3 we reconstruct the spectral density S(ω) for this optimal choice of T and show the fidelity of this estimate as a function of   Fig. 2. The optimal filter operation time depends on Γ and on the protocol, and was found to be in the range [5,25] for the AS protocol and [2,5] for the FO protocol. the dephasing Γ. Finally, the reconstructed functional shape of the power spectral density S(ω) obtained by each of the sensing protocols is shown in Fig. 4 for Γ = 0 and in Fig. 5 for Γ = 0.4, respectively. Note that if the measurements are affected by noise, for the FO protocol the estimation accuracy can be improved just by omitting in the expansion of Eq. (14) the terms with the smallest eigenvalues λ k of the matrix A. In the simulations, depending on T we have used between 10-20 terms out of K = 20 (the larger T, the less eigenvalues are omitted) in order to minimize the effect of the noise.
Application: NV-centers in diamond. The above examples are calculated with dimensionless variables and the results hold with suitable scaling of the parameters. As an applicative example, we can consider a nitrogen-vacancy defect in diamond. Typical dephasing rates are of the order of 1/Γ = 100 μs 42 . By fixing the dephasing rate, the units of all other parameters are defined so that we can now interpret Fig. 5 as an example with NV-centers in diamond with parameters 1/Γ = 100 μs, Δp max = 0.01, and ω ∈ [0, 2π × 0.4 kHz]. By taking such values for the noise parameters, the corresponding operation time is, respectively, equal to T = 80 μs for the FO protocol and T = 200 μs for the AS protocol.
Fighting spectral leakage via entangled multi-qubit probes. In this section, we introduce estimation protocols employing entangled multi-qubit probes and show how they can be used to fight spectral leakage in the sensing protocol. Spectral leakage is a known source of error that arises from the fact that the support of the filter functions is larger than the bandwidth they can analyse 29,34 . As a consequence, spectral contributions of the signal S(ω) outside the analysed bandwidth, but within the support of the filters, can disturb the reconstruction of the signal also within the analysed bandwidth.
The key idea of our entangled multi-qubit probes is that instead of applying the initial π/2 pulses, we prepare an N qubit system in the GHZ state Then, during the main part of the protocol, we apply single qubit pulses to each qubit individually. The pulse modulation function y(t) is then generalized to the single qubit modulation function y i (t) (i = 1, …, N), where a switch of sign of y j (t) indicates the application of a single qubit π-pulse on qubit j. The final filter function for N qubits then reads We can then use these multi-qubit probe filters for the FO protocol. For more details we refer to the Methods. We test the multi-qubit filters by reconstructing in a range ω ∈ [0, ω c = 10] with filters designed for that range. The first two terms are as in the previous section.
The last term has its main contribution outside the analysed range and thus potentially leads to spectral leakage. We choose its amplitude to be large enough to effectively disturb the measurement of the spectrum in the range ω ∈ [0, ω c = 10] using the methods presented in the previous section. Figure 6 shows how the estimation fidelity (FO protocol) scales with the number of qubits in the multi-qubit probe. Note that this is done only for the FO protocol since the AS protocol does not allow for a straightforward extension to multiple qubits as it relies on Dirac-Delta-shaped filters. Figure 7 gives an insight to the dependency of the estimation fidelity on the number of qubits: The peak of the power spectral density outside the filter range causes a distortion of the reconstruction when using only one qubit (spectral leakage). The reason is that the pulse modulation function introduces also higher frequencies (i.e. outside the analysed bandwidth or filter range) due to the stepwise modulation. A multi-qubit probe allows for intermediate steps and thus suppresses these high-frequency contributions in the filters. As a consequence, the peak outside the filter range does not significantly disturb the reconstruction within the analysed range. This is true already for only 3-5 qubits.
Application: Trapped Ions. The essential ingredient in this section is the creation of the GHZ-states. These states can be conveniently generated in linear ion traps via a Mølmer-Sørensen-gate operation 53 . Each ion represents one qubit and the gate operation generates a GHZ-state on the spin states of the ions. While the gate is very robust with respect to the state of the vibrational mode, still the states are not perfectly generated, especially for more than two qubits. In our model, we encode this state-preparation error in the value of Δp max . Inspired by the experimental values in 53 we choose Δp max = 0.01, 0.02, 0.03, 0.04, 0.1 for N = 1, 2, 3, 4, 6 qubits (note that we omit N = 5, since in the reference no value for the fidelity with N = 5 is presented), respectively and a dephasing rate Γ = 0.01 ms −1 . We simulate the measurement of the signal in a range ω ∈ [0, ω c ] with filters designed for that range. The frequency parameters are ω 1 = 2π × 1 kHz, ω 2 = 2π × 3 kHz, ω 3 = 2π × 10 kHz, and ω c = 2π × 5 kHz. Figure 8 shows the resulting fidelity of reconstruction. The optimal filter operation time turned out to be T = 10 ms for N = 1 qubit and T = 4 ms for N ≥ 2 qubits. Smaller filter operation times do not produce the frequency range accurately enough. Larger filter operation times are affected by dephasing, especially collective dephasing when more than one qubit is used. The highest fidelity is . This value could be improved by correcting for deterministic errors in the state preparation, reducing the state preparation error, or reducing the dephasing rate. By doing so, the optimum could be reached for N > 2.
Optimal control filters via Fisher information. By measuring the survival probability at the end of the filter sequence, we can indeed determine the decoherence function χ(t) and the following functional   Note that we can choose the cut-off frequency ω c in such a way that the support of the signal is not truncated, c . Alternatively, we can add an additional term to the control objective ξ imposing a constraint on the optimized filters such that they have a vanishing contribution for ω > ω c .
As shown in the Methods section, this optimization procedure corresponds to the maximization of the Fisher Information associated with the sensitivity of the measurement in the direction (within its functional space 50 ) of S(ω), and thus to the maximum sensitivity of the filter with respect to the signal. A similar approach for single-parameter estimation has allowed to determine the characteristic width of the spectrum with maximum sensitivity 47 .
We test the optimization numerically by trying to design a filter that maximizes the functional ξ for two examples of a given spectral density S(ω): we fix S(ω) and perform the optimization of the control pulses via the DCRAB algorithm 54 that straightforwardly allows to include constraints on the pulse bandwidth (or ω c ) as well as the additional constraint reflecting the use of π-pulses for single or multi-qubit probes. Figure 9 shows the fidelity ξ S / c for the Lorentzian power spectral density as a function of the number N of qubits employed for the filtering. The green horizontal line shows the near-unity fidelity obtained in the ideal limit of a continuous pulse modulation (an infinite number of qubits or a continuous modulation of Ω(t) -see Methods). The inset shows how the fidelity scales with the pulse duration T for 1 qubit (yellow) and 4 qubits (blue). Note that in both cases the fidelity peaks at T = 5. We choose this optimal value for the time T also for all other optimal-control filter (OCF) results presented hereafter. However, in general this optimal value of T depends on the desired frequency resolution: if the signal is less smooth (as a function of ω) than in our example, the optimal value of T will be larger; if it is smoother, the optimal T will be smaller. Figure 10 shows the corresponding approximations of S(ω) obtained by 1 (yellow) and 6 (blue) qubits, respectively. The green line shows the result in the limit of a continuous pulse modulation. The insets show the shape of the modulation functions for the two discrete cases. Figure 11 shows corresponding the results for an example of a double Lorentzian power spectral density Fast detection of time-dependent spectra. We now consider a complex dynamical system that can exhibit two different competing processes, each leading to a different spectral behaviour of an external field (signal). We further assume that the system is not stable but rather oscillating between the two regimes and model the composite signal as  Here, S 1 (ω) and S 2 (ω) are the signals corresponding to the two regimes and S 1,0 , S 2,0 are normalization constants. The time-dependent coefficients s 1 (t) and s 2 (t) model the oscillation between the two spectral components. Figure 12 shows the signal S(ω, t) for different time instances corresponding to the value of the coefficients s 1 (t) = 1, s 2 (t) = 0 (black), s 1 (t) = 0, s 2 (t) = 1 (red), and s 1 (t) = 0.5, s 2 (t) = 0.5 (green).  In this section, we want to study how these time-dependent coefficients can be measured by the methods developed in this paper, under the assumption that we know the two single components S 1 (ω) and S 2 (ω) (e.g. by measuring them in a stable configuration). We then analyse two approaches. In the first one (FO) we apply a sequence of basis filter functions and reconstruct the signal by orthogonalization. The coefficients are obtained by comparing the reconstructed signal with the ansatz of Eq. (26). In the simulations we apply filters of T = 5 duration and apply K = 10 basis filter functions, leading to a total duration of T c = 50 for the measurement of one time instance of the coefficients s 1 and s 2 (we assume an instantaneous population measurement). We then repeat this procedure until we cover the time span we want to investigate. In the second approach, we use the OCF protocol discussed above. We first apply the filter optimized for S 1 (ω) and then the one for S 2 (ω). Later, we use the measured overlap with the signal to calculate s 1 and s 2 . Again, we use filter operation times of T = 5, such that total duration is T c = 10 for one time instance. Thus it is 5 times faster compared to the first approach. In the simulations we test it for 1 and for 6 qubits.
We show the results for two different values of ω osc (0.004π and 0.01π) and analyse the time interval t ∈ [0, 500]. The first approach can thus sample 10 points, while the second approach 50 points. Figure 13 (for ω osc = 0.004π) and Fig. 14 (for ω osc = 0.01π) show the oscillation of the coefficient s 2 (t) as measured via the two approaches. For the slower oscillation, the sampling rate of the first scheme (FO) is high enough to resolve the oscillation, and both schemes work well, while the second approach has a clear advantage when operated with 6 qubits. For the faster oscillation, however, the sampling rate of the first approach (FO) is too small, and the measured points do not correspond to the signal. The second approach (OCF), instead still works very well.

Discussion
We have introduced new concepts in quantum sensing of the power spectral density of a signal, able to considerably speed up the sensing protocols and make them more robust against noise affecting the probe and the signal. In detail, we have introduced the concept of filter orthogonalization allowing us to use more arbitrary shapes of the filter functions, that might have future applications in the detection of multiple baths by distinguishing different single-bath contributions. Then, we have studied a new type of error model that includes probe dephasing and finite detector precision and thus shown the robustness of the protocols in realistic experimental situations. We have also included statistical noise and analyzed its effect on the ultimate precision in terms of the Fisher information. In this regard, we have expanded the approach of 47 , where the signal is parameterized and the precision of the parameter measurement is optimized: By means of the concept of the Fisher Information Operator we provide a tool to optimize the detection of an arbitrary functional shape of the signal. We have implemented this approach by optimizing the filter overlap with a given signal; this extends optimal control theory to filter design in quantum noise spectroscopy. As an application we have demonstrated that this new method provides a fast way to measure the contribution of two components of a composite signal. We have given two explicit examples of possible experimental realizations, both for single qubits (with NV-centers in diamond) and entangled multi-qubits (with trapped ions), using parameters from existing experiments, showing thus that an experimental realization of our new algorithms is feasible and advantageous with current technology and set-ups. Finally, we have shown how multi-qubit probes and entanglement offer similar advantages as more complex control of the probe 49 , and we believe that combining the two in a future work could provide new impulses to the study of networks of entangled quantum probes and the detection of both temporal and spatial resolutions 55,56 .

Methods
Estimation of the decoherence function. The presented derivation follows roughly 20-22 . We study a twolevel system (with levels |0〉 and |1〉) governed by the Hamiltonian where H c and H S are respectively the control and signal Hamiltonian, Ω(t) the applied control field, while E(t) a stochastic field stemming from an interaction with the environment. Moreover, in Eq. (28) σ x = |0〉 〈1| + |1〉 〈0| and σ z = |1〉 〈1| − |0〉 〈0| are Pauli operators. The system is initially prepared in the state |0〉, which is also the state that we can measure. An initial π/2-pulse brings the population into the state |↑ = , ω osc = 0.004π) and as reconstructed by the two approaches FO (for 1 qubit) and OCF (for 1 and 6 qubits). , ω osc = 0.01π) and as reconstructed by the two approaches FO (for 1 qubit) and OCF (for 1 and 6 qubits). 0 1 2 In this case, the pulse modulation function takes the discrete values y(t) ∈ {−2, 0, 2}. Finally, to generalize this estimation procedure to a multi-qubits framework, we consider the GHZ state with modulation functions y i (t) for each of the single qubits and y(t) ∈ {−N, 2 − N, …, N − 2, N}. The final filter function for N qubits, then, reads rely on performing measurements of the survival probability p(t) for different filter functions. The estimate ω  S( ) will generally deviate from S(ω) and we want to analyse how large can be such a deviation in a given direction ω S( ) in function space. If we assume that the main source of error is the statistical noise in the measurement of the survival probability and that on each realization of the measurement we obtain binary outcomes (survival or not), we can estimate the error ε ω S( ) (with ε being a constant quantifying the deviation of ω  S( ) from S(ω) in the direction ω S( )) by means of the Fisher Information Operator (FIO) as given in 50 . The latter is a straightforward generalization of the Fisher Information matrix to an (infinite dimensional) function space. The single FIO F k , associated with the measurement of the k-th probability p k after the application of the filter function F k , can be obtained by computing the derivative of p k with respect to χ k (t) (decoherence function of the probability p k ), i.e.  formally dividing by δS(ω) yields the functional derivative δχ k /δS, which is an element of the dual space of the tangent space of the spectral density functions, thus a linear mapping from the admissible changes δS(ω) to a real number δχ k . We have expressed this fact by means of the ket notation 〈⋅| as above.
Then, by using the chain rule, the FIO F k is equal to Finally, if we perform a sequence of measurements of p(t) by using several filter functions F k , k = 1, …, K, then the overall FIO F (for the sequence of measurements) is additive, since the measurements are statistically independent, and, as a result, we get The rank of this FIO F corresponds to the number of linear independent vectors ⟨ S k | χ ∂ ∂ , which is equal to the number of linear independent filter functions F k adopted for the estimation of S(ω). Thus, for an adequate choice of the pulse sequences (i.e. of the filter functions) the rank of the FIO equals to K. This means, that using such K filter functions we can determine the behaviour of the spectrum in a K-dimensional subspace of the function space of the spectral densities. This is a generalization of the common approach (e.g. for the AS protocol), where K filter functions are used to determine the value of the spectrum at K discrete values of the frequency. The behaviour of the spectrum outside of this subspace (or, in the discrete case, for any other value of the frequency) is completely inaccessible by the given choice of filter functions.
Cramér-Rao bound. A bound for the error ε is given by the Cramér-Rao bound associated with the Fisher Information for the parameter estimation of ε. As shown in 50 , this Fisher information is obtained from the FIO as where S( ) ω is, as before, a given direction in the function space of the spectral densities. The Cramér-Rao bound for ε is then given by ε ≥ . F

(54)
S Note that this bound is finite iff ω S( ) has finite overlap with at least one filter function F k (ω). This corresponds to the fact that the application of K (linearly independent) filters determines the signal S(ω) in a K-dimensional subspace of the full function space.