Real-time adaptive estimation of decoherence timescales for a single qubit

Characterising the time over which quantum coherence survives is critical for any implementation of quantum bits, memories and sensors. The usual method for determining a quantum system's decoherence rate involves a suite of experiments probing the entire expected range of this parameter, and extracting the resulting estimation in post-processing. Here we present an adaptive multi-parameter Bayesian approach, based on a simple analytical update rule, to estimate the key decoherence timescales ($T_1$, $T_2^*$ and $T_2$) and the corresponding decay exponent of a quantum system in real time, using information gained in preceding experiments. This approach reduces the time required to reach a given uncertainty by a factor up to an order of magnitude, depending on the specific experiment, compared to the standard protocol of curve fitting. A further speed-up of a factor $\sim 2$ can be realised by performing our optimisation with respect to sensitivity as opposed to variance.


INTRODUCTION
Decoherence, resulting from the interaction of a quantum system with its environment, is a key performance indicator for qubits in quantum technologies [1] including quantum communication, computation and sensing.Decoherence timescales determine the storage time for quantum memories and quantum repeaters, a crucial metric for quantum communication networks [2][3][4][5].Rapid benchmarking of decoherence timescales in platforms such as superconducting qubits [6,7] or silicon spin qubits [8,9], is a critical validation and quality assurance step for the development of large-scale quantum computing architectures, and has the potential to improve error correction protocols efficiently close to faulttolerance thresholds.In quantum sensing, the role of decoherence is two-fold.On one side, decoherence sets the ultimate performance limit of the sensors [10].On the other hand, decoherence itself can be the quantity measured by a quantum sensor, as it provides information about the environment.An example of this is relaxometry, where the rate at which a polarised quantum sensor reaches the thermal equilibrium configuration gives information about different physical processes in the environment [11][12][13][14].
Decoherence rates can be measured by preparing the system into a known quantum state and probing it at varying time delays to determine the probability of decay from its initial state.The standard protocol for decoherence estimation involves a series of measurements with * c.bonato@hw.ac.uk time delays set over a pre-determined range, reflecting the expected value of the decoherence rate, and fitting of the result to a decay function.As the range of time delays is determined in advance, some of the measurements will provide little information on the decoherence of the system, since the time delay is either much shorter than the true decoherence rate, resulting in no decay, or much longer, resulting in complete decay.
Here we introduce an real-time adaptive protocol to measure decoherence timescales T 1 , T * 2 and T 2 for a single qubit [15], respectively corresponding to relaxation, dephasing, and echo decay time [1], together with the coherence decay exponent β.While the proposed algorithms are very general and can be applied to any quantum architecture, our experiments are implemented on a single spin qubit associated with a nitrogen-vacancy (NV) centre in diamond.
Adaptive techniques have been shown to be central to progress across a broad range of quantum technologies [16].Early work in this field involved the implementation of adaptive quantum phase estimation algorithms on photonic systems [17], later extended to frequency estimation with applications to (static) DC magnetometry with single electron spins [18][19][20].Alternative adaptive protocols for the estimation of static magnetic fields are based on sequential Bayesian experiment design [21,22] and ad-hoc heuristics [23], later applied to the characterisation of a single nuclear spin [24].Real-time adaptation of experimental settings has also been shown to be advantageous when measuring spin relaxation [25] or tracking the magnetic resonance of a single electron spin in real-time [26][27][28][29].Furthermore, adaptive techniques have been investigated to enhance photonic quantum sen-sors [30][31][32], as a control tool for quantum state preparation [33] and to extend quantum coherence of a qubit by manipulating the environment [34][35][36].
Despite these pioneering experiments, several important methodological questions still remain open.A priority concern is that adaptive protocols introduce an overhead, given by the time required to compute settings on the fly for the next iteration.It is crucial to minimise this computation time, since it can slow the protocol down to the point that the overhead can reverse the gain in measurement speed compared to a simple parameter sweep.This has not been considered in many cases, in particular where algorithms were investigated through computer simulations [21,22,27,37] or as off-line processing of pre-existing experimental data [23].While the optimisation of complex utility functions can possibly deliver the best theoretical results, this could be practically less advantageous than near-optimal approaches with very fast update rules in minimising total measurement duration.A second issue is related to the fact that, for multi-parameter Hamiltonian estimation, standard approaches such as the maximisation of Fisher information can fail, as the Fisher information matrix becomes singular when controlling the evolution time [38].This has stimulated researchers to find ad-hoc heuristics, for example, the particle guess heuristic [23,24,38] for the estimation of Hamiltonian terms; these heuristics, however, do not necessarily work beyond Hamiltonian estimation.A third question is related to what quantity should be optimised.Previous work has targeted the minimisation of the variance of the probability distribution for the quantity of interest [24,39].While this is clear when all measurements feature the same duration, the answer is less straightforward when adapting the probing time.If two measurements with different probing times result in a similar variance, the protocol should prefer the shorter one, minimising the overall sensing time.
Here we address these open questions, presenting theoretical and experimental data about the adaptive estimation of decoherence for a single qubit, using NV centres as a case study.Compared to other recent investigations of adaptive protocols [23][24][25], our experiments utilize a very simple analytical update rule based on the concept of Fisher information and the Cramér-Rao bound.By exploiting state-of-the-art fast electronics, we experimentally perform the real-time processing in than 50µs, an order of magnitude shorter than previous real-time experiments [24], negligible compared to the duration of each measurement.Such a short timescale makes our approach useful for qubits where fast single-shot readout is available such as trapped ions [40], superconducting qubits [41] and several types of spin qubits [42][43][44][45], and could be further shortened in future work by implementing the protocols on field-programmable gate array (FPGA) hardware.
In the case of multi-parameter estimation, previous work on Hamiltonian estimation had pointed out that the Cramér-Rao bound cannot be used in the optimisation as the Fisher information matrix is singular and cannot be inverted [38].Here we address this issue by utilising multiple probing times, showing that the Fisher information matrix can be inverted and that the corresponding adaptive scheme provides better performance than nonadaptive approaches.Finally, we discuss what quantity needs to be targeted to achieve the best sensor performance, experimentally demonstrating the superiority of optimizing sensitivity, defined as variance multiplied by time, over optimizing variance.As a figure of merit, sensitivity encourages faster measurements.
Our work tackles these general questions using the characterisation of decoherence as a test case.While adaptive approaches have been investigated in the case of phase and frequency estimation [17][18][19][20][21][22][23][24], also in relation to Hamiltonian learning [23], the case of decoherence is much less explored, with only one work targeting the estimation of the relaxation timescale T 1 [25].Here we provide the first complete characterisation of the three decoherence timescales typically used in experiments (T 1 , T * 2 and T 2 ), together with the decoherence decay exponent β.

II. THEORY
Decoherence and relaxation are processes induced by the interaction of a qubit with its environment, leading to random transitions between states or random phase accumulation during the evolution of the qubit.These processes are typically estimated by preparing a quantum state and tracking the probability of still measuring the initial state over time, which can be captured by a functional form [10] Although the noise processes induced by interaction with the environment can be complex, χ(t) can often be approximated by a simple power law: where T χ and β depend on the specific noise process [1].
For white noise, the decay is exponential with β = 1.For a generic 1/f q decay, relevant for example for superconducting qubits, with a noise spectral density as ∝ ω q , χ(t) scales as χ(t) ∝ (t/T χ ) 1+q [46].In the case of a single electronic spin dipolarly coupled to a diluted bath of nuclear spins, the decay exponents have been thoroughly investigated, with analytical solutions available for different parameter regimes [47].If the intra-bath coupling can be neglected, the free induction decay of a single spin is approximately Gaussian (β = 2) [48,49].The Hahn echo decay exponent T 2 can vary, typically between β ∼ 1.5 − 4 depending on the specific bath parameters and applied static magnetic field [47,50].(a) Schematic of the real-time adaptive protocol demonstrated in this work, using the electronic spin associated with a nitrogen-vacancy (NV) centre in diamond.An Arbitrary Waveform Generator (AWG) is used to generate pulses for manipulation of the spin qubit.The spin state is then optically measured and the detected photon count rate is used by the microcontroller to estimate the value of the decoherence timescale Tχ (and the decay exponent β) using Bayesian inference.The microcontroller also computes the optimal probing time τ for the next measurement, passing the value to the AWG, which builds the next estimation sequence accordingly.(b) In our experiment, a single interrogation of the qubit does not provide sufficient information to discriminate qubit state.Hence R measurements are performed, and the resulting number r of detected photons are used to update p(Tχ) through Bayes' rule and to compute τ for the next experiment.(c) Example of an experimental adaptive estimation sequence, with p(T * 2 ) shown for each measurement epoch (here the decay exponent is known a-priori as β = 2).The probability p(T * 2 ), initially uniform, converges towards a narrow peak as more measurement outcomes are accumulated.The experimental adaptively-optimised values for τ are shown on the bottom plot.
In Sec.II A, we assume the decay exponent β to be known, and we only focus on estimating the decay time T χ .This is a practically-relevant situation in cases where the nature of the bath is well-understood and the decay exponent β is known, at least approximately, a priori.We then extend our analysis to the simultaneous estimation of T χ and β in Sec.II B.
Fig. 1(a) sketches the operation of the real-time adaptive sensing system developed in our study.We have utilised the electronic spin associated with a nitrogenvacancy (NV) center in diamond as the qubit, which is initialised and readout by optical pulses.The qubit state is manipulated by microwave (MW) pulses, created in real-time by an Arbitrary Waveform Generator (AWG) based on an external digital input.After the application of a pulse sequence, the qubit is optically readout, with the spin state information enconded in the number photoluminescence photon counts during optical illumination.The core of our adaptive system is a real-time microcontroller, which uses the detected photon count rate to estimate the values of the decoherence timescale (T χ ) and the decay exponent (β) via Bayesian inference.As shown in the inset, the probability distribution starts out as uniformly flat, but begins to converge around the true value after a few iterations.Based on the estimated value in the current iteration, the microcontroller computes the optimal probing time (τ ) for the subsequent measurement and communicates this value to the AWG, which then constructs the next estimation sequence accordingly.This cycle repeats for several iterations until a desired level of error in the estimation of the target quantity is reached.Fig. 1(b) shows the flow of the experimental estimation sequence.For our experiments, a single measurement of the qubit lacks the information required for discriminating its state effectively.Therefore, we conduct R measurements, to obtain r detected photons, enough to discriminate the spin state.Such counts are then utilised to update the probability distribution p(T χ ) using Bayes' rule.After the Bayesian update, updated probability distribution is used to compute the optimal settings and provide feedback for the subsequent measurements.Fig. 1(c) shows an example of experimental estimation of T * 2 , performed by an adaptive Ramsey experiment, plotted as the evolution of p(T * 2 ) for increasing estimation epochs.In the beginning, p(T * 2 ) is a uniform distribution in the range 0-8 µs, which then converges to a singly-peaked distribution after more and more measurement outcomes are processed.In the case of an NV centre in a highpurity diamond, the decay is expected to be Gaussian (β = 2) [49].As described later in Sec.II A, the optimal adaptive rule for this case is to choose the probing time as τ opt ∼ 0.89 • T * 2 (see Eq. 20, where T * 2 is the current estimate of T * 2 computed from the probability distribution p(T * 2 ).The chosen values for τ are shown on the bottom plot, illustrating how they converge very fast to the optimal value τ opt ∼ 0.89 • (T * 2 ) true ∼ 2.23 µs. .

A. Adaptive Bayesian estimation
We utilize Bayesian inference, exploiting Bayes' theorem to update knowledge about the decoherence time T χ and decay exponent β in the light of a set of new measurement outcomes denoted by ⃗ m = {m 1 ,m 2 , . . .}. Thanks to its flexibility in accounting for experimental imperfections and for integrating real-time adaptation of the experimental setting while remaining easy to interpret mathematically, the Bayesian framework [16,51] has been widely applied in quantum technology, from sensing [19,[23][24][25], to the tuning of quantum circuits [52,53], and model learning [54].In this section, we will restrict the discussion to the characterisation of the decoherence time T χ ; the extension to a multi-parameter case, with the simultaneous estimation of T χ and β, will be presented in Sec.II B.
For each binary measurement outcome m n (m n = 0, 1), the probability distribution of T χ , which represents our knowledge about T χ , is updated as where m 1:n = {m 1 , . . ., m n }.Here, P (T χ |m 1:(n−1) ) is the posterior probability after (n−1)-th update and proceeds to serve as prior distribution at the n-th iteration, and P (m n |T χ ) is the likelihood function Note that this likelihood depends on τ (which we will adjust later) but this dependency is omitted in P (m|T χ ) to simplify the notation.Our approach to adaptive estimation is to derive a simple expression for the optimal parameter settings, that can be computed in real-time by an analytical formula without adding much extra computation time to the sensing process.
A conventional approach to updating τ adaptively would be to use the information gain as criterion [21,52,55].However, this involves integrals (with respect to T χ ) requiring numerical evaluation and an associated significant computational overhead.Here, we instead employ an approximation of the Bayesian information matrix (BIM) (a 1 × 1 matrix, in the case of a single parameter) [56] which links to the classical Fisher information [57].While computing the BIM also requires the computation of an integral (more precisely an expectation) with respect to T χ , this can now be easily approximated as explained in Appendix A.
The Cramér-Rao lower bound (CRLB) of T χ represents the minimum reachable variance for any (unbiased) estimator of T χ , and is inversely proportional to the Fisher information F .Thus, maximizing F with respect to the control parameter τ is expected to improve our estimate of T χ .
As shown in Appendix A, we examine a Bayesian form of the CRLB, computing the corresponding Fisher information F B can be computed as: where Tχ is a point estimate of T χ before each measurement.
While we are unable to maximize FB (τ ) analytically, approximate solutions exist.We found that the heuristic leads to satisfactory results, where ξ is a parameter that depends on β.Some numerically-computed values for ξ are listed in Table I, for some common values of β.
The Fisher information F as a function of the ground truth value for T χ = T * 2 and the probing time τ , from Eq. ( 17) is plotted in Fig. 2(a).F (τ, T * 2 ) is normalised by its maximum with respect to τ , for each value of T * 2 .The plot shows clearly that the maximum of F (τ ) has a linear dependence on T * 2 , following Eq.( 17) (shown as the red dashed line).

B. Multi-parameter estimation
In many practical situations, it is important to learn both the decoherence timescale, T χ , and the noise exponent, β, as the latter provides useful information about the nature of the qubit environment., assuming that single shot qubit readout is available: random choice of τ (orange), choice of τ obtained maximizing the Fisher information F in Eq. ( 17) (blue) and choice of τ obtained maximizing the rescaled Fisher information FT in Eq. ( 9) (red).The x-axis shows the total probing time (cumulative over epochs), while the y-axis shows the uncertainty, defined as the root mean squared error (RMSE) from the ground truth value.The solid black line corresponds to the theoretical limit set by the Cramér-Rao lower bound (CRLB) for the Fisher information in Eq. ( 17).Adaptive protocols outperform the non-adaptive (random) protocol.(c) Simulations for the same protocols as in (b), for the case where single shot readout is not available (using experimental values p cl (|0⟩) = 0.0187 and p cl (|1⟩) = 0.0148, and R = 50000).The solid black line corresponds to the limit set by the CRLB for the Fisher information in Eq. ( 17).The dashed black line shows the CRLB limit for the Fisher information FE in Eq. (10).Directly extending the approach for a single parameter estimate is unfruitful as the determinant of the BIM is zero, as reported previously [38].The intuitive reason for this is that the determination of two parameters from measurements in one single setting creates correlations between the two parameters, resulting in a singular Fisher information matrix.
We address this issue by using two consecutive measurement results, m and n, taken at times τ 0 and τ 1 .In this case, a non-zero determinant can be found by maximising the BIM on p(m|T χ , β) × p(n|T χ , β).Assuming binary measurement outcomes, the determinant of this matrix is given by: det FB = We did not find a way of calculating the maxima of this function analytically.Instead, here it is approximated numerically to determine the update heuristic.Starting by numerically estimating the value of τ 1 /T χ , which maximises Eq. ( 7) across a range of values of τ 0 /T χ and β, we can fit a piece-wise linear approximation for the best choice of τ 1 that works well for 1 < β < 5.This approximation is given by: C. Optimizing sensitivity The approach we introduced in the previous section aims at maximizing the BIM as a proxy for minimizing the uncertainty of T χ .However, this criterion does not account for the overall measurement duration.For example, if two measurements (obtained from two different values of τ ) lead to similar BIMs, the shorter one of the two should be more favourable, as our goal is to reach the smallest possible uncertainty in the shortest time.We therefore also consider an alternative approach which, instead of targeting the minimisation of the mean-squared error (MSE), is designed to improve the sensitivity, defined as η 2 = M SE • τ [37].This leads to a new criterion that is the BIM rescaled by the probing time: Again, numerically solving for (τ / Tχ ), we can find the optimal value τ opt = ξ(F T )• Tχ .The bottom row in Table I shows different values for the multiplication factor ξ for different β values.We are unable to find a maximum for the re-scaled Fisher information when β = 1.

D. Theoretical limits
The fundamental performance limit for decoherence timescale learning is given by the CRLB for the Fisher information in Eq. (17).Simulations for the performance of our proposed algorithm for perfect single-shot readout are presented in Fig. 2(b), comparing three strategies to learn T * 2 : random choice of τ (orange curve), choice of τ obtained maximizing the Fisher information F in Eq. ( 17) (i.e.minimizing the lower bound on the variance, blue curve) and choice of τ obtained maximizing the rescaled Fisher information F T in Eq. ( 9) (i.e.minimizing the 'sensitivity', red curve).The solid black line corresponds to the theoretical limit set by the Cramér-Rao lower bound (CRLB) for the Fisher information expression in Eq. ( 17).The curves represent the performance averaged over 500 repetitions of each protocol.The simulations show that both adaptive approaches outperform the non-adaptive (random τ ) protocol (orange solid line), with the adaptive protocol maximizing F T (red line) performing better than the maximisation of F (blue line).Numerical simulations in Fig. 2(b) confirm that the estimator is asymptotically unbiased, as the root-meansquare-error (RMSE) becomes smaller and smaller with increasing sensing time, approaching the fundamental estimation limit set by the Cramer-Rao bound.
To disentangle the effect of sub-optimal readout from the performance of the proposed estimation algorithm, we can compute the best performance achievable for a given readout strategy.This can be computed from the Cramér-Rao bound, considering the classical Fisher information F E for the experimental likelihood given by Eq. (12): Once again, the ultimate lower bound on the variance obtainable in the limit of many identical repetitions (epochs) is proportional to the inverse F E .To find the optimal probing time, given experimental values for α and V , we, therefore, minimize 1/F E , displayed as the solid black line in Fig. 2(b) and 2(c).

III. EXPERIMENTAL IMPLEMENTATION A. Setup and estimation algorithm
We demonstrate the proposed protocols using the single electron spin associated with a nitrogen-vacancy (NV) centre in diamond [58,59], created by a laser writing method [60,61].The NV electron spin is polarised and measured optically, even at room temperature; it can be further controlled by microwave pulses.
Whilst optical measurement at room temperature is available, it is not possible to obtain binary single shot read-out, unlike at cryogenic temperatures [43].The absence of a single shot readout can be circumvented by repeating each measurement sequence R times.At the n-th iteration, the number r n of detected photons is then used to update the distribution of T χ as The likelihood is a binomial distribution (as a result of R independent and identical Bernoulli experiments).The probability of detection of a photon click in a single repetition is given by This probability depends on V and α, which can be computed from the experimental probabilities to detect a click when the qubit is in and α = p cl (|0⟩) + p cl (|1⟩) /2.However, since the probability of detecting a photon in one measurement cycle is very small and R is large, p(r n |T χ , R) can be well approximated by a Gaussian distribution [20,37] Note that as the adaptive update rule in Eq. ( 20) does not depend on the outcome but only on the current estimate of T χ from p(T χ |r 1:n ), it can be used in the absence of single shot readout.
The experimental setup is sketched in Fig. 3(a) (more details in Appendix B).A real-time microcontroller performs the update of the probability distribution of T χ via a particle filter and chooses the optimal probing time τ , in a time-scale of about 50 µs.The value of τ is then passed to an arbitrary waveform generator (AWG) which, in real-time, generates the appropriate spin manipulation pulse sequence (including laser and microwave pulses).A computational latency of 50 µs is a factor 20 faster than previous real-time quantum sensing experiments at room-temperature [24], and much shorter than even the shortest measurement timescales in our experiments (tens of milliseconds, for τ ∼ 1 µs and R = 10 4 ).

Algorithm 1 Adaptive estimation algorithm
Input: p0(x): prior probability distribution for x = Tχ; K: number of particles; N : number of epochs; aLW : Liu-West re-sampling parameter; tRS: re-sampling threshold 1: procedure AdaptiveEstimation(n, p0, N , aLW , tRS) ri ← R j=1 mj 10: {ω k } ← {1/K} 15: return Tχ 16: The estimation algorithm is detailed in Algorithm 1.We first choose a discretisation of the distribution of T χ , described by a number K of particles {x k } distributed according to the prior probability p 0 (T χ ) (uniform in our case).We initially set the weights of each particle to 1/K, and K = 100.At each iteration, we compute the mean Tχ of the distribution (line 5) and use it to set the next probing time to τ = ξ • Tχ (line 6).We then perform the selected experiment R times, detecting r photons.The next step updates the probability distribution p(T χ ) according to Bayes rule (line 11), normalizing the distribution (line 12).At this point, if the distribution features large areas with small weights (described by the condition in line 13), is re-sampled according to the Liu-West algorithm [62].We compute the variance σ 2 and we then sample new particles (line 25) from a Gaussian distribution with variance σ 2 and mean: where a LW is the Liu-West parameter, which determines how much the new sampling preserves the original {x k } and how much it reflects the properties (mean µ) of current p(T χ ).Resampling is performed in parallel to data acquisition (see Fig. 3), so that it does not add any additional delay to the overall measurement time.
In the experiments described in Sec.III B we compare the adaptive protocols with a non-adaptive alternative The RF switches select the first shaded range, with a difference in photoluminescence intensity between the two spin states, as signal (directed to counter-A) and the second range, with no difference, as background (directed to counter-B) to monitor drifts.(b) Flow-chart representing the control flow of the adaptive experiments.The control flow of the microcontroller is agnostic to the actual measurement being performed, which is activated by trigger signals.Resampling, if necessary, is performed in parallel to data acquisition (i.e. between the triggering of the AWG with counters starting and the counters stopping) so that no overhead is added to the procedure.
where the probing time τ is chosen randomly within an expected range, known a-priori.The choice of this range affects the performance of the protocols, as adaptive techniques are more effective when there is a large uncertainty in the parameter to be estimated (i.e. for high dynamic range estimation).If the parameter is already known with a good approximation, then settings can already be optimised a-priori, without the need for real-time adaptation.In the experiments presented in Sec.III B we will select the a-priori parameter range based on typical ranges of variation of decoherence timescales for NV centres in diamond at room temperature.

B. Experimental results
Fig. 4 compares the performance of adaptive and nonadaptive protocols to estimate the static dephasing T * 2 timescale more extensively.We compare three different values of readout repetitions R (R = 10 6 , R = 10 5 , R = 10 4 ).Given our experimentally measured values of p cl (|0⟩) = 0.0186 ± 0.0017, and p cl (|0⟩) = 0.0148 ± 0.0016 photons detected per readout for the two qubit states {|0⟩ , |1⟩}, the three different numbers of readout repetitions R correspond, respectively, to mean photon num-bers ⟨n⟩ ∼ 16700, ⟨n⟩ ∼ 1670 and ⟨n⟩ ∼ 167.We implement two forms of non-adaptive protocols, one in which τ is chosen at random within the given range 0 − 8 µs (based on typical experimentally observed ranges), and one in which τ is swept uniformly across the range.For a fair comparison, the total number of measurements N • R (with N being the number of epochs) was kept constant for different R, in order to keep to total measurement time fixed.All curves are averaged over 110 repetitions of the whole estimation sequence in order to obtain the mean performance and the 95% confidence interval.Fig. 4(a) highlights the comparison of measurement uncertainty versus total measurement time for R = 10 6 repetitions per experiment, starting from a flat prior corresponding to the same starting point for all three curves.The evolution of the uncertainty is visibly distinct for the three protocols.The adaptive protocol starts learning about the unknown parameter within just a few epochs, leading to a quick drop in the uncertainty compared to the non-adaptive protocols.For the random choice of delay τ the learning is intermediate in performance between the adaptive protocol and the non-adaptive sweep.The sweep protocol performs particularly badly in the first few estimation epochs, as it performs measurements with probing times much shorter than the decoherence  rate, resulting in very little information gained.Once it sweeps across probing times closer to the true decoherence rate, uncertainty is reduced much faster.The protocol with random probing times, on the other hand, samples from the whole expected range uniformly, so that it has a higher probability of retrieving significant information even in the first epochs.The adaptive protocol outperforms both non-adaptive ones, as it quickly learns the optimal probing time to retrieve the most information about the decoherence rate.Comparison between the three experimental curves highlights that, for example, choosing an uncertainty of 1 µs, the adaptive scheme is about 10 times faster than both the non-adaptive schemes.Fig. 4(b) shows the comparison for R = 10 5 , with a similar trend, and the learning achieved by the adaptive protocol is about 4 times faster than the random.The sweep protocol performs a factor of 20 times worse.A similar advantage of the adaptive protocol is also shown for R = 10 4 [see Fig. 4(c)].
In Fig. 4(d), a comparison of the performance of the adaptive scheme for different R values, shows that R = 10 4 (corresponding to ⟨n⟩ ∼ 167 photons detected on average) achieves the minimum uncertainty obtained for any measured probing time, highlighting that more frequent adaptive updates are beneficial.Decreasing R even further, in the Gaussian approximation in Eq. ( 14), does not improve performance further.Therefore we fix R = 10 4 for the experiments in the rest of the paper.We also focus only on the version of the non-adaptive protocol with a random choice of τ , as it performs better than parameter sweeping.Each measurement is performed by preparing the qubit in (|0⟩ + |1⟩) /2 and detecting the overlap with the initial state after a time 2 • τ , with a spin-flip at time τ to re-phase the state evolution cancelling the effect of static magnetic uncertainty.The parameter range is [5 × 10 −6 s, 1 × 10 −3 s] and the ground truth is 35 µs.The standard deviation is averaged over 40 repetitions of the protocol.In all plots, shaded regions correspond to the 95% confidence interval.The confidence interval is much wider for T1 estimation as the number of protocol repetitions is quite small (due to the long time required given the ms-scale delays).
Results in Fig. 4 show that the sweep protocol performs poorly, due to very little information gain at the start.Therefore in the subsequent measurements (Fig. 5 and Fig. 6) we restricted ourselves to a comparison between the random (non-adaptive) and adaptive schemes.In Fig. 5, we compare the adaptive and non-adaptive protocols for the estimation of T 1 (relaxation) and T 2 (echo decay) timescales.In Fig. 5(top) plot we find that for T 1 estimation the adaptive protocol performs ∼ 2 times better than non-adaptive random scheme.This performance is limited by experimental uncertainties, as the long delay time combined with R = 10 4 repetitions leads to greater fluctuations over time as opposed to the other estimation sequences (T 2 , T * 2 ) presented here, where the targeted decoherence timescale in µs regime.In Fig. 5(bottom) plot we see that the adaptive protocol outperforms the nonadaptive by significant margin.Such a remarkable gain is achieved as our knowledge about T 2 for an NV centre in buld diamond is typically more uncertain apriori than the value of T 1 .In both cases, the adaptive strategy greatly outperforms the non-adaptive protocol.
We verify the performance of multi-parameter estimation [see Eq. ( 8)] using Ramsey measurement (T χ = T * 2 ).The choice of optimal sensing time τ is based on using the piecewise linear approximation, discussed in Appendix B. In Fig. 6, we compare the adaptive protocol against random measurement scheme, plotting the estimation error of both T * 2 and β versus total probing time, demonstrating a clear advantage for the adaptive algorithm.The gain observed in the T * 2 plot is ∼ 10 times for adaptive estimation over the random non-adaptive scheme.While β estimation plot shows a gain of ∼ 20 times for the adaptive scheme.For this experiment the number of particles (K) over which the Bayesian update is performed was increased to 500 to better sample the two parameter space.Fig. 7 compares adaptive estimation of the dephasing time T * 2 when maximizing the BIM (related to F ) or the BIM rescaled by the probing time (related to F T ).As discussed, this results in different multiplication factor ξ for the adaptive choice of τ rule as given in Table I.Fig. 7 shows that maximisation of F T is a factor ∼2 faster than maximisation of F , as a function of the total probing time, as expected from Sec. II C. It is worth noting that both protocols shown in Fig. 7 eventually reach the same scaling as the Fisher Information, but with different offsets.This experimental result is consistent with the theoretical prediction from Fig. 2. Similarly, the experimental results shown earlier display strong correspondence with theoretical predictions, providing a robust base for the improved performance of adaptive protocols.
Generally, while all experimental datasets reproduce the predicted trends, they do not appear to reach the theoretical limits.We believe this is due to residual experimental imperfections or fluctuations not accounted for by the likelihood function we use in the Bayesian inference process, which could be incorporated through data-driven or "gray-box" approaches [63].Experimental comparison between maximizing F and FT .An adaptive choice computed maximizing the Fisher information rescaled by the total probing time [see Eq. ( 9)] achieves a smaller uncertainty for a given probing time than simply maximizing the Fisher information.The experimental estimations of T * 2 (β = 2) are performed with R = 10 4 , with parameter range T * 2 ∈ [1 × 10 −7 s, 8 × 10 −6 s], ground truth 2.5 µs and averaged over 35 repetitions of the whole protocol.The dashed black line corresponds to the theoretical limit given by FE in Eq. (10).Both experimental curves match the scaling dictated by the CRLB for FE.

IV. DISCUSSION AND OUTLOOK
Using an adaptive Bayesian approach based on highspeed electronics, we have compared approaches to estimating decoherence timescales and decay exponents, showing that real-time adaptation of the probing time τ significantly outperforms non-adaptive approaches.We have also investigated different adaptation heuristics, demonstrating an advantage in sensitivity optimisation in situations where time is the resource to be minimised.
Real-time adaptation of experimental settings is most valuable when there is a large uncertainty in the parameter to be estimated.In this case, pre-optimising the experiment offline is not trivial, as the optimal settings typically depend on the unknown value of the parameter we are measuring.On the other hand, real-time adaptation is less important when the parameter to be estimated varies on a very small range, as in this case the experimental settings can be optimised offline based on the strong available prior information.In this paper, we have chosen reasonable ranges for the decoherence timescales and exponents we estimate, based on our own experience in the field and the variability observed in the literature.
The techniques demonstrated here will directly improve the performance of quantum relaxometry, where measuring the decoherence of a qubit is used to extract useful information from noise in the environment.An example is the use of spins in nano-diamonds to measure the concentration of radicals, exploiting the fact that they induce magnetic noise that shortens the sensor spin T 1 timescale.This has been applied, for example, to probe specific chemical reactions [64] or the concentration of radical oxygen species inside living cells [13,65].The time required to take a measurement limits the bandwidth over which fast signal variations can be tracked, a very important parameter when dealing for example with biological processes in living cells.
Our technique can be readily extended to more complex experiments and pulse sequences, for example, double electron-electron resonance (DEER) or electronnuclear double resonance (ENDOR) where measurements on electronic spins are used to infer the dynamics and decoherence of other electronic or nuclear spin.DEER experiments on single NV centres close to the diamond surface, for example, reveal changes in decay exponent that elucidate the physics of dark spins on the diamond surface [66].
While our experiments are performed using a single electronic spin associated with an NV centre in diamond, our approach and our analysis are very general and can be readily applied to any qubit system and to different applications.For example, our approach could be used for fast characterisation of quantum memories, or multiqubit quantum processors.
In summary, we have shown that real-time adaptive estimation of decoherence, through a simple analytical optimisation heuristic, leads to a considerable speed-up compared to non-adaptive schemes.Extending this approach to more complex settings, such as towards the detection of multiple individual nuclear spins and nanoscale magnetic resonance [67][68][69][70] , holds the promise to deliver even greater advantages.

APPENDIX A: DERIVATION OF ADAPTIVE RULES
As discussed in Sec.II A, we optimize the settings for the next measurement by considering an approximation of the Bayesian information matrix (BIM) (a 1×1 matrix, in the case of a single parameter) [56] which links to the classical Fisher information [57].
For a binary outcome m = 0 or 1, the classical Fisher information for the decoherence timescale T χ is where p(m|T χ ) is the probability to detect outcome m given a decoherence timescale T χ described by Eq. ( 4).The classical Fisher information F is then given by: The Fisher information F is inversely proportional to the Cramér-Rao lower bound (CRLB) of T χ , which represents the minimum reachable variance for any (unbiased) estimator of T χ .Thus, maximizing F with respect to τ is expected to improve our estimate of T χ .The value of τ that maximizes Eq. ( 17) i) needs to be approximated numerically and more importantly ii) depends on T χ which is unknown [this dependency is not explicit in the notation F (τ ) but visible in Eq. ( 17)].In this context, where prior information is known about T χ from previous measurements [see Eq. ( 3)], the Bayesian CRLB [56] is a more adapted criterion as it averages the Fisher information over the high-density region of the prior distribution of T χ .In the single-parameter setting considered here, the BIM is given (see [56]) by where E Tχ [•] denotes the statistical expectation with respect to P (T χ |m 1:(n−1) ).The Bayesian CRLB (BCRLB) is the inverse of F B (τ ) and bounds the performance [mean squared error (MSE)] of estimators of T χ [56].Instead of maximizing the information gain or minimizing the MSE, we use the BIM as a proxy reward function to be maximised when adapting τ .First, it is worth noting that the second term in the expectation in Eq. ( 18) does not depend on τ .Thus, it it sufficient to maximize E Tχ [F (τ )].
To overcome this intractable integral, we use where Tχ is a point estimate of T χ .This corresponds to replacing, in the integral E Tχ [F (τ )], the distribution P (T χ |m 1:(n−1) ) by a Dirac Delta function centered at Tχ .
If the prior distribution P (T χ |m 1:(n−1) ) is concentrated around Tχ , one option is to choose Tχ as the most likely value of T χ [15].Here we used instead the prior mean as it can be more robust when P (T χ |m 1:(n−1) ) has a large variance and/or is skewed.Note that here we do not specify how P (T χ |m 1:(n−1) ) and its moments are computed.This is detailed in Sec.III A.
As discussed in Sec.II A, we find an approximate maximum of Eq. ( 19) as , where ξ is a parameter that depends on β.

APPENDIX B: EXPERIMENTAL SETUP
Sample.The sample consists of an electronic-grade CVD diamond plate (Element Six).Isolated NV centres with good spin coherence are created by the laserwriting of vacancies at an array of sites in the diamond followed by thermal annealing to form NV centres with background nitrogen impurities.
Optics.The system is based on a custom-made confocal microscope.The sample, mounted on a custom PCB on a Newport stage, is imaged through an oil immersion objective (Olympus RMS100X-PFO).Photoluminescence is excited by a CW 532nm diode-pumped solid state laser (CNI MLL-U-532-200mW), pulsed by an acousto-optic modulator (Isomet 1250C-829A) in doublepass configuration.The laser beam is scanned across the sample by a galvo mirror pair (Thorlabs GVS012/M), using a 4f optical system.The photoluminescence is collected through a dichroic mirror (Semrock LP02-633RU-25) and detected by a fibre-coupled single photon avalanche photodiode (Excelitas SPCM-AQRH-15-FC).The optical setup is mounted inside a styrofoam box and temperature stabilised by a Peltier element controlled by a PID loop (Meerstetter TEC-1092), down to 10 mK rms.
Magnetic field.A magnetic field of 460 gauss, aligned with the NV axis, is applied by a permanent SmCo magnet mounted on a 3D motorised translation stage (Standa 8MT173-30DCE2-XYZ).The magnetic field is selected at the NV centre excited-state level anticrossing to polarize the 14 N nuclear spin, so that no hyperfine lines are present in the Ramsey experiments.The magnet is encased in a 3D-printed holder, and thermally stabilised by a Peltier element controlled by a PID loop (Meerstetter TEC-1091), down to 10 mK rms.
Microwaves.Microwaves are generated by singlesideband modulation of a tone at 1.53 GHz, from a vector source (RohdeSchwarz SMBV100A), with a lowfrequency signal (40 MHz) from a Zurich Instruments HDAWG4 arbitrary waveform generator (2.4 GSa/s, 16 bits vertical resolution, 750 MHz signal bandwidth).The modulated microwave signal is amplified (Amplifier Research 15S1G6, 0.7-6 GHz, 15 W) and directed to a wire across the sample, through a circulator (Aaren Technology 11B-GX017027G-AF), to create an RF magnetic field to drive the (m s = 0 ⇐⇒ m s = −1) electron spin resonance.
Adaptive electronics.The real-time experiment optimisation is handled by a hard real-time micro-controller (AdWin Pro II, Jäger Computergesteuerte Messtechnik).Electrical pulses from the photon detectors are sent to a pair of switches (Minicircuit ZASW-2-50DRA+), that redirect the window A ([t 0 , t 0 + 270 ns]) to counter-1 in the AdWin and the window B ([t 0 + 4 µs, t 0 + 4 µs + 400ns]) to counter-2 (t 0 is the start time of the excitation laser pulse).Counts in window-A correspond to photon counts that discriminate the spin signal [see Fig. 3(a) inset], since m s = 0 and m s = ±1 exhibit different photoluminescence intensities.Counts in window-B are not spin-dependent and are used to monitor the system's stability.
The microcontroller reads the values in the two counters and uses the signal counts in Window-A to perform a Bayesian update and to select the optimal probing time τ for the next measurement.The value f τ is passed, as an 8-bit integer to the HDAWG through its digital IO port.For each batch of R measurements, the HDAWG reads the value from the DIO port and constructs the next pulse sequence accordingly.The HDAWG generates control pulses for the acousto-optic modulator, singlesideband modulation of the microwave signal and to control the switch to re-direct electrical pulses from the detector to the microcontroller counters.
Details of micro-controller (MCU) operation flow.The operation of the MCU is described by the flowchart in Fig. 3(b).At each iteration, the MCU triggers the AWG, which delivers all the control pulse sequences, and starts the counters.If resampling is needed, it is started at this point and carried out in parallel to data acquisition so it does not add any additional overhead.Once the MCU receives a trigger back from the AWG, signalling that all pulse sequences have been executed, it stops and reads the counters, which have accumulated the total count rate for all R repetitions.The MCU then performs the Bayesian update of p(T χ ) and selects the next τ , sending the value to the AWG to compile the next qubit control sequence.
The timing of the experiment was benchmarked with respect to the internal clock of the microcontroller, finding T µs ∼ 0.255 • n, where n is the number of particles discretising T χ , i.e. about 50 µs for n = 200 particles.
Details of multiparameter approximate adaptive rule.As we could not find an analytical solution to create an update rule, instead we evaluate the equation numerically across our region of interest and fit an approximation.The results are shown in Fig. 8.We numerically evaluate Eq. ( 7), modified to include terms to quantify measurement readout occurring over R repetitions in the absence of single shot readout.We aim to find the simplest approximation to replicate its result.Observation of the steep jump in value just below τ 0 = 1 suggested the use of a piecewise polynomial fit in τ 0 and β.Using the polyfit function of the numpy package we find optimal parameters for linear fits in two separate regions.In between these regions, a heuristic choice is made to pick a measurement time 0.7 × T * 2 .Simulation of experiments using this piecewise linear approximation rule showed similar performance to the use of the full numerical approximation. .The accuracy of the numerical approximation (left) is limited by the 64-bit resolution of the floating point used in the calculation, resulting in a region in the upper right of the numerical approximation plot that was manually set to zero.The piece-wise linear approximation given by Eq. ( 8) is plotted in the centre, and the absolute difference between these two approximations is given in the right-hand side plot.

Figure 1 .
Figure 1.Real-time adaptive feedback.(a)Schematic of the real-time adaptive protocol demonstrated in this work, using the electronic spin associated with a nitrogen-vacancy (NV) centre in diamond.An Arbitrary Waveform Generator (AWG) is used to generate pulses for manipulation of the spin qubit.The spin state is then optically measured and the detected photon count rate is used by the microcontroller to estimate the value of the decoherence timescale Tχ (and the decay exponent β) using Bayesian inference.The microcontroller also computes the optimal probing time τ for the next measurement, passing the value to the AWG, which builds the next estimation sequence accordingly.(b) In our experiment, a single interrogation of the qubit does not provide sufficient information to discriminate qubit state.Hence R measurements are performed, and the resulting number r of detected photons are used to update p(Tχ) through Bayes' rule and to compute τ for the next experiment.(c) Example of an experimental adaptive estimation sequence, with p(T *2 ) shown for each measurement epoch (here the decay exponent is known a-priori as β = 2).The probability p(T * 2 ), initially uniform, converges towards a narrow peak as more measurement outcomes are accumulated.The experimental adaptively-optimised values for τ are shown on the bottom plot.

Figure 2 .
Figure 2. Numerical simulations for T * 2 estimation.(a) Fisher information F as a function of the ground truth value for Tχ = T * 2 and the probing time τ , from Eq. (17).The plot shows that the maximum of F (τ ) has a linear dependence on T * 2 , following Eq.(17) (shown as the red dashed line).(b) Comparison between different strategies to learn T * 2, assuming that single shot qubit readout is available: random choice of τ (orange), choice of τ obtained maximizing the Fisher information F in Eq. (17) (blue) and choice of τ obtained maximizing the rescaled Fisher information FT in Eq. (9) (red).The x-axis shows the total probing time (cumulative over epochs), while the y-axis shows the uncertainty, defined as the root mean squared error (RMSE) from the ground truth value.The solid black line corresponds to the theoretical limit set by the Cramér-Rao lower bound (CRLB) for the Fisher information in Eq. (17).Adaptive protocols outperform the non-adaptive (random) protocol.(c) Simulations for the same protocols as in (b), for the case where single shot readout is not available (using experimental values p cl (|0⟩) = 0.0187 and p cl (|1⟩) = 0.0148, and R = 50000).The solid black line corresponds to the limit set by the CRLB for the Fisher information in Eq. (17).The dashed black line shows the CRLB limit for the Fisher information FE in Eq. (10).

Figure 4 . 2 .
Figure 4. Experimental comparison between the adaptive and non-adaptive estimation of the dephasing time T * 2 .The dephasing time T * 2 is measured by a Ramsey experiment.The spin state, initially polarised by a laser pulse in |0⟩, is prepared in an initial superposition state (|0⟩ + |1⟩)/2 with a microwave π/2 pulse.A second π/2 pulse, after a delay τ , converts the phase information into the population of the |0⟩ state, which is then readout optically.Our adaptive protocol (blue line) outperforms non-adaptive protocols, either sweeping τ on a pre-determined range (green line), or randomly picking τ (orange line), for a different number of readout repetitions: R = 10 6 (a), R = 10 5 (b) and R = 10 4 (c).In each sub-plot, we show the root-mean-squared error (RMSE), averaged over 110 experimental protocol repetitions, as a function of the total probing time in seconds.(d) Comparison of the performance of the adaptive protocol for different values of R, with the same experimental data as in (a), (b), and (c) plotted together for ease of comparison.For all plots the parameter range is T f 2 * ∈ [1 × 10 −7 s, 8 × 10 −6 s], the ground truth is 2.5 µs and the shaded areas correspond to 95% confidence intervals.

Figure 5 .
Figure 5. Experimental comparison of adaptive and non-adaptive estimation for T1 and T2.In both cases the adaptive choice of τ is performed maximizing F , and the number of readout repetitions is R = 10 4 .(a) Standard deviation as a function of probing time in the estimation of spin relaxation T1 (β = 1).Each measurement is performed by initializing the electron spin in |0⟩ and detecting the probability of being in |0⟩ after a delay τ .The parameter range is [1 × 10 −4 s, 8 × 10 −3 s] and the ground truth is 1.45 ms.The standard deviation is averaged over 10 repetitions of the protocol.(b) Standard deviation as a function of probing time in the estimation of the Hahn echo decay time T2 (β = 3/2).Each measurement is performed by preparing the qubit in (|0⟩ + |1⟩) /2 and detecting the overlap with the initial state after a time 2 • τ , with a spin-flip at time τ to re-phase the state evolution cancelling the effect of static magnetic uncertainty.The parameter range is [5 × 10 −6 s, 1 × 10 −3 s] and the ground truth is 35 µs.The standard deviation is averaged over 40 repetitions of the protocol.In all plots, shaded regions correspond to the 95% confidence interval.The confidence interval is much wider for T1 estimation as the number of protocol repetitions is quite small (due to the long time required given the ms-scale delays).

Figure 8 .
Figure 8. Piece-wise linear approximation accuracy.Plots demonstrate the optimal choice of next measurement time (τ1) vs the current estimate of coherence decay exponent β and previous measurement time τ0.Both measurement times are normalised to the current estimate of T *2 .The accuracy of the numerical approximation (left) is limited by the 64-bit resolution of the floating point used in the calculation, resulting in a region in the upper right of the numerical approximation plot that was manually set to zero.The piece-wise linear approximation given by Eq. (8) is plotted in the centre, and the absolute difference between these two approximations is given in the right-hand side plot.

Table I .
Choice (20)ptimal probing time.The optimal probing time τopt is chosen according to Eq.(20), with factor ξ. Numerically computed values of ξ for different decay exponents β are shown in the table, computed by maximizing either the classical Fisher information F [see Eq. (17)] or the rescaled Fisher information FT [see Eq. (9)].No maximum of FT can be found for exponential decay, when β = 1.