Optimal control of a quantum sensor: A fast algorithm based on an analytic solution

Quantum sensors can show unprecedented sensitivities, provided they are controlled in a very specific, optimal way. Here, we consider a spin sensor of time-varying fields in the presence of dephasing noise, and we show that the problem of finding the pulsed control field that optimizes the sensitivity (i.e., the smallest detectable signal) can be mapped to the determination of the ground state of a spin chain. We find an approximate but analytic solution of this problem, which provides a lower bound for the sensitivity and a pulsed control very close to optimal, which we further use as initial guess for realizing a fast simulated annealing algorithm. We experimentally demonstrate the sensitivity improvement for a spin-qubit magnetometer based on a nitrogen-vacancy center in diamond.


Introduction
Quantum systems are notoriously sensitive to external influences.This sensitivity is a core element in the development of quantum technologies, as is the case of quantum sensing, which takes advantage of quantum coherence to detect weak or nanoscale signals.Quantum sensing devices can in principle attain precision, accuracy, and repeatability reaching fundamental limits [1,2].However, the extreme sensitivity to external perturbations also causes the quantum sensor to couple with detrimental noise sources that induce decoherence, therefore limiting the interaction time with the target signal.
Here, we introduce a method to find optimal control protocols [3] for ac quantum sensing in the presence of dephasing noise.The sensitivity, i.e., the smallest detectable signal, is the figure of merit of the optimization, provided that the spectra of the target signal and of the noise are known.Such optimization problem is in general a complex classical problem.Our method, that draws an analogy between pulsed dynamical decoupling (DD) protocols [4][5][6][7][8] and spin glass systems [9], maximizes the phase acquired by the quantum sensor due to the target ac field while minimizing the noise detrimental effect.The optimal control fields yield an improved sensitivity with respect to commonly used protocols, as we experimentally demonstrate using a spin-qubit magnetometer based on a Nitrogen-Vacancy (NV) center in diamond [10][11][12][13][14].
More in detail, we find that the problem of optimizing the control protocol for our quantum sensor (i.e., optimizing its sensitivity) is homologous to that of finding the ground state of a classical Ising spin Hamiltonian, as depicted in Fig. 1.The control π-pulse times correspond to the locations in the chain of the domain walls.The couplings between the model spins, which encode the noise autocorrelation, are of both signs, and this is customary in optimization problems.The antiferromagnetic couplings capture the frustration between the different terms in the Hamiltonian, which then prima facie is that of a spin-glass model-which does not mean that there is a spin-glass phase at low temperature (see later).
The study of optimization problems in statistical physics is a large field of research in disordered systems, with far-reaching connections to the physics of spin glasses [15,16] and other frustrated, classical and quantum models [17][18][19][20][21][22][23].Optimization problems in quantum control can show some degree of frustration, with terms that compete in a similar way in Figure 1: (a) A single spin sensor is used to detect an AC target magnetic field b(t ).(b) An optimal control field applied to the spin sensor increases its coherence, hence improving its sensitivity.(c) The difficult problem of finding an optimal control sequence can be mapped (upon time discretization) into a problem of finding the ground state of a virtual spin chain.which ferromagnetic and anti-ferromagnetic bonds compete in spin glasses [24].We find, however, that in the specific case of the optimal control of a qubit sensor, by trading the Ising 2 spins for the continuous spins of a spherical model (SM) [25,26] one gets rid of frustration altogether, and the model shows little signs of competing equilibria at low temperature, typical of replica-symmetry-broken phases [9,27].Since the ground state of the spherical model can be found analytically if the spectra of the signal and of the noise are known, we obtain from this both a lower bound for the sensitivity 1 , and a quasi-optimal controlled pulsed field.This quasi-optimal sequence can then be fed to a simulated annealing (SA) algorithm [28][29][30][31], in order to find the optimal one with little computational effort.Such annealed sequences show, in agreement with the experiments, very good sensitivities (only about 20% worse than the bound).Our method is thus superior to standard DD protocols as Carr-Purcell (CP), or minimal generalizations of the latter to colored signals, as we discuss in detail.Finally, to show the unparalleled performance of the algorithm, which can open the door to real-time optimization in sensing, we run it on a Raspberry Pi microcomputer, where it takes milliseconds to find the optimal solutions.

Optimized dynamical decoupling for sensing
We consider a single spin-qubit sensor of time-varying magnetic fields, in the presence of dephasing noise.This quantum sensing task can be described as a compromise between spin phase accumulation due to the external target field to be measured b(t ) ≡ bh(t ), and refocusing of the non-Markovian noise, obtained via dynamical decoupling (DD) protocols [4][5][6][7][8].Above, b is the magnetic field strength to be detected, and h(t ) a known, dimensionless function specifying its time dependence.Notice that optimizing the sensor's response to a target field with a time dependence h(t ) that is known beforehand has several applications, e.g., sensing spins ensembles or spins in molecules [1,32,33], or sensing weak signals coming from biological samples-such as action potentials [34].Moreover, our optimization method can be part of a protocol where the time-dependence is first (sub-optimally) determined and then then the amplitude detected optimally.
As in Hahn's echo [35][36][37], a DD sequence is implemented by applying a series of n π-pulses that act as time reversal for the phase acquired by the qubit during its free evolution, and can be described by a modulation function y : [0, T ] ∋ t → {−1, 1} (see Fig. 1b).The DD sequence is embedded within a Ramsey interferometer, hence the qubit coherence is mapped onto the probability of the qubit to stay its initial state |0〉: Here, ϕ is the phase acquired by the qubit during the sensing time T : with γ the coupling to the field (e.g., the electronic gyromagnetic ratio of the spin sensor).
The noise-induced decoherence function is the convolution between the noise spectral density (NSD) S(ω) and the filter function Y (T, ω) = iω T 0 d t e −iωt y(t ).Note that we neglect the effect of the target field on the noise source [38] and we assume the noise to be a stationary Gaussian process.
Dynamical decoupling is a very versatile control technique, with a virtually infinite space of degrees of freedom spanned by all the possible distributions of π pulses, even at finite sensing time T .One of the most common DD sequences is the Carr-Purcell (CP) sequence [36,37], formed by a set of equidistant pulses.Non-equidistant sequences have been proposed and experimentally tested, e.g. in Refs.[7,[39][40][41][42][43].Each of these sequences has internal degrees of freedom, that can be tuned to increase the sensing capabilities for specific target fields.Another example is what we call the "generalized Carr-Purcell" (gCP) protocol, in which π pulses are applied when the signal b(t ) changes sign, i.e. in correspondence to its zeros.All these DD sequences are already optimal for simple target fields that happen to be well off-resonance with the noise, both for the fundamental frequencies and for all the higher harmonics.However, as the complexity of the target field increases, it increases also the difficulty to find a pulse sequence that successfully filters out the noise components, while still maintaining the sensitivity to the target field.
A possible approach is to use an optimization algorithm, to find a π-pulse sequence that optimizes a desired figure of merit, for example the sensitivity, i.e. the smallest detectable signal.Indeed, in a slope detection protocol [1] (suited for DD and Ramsey experiments), the sensitivity η is the minimum variation of b that can be measured, and it is defined as [1,44] (see also Appendix A for a derivation).The sensitivity η is therefore a compromise between noise cancellation (minimizing χ(T )) and target ac field encoding (maximizing ϕ(T, b)), and it is hard to optimize.This concept was proposed and demonstrated experimentally for an NV center used as a quantum magnetometer [45].Despite the achieved improvements, the computational complexity of the above optimization problem limited its applicability.
In our approach, instead, we recast the cost function η as the Hamiltonian of a classical Ising spin system.In this way, the continuous optimization problem for the minimization of the sensitivity of a NV-center magnetometer is re-interpreted as a discrete energy minimization problem.Specifically, we define the new cost function to be the (dimensionless) logarithmic sensitivity and we show in Sec. 3 that upon time discretization ε becomes an Ising Hamiltonian, albeit with sign-alternating, long-range interactions and a peculiar logarithmic field-spin coupling.Before doing that, however, we show how the problem can be tackled in continuous time, and by means of a reasonable approximation.

A variational approach
Our task is to find the optimal function y(t ) which minimizes the sensitivity η, Eq. ( 4), or the logarithmic sensitivity ε, Eq. ( 5).First of all, we anticipate why simple choices for y(t ) do not yield good results for generic sensing tasks.Looking at Eqs. ( 2)-( 4), one understands that the minimum detectable signal η is determined by a competition of the signal, through ϕ, and the noise, through χ.Commonly used DD protocols, as CP sequences, focus only on the properties of the signal, trying to amplify it irrespective of the noise (or assuming a zero-to-low-frequency noise).So, either using a CP sequence to amplify one frequency the signal is composed of, or taking y(t ) ∝ h(t ) to mimic as close possible the signal (the strategy dubbed gCP above), fail when the noise and the signal share common frequencies.Nevertheless, with the procedure outlined below, we show how it is possible to "orthogonalize" the DD sequence with respect to the noise to minimize the overlap χ, while keeping it "parallel" to the signal to maximize ϕ.In passing, we obtain useful analytical results that allow us to assess the performance of our method.Let us rewrite ε as [see Eqs. ( 2),( 3) and ( 5)] with the noise autocorrelation function, which depends only on the difference t − t ′ by stationarity of the noise.Notice also that J is a positive operator even though J(t ′ , t ) can take up any values in .Then, in order to find y(t ) that minimizes ε, we start by imposing the constraint y(t ) 2 = 1 for all t via a continuous set of Lagrange multipliers, i.e. via a function λ(t ): We need to find the stationary point of F [y] with respect to y(t ) and λ(t ).Formally, the saddle point equations are One can see that the extreme with respect to λ simply gives the constraint.The formal solution of the above equations is where λ stands for the diagonal operator λ(t )δ(t − t ′ ), and Above, with the notation 1 J+λ t ,t ′ we indicate the kernel of the operator inverse in the space of linear operators acting on square-integrable functions, and the quantity D can be interpreted as a self-consistent normalization for y(t ).By plugging Eq. ( 11) in Eq. ( 8), one can express the cost function at the saddle as The last expression is a function of λ(t ) only and one can, in principle, find its saddle point and substitute it in Eq. ( 11) to obtain the optimum DD sequence.Short of solving exactly the model in Eq. ( 8), we can get good results to guide the experiment by simplifying the space in which we are searching for the minimum.We can do this in two ways: either we keep y(t ) defined on (i.e.we keep the time continuum) and we give more structure to λ(t ), or we discretize time and enforce the constraint y(t )2 = 1 exactly (therefore getting rid of λ).These two approaches will be implemented in the following.

Spherical approximation
In order to make progress, we substitute for the moment the constraint y(t ) 2 = 1, for all t , with the constraint 1 This is equivalent to finding the stationary point of F [y, λ], Eq. ( 8), in the subspace in which λ(t ) ≡ λ ≡ const.Inspired by the physics of spin glasses [25,26], we call the resulting approximation spherical model (SM) 2 .Spherical models are often good mean field models of spin glasses and of their dynamics [25,26,46], and this case will prove to be of similar nature despite the unusual logarithmic field coupling term.By setting λ(t ) = λ we obtain a function of the single parameter The continuous spins s i ∈ (blue line) can be converted into Ising spins s i = ±1, necessary for the π-pulses, by using the sign function (orange line): this step corresponds to the projection (arrow) in panel (a).The sensitivity can be improved further with a few iterations of SA to get a close-by sequence (dashed black line).In this example, the trichromatic target signal and the noise are equal to the ones used in the experiment (see text).
where J + λ is the operator with integral kernel J(t ′ , t ) + λδ(t ′ − t ), as above.Minimizing with respect to λ ∈ , one finds η SM = e ε SM /γ T .This is a lower bound for the sensitivity, η > η SM , because the minimum of ε SM is found by searching for a y(t ) over a larger space of functions (Eq.( 15) is weaker than constraint y(t ) 2 = 1), as shown schematically in Fig. 2.
In principle the bound is not sharp, however it provides a quick and accurate measure of the goodness of our results.Moreover, we have found by experience that it is in practice pretty close to being sharp and that it can hardly be improved analytically by adding more freedom to the function λ(t ) beyond the constant λ(t ) = λ.For example the test function giving a two-parameters space (λ 1 , λ 2 ) for minimization, gives at most a few percent increase on the bound on η.We therefore use it as if it were sharp.One can define for any DD sequence the dimensionless quantity η SM /η < 1.We will see in the next section how different approximate solutions give different values of this quantity.Moreover, we will see how the solution of the SM, although not per se a DD sequence, can function as a starting point for finding an optimal DD sequence.

Time discretization and Simulated Annealing
Let us now focus on the second method: time discretization.We discretize the sensing time T into small time intervals ∆t , to obtain a sequence of times t i = i∆t with i ∈ 1, ..., N = T /∆t .The interval ∆t is the smallest time we allow the π-pulses of the DD sequence to be separated by.Apart from the physical limit given by the experimental apparatus, which sets a minimum ∆t , one does not expect to need in the optimal solution π-pulses separated by much less than the minimum period of h(t ), if it exists (the spectrum of h(t ) can extend up to infinite frequency).The modulation function at each of these times is y(t i ) = ±1, which dictates the sign of the phase acquired by the spin qubit during the time interval [t i − ∆t ,t i ].We can therefore write the modulation function as where s i = ±1, and as before χ J i j s i s j (19) where represents the interaction with a normalized target ac field, and represents the interaction with the detrimental noise.We can now express the new cost function as this closely resembles the Hamiltonian of the Ising spin glass problem for a set of N spins s i .
The ground state for this Hamiltonian can be used to obtain a modulation function, therefore a DD sequence, that minimizes the sensitivity η.
At first sight, minimizing ε in Eq. ( 22) on the hypercube {s i } ∈ {−1, 1} N seems a difficult problem, since the couplings J i j can be of both signs.Therefore, one is tempted to use a simulated annealing (SA) minimization algorithm [28][29][30] to find the minimum of the energy ε.However, the starting configuration strongly affects the performance of SA, both its final value and, at least as importantly, the time to reach it.With this in mind we turn to the SM solved in the previous section but with our discretized time, in terms of which the spherical constraint reads N i=1 y 2 i = N.In the discretized form, the solution of the SM is (see Eq. ( 11)) Above, we introduced the Fourier transform of the signal term hk = 1 N j e −i 2πk N j h j , and of the noise term Jk = 1 N j e −i 2πk N j J i,i− j : indeed, since λ(t ) is constant and J i j depends only on the difference i − j , the matrix J + λ is diagonal in Fourier space 3 .The value of λ is chosen to enforce the spherical constraint, and , see Eq. (12).One can notice that in Fourier space the optimal solution is aligned with the field, and orthogonal to the noise.
An example solution is shown in Fig. 2. The values of y i do not form a sequence of ±1, but the solution is reasonably close to the minimum of the exact functional, Eq. ( 22), over the hypercube {−1, 1} N .We can now use the solution in Eq. ( 23) as a starting point to find the optimal sequence.To do so, we first define s i = sign(y i ) ∈ {−1, 1} and then run few steps of SA moving only the domain walls, i.e. flipping only spins which are on a sign change: s i = −s i+1 .The π-pulse sequence is, as before, the sequence of times where the spins change sign (the position of the domain walls in the spin chain).
We test our procedure on an ensemble of test cases constructed as follows.The signal is a superposition of monochromatic waves h(t ) = N freq n=1 A n cos(ω n t + φ n ): we fix N freq = 7 and extract uniformly random frequencies in the interval [0, 1] MHz, uniformly random phases φ n , and uniformly random amplitudes A n s.t.
The noise spectrum is instead a gaussian centered at 0.4316 MHz, and with standard deviation 0.016 MHz: thus, it is close to (but a little bit stronger than) the experimentally relevant situation discussed in the next section.Interestingly, there is a strong overlap between the signal and noise frequency bandwidths.
First, we use the generalized Carr-Purcell (gCP) protocol introduced above.This procedure is simple but not very effective: on average, it returns a sensitivity between 1.5 to 3 times as large as the optimal one, monotonically increasing with the time of the sampling (see Fig. 3a).The sensitivity degradation with time is caused by the fact that the gCP sequences do not take into account the dephasing noise.Hence, as time increases the accumulation of noise by the sensor degrades its sensitivity.Second, we use the solution of the SM, that is, s i = sign(y i ), as DD sequence: this gives a better solution, as the sequence attempts to partially filter out the noise, but it is still not optimal.The best results, however, are obtained by running a fixed number of steps of SA starting from either a random DD sequence (SA, more on this below), from the gCP DD sequence (gCP+SA), or from the sign(SM) DD sequence (sign(SM)+SA).All these three cases perform the best because the SA algorithm is able to find good local minima of the optimization landscape.As it is seen in Fig. 3a, the sign(SM)+SA sequence gives the overall best result, with a solution close to the upper bound given by the SM itself (before projecting on the hypercube).It is important to stress that, although the ratio η SM /η for sign(SM)+SA is close to be constant as a function of time, eventually the sensor will not be able to detect any signals due to decoherence beyond dephasing (not considered in our model), e.g., T is limited by the spin-lattice relaxation time T 1 .For example, at room temperature T 1 ≃ 1 ms for NV spin sensors.
We can understand the comparative performance of the different control protocols as follows: The gCP case performs the worst because it ignores the noise, leading to dephasing as time increases.This effect is mitigated by the sign(SM) case, that accounts for some effect of the noise.The SA cases perform the best since by construction they represent good local minima of the optimization landscape, found by the numerical sampling procedure.
We can finally compare sign(SM)+SA with the "unbiased" SA optimization, which starts at infinite temperature from a uniformly random sequence of s i = ±1.Not only it is outperformed by sign(SM)+SA (especially in convergence time), but its solutions typically require more π pulses than needed (see the following).In this case, to reduce the number of π pulses it is possible to introduce by hand a ferromagnetic coupling term in the Hamiltonian: with K > 0 to be tuned.One can see in Fig. 3b that the best sensitivity is however still obtained with the combination of the SM solution and SA optimization.Additionally, from Fig. 3b one can also understand that the optimal solution represents the best trade-off between number of π pulses (which the experimenter would like to maintain low) and sensitivity.
Figure 3: (a) Comparison of performances, over a broad ensemble of parameters, for the DD sequences discussed in the main text: generalized Carr-Purcell (gCP), spherical model projected with the sign function onto the hypercube (sign(SM)), simulated annealing (SA), and SA optimization starting from gCP and sign(SM).One can see that the best results are obtained for the SA optimization guided by the SM solution.The data refer to the ensemble of random test signals described in the main text: the dots are the average values, and the shaded area represents the 20-80 percentile of the distribution of results.The discretization interval is ∆t = 0.1 µs.(b) Single instance of a random signal, corresponding to T = 100 µs.We show the average sensitivity and the number of π pulses (the error bars correspond to one standard deviation over the ensemble of annealing realizations) of the solutions coming from the unbiased SA, i.e. starting from infinite temperature (purple to green circles), and from the SA guided by the SM solution (red square).The unbiased SA needs a ferromagnetic term ∝ K , see Eq. (24), with K to be optimized over, in order to keep under control the number of π pulses.From this plot, one learns that first, the optimal solution represents also the best trade-off between number of π pulses and sensitivity, and second, that the SA optimization guided from the SM performs better, and with less fluctuations.Here, each unbiased SA procedure uses 10 5 Monte Carlo steps and a power-law temperature ramp, while only 10 3 steps are needed for the SA from the SM solution.
To conclude, we stress that our optimization procedure is very fast, if compared to standard, general-purpose routines.In particular, we were able to run our codes on a Raspberry Pi microcomputer, where the single instance takes ∼ 0.5 s for the unbiased SA algorithm, and ∼ 0.02 s for the solution of the SM and subsequent annealing (using N = 500 spins).Taking in consideration that few instances of the sign(SM)+SA protocol are sufficient to obtain a good result, while the optimization over the parameter K requires hundreds, if not thousands, of separate SA runs, the gain provided by our method becomes apparent.This fact also opens the door to the miniaturization of the control electronics, in view of possible technological applications of quantum sensing.

Experiment
While our method is general and applicable to any spin-qubit sensor, we exemplify it through experiments with a single NV center in bulk diamond with naturally abundant 13 C nuclear spins, at room temperature.The ground state electron spin of the NV center can be initialized and measured by exploiting spin-dependent fluorescence, and can be coherently manipulated by microwaves [14].We consider the two ground-state spin levels, m S = 0 and m S = +1, to form the computational basis of the qubit sensor {|0〉, |1〉} (see Appendix B).The main source of noise for the NV spin qubit derives from the collective effect of 13 C impurities randomly oriented in the diamond lattice.
In the presence of a relatively high bias field (≳ 150 G), the collective effect of the nuclear spin bath on the NV spin is effectively described as a classical stochastic field, with gaussian noise spectral density (NSD) centered at the 13 C Larmor frequency ν L = γB [48,49], where γ is the 13 C gyromagnetic ratio.We preliminarily characterize the NSD of the NV spin sensor as in Ref. [49].The direct coupling between the target field and the nuclear spins is negligible due to the small nuclear magnetic moment [38], and the indirect coupling via the NV electronic spin is also negligible due to the presence of the strong bias field [49].Therefore, the NV spin dynamics is well described by Eq. (1).For the experiments we present throughout this article we used a bias magnetic field of 403.As a test case for our optimal control method versus standard control, we consider a threechromatic target signal, with h(t ) = +1 i=−1 A i cos(2πν i t ), where ν i = {0.1150,0.2125, 0.1450} MHz are the frequency components, and A i = {0.288,0.335, 0.377} are the relative amplitudes, respectively for i = −1, 0, +1.
In Fig. 4(a) we show the NV spin coherence P(nτ, b) under Carr-Purcell (CP)-type DD control, formed by n pulses with uniform interpulse spacing τ = T /n, as a function of τ.The value of b at the position of the NV defect inside the diamond is obtained from minimizing the squared residuals between experiment (gray bullets) and simulation (gray line), for which b is the only free parameter (see Appendix B.1 for more details).
The CP pulse sequence acts as a quasi-monochromatic filter centered at 1/τ, so that a single component of b(t ) can be sensed in each experimental realization.As a consequence, P(nτ, b) in Fig. 4(a) shows collapses occurring at τ ∼ 1/2ν i .Notice that the collapse corresponding to the frequency component ν +1 (τ ≃ 3.448 µs) cannot be resolved from noise since the first harmonic of the filter function roughly coincides with the NSD peak (ν +1 ≃ ν L /3) [Fig.4(b)].To detect the three components of the target signal and filter out the NSD, we need an optimized sequence.We thus apply the optimization algorithm detailed before to solve this experimental sensing problem.
In order to confirm the theoretical prediction on how the optimized DD sequence can outperform the standard control in terms of sensitivity, we performed measurements of the sensitivity itself.Specifically we used three different CP sequences, each with time between pulses τ = 1 2ν i , for i = −1, 0, +1.Having a previous knowledge of the NSD allows us to predict the sensitivity of the the spin sensor using equations ( 3), (2), and (4), for any given DD sequence, and for any target AC signal b(t ).In Fig. 5a we show the estimated values for the inverse of the sensitivity as a function of the sensing time T = nτ.Since τ = 1 2ν i is fixed for each of the CP sequences, the variation of T corresponds to a variation of the number of pulses n.Notice how for τ = 1 2ν +1 , the inverse of the sensitivity rapidly goes to zero.The estimated inverse sensitivity for the optimized sequence sign(SM)+SA is also shown in Fig. 5a.The inverse sensitivity increases as a function of T , although we expect it to decrease at longer times due to decoherence.In particular we know that for NV spin qubits the spinlattice relaxation time T 1 ultimately limits the sensing time T .However, even at shorter times T < T 1 the sensitivity could be limited by other experimental factors, the most probable one being π-pulse imperfections.
In the experiment, we measure P(T, b) as a function of the field amplitude b at a fixed   1) and ( 2)), and therefore we can obtain the values of η using Eq. ( 4).The sensitivity measured experimentally shows an excellent agreement with the expected simulated values (see Fig 5a).See Appendix C for two additional test cases: one for a monochromatic target signal such that the fifth harmonic of the NSD coincides with the frequency of the target signal; and one for a target signal with seven frequency components, all close to the NSD peak.These two cases confirm the results of the experiments shown in the main text.
The sensitivity reported in Fig. 5a is the result of the optimized DD control sequence assuming that the measurement of b follows a slope detection protocol [see Sec.IV.E in Ref. [1]].Therefore, to estimate the sensitivity it is sufficient to obtain the maximum slope of P(T, b) as a function of b.In order to achieve a precise estimation of 1/η, we measured more than one period of P(T, b) [Fig.5b].Notice that the same DD control sequence would be used to measure the value of b.One just needs to measure P(T, b), in the range around P(T, b) ≃ 0.5.Outwardly, an improved sensitivity thus comes at the cost of a smaller bandwidth.Notice though that the bandwidth can be effectively extended by changing the range around which P(T, b) ≃ 0.5, which is achieved by tuning the phase of the oscillations of P(T, b) as a function of b (or, in practice, the phase of the final π/2 pulse).

Conclusion
We have shown that the problem of finding an optimal solution to quantum control a single spin system for quantum sensing can be solved by first finding the ground state of a solvable spherical model of classical spins, and then using this as a starting point for a simulated annealing algorithm.In this way, the optimization algorithm is able to find a control sequence that shows a significant improvement to the sensitivity with respect to standard control sequences.In addition, from the spherical model we found a theoretical bound on the sensitivity.Although the spherical model can be mapped to a control sequence that gives relatively good results, state, imposing a zero-field-splitting of D g ≃ 2.87 GHz.An external bias field, aligned with the spin quantization axis, removes the degeneracy between the m S = ±1 states, allowing to individually address the m S = 0 ↔ m S = +1 transition using on-resonance microwave radiation.By using microwave pulses with a appropriate duration, amplitude and phase, it is possible to apply any kind of gate to the single two level system.Therefore, the two level system formed by the m S = 0 (|0〉) and m S = +1 (|1〉) states fulfills the requirements to be used as a qubit based magnetometer.

B.1 Characterization of the amplitude of the target signal
The target signal is delivered via a signal radio-frequency (RF) generator connected to the same wire, placed close to the diamond, that delivers the MW control field.We can control the amplitude of the target field by changing the output amplitude of the RF generator.However, the absolute value of the amplitude of the target field b has to be characterized in order to take into account the attenuation of the circuit, the emission efficacy of the wire (which depends on the RF frequency) and the distance between the wire and the NV defect.To achieve such characterization, as explained in the main text, we measure the spin dynamics for a CP sequence as a function of the sequence interpulse time, and we compare with the simulation to minimize the residuals using b as the only free parameter.By performing this measurements for different values of the RF generator output amplitude a RF , we can extract a relation between a RF (in [Vpp]) and the amplitude of the target magnetic field b (in [T]).

C Additional test cases
In order to reinforce our results, we repeated the analysis presented in the main text for two different target signals.A monochromatic target signal that coincides with one of the NSD harmonics, and a 7-chromatic target signal that accentuates the difference between the generalized CP and the optimal solution.

C.1 Second test case: Monochromatic target signal
If we want to detect a monochromatic target signal b(t ), in most cases a Carr-Purcell CP sequence of equidistant pulses is the best way to increase the sensor's response to that target signal and filter out the noise.This is due to the quasi-monochromatic filter function associated with a CP sequence.Assuming that τ is the time between pulses, the filter function shows a peak centered at ω/2π = 1 2τ .However, the filter function is not exactly monochromatic, it shows harmonics at ω/2π = 1 2(2ℓ+1)τ , with ℓ ∈ {1, 2, ...}.Therefore, if the frequency associated with b(t ) is close to ω L /(2ℓ + 1), then a CP sequence will amplify the effect of both, the target signal and the noise, leading to not-optimal sensitivities.
Here we used the optimization algorithm described in the main text in order to obtain optimal sequences for this problem.In particular, we explored the case of a monochromatic signal with frequency ν mono = 39.29 kHz, which is close enough to ν L /11 so that the 5-th harmonic of the CP sequence coincides with the noise components.We used the same NSD S(ω) as in the three-chromatic case Similarly to the case detailed in the main text, the optimal sequences improve the sensitivity of the quantum sensor, resulting in some cases to an inverse sensitivity that is close to a twice the one from the CP sequence.In the Figure 6: Results for the case of a monochromatic target signal.(a) Probability to remain in the state |0〉 as a function of b, for fixed sensing times T for an optimal DD sequence (blue), and for a CP sequence (orange).The values of the sensing time and of the number of pulses for both sequences are shown as titles of the plots.A cosine function is fitted (solid lines) to the experimental data (bullets with errorbars) in order to obtain 1/η (see main text).(b) Inverse sensitivity as a function of the sensing T .Blue data corresponds to the optimized sequences obtained with simulated annealing (SA).Orange data corresponds to the CP sequences with τ = 12.726 µs.We found a good agreement between the predicted values (dotted lines) and the experimental values (bullets with errorbars).monochromatic case explored here, the sensitivity gets worse when increasing the sensing time beyond 100 µs.Instead the optimal solutions are able to improve the sensitivity even for times T > 300 µs.For T ≃ 100 µs, and longer sensing times, the optimized sequences achieve higher values of 1/η than the maximum value achieved by a CP sequence.

C.2 Third test case: 7-chromatic target signal
We have explored the case of a target signal with 7 frequency components, as specified in Fig. 7  (a-b).As in the main text, we used the optimization algorithm either to find the approximated spherical solution, or the solution using simulated annealing (SA) in order to minimize the sensitivity.The predicted values of the inverse sensitivity, together with their experimental values are shown in Fig. 7(c).Similarly to the previous test cases, the optimal sequences improve the sensitivity of our quantum sensor.In this case, the sensitivity obtained with the optimal solutions almost 1/2, and 1/3 with respect to the generalized CP (gCP) sequence for T = 80 µs, and T = 160 µs, respectively.

Figure 2 :
Figure 2: (a) Sketch of the spherical model.The N-dimensional, hyper-spherical surface (in blue) strictly contains the hypercube {+1, −1} N (black dots), each point of which encodes the configuration of classical spins in the Ising model.Therefore, the solution of the spherical model (red square) is, in general, not a point on the hypercube, but it can be projected (arrow) onto the latter, giving a good value for the sensitivity.(b) Comparison between solutions for 200 spins (T = 32 µs, ∆t = 0.16 µs).The continuous spins s i ∈ (blue line) can be converted into Ising spins s i = ±1, necessary for the π-pulses, by using the sign function (orange line): this step corresponds to the projection (arrow) in panel (a).The sensitivity can be improved further with a few iterations of SA to get a close-by sequence (dashed black line).In this example, the trichromatic target signal and the noise are equal to the ones used in the experiment (see text).
[a,b] is the characteristic function of an interval [a, b].Writing the modulation function in this way allows us to recast Eqs.(2) and (3) respectively as ϕ(T, b) = T γb

Figure 4 :
Figure 4: (a) Dynamics of the NV spin qubit under a DD sequence with n = 16 equidistant pulses for a trichromatic signal (see text).The NV spin coherence is mapped onto the probability of the NV spin to be in the state |0〉), P(nτ, b).Gray bullets: experimental data.Black dashed line: simulated spin coherence in the presence of noise, without any external target signals.Orange, red, and purple dashed lines: simulated spin coherence in the presence of monochromatic target fields with ω 1 , ω 2 , and ω 3 , respectively, with no noise.Gray solid line: simulated data combining all of the above using Eq.(1).Residuals between gray experimental data and gray solid line are shown in the bottom plot.(b) NSD given by the nuclear spin environment of the NV sensor (black line); fast Fourier transform (FFT) of the target signal h(t ) (gray line).Vertical dotted lines: frequency components of the target signal, and center of the NSD.Orange, purple, and red lines: filter function for a CP sequence with T = n 2ν i , for i = −1, 0, and +1, respectively.Blue line: filter function of the optimized sequence.Inset: examples of time distribution of π pulses.

Figure 5 :
Figure 5: (a) Experimental values of the inverse sensitivity for the optimized sequences (blue circles; for ∆t = 160 ns), and for the CP sequences (orange, red, and purple triangles).The predicted values of 1/η are represented by dotted lines.Black dashed line: theoretical upper bound of 1/η, obtained from the solution of the spherical model in the continuum limit.Solid black line: predicted values of 1/η for the gCP sequence.Inset: Ratio η SM /η for the sign(SM)+SA and for the gCP sequences (blue and black, respectively).(b) P(T, b) as a function of the amplitude of b, at T ≃ 152 µs.Same color code as in (a).Lines are a cosine fit (see text).
. The experimental values of 1/η are obtained from the measurement of P(T, b) as a function of b.The results of P(T, b) for one value of the sensing time T are shown in Fig. 6(a).The predicted values of the inverse sensitivity, together with their experimental values are shown in Fig. 6(b).

Figure 7 :
Figure 7: Results for the case of a target signal with seven frequency components.(a) Table to indicate the amplitude, frequency and phase of each component of the target signal h(t ) = 6 i=0 A i cos(2πν i t + c i ).(b) Fast Fourier transform (FFT) of the target signal.(c) Inverse sensitivity for T = 80 µs and 160 µs.The predicted values (squares) and the experimental values (bullets with errorbars) show that the sequences obtained from the spherical solution (Sph.) or from the simulated annealing solution (SA) result in an improved sensitivity with respect to the generalized CP (gCP) sequences.