Resource-efficient adaptive Bayesian tracking of magnetic fields with a quantum sensor
RResource-efficient adaptive Bayesian tracking ofmagnetic fields with a quantum sensor
K. Craigie, E. M. Gauger, Y. Altmann and C. Bonato
School of Engineering and Physical Sciences, SUPA, Heriot-Watt University,Edinburgh, EH14 4AS, UKE-mail: [email protected], [email protected]
Abstract.
By addressing single electron spins through Ramsey experiments,nitrogen-vacancy centres can act as high-resolution sensors of magnetic field. Inapplications where the magnetic field may be changing rapidly, total sensing timeis crucial and must be minimised. Bayesian estimation and adaptive experimentoptimisation protocols work by computing the probability distribution of the magneticeld based on measurement outcomes and, by computing acquisition settings for thenext measurement. These protocols can speed up the sensing process by reducingthe number of measurements required. However, the computations feeding into thenext iteration measurement settings must be performed quickly enough to allow real-time updates. This paper addresses the issue of computational speed by implementingan approximated Bayesian estimation technique, where probability distributions areapproximated by a superposition of Gaussian functions. Given that only threeparameters are required to fully describe a Gaussian, we find that the magneticfield probability distribution can typically be described by fewer than ten numbers,achieving a reduction in the number of operations by factor 20 compared to existingapproaches, allowing for faster processing.
Keywords : Nitrogen-vacancy centre, Quantum metrology, Magnetic field tracking,Bayesian filtering
1. Introduction
Control and measurement of individual electron spins, achieved in the last two decades,enables highly sensitive measurements of magnetic fields [1–3], with spatial resolutionon the order of tens of nanometres [4]. This is typically achieved through a defect indiamond, the nitrogen-vacancy (NV) centre, which allows spin readout even at roomtemperature through Optically Detected Magnetic Resonance (ODMR) and long spincoherence time [5, 6].NV centres can also be used as nanoscale sensors for temperature [7], strain [8,9] andelectric eld [10, 11]. Their sensing capabilities have been applied to studying nanoscalemagnetic phenomena in materials and biological processes [4,12]. For example, biologicalcompatibility for nanodiamonds containing NV centres allows monitoring of nanoscale, a r X i v : . [ qu a n t - ph ] A ug esource-efficient Bayesian protocol for quantum magnetic field tracking
2. Background
A magnetic field applied to the electron spin induces a Zeeman splitting of the energylevels, which can be measured by a Ramsey experiment [28]. In a Ramsey experiment,an equal spin superposition freely evolves under the applied magnetic field B , so thatthe spin eigenvalues acquire a relative phase, corresponding to a rotation at the Larmorfrequency. The probability for outcome µ ∈ {
0; 1 } given a Larmor frequency f B (corre-sponding to a magnetic field B = f B /γ , where γ is the gyromagnetic ratio, is P ( µ | f B , θ, τ ) = 1 + e iµπ cos (2 πτ f B + θ )2 . (1)where τ is the sensing time and θ the rotation angle of the measurement basis (which iscontrolled by the phase of the second π/ esource-efficient Bayesian protocol for quantum magnetic field tracking P ( x | y, z ) denotes the distribution of x , conditioned on the value of ( y, z ),ie, the distribution of x given y and z .In this paper, we consider the problem of tracking a magnetic field (through f B )that fluctuates on timescales longer than a single measurement time, but potentiallyshorter than a large number of repetitions of the measurement time. We assume thefluctuations to be described by a Wiener process yielding a sequence of Larmor frequencyvalues that are generated as f ( t + δt ) B = f ( t ) B + κ d W ( t ) . (2)The diffusion coefficient, κ , is a measure of the magnetic field rate of change andd W ( t ) is an innitesimal Wiener increment during a time inverval δt . This process issimulated by discretising the time axis to intervals of length τ min (minimum sensingtime), and generating a normal distribution with variance τ min . In this way, a changingLarmor frequency signal is generated that will act as the ground truth in trackingsimulations. More generally for any small time interval δt , equation 2, leads to thefollowing Gaussian random walk distribution P (cid:18) f ( t + δt ) B (cid:12)(cid:12)(cid:12)(cid:12) f tB , κ (cid:19) ∝ exp − (cid:16) f ( t + δt ) B − f ( t ) B (cid:17) δtκ , (3)which will be used in the tracking algorithm described in the next section. Bayesian online tracking of f B consists of approximating in real time of the probabilitydistribution f B , which can in turn allow us to optimise the experimental settingsthrough θ and τ in equation 1. After the n -th measurement, the (posterior) distribution P ( f ( t n ) B | µ ( t n ) , θ n , τ n ), of f ( t n ) B is updated using Bayes’ rule P ( f ( t n ) B | µ ( t n ) , θ n , τ n ) ∝ P ( µ ( t n ) | f ( t ) B , θ n , τ n ) P ( f ( t n ) B | µ ( t n − ) , θ n , τ n ) , (4)where µ ( t n ) = (cid:110) µ ( t ) , . . . , µ ( t n ) (cid:111) . In equation 4, P ( f ( t n ) B | µ ( t n − ) , θ n , τ n ) acts as a priordistribution and can be obtained via P ( f ( t n ) B | µ ( t n − ) , θ n , τ n ) = (cid:90) P (cid:16) f ( t ) B (cid:12)(cid:12)(cid:12) f ( t n − ) B , κ ) P ( f ( t n − ) B | µ ( t n − ) , θ n , τ n )d f ( t n − ) B , (5)with P ( f ( t n − ) B | µ ( t n − ) , θ n , τ n ) obtained via the Bayesian update after the ( n − θ n , τ n ), are described in more detail in Sections3.5.2 and 3.5.1.Figure 1 illustrates the main principle of Bayesian adaptive tracking algorithms.However, as discussed in the introduction, exact Bayesian inference based on equations4 and 5 is not tractable because of the shape of the likelihood in equation 1, whichmakes the integral in equation 5 not computationally tractable efficiently. While it istheoretically possible to use particle filters [25], their computational complexity make esource-efficient Bayesian protocol for quantum magnetic field tracking Initialise PrecessionBayesian Update (equation 4)
Prediction (equation 5)
Read-out θ n τ n μ (t n ) Update probability distributionRamsey Measurement
Figure 1.
Graphical representation of adaptive protocol, showing flow of information.The measurement outcome, μ , from the Ramsey experiment is fed into the Bayesianupdate, followed by a prediction step. From this computation, the adaptive phase, θ n ,and sensing time, τ n , are determined. These values are fed into the following Ramseyexperiment. them less attractive than approximate methods using Gaussian mixture approximations,as proposed here.For completeness, a simplified example of the tracking protocol, similar to thatdescribed in [22], is provided in figure 2, where the adaptive phase and sensing time arechosen to completely eliminate one of the two peaks in panel (d).
3. Method
In this article, we aim to reduce computational time and memory requirements byapproximating P ( f ( t n ) B | µ ( t n ) , θ n , τ n ) as a finite sum of Gaussian distributions. This isperformed by approximating the likelihood P ( µ ( t n ) | f ( t n ) B , θ n , τ n ) by a weighted sum ofGaussian functions (with respect to f ( t n ) B ). Figure 2a shows how the actual cosine-shapedlikelihood function is approximated by two shifted Gaussian peaks. As is apparent, theapproximation is relatively accurate at the top of the peaks, but is rather poor atthe bottom due to the Gaussian tails. In the following, we address the dual questionof how well a magnetic field can be tracked despite this simplification, and examinethe associated advantage of doing so in terms of the algorithm’s reduction in memoryrequirements and computation time. We begin by approximating the initial likelihood in equation 1 as: P ( µ | f B , θ n , τ n ) ≈ N G (cid:88) l =0 Ae − ( fB − al )22 σ a . (6)This approximation can be computed by a Taylor expansion of equation 1 at 2 πτ f B + θ n + µπ ≈ πl and re-writing as a Taylor expansion for an exponential. esource-efficient Bayesian protocol for quantum magnetic field tracking
20 0 20 f B [ MHz ] P ( f B )
20 0 20 f B [ MHz ] P ( f B )
20 0 20 f B [ MHz ] P ( f B )
20 0 20 f B [ MHz ] P ( f B )
20 0 20 f B [ MHz ] P ( f B ) (a)(b) (c)
20 0 20 f B [ MHz ] P ( f B ) (d)
20 0 20 f B [ MHz ] P ( f B ) (e)(f) (g) Non-Gaussian P(f B )Gaussian P(f B )P( μ = 1|f B )P( μ = 0|f B ) Figure 2.
Tracking steps after a change in magnetic field. (a) likelihood function withGaussian approximation, (b) prior and likelihood function for measurement n=1, (c)posterior for n=1, μ = 1, with Gaussian approximation, (d) prior and likelihood functionfor measurement n=2, (e) posterior for n=2, μ = 1, with Gaussian approximation, (f)prior and likelihood function for measurement n=1, (g) posterior for n=3, μ = 1, withGaussian approximation. πτ n f B + θ n + µπ )2 ≈ − (2 πτ f B + θ n + µπ ) ≈ exp − ( πτ n f B + θn + µπ ) (7)For the initial likelihood and under this approximation, we let all Gaussians have thesame amplitude A = 1 (the normalisation is handled) and width σ a = 1 / ( √ πτ ),meaning they only differ in their centres a l as follows: a l = 2 πl + πµ + θ n πτ . (8) esource-efficient Bayesian protocol for quantum magnetic field tracking τ ). This implies that the required number N G of Gaussiansfor the approximation is given by 2 N + 1, since this is the maximum sensing coefficient(plus one additional one to cater for edge peaks which are only partially visible). Once the likelihood is approximated by a finite sum of Gaussians, the Bayesian updatein equation 4 becomes tractable and the posterior distribution reduced to a productof two mixtures of Gaussians (one arising from the prior distribution and one from thelikelihood function). Defining
A, a, σ a and B, b, σ b as, respectively, the amplitude, centreand standard deviation of two Gaussians, their product will be another Gaussian thatis fully characterised by the parameters C, c, σ c with σ c = (cid:118)(cid:117)(cid:117)(cid:116) σ a σ b σ a + σ b , (9) c = aσ b + bσ a σ a + σ b , (10) C = ABe ( aσ b + bσ a )2 / ( σ a + σ b ) − ( a σ b + b σ a )(2 σ aσ b ) . (11)Note that after the initial assignment, the amplitudes of different Gaussians will nolonger generally be identical, requiring calculation of C above. As described in equation 5, the (predictive) probability distribution of f ( t n ) B at time t n = t n − + δt n can be found as the convolution of the probability distribution attime t n − with a zero-mean Gaussian with variance δt n κ . The time elapsed betweenmeasurements is given by δt n = τ n + t oh , where τ n is the sensing time and t oh is theoverhead time. This overhead time is added to account for the time it takes to physicallycarry out measurements in a lab. Naturally, this value varies depending on the preciseexperimental equipment and set-up. The distributions before and after prediction aretherefore: P ( f ( t n − ) B | µ ( t n − ) , θ n , τ n ) = (cid:88) l C l e − (cid:16) f ( tn − B − cl (cid:17) σ l , (12) P ( f ( t n ) B | µ ( t n − ) , θ n , τ n ) = (cid:88) l C l σ l (cid:113) σ l + κ dt e − ( f ( tn ) B − cl ) σ l + κ δtn ) , (13) esource-efficient Bayesian protocol for quantum magnetic field tracking σ (cid:48) = √ σ + κ δt and D = Cσσ (cid:48) . It is important to note that, using mixtures of Gaussians to approximate the likelihoodand the prior distribution (say with m and n terms, respectively), the number of Gaus-sians ( m × n ) in the resulting posterior distribution keeps on increasing over time. Toprevent this, we introduce an automatic pruning step based on amplitude thresholding.Any Gaussian with amplitude smaller than the predefined threshold is discarded.Additionally, we combine Gaussians that are very similar to each other, as defined bytheir Kullback-Leibler (KL) divergence [29]: KL ( g , g ) = log σ σ + σ + ( a − a ) σ − , (14)where g and g are the two Gaussians being compared and a , a are their respectivecentral positions. When the KL divergence for a pair of Gaussians is below a giventhreshold, the pair is merged into one single Gaussian whose mean and variance areobtained by averaging the parameters of the two original Gaussians and whose amplitudeis the sum of the original amplitudes. We found that the KL divergence threshold KL th = 0 .
001 for merging and the amplitude threshold A th = 0 .
04 for pruning workwell in practice within the parameter range we have studied.The proposed tracking scheme might temporarily lose the track of f B and this yieldsall Gaussians having small amplitudes. If no amplitude exceeds the threshold, all theprevious amplitudes and positions are kept unchanged and the variances are doubled.By simply broadening the previous distribution when the tracking is lost, the protocolpicks up the signal again after a few iterations. The proposed protocol is divided into phases, namely the initial adaptive sensing,followed by the adaptive tracking phase. The sensing portion is required to obtaina starting Larmor frequency value, which can then be tracked. Experimentally, thiscould also be obtained from locating the initial Larmor frequency by observing thespin resonance signal while sweeping the spin drive frequency. During sensing, themeasurement time is not chosen adaptively, but in a predetermined sequence to narrowin on the Larmor frequency. Repeating each sensing time is required to minimise errorsin sensing (see line 4 of Algorithm 1). However, the goal of adaptive sensing is to usethe fewest Ramsey measurements and so a balance must be achieved. The number ofRamsey measurements performed at each sensing time is determined by integers G andF, which have been selected in accordance with previous work on adaptive sensing [24].Due to this previous work studying in detail the initial sensing process, as well as theGaussian approximation implementation providing an advantage only for tracking, thisphase is not of primary interest for the current work. Algorithm 1 gives an overview esource-efficient Bayesian protocol for quantum magnetic field tracking f B typically evolvesfrom a uniform distribution (prior to any measurements) to an almost unimodal densitywhich can be well approximated by a single Gaussian distribution. Algorithm 1 : Adaptive sensing overview. Variables: adaptive phase ( θ n ),measurement outcome ( µ ( t n ) ), sensing time ( τ n ), minimum sensing time ( τ min ), numberof sensing times (N). The remaining variables, M n , G and F, define the number ofrepetitions for each sensing time. for n= 0 to N-1 do τ n = 2 n τ min choose θ n M n = G + F ( n − for m = 1 to M n : do µ ( t n ) = Ramsey ( θ = θ n , τ = τ n ) Bayesian update ( µ = µ ( t n ) , θ = θ ,n τ = τ n ) end for end for In contrast to the first phase, during the second phase of the protocol, the sensingtime is chosen before each Ramsey measurement, based on the level of uncertainty ofthe current probability distribution of f B . Algorithm 2 illustrates the sequence of stepsinvolved in adaptive tracking. In both tracking and sensing, the Gaussian approximationonly affects steps ”choose θ n ,” ”Bayesian update” and ”Prediction”. The other stepsdo not use P ( f ( t ) B | µ ( t ) , θ n , τ n ) and thus remain unaffected. Algorithm 2 : Adaptive tracking overview. Variables: adaptive phase ( θ n ),measurement outcome ( µ ( t n ) ), sensing time ( τ n ). for total tracking time interval do choose τ n choose θ n µ ( t n ) = Ramsey ( θ = θ n , τ = τ n ) Bayesian update ( µ = µ ( t n ) , θ = θ n , τ = τ n ) P rediction ( κ ) end for In this work, to determine whether the optimum sensingtime has been chosen, we use the figure of merit proposed in [22] and computed froman estimate of the standard deviation of the posterior distribution of f B (as detailed inequation 21 of [22]). This figure of merit is associated with a threshold above which thesensing time is judged insufficiently accurate and reduced by a factor two for the nextmeasurement. This procedure is repeated until the threshold condition is met. If the esource-efficient Bayesian protocol for quantum magnetic field tracking In order to maximise the information from each measurement,we adaptively set the angle of the measurement basis as [20]: θ n = 12 arg { p t n } , (15)where p t n is the prior probability distribution in Fourier space, p k for k = 2 t n , where t n is the sensing time coefficient. In terms of Gaussians, p t n is computed as: p t n = (cid:88) l √ πC l σ l e − π (2 t n ) τ σ l + l π (2 t n ) τc l (16)
4. Results
To test the performance of the Gaussian-approximation tracking protocol, we performednumerical simulations. We assume perfect spin readout fidelity and T ∗ = 100 µ s.Other constant parameters were the amplitude pruning threshold A th = 0 .
04, mergingthreshold KL th = 0 .
001 as well as G = 5 and F = 3. The time interval over which thesignal was tracked was set to 50 ms for direct comparisons and varied during statisticalcomparisons, so that each run would contain 1000 Ramsey measurements. Overheadtime, the extra experimental time taken to perform a Ramsey measurement, and κ , theprediction coefficient, were varied to examine the robustness of the protocol.The performance is assessed with the mean squared error (MSE) (cid:15) = 1 T (cid:90) T | f b − f estb | , (17)where f est is the estimated frequency and f B the true frequency. Figure 3 illustrates asuccessful tracking run, characterised by a typical MSE value of around 0.08 MHz/ms.In addition, to assess the improvement in terms of computational cost, we comparedthe average number of parameters used in the discretisation of the distribution of f B .For the non-approximate implementation, this is the number of equally-space frequencypoints in the discretisation. For the Gaussian case, this is simply the number of Gaussianparameters i.e. 3 × the number of Gaussians employed for approximating the probabilitydistribution. As suggested by figure 2, without factoring computation time, the reliability of thetracking is generally degraded when using Gaussian tracking. In figure 4, panel a, 50 msof the exact same magnetic field fluctuations were tracked 200 times using both methods.Neither are perfect but it is clear that the first histogram bin contains slightly fewerGaussian runs (88.5%) than non-Gaussian runs (93%). Any run in this first bin hassuccessfully tracked f B without losing track of the Larmor frequency in a major way. esource-efficient Bayesian protocol for quantum magnetic field tracking f B [ M H z ] f B [ M H z ] Estimated frequencyTrue frequency
Time [ ms ] f B [ M H z ] P ( f B ) Estimated frequencyTrue frequency
Time [ ms ] Figure 3. (a) example Gaussian tracking run, with an MSE of 0.078 MHz/ms, κ = 10MHz H z / , t oh = 10 µs . Bottom subplot shows the difference between the estimationand ground truth, green colour indicates the probability distribution calculated by theBayesian update We chose to vary overhead time and prediction coefficient, κ , as they had been usedto benchmark previous comparisons of adaptive versus non-adaptive schemes [22]. Inthis way, we could test the robustness of the methods, as these are also not parameterswe have control over in an experiment. From subplots 4b and c, it can be seen thatthe mean squared error tends to increase as these variables increase. We also see thatthe Gaussian method is overall more likely to break down at larger values of κ andoverhead time, than the smaller values. This is because, at larger values of both thesevariables, we can fit in fewer measurements per change of magnetic field. We have lesstime to narrow in on the correct frequency, and this puts the Gaussian method at a esource-efficient Bayesian protocol for quantum magnetic field tracking MSE [
MHz / ms ] N u m b e r o f r un s Non-GaussianGaussian A v g . N u m b e r s
10 20 30 40 50 [ MHz Hz ] M S E [ M H z / m s ] (a) (b)(c) A v g . N u m b e r s Overhead time [ s ] M S E [ M H z / m s ] Figure 4. (a) direct comparison of 200 tracking runs on the same set of ’true’ Larmorfrequency values, for both Gaussian approximated and non-approximate algorithms. Asingle run is described in figure 3 (b) sweeping prediction coefficient, κ , and performinga statistical comparison between the two methods. Each data point is 1000 trackingruns, overhead time fixed at 10 µs (c) similarly sweeping overhead time, with κ set to10 MHz H z / . Again, each data point is 1000 runs. disadvantage, since we cannot factor in the potential computational time reduction inthe simulation.The average numbers is a reflection of the Bayesian update computation time.In the case of the original method, this is mostly in the thousands. Note that thesenumbers are not constant as the method discretises of the Larmor frequency probabilitydistribution with a grid whose resolution is proportional to the number of sensing times N . Moreover, N varies depending on κ and t oh since the maximum sensing time isoptimised according to equation 14 of [22]. For the Gaussian method though, the averagenumber of parameters is around eight or nine for this length of tracking run, suggestinga substantial reduction in computational complexity and thus a gain in computationaltime. The average number of parameters is dragged up by the initial sensing, whichinitially requires hundreds of Gaussian parameters to describe it. The tracking usesmostly three, sometimes six, and rarely nine or more parameters, indicating that asingle Gaussian peak, with only occasionally a second or third, is largely sufficientand delivers adequate performance. This translated to a reduction in the number ofoperations required to track a changing magnetic field by a factor of 20. Therefore, onewould expect this Gaussian-approximated computation to be 20 times faster than the esource-efficient Bayesian protocol for quantum magnetic field tracking
5. Conclusion
In this paper, we have simulated and tracked a fluctuating magnetic field via Ramseyexperiments on an NV centre. Ramsey measurements were optimised through anadaptive Bayesian update protocol, with Gaussian approximation of all probabilitydistributions. A comparison of the Gaussian-approximated with the original protocolrevealed that the approximation yielded potentially significant increases in computationspeed, with the number of operations involved in calculating the Bayesian update andadaptive measurement inputs decreased by factor 20. As expected, the approximationperformed slightly poorer at tracking than the non-approximated protocol. However,under the simulation parameters detailed here, the fail rate of the tracking increasedonly by around 5%.This protocol could find applications in sensing settings where one needs to tracka fluctuating signal with a statistics that is, at least approximately, known in advance.For example, one could use an NV centre in a nanodiamond to monitor temperature driftinside a living cell [15,16]. Measuring temperature is an integral part of studying energymetabolism [30] or developmental processes [31,32]. One other issue with nanodiamonds,is that while moving in a fluid medium, they rotate considerably, often very rapidly. Ourprotocol could be extended to track this rotation with minimal resource consumption.Another possible application is in experiments with levitated nano-diamond, whereour technique could be used for fast tracking of rotation. In these experiments, thenanodiamond containing NV centres is held in place translationally using ion traps[33,34], optical traps [35] or magnetic traps [36]. The librational (rotational) frequenciesof trapped nanodiamond vary from 100s [37] of Hz to 1GHz [38]. For the lowerfrequencies, in which tracking resolutions of 1ms per data point are suitable, our methodcould be directly applied. However, the method could also be used in conjunction withac magnetomety [39, 40], using spin-echo instead of Ramsey measurements, to achievetracking of faster, periodic librations.Tracking provides the information required for realigning the nanodiamondorientation, for whichever feedback mechanism the traps use. This feedback mechanismcould potentially be 3D Helmholtz coils in the case of ion traps or optical traps. Fora magnetic trap, it has been proposed [41] that the diamond orientation be confinedwith an electrode using the dielectric force on the non-spherical diamond. Though thismethod does not conventionally require feedback, tracking may nonetheless prove usefulin testing the confinement.Finally, the Gaussian-approximation described here could be applied to tracka quantum signal, such as the magnetic field arising from a bath of nuclear spinssurrounding a central electron spin. Previous theoretical work has shown that, byadaptively tracking the fluctuating nuclear magnetic field and narrowing its distributionthrough the back-action of the quantum measurement process, one can considerably esource-efficient Bayesian protocol for quantum magnetic field tracking
Acknowledgments
The authors would like to thank Eleanor Scerri for helpful discussions. This projectis supported by the Engineering and Physical Sciences Research council throughgrant EP/S000550/1, and by a Weizmann-UK Joint Research Program grant. K.C.acknowledges studentship funding from EPSRC under grant no. EP/L015110/1. Y.A. issupported by the Royal Academy of Engineering under the Research Fellowship schemeRF201617/16/31.
References [1] Degen C L, Reinhard F and Cappellaro P 2017 Rev. Mod. Phys.
19, 7 esource-efficient Bayesian protocol for quantum magnetic field tracking [19] Dinani H T, Berry D W, Gonzalez R, Maze J R and Bonato C 2019 Phys. Rev. B Preprint )[22] Bonato C and Berry D W 2017 Phys Rev A ISSN 20452322[27] Yan F F, Yi A L, Wang J F, Li Q, Yu P, Zhang J X, Gali A, Wang Y, Xu J S, Ou X, Li C F andGuo G C 2020 npj Quantum Information ISSN 080802[41] Pedernales J S, Morley G W and Plenio M B 2019 Motional dynamical decoupling for matter-waveinterferometry (
Preprint arXiv:1906.00835 )[42] Scerri E, Gauger E M and Bonato C 2020 New J.Phys.3,22