Bayesian-Based Hybrid Method for Rapid Optimization of NV Center Sensors

NV centers are among the most promising platforms in the field of quantum sensing. Magnetometry based on NV centers, especially, has achieved concrete development in areas of biomedicine and medical diagnostics. Improving the sensitivity of NV center sensors under wide inhomogeneous broadening and fieldamplitude drift is a crucial issue of continuous concern that relies on the coherent control of NV centers with high average fidelity. Quantum optimal control (QOC) methods provide access to this target; nevertheless, the high time consumption of current methods due to the large number of needful sample points as well as the complexity of the parameter space has hindered their usability. In this paper, we propose the Bayesian estimation phase-modulated (B-PM) method to tackle this problem. In the case of the state transforming of an NV center ensemble, the B-PM method reduced the time consumption by more than 90% compared with the conventional standard Fourier basis (SFB) method while increasing the average fidelity from 0.894 to 0.905. In the AC magnetometry scenario, the optimized control pulse obtained with the B-PM method achieved an eight-fold extension of coherence time T2 compared with the rectangular π pulse. Similar application can be made in other sensing situations. As a general algorithm, the B-PM method can be further extended to the open- and closed-loop optimization of complex systems based on a variety of quantum platforms.


I. INTRODUCTION
Nitrogen-vacancy (NV) center in diamond shows its bright prospect in quantum sensing of magnetic field [1][2][3], electric field [4], temperature [5,6] and strain [7].Among these areas, researches on NV center based ultrasensitive magnetometry have achieved fast development [8,9] and are on the road to practical and commercial applications in biomedicine and diagnostics [10][11][12][13].The inhomogeneous broadening due to ambient nuclear spins and external bias field is one of the main obstacle to further improving the sensitivity of the sensor.To alleviate the problem, dynamic decoupling (DD) is widely applied in sensing strategies based on NV center [14][15][16], and varies of optimal method are used to increase the performance of DD sequence [17][18][19], since the inhomogeneous broadening will also damages the fidelity of control pulses with finite pulse length that construct the DD series.Adiabatic strategies can significantly prolong the decoherence time  2 and sensitivity compared to the conventional flat  pulses [18], while smooth shaped pulse design based on numerical optimization [19] can adapt to a wider pulse length range where adiabatic condition is not satisfied.One major drawback of numerical optimal design is its time efficiency.On one hand, multiple sample points need to be measured to give a accurate value of the objective function denoting the average fidelity over the frequency and field amplitude broadening, such the processing time of calling the objective function once is prolonged.On the other hand, the total calling times of the objective function will increase along with the number of parameters, while in most cases a larger number of parameter (over 10) is needful to guarantee a fair fidelity.A rapid while well-performed * tianjiazhao@tyut.edu.cnoptimization method is anticipated to improve the usability of such control strategies.
In this work we propose the Bayesian estimation phasemodulated (B-PM) method to overcome the time consumption problem.Unlike the developed Bayesian methods [20,21], the B-PM method grafts the Bayesian estimation model onto the direct search method, such circumvent the complex process to calculate the acquisition functions [22].Further taking advantage to the phase-modulated method, the B-PM method makes itself a efficient hybrid optimization method for robust quantum control against noise.The objective function of optimization process under consideration is the average value of multiple function with different detuning and amplitude value of the control filed.The value of these multiple function can specifically refer to the fidelity between the final and target states, gate operators and so on.Using Bayesian based estimation model, we make an accurate prediction of the average fidelity based on a small number of sample points, such the computation time of calling the objective function once is reduced.In addition, adopting phase-modulated base allow the control field comprising multiple frequency components with less parameters, which lead to a significant decrease of the necessary total calling number of the objective function to find the local optimal results.We further verify the computation time of estimation process is negligible compared with the sample measuring process.Overall the B-PM method decrease the total time consumed during the entire optimization process.We firstly apply the B-PM method to the state flipping of NV center ensemble.Comparing with the conventional standard Fourier bases (SFB) method without estimation, the B-PM method increase the average fidelity from 0.894 to 0.905 with only 9.3% total time consumption.When applied to the sensing strategy of AC magnetic signal, the B-PM shaped pulse prolong the decoherence time  2 by 8 times compared to the conventional arXiv:2302.08410v1[quant-ph] 16 Feb 2023 rectangular  pulse.The B-PM method can be extended to DC sensing strategy [23], and the open-and closed-loop optimization process for open systems [24,25] and many-body systems [26][27][28].

A. Optimal control model of NV center ensemble
The nitrogen-vacancy center (NV center) in diamond [29] has an triplet ground state with electron spin  = 1 and zerofield splitting  = 2 × 2.88 GHz.The energy gap between ground and excited state is 1.945 ev (637 nm), and due to uneven radiation probability from excited state |  = ±1 and |  = 0 to the radiationless metastable intermediate state, NV center can be optically initialized and read out.The structure of NV center and energy level scheme is shown in Fig. 1 (a) and (b).The Hamiltonian of the triplet ground state of NV center can be expressed as where  is the zero-field splitting, B is the magnetic vector,  NV = 2 ×2.8 MHz/G is the gyromagnetic ratio of NV center,  ele is the electric interaction term with coupling coefficient ∼ Hz/(V/m) and  hf is the hyperfine coupling term express interactions between NV center and ambient nuclear spins.Under the magnetic field along  direction B =   and neglecting the electric coupling term, the Hamiltonian of the two-level subspace formed by the   = 0 ground state | 0 and the   = −1 (or   = 1) ground state | 1 can be written as One simple but requisite control target is to flip all spins in the NV ensemble from one state to another with high average fidelity.Consider an ensemble of two-level system described by equation 2 which is controlled by a time-dependent field (), the Hamiltonian of each spin in the ensemble can be represented as with  0 =  +   the unperturbed energy gap between | 0 and | 1 ,  the normally distributed detuning term as a result of inhomogeneous local ambient among the ensemble, and  the amplitude drift factor varies with different control trials.The probability density of  and  can be noted as and respectively, where we take the full width at half maximum (FWHM) of () as 2 × 26.5MHz, corresponding to a dephasing time of  * 2 ≈ 20ns [18], and FWHM of () as 0.5.Figure 1 (c-e) depict the fidelity as a function of detuning  and amplitude drift  of the rectangular control field () = / = 2 × 10 MHz, where  is the evolution time, (, , ) is the final state, and the initial state and the target state is taken as (, , 0) =  0 and   =  1 respectively.Based on Figure 1 (c) we take the sample range as  ∈ 2 × [−10, 10] MHz and  ∈ [0.5, 1.5], as is showed in Figure 1 (f), and the optimization target is to improve the average fidelity among this range.
Taking into account the inhomogeneous detuning as well as the amplitude drift, the explicit form of the average fidelity can be represented as For fixed ,  and , (, , ) is functional of control field (), therefore is the function of parameters  that construct ().For example, () can be constructed based on the standard Fourier basis as with  = {, , }, or on the phase modulated Fourier basis [30] as with  = {, , }.The optimization target, therefore, is to find the parameters  that gives the highest value of F while the maximum amplitude of () is limited to  max .We concisely describe this optimization model as Plenty of optimization methods can be used to solve this problem, with  being the optimization parameters and F being the objective function.Ideally, the value of F should be calculated by averaging the fidelity of a large number of samples taken from the whole ensemble with several repeated control processes.In practice, this time-consuming process is usually replace by taken the weighted average of several equidistant sampling among certain range and can be represented as wherein the distribution of  and  is assumed to be known, ,  are the number of different values of  and  respectively, is the normalization constant.Even so, this process still consume a substantial amount of time.Below we will show how the Bayesian-based estimation method can be used in the optimization process and combined with the PM method to efficiently reduce the total execution time.

B. The estimation model
One essential task of the optimal process is to build an computationally cheap estimation model, also known as the metamodel [31][32][33], of the expensive function.Bayesian statics and analysis methods provide a powerful tool to complete this task.From the viewpoint of Bayesian statistics, any event can be endowed with a probabilistic property, whether or not it is a stochastic event.So the undetermined parameter  in a static model can be taken as a realization of stochastic process  subject to a prior distribution (), which represents the available knowledge or simply beliefs about  before any observation of samples.This prior knowledge need to be updated with the information in observed sample data , to give the final description about  represent by the posterior distribution (|).This relationship is explicitly formulated by the Bayesian theorem, wherein the likelihood function L (|) = (|) is defined as the conditional probability distribution of the given parameters of the data [34].
We limit the considered target true functions to be deterministic, which means all trials with the same input parameters gives identical function value.Follow the principle of Bayesian statistic, the deterministic response function () of a -dimension variable  = { 1 ,  2 , • • • ,   } can be treated as a realization of stochastic process  () [35]: wherein  ℎ () is function of ,  ℎ are unknown coefficients to be estimated,  () is the error term, which is a random process with zero mean and covariance where   and   are two sets of variables,  is the process variance and Corr   ,   represents the correlation.Equation ( 13) can be regarded as the Bayesian prior on the true response function, in which the right hand side resembles the form of linear regression model, except that the errors of different  are correlated rather then independent with each other.We consider the common and simple case where the stochastic process is stationary, and the correlation takes the specific form of [35] Corr with where  ℎ and  ℎ are parameters that need to be determined.This specific form corresponding to a product of Ornstein-Uhlenbeck process at  ℎ = 1, which is widely used in physical science models.Furthermore, with the correlation form of Equation ( 15) and ( 16), the regression parts ℎ  ℎ  ℎ () in the stochastic process model can be simply replaced by a constant  without undermining the predictive performance [36].Hereafter we use the simplified stochastic process model A typical methods for analyzing such stationary Gaussian process model is the Kriging method [37], which is flexible and robust for global estimation [38] and is widely applied in fields of spatial statistics [37,39], design of computer experiments [35,40] and Bayesian optimization [36].Besides the ordinary Kriging method, various analysis methods based on Kriging method as well as Bayesian approach are also well studied for Gaussian process with non-stationary correlation and other kinds of stochastic process [38,[41][42][43][44].Here we adopt the best linear unbiased predictor (BLUP) of Equation (17) given by the Kriging method to form the estimation model.Supposing the response function has been evaluated at The BLUP of () can be represented as [36] where R is  ×  correlation matrix with entries is the vector of correlation between the errors at sample and untried input , is the generalized least-squares estimation of .
Selecting the parameters in Equation ( 16) properly we can expect a predict function that describing the true function () well using only the limited number of sample points .One useful method to determine these parameters is the maximum likelihood estimation (MLE).For the Gaussian process we take here, the likelihood function of the samples are [36] 1 where |R| represents the determinant of R. The analytical solutions of  and  maximizing Equation (19) are which is exactly the generalized least-squares estimation of , and Such we only need to find the values of  ℎ and  ℎ that give the maximum value of Equation (19), which can be completed by any optimization method at hand.Corresponding to the specific model of two-level NV electron spin system, once the control field () and control time  are fixed, the target function  (, ) is the deterministic response function of two-dimension variable  = {, } , and can be approximately estimated by the predict function f () which takes the same form of ŷ() in Equation (18).Specifically, we write the Hamiltonian in the interaction picture with respect to  0 2   as where   () is the time-depended Hamiltonian.Neglecting counter-rotating terms in a rotating-wave approximation,   () can be expressed as where Ω  () and Ω  () take the form for SFB method, and the form for PM method.The target response function can be expressed as where is the evolution operator and T is the time-order operator.So far we have established the Bayesian based estimation model of the final state fidelity under certain frequency detuning  and amplitude bias ratio .With this usable tool, we could go further to construct a high efficiency optimization method to explore the optimal control field of ensemble system.

C. The hybrid optimization method
Now we introduce the hybrid B-PM optimization method for solving the problem described by Equation (10).As represented in Figure 2, the whole optimization process can be divided into three parts, initially construct the available predict function, then search for the optimal parameter  using the predict function as the objective function, lastly calculate the true value of objective function in Equation (11) with  ×  = 2500 sample points.The first segment is the essential part of this method, and the following strategies are adopted to ensure the convergence and improve the optimization results.1. Sample position: In the sample selection step, instead of completely randomly picking the sample points, we add random bias on even-distributed coordinate values to give the randomly yet uniformly distributed sample position.This choice improves the estimation performance especially in cases with small sample number.

Model parameter:
The model parameter  and  are selected based on a random initial filed  ini (), and is fixed during one optimization process.This tactic ensures the monodromy of the objective function and the convergence of the subsequent search process, reduces the computation time spend on building the predict function, and doesn't damage the estimation accuracy.Detailed verification are showed in Figure 4 of next section.
3. Cross valiation: After constructing a predict function, an cross valiation is applied to eliminate the lowperforming ones.Each function value of the  sample points () is successively regarded as unknown quantity and is predicted based on the model parameter ,  and other  − 1 samples, given a corresponding predict vector  pred ().Then a linear fitting is performed on those points located at (),  pred () and we use the slope  fit as a criterion.In principle, a  fit far away from 1 indicates a bad performance of the predictor, and by examining the fitted slop of bad predict models occurred in the optimization process we found bad models corresponding to small values of  fit .Therefore we set  fit > 0.6 as the threshold.If this condition is fulfilled, the predict function would be used as the objective function in the following direct search process.Otherwise a new model need to be build from the very beginning.
After the valid predictor is constructed, it is used as the objective function of the following direct search process, where we using the phase-modulated to further improve the efficiency.The initial parameter  ini is taken in the range with Ω max being the maximum field amplitude and  the evolution time.The SFB method is also investigated for comparison, with the initial parameter For both methods, we take Ω max = 2 × 10 MHz and  = 100 ns, and the total field amplitude is constrained as √︃ Ω 2  () + Ω 2  () ≤ Ω max , the frequency parameters   ,   and   are constrained in range [0, 5 × 2/], and the phase parameters   and   are constrained in range [0, 2].

D. AC magnetometry
Magnetometry is one of the most promising field that NV center can play an important part with its outstanding characters such as fine biological compatibility, normal pressure and temperature operating conditions, and relatively long spin lifetimes.Pulsed and continuous dynamical decoupling (DD) has been applied in coherent control and magnetic sensing of NV center.In practice, large inhomogeneous broadening of NV center ensemble will decrease the controllability of DD sequence, and ultimately hinder the sensitivity in magnetometry.Adiabatic strategy has been applied to tackle such problem and achieve good results.When the frequency of signal is not low enough to fulfill the adiabatic condition, finding optimal field with limited pulse length to achieve certain gate operation become a key approach to improve the sensitivity.Here we use the B-PM method to optimize the control field and gives robust Pauli-X and Pauli-Y gate under inhomogeneous broadening, which composes the XY-8 DD sequence in the AC magnetometry process, see the sketch in Figure 7 (a).
During the pulse period, the control Hamiltonian  c in Equation ( 22) takes the form of for Pauli-X and for Pauli-Y gate, with sin    .
(32) For gate optimization the objective function takes the form where is the gate fidelity of fixed  and ,  Tar is the target quantum gate,  , is the actual operator,takes the same form as in Equation ( 29), N ,  (  ) and  (  ) take the same form and value as in Equation (11).We apply the B-PM method to search for the robust control field with the amplitude constraint √︃ Ω 2 1 () + Ω 2 2 () Ω max , and obtain an optimized shaped pulse which constructs the XY-8 DD sequence in the AC magnetometry strategy.
As is showed in Figure 7 (a), the AC magnetic signal to be measured oscillates at frequency  s , and the control pulses is applied on time nodes where the AC signal changes its direction.Due to the power limitation of the microwave field, each pulse possesses a finite pulse length  pulse , and its value follows the relationship   = /  pulse +  pulse , where  pulse is the time separation between two adjacent pulses.The amplitude of the AC magnetic signal  ac can be read out from the Ramsey oscillation of the NV center sensor.At the beginning of the measurement, the NV centers in ensemble are initialized into state |0 , flipped by a /2 pulse and begin to rotate around   axis.After a time period , one 3/2 pulse is applied to the sensors and their population on |0 is measured.In ideal condition with instantaneous gate pulses, the population can be expressed as  0 = (1 + cos(2()) ()) /2, with the phase term () ≡ ∫  0  |cos ( s  )|  , and the exponential decay term  () induced by the inhomogeneous broadening and dynamic noise.At coherent time  =  2 , the decay term goes down to 1/.In practice, the finite pulse length leads to a deviation of the real phase term from ().
Specifically, the total Hamiltonian in the rotating wave approximation is where a time-dependent dynamic noise term () is considered, formed by the Ornstein-Uhlenbeck process where  and  is the correlation time and the diffusion constant of the noise respectively, valued as  = 20 s and √︁ /2 = 2 × 50KHz [18], and  ∼ N (0, 1).The controlled Hamiltonian  c only take a non-zero value during the pulse duration  pulse .We take the frequency of the signal to be   = /  pulse +  pulse = 2.5 MHz, where  pulse is the time separation between two adjacent pulses.We take the maximum limit of field amplitude as Ω max = 2 × 10 MHz, so for rectangular pulse,  pulse = 50 ns, during which   = Ω max   for Pauli-X gate and   = Ω max   for Pauli-Y gate, and the time separation is  pulse = 350 ns.For optimization pulses we set  pulse = 100 ns,  pulse = 300 ns.

A. Feasibility of the estimation model
We first demonstrate the feasibility of the estimation model based on its estimation accuracy and time efficiency.Figure 3 shows its performance on predicting the fidelity of final state under certain frequency detuning  and amplitude drift factor .For a fixed control field showed in Figure 3 (a), the value of  (, ) varies with  and , which we show by a 2-D color distribution image in Figure 3 (b), where 50 × 50 sample points are taken to ensure it's resolution.If the available sample number decreases to 9 and 16 respectively, the resolution ratio drops clearly, as shown in Figure 3 (c) and (d).However, based on the same number of sample points, the image depicted by 2500 prediction values f () in Figure 3 (e) and (f) shows an obviously improvement in the resolution ratio and the accuracy compared to their counterpart in Figure 3 (c) and  (d) respectively.By applying the estimation method, more information about how the target function distributed in the parameter space is extracted based on the knowledge provided by the sample points.
Another concern is how much time would this prediction process takes, or if it has advantage on computing time and cost.Here we use the average computing time of 100 prediction process to give a general description.The control field in each process takes the form of () =  cos  0  +   sin () , with evolution time  = 100 ns, and the parameter ,  and  are randomly taken from the range  ∈ [0, 10 × 2] MHz,  ∈ [0, 2/] and  ∈ [0, 2/].Each prediction process using 16 sample point, and the mission is to calculate the average fidelity (see Equation (11) ) in range  ∈ 2 × [−10, 10] MHz

B. Optimization efficiency of the B-PM method
The results in Figure 4 has demonstrated using fixed estimation model, the total computation time is decided by the calling time of the function  (, ). Figure 5 shows a comparison between B-PM method, PM method, Bayesian estimation SFB (B-SFB) method and SFB method in terms of final value of F obj and calling times of function  (, ) during the search process.Results with  D = 1 are showed in Figure 5 (a) and (c), results with  D = 1 are showed in Figure 5 (b) and  (d).Each results are based on 100 trials with random initial parameters.On the whole, the B-PM method gives result of F obj = 0.905 with average calling times 1252, while the SFB method gives result of F obj = 0.894 with average calling times 13479, such the B-PM reaches a improved average fidelity using only 9.3% computation resource compared with SFB method.
Figure 6 shows the detail optimization results of B-PM method with  D = 1,  = 9 and SFB method with  D = 2,  ×  = 16.Figure 6 (a) and (d) shows the final objective function values of 100 trials with random initial parameters , in which we mark the value used in the optimization process by black points, and their corresponding more accurate value calculated by  ×  = 2500 true sample points by red points.As is showed in Figure 6 (a), among 100 trials of B-PM method, 42 results give F obj ≥ 0.9 and 86 results give F obj ≥ 0.87, while for SFB method showed in Figure 6

C. Sensitivity improvement in AC magnetometry
The optimal control field of Pauli-X and Pauli-Y gate given by the B-PM method are showed in Fig.

IV. DISCUSSION
Our work shows the Bayesian based estimation model can be effectively applied to estimating the fidelity of the state transformation and the gate construction of NV center under wide inhomogeneous broadening noise.The estimation accuracy as well as the time efficiency make it an ideal tool to constitute a new practical optimization method, which we denote as Bayesian estimation phase-modulated (B-PM) method.B-PM method can be combined with various of search algorithms, i.e., the direct search algorithm [45,46], the genetic algorithm [47] and the gradient based algorithm [48][49][50][51].
Here we adopt the widely used Nelder-Mead direct search method [52], related optimization tools is easily available on most of common programming platforms.The Nelder-Mead method seeks the best parameter point by successively construct new points to replace the worst point following a heuristic simplex approach.Consequently the objective function need to be frequently called to rank the points during the searching process, which consume the majority of the total processing time.The hybrid Bayesian estimation phasemodulated (B-PM) method reduce the processing time from two aspects.On the one hand, by setting the predict function as the objective function, the time needed to compute the objective function value once is reduced, as showed in Figure 4. On the other hand, the phase-modulated method can significantly decrease the total number that objective function are called, by introducing a simpler landscape of the parameter space.Taken these two factors together, the advantage of the B-PM method becomes obvious compared to the commonly used standard Fourier base (SFB) method.To be specific, compared to the best result given by the SFB method, the B-PM method raises the optimal value of objective function from 0.894 to 0.905, with more than 90% decrease of the average time consump- tion.Moreover, among 100 results begin from random initial parameters, only 1 results of SFB ( D = 2,  ×  = 16) method give F obj ≥ 0.89, while 42 results of B-PM ( D = 1,  = 9) method give F obj ≥ 0.9, indicating the B-PM method is much more robust to reach a fine result when the amount of total trials is limited.When making a comparison between B-PM and PM method, people may query that the improvement induced by Bayesian estimation model is not obvious.This is true for the current system studied in this article, and needs to be tested in more physical models and systems.However, we stress here that based on the results showed by Figure 5, the B-PM method guarantees a higher value of objective function when the computation resources are on the same scale compared with PM method.This is crucial when the computation resource is severely stressed, and when the control performance is sensitive to the precise value of the objective function.
We further display the utility of B-PM method in NV center based AC magnetometry, especially in high frequency cases where the adiabatic strategy is not viable.The numerical simulation result shows the optimized pulse achieves a four-fold extension of the coherence time  2 and a two fold improvement of the sensitivity compared to the conventional rectangular  pulses with the same maximum amplitude.Similar optimization strategy can applied to other sensing cases, e.g., the spin bath driving [23,53,54] process for DC magnetometry sensing.
One important potential application area of B-PM method is the closed-loop optimization [26][27][28] for complex system such as many-body [55]and many electron systems [56].In these cases numerical simulations fail to describe the system accurately due to a lack of information of the complex system, and the experimental results are directly used as the value of objective function in optimization process.The advantage of B-PM method in reducing the requisite total sample number during the process will be more prominent in such conditions since the time needed to experimentally obtain one sample data is generally much longer compared to the simulation process.Besides NV center, the B-PM method is also expected be applied to optimization of other prevailing quantum platforms such as trapped ions [57][58][59], cold atoms and superconducting qubits [60,61].

FIG. 1 .
FIG. 1.(a) Schematic of NV center in diamond lattice.(b) Schematic diagram of the energy level of NV center.(c) Fidelity of state flip of NV center using rectangular control field () = 10 MHz.(d) The probability density of frequency detuning .(e) The probability density of amplitude drift factor . (f) The sampling range of numerical simulating when computing the value of average fidelity F .

FIG. 4 .
FIG. 4. (a) Average computation time of the objective function as a function of calculating sample number  × . 100 processes under random control fields are calculated, and the estimation model parameter is updated in each process when computing Fobj .Insert: The zoom-in graph for  ×  ≤ 100.(b) Average function value deviation as a function of calculating sample number  × .Based on the same data with (a).The deviation is calculated by subtracting the values of  (, ) with  ×  = 2500, i.e., |F obj − F obj(  × =2500) | for true value based process and | Fobj − F obj(  × =2500) | for predictor based process.(c) Average computation time of the objective function as a function of calculating sample number  × , while the estimation model parameter is fixed in all the 100 random processes.Insert: The zoom-in graph for  ×  ≤ 100.(d) Average function value deviation as a function of calculating sample number  × .Based on the same data with (c).The control fields take the form () =  cos  0  +   sin () , with the evolution time  = 100 ns, and ,  and  randomly taken from the range  ∈ [0, 10 × 2] MHz,  ∈ [0, 2/] and  ∈ [0, 2/].The sample number used in the estimation model is  = 16.

Figure 4 (
b) and (d) shows the estimation accuracy of Fobj and F obj with different calculated sample number  × . Figure 4 (b) uses the same data with Figure 4 (a), Figure 4 (d) uses the same data with Figure 4 (c).We use the value of F obj with  ×  = 2500 as the benchmark to calibrate the accuracy, denoted as F obj(  × =2500) .The estimation accuracy is expressed as the value deviation from F obj(  × =2500) , that is | Fobj − F obj(  × =2500) | for estimation objective func-

N D = 1 N D = 2 N D = 2 FIG. 5 .
FIG. 5. (a) The optimized fidelity of B-PM, PM, B-SFB and SFB methods with parameter sets number  D = 1.(b) The optimized fidelity of B-PM, PM, B-SFB and SFB methods with parameter sets number  D = 2. (c, d) Average calling times of function  (, ) in Equation (6)) that gives the results in (a) and (b) respectively.All results are based on 100 trials with random initial parameters.The total evolution time is taken as  = 100 ns, the maximal field amplitude is bounded as max |()| Ω max = 2 × 10 MHz.
(b), only 3 results out of the 100 trials give F obj ≥ 0.87. Figure 6 (b) and (e) display the shape of optimized control field in the interaction picture, Figure 6 (c) and (f) display the fidelity distribution in region  ∈ 2 × [−10, 10] MHz and  ∈ [0.5, 1.5].
(a), with  D = 2,  = {0.0583,0.0046},  = {0.0844,0.1493} and  = {0.0307,0.0413}.Figure 7 (b) shows the population of the sensor in |0 , measured at every terminal of single XY-8 period, such the time interval between two data points is 3.2 s.1000 measurement is made under different  obey Gaussian distribution with zero mean value and FWHM = 2 × 26.5 MHz.Because of inhomogeneous broaden of  as well as the dynamic noise (), the population of sensor decays with time, featured by  2 ≈ 180 s for rectangular XY-8 pulse and  2 ≈ 1500 s for PM XY-8 pulse respectively.