Machine optimized quantum metrology of concurrent entanglement generation and sensing

Entanglement is one of the key ingredients for enhancing the measurement precision of quantum sensors. Generally, there is a trade-off between state preparation and sensing within a limited coherence time. To fully exploit temporal resources, concurrent entanglement generation and sensing with designed sequence of rotations are proposed. Based on twist-and-turn dynamics, modulated rotations along only one axis may be sufficient to drive the state to the optimal one for tiny estimated parameter. However, when the estimated parameter is not tiny, it may impact the evolved state and hence degrade the final measurement precision. Here, we introduce another modulated rotations along different axis and find out the optimal control sequences by means of machine optimization. The optimal measurement precision bounds become independent on the estimated parameter, which improves the dynamic range of the machine designed sensors. Particularly, by optimizing the interaction strength for different particle number and the time-modulated rotations along two different axes via machine optimization, the Heisenberg-limited precision scaling can be attained. Our work points out a way for designing optimized quantum-enhanced metrology protocols, which is promising for developing practical quantum sensors.

A general quantum parameter estimation procedure has three stages [39,40]: state preparation, sensing, and readout measurement. A quantum state is prepared and then evolves under a parameter-dependent time-evolution. Finally certain observable on the evolved state is measured and used for extracting the parameter. It is fit to regard the quantum parameter estimation procedure as an optimal control problem, where the optimal control protocol can be found by maximizing the performance metric for given resources [41,42]. To achieve the best measurement precision, all stages need to be optimized. Recently, the machine optimization techniques are applied with success to design optimal state preparation [42][43][44][45], sensing [46][47][48] or measurement [45,49] separately for better estimation precision.
For the traditional three-stage schemes, there is a trade-off between preparation time T pre and sensing time T sen because of the limited coherence time in realistic experiments [50,51]. Longer preparation time T pre could introduce more metrologically-useful entanglement while longer sensing times T sen could accumulate a better phase for estimation. In order to maximize the benefits of entanglement for a fixed total duration T = T pre + T sen , schemes on concurrent entanglement generation and sensing are proposed [52,53]. Combining these two stages can devote all the time for parameter-encoding and yield better measurement precision with the same temporal resource T. For example, a chaotic quantum sensing scheme based on the continuous harmonic driving is proposed in Bose-Josephson system [54]. By applying nonlinear control pulses to the dynamics of the quantum sensor, quantum enhancements in measurement precision can be achieved [55,56]. In addition, by combining the OAT or TNT dynamics with a machine-designed rotation time-sequence, a higher sensitivity can be achieved [53].
For the schemes of concurrent entanglement generation and sensing, the value of the estimated parameter will have an influence on the dynamics and hence affect the final measurement precision. Take the TNT dynamics for example, the phase accumulation related to the estimated parameter ω can be equivalent to the rotation around z axis on the generalized Bloch sphere [1]. When ω ≈ 0, the effect of the rotation around z axis is tiny. Then, starting from an initial SCS, only the time-dependent modulation around x axis via machine optimization is enough to achieve the high-precision measurement [53]. Differently, when |ω| 0, the effect of the rotation around z axis plays a role and breaks the parity-symmetry [57]. Under this influence, only the time-dependent modulation around x axis may not be sufficient and the evolved state may be distorted far from the optimal one and hence degrade the measurement precision. For practical atomic clocks [6] or magnetometers [15], ω is unknown and one need to estimate ω with high-precision in a wide range [4]. Therefore, can we find the optimal control to compensate the influence of the estimated parameter and achieve high-precision measurement with a high dynamic range?
In this article, we investigate the machine optimization of concurrent entanglement generation and sensing via many-body quantum interferometry with TNT dynamics. By introducing additional variational rotation along other direction, we can find the optimal control sequences to cancel out the influence of the estimated parameter and obtain a higher dynamic range. Our article is structured as follows. In section 2, we will briefly introduce the model of interest and some basic concepts in quantum metrology. In section 3, we will first introduce the machine optimization with single variable and then show the influences under different estimated parameters. In section 4, we will present the multi-variable machine optimization scheme, and show how to eliminate the influences of the estimated parameter on the final measurement precision. By introducing another time-dependent variable and maximizing the quantum and classical Fisher information (CFI) respectively, we obtain the measurement precision bounds and the corresponding optimal control sequences. The measurement precision scalings versus the total particle number are also carefully analyzed. Further, the performance of our scheme under imperfect detection is discussed and the CFI in the presence of detection noise is analyzed. Finally, we will give a brief summary and discussion in section 5.

Model
In the pseudospin picture, an ensemble of N two-level atoms can be regraded as N identical spin-1 2 particles. Such a system can be conveniently described by a collective spin with spin length J = N 2 and displayed on a generalized Bloch sphere. The collective spin operatorĴ contains three components, J x = 1 2 (â †b +b †â ),Ĵ y = i 2 (â †b −b †â ),Ĵ z = 1 2 (b †b −â †â ), whereâ andb are the annihilation operators for particles in level |a and |b , respectively. These collective spin operators obey the general angular momentum commutation relations where ijk is the Levi-Civita symbol. A system state can be expressed by |Ψ = N/2 m=−N/2 C m |m , where |m is the Dicke basis denoting N/2 − m particles in |a and N/2 + m particles in |b . For convenience, we set = 1 below.
For concurrent entanglement generation and sensing, the general TNT Hamiltonian in atomic ensemble can be described asĤ The first term is the twisting dynamics with twisting strength χ and generates the entanglement between particles. The second term is the turning dynamics and Ω x (t) is the modulated rotations control aroundĴ x axis. The last term represents the sensing procedure and leads to a rotation aroundĴ z during the dynamics, which is related to the estimated parameter ω. For a duration time T, an initial pure state |Ψ in evolves into an output state |Ψ out (ω, T) under the Hamiltonian (2). One can perform certain measurementX on the output state and extract the information of ω. The estimation precision of ω is bounded by the Cramér-Rao bound (CRB), Δω Δω CRB ≡ 1/vF C , where v is the number of repeats of the estimation procedure, and is the CFI, which is a measure of the ability for extracting the parameter ω from the measurement data of the final states. P(α|ω) is the conditional probability of the measurementX with outcome α given the parameter ω. The CFI can be maximized by going through all possible generalized quantum measurements and strategies. The corresponding measurement precision is bounded by the quantum Cramér-Rao bound (QCRB), Δω Δω CRB Δω QCRB ≡ 1/vF Q . Here, the quantum Fisher information (QFI) is expressed as which is a function only dependent on the final state |Ψ out (ω, T) and its derivative |∂ ω Ψ out = ∂|Ψ out (ω, T) /∂ω with respect to ω. In the following, we select QFI and CFI as the objective function to obtain the optimal control by machine optimization under different conditions.

Machine optimization with single variable
In this section, we first revisit part of the results in reference [53], where the pioneer example of a machine-designed sensor is presented. For the Hamiltonian (2) with duration time T, by dividing T into n equal segments, the time-dependent rotation controls can be parametrized as with Λ x (t) a piecewise step function where in each segment can be varied individually (k = 1, 2, . . . , n). Thus, Ω x (t) involves n variational parameters. Starting from an initial SCS, for fixed χT and ω, one can obtain the maximal QFI F Q (ω, T) by optimizing Ω x (t). With optimization algorithms (see appendix B for details), one can find the optimal controls of Λ x (t) that maximizes the QFI. This protocol has recently been demonstrated that the QFI can be greatly improved compared with the traditional OAT and TNT schemes [53].
For an example with N = 100, χT = 0.1 and ω = 0, the QFI can be achieved approximately 29N with N the QFI of SQL. This shows the ability for sensitive estimation for tiny parameter in the vicinity of ω = 0. For better visualization, we display the Husimi distributions of the evolved states obtained via machine optimization of Ω x (t) in figure 1. The column c in figure 1 is the Husimi distributions of the evolved states at time t = 0, 0.5T, T with the special case of ω = 0. The time-dependent rotation aroundĴ x speed up the twisting dynamics. The state evolves from an initial SCS to squeezed, over-squeezed, and fragmented states sequentially. In the end, the quantum state becomes 'ribbon-shaped', localized at the north and south poles of the generalized Bloch sphere.
For general case with |ω| > 0, we find that the maximal QFI decreases when the estimated parameter is not tiny. Pictorially, there is a non-negligible rotation aroundĴ z during the evolution in the presence of large ω. Under the joint action ofĴ x andĴ z , the evolved state no longer preserves its parity symmetry and deviate from the optimal one. In figure 1, the remaining columns (a), (b), (d) and(e) represent the Husimi distributions of the evolved states with ω = −π/2, −π/4, π/4, π/2, respectively. The corresponding optimal QFI dramatically decreases as |ω| gets larger from 0. The changes of QFI with ω are shown with red dotted line in figure 2. When |ω| > 0, the optimal QFI is smaller than the special case of ω = 0.

Machine optimization with multiple variables
In practice, ω can be large and only optimizing one rotation aroundĴ x may not be the optimal scheme when |ω| > 0. To compensate the influence of ω, we propose a ω-independent machine optimization scheme with multiple variables.

Optimization of Ω x (t) and Ω y (t)
Besides Ω x (t), we add another modulated controls Ω y (t) aroundĴ y . The Hamiltonian we consider becomeŝ where Ω x (t) is parametrized as equations (5) and (6), and Ω y (t) is parametrized the same way as Ω x (t), i.e., and In this case, the time-dependent variables Ω x (t) and Ω y (t) totally involve 2n variational parameters. Still, we start from an initial SCS, and let the state evolve under the Hamiltonian (7) for T. Similarly, we set the time segments n = 10, the particle number N = 100 and χT = 0.1. A set of initial values of Λ x and Λ y are chosen randomly and vary individually in each segment. Through the method of machine optimization, we can obtain the optimum Λ x and Λ y that maximizing the QFI. The machine optimized QFI as a function of ω is shown with red solid lin in figure 2.
With the additional modulated controls aroundĴ y , we can find that the optimal QFI versus ω almost remains unchanged. For ω = 0, the optimal QFI is consistent with the one obtained via machine optimization of Ω x (t). As ω deviates from zero, the optimal QFI with the additional modulated controls aroundĴ y remains unchanged while it decays quickly with ω in single variable machine optimization scheme. The bi-modulation of Ω x (t) and Ω y (t) can perfectly compensate the influence of the term ωĴ z , and the final state at time T keeps nearly the same structure as the one in figure 1(c3). Figure 3 displays the evolved states obtained via machine optimization of both Ω x (t) and Ω y (t) at time t = 0, 0.5T, T for different ω on the generalized Bloch sphere. Unlike the machine optimization with only Ω x (t), the final evolved states recover 'ribbon-shaped' for different ω. The difference between these spin states is the center of the distribution, owing to the different rotation rates ω aroundĴ z . In other words, the bi-modulation aroundĴ x andĴ y makes the final evolved states close to the optimal one and keeps the QFI nearly unchanged while ωĴ z plays a role.
As a reference, the optimization result of Λ x (t) obtained from the machine optimization with only Ω x (t) in the case of ω = 0 is shown in the first column of figure 4. In the second column of figure 4, the optimization results of Λ x (t) and Λ y (t) for ω = 0 are shown. Λ x (t) are almost the same with the reference sequence and the modulation amplitudes of Λ y (t) are very close to zero, which is consistent with the result of only modulating Λ x (t). As |ω| > 0 gradually increases, the optimization sequences of Λ x (t) still have a similar shape to the reference sequence, while the modulation aroundĴ y begins to work. The modulated amplitudes of Λ y (t) increase obviously.
Up till now, we select the QFI as the objective function for optimization. The QFI sets a mathematically ultimate precision bound, which is only relevant to the final evolved state and has nothing to do with which observable is selected. In practice, one should select certain physical observable for measurement, and the CFI is more practical than the QFI. Thus, instead of QFI, we change the objective function to the CFI. Specifically, we choose the half of population differenceĴ z and use equation (3) to calculate the corresponding CFI. Following the same optimization procedure, we obtain the results for different ω, see figure 2. Similarly, we find that the optimal CFI decreases as |ω| > 0 gets larger when only Ω x (t) is optimized. While using the machine optimization with double variables Ω x (t) and Ω y (t), the influence of ω can be eliminated with the optimal CFI almost independent on ω itself. Although the optimal CFI is a bit smaller than the QFI, it still surpasses the SQL.

Pre-estimation on the unknown parameter
Our optimization scheme is based on the precise knowledge of the parameter to be estimated. However, the parameter to be estimated is unknown and it may leads to the infeasibility of our optimization in real parameter estimation process. In order to remedy this defect, a rough pre-estimation of the unknown parameters can be introduced [58,59]. We replace the true value with an estimated value by a pre-estimation process (see appendix A for details). To be specific, an SCS undergo a complete parameter estimation process and we get a first estimated value ω est . Since there is no entanglement, the best achievable estimation precision δω is equal to 1/ √ NT and the true value is in the range of ω est ± 1 δω . Now the selected objective function is F Q (ω est , T) and the corresponding optimal control parameters depend on ω est . We use this set of parameters for estimation and further explore how the deviation of the true value affects the estimation precision. For comparison, we normalize the deviation of ω with ω nor = (ω − ω est )/δω and define the decay of QFI (in dB) as follow: F Q (ω nor ; ω est , T) means that QFI versus ω nor for fixed ω est and T. According to figure 5, the change of decay is a parabola, symmetric about the ideal estimated point ω = ω est . Decay grows as ω nor deviates zero and the maximal decay is at the edge of the uncertainty range. The maximal decay is small and it decreases as total particle number grows. The inset shows that the maximal decay remain unchanged for different pre-estimated values ω est . That is, our optimization scheme after a pre-estimation is valid in a high dynamic range of the estimated parameter. Therefore, a rough estimation on ω is sufficient to find out the optimal Ω x (t) and Ω y (t).

Optimization of χ, Ω x (t) and Ω y (t)
After obtaining the optimal QFI with fixed particle number, we can further explore the precision scaling versus particle number. For an estimated parameter ω = 0.1, we first fix χT = 0.1, obtain the optimal QFI via machine optimization with Ω x (t) and Ω y (t) under different particle number and calculate the corresponding ultimate precision bound F −1/2 Q . In figure 6, the precision scaling of the optimization results are indicated by red dots. Unfortunately, we cannot get a well-defined precision scaling with respect to N. When N is small, the precision scaling is only a bit better than the SQL, while as N gradually increases the precision scaling approach the HL.
The reason why the precision scaling transits from SQL to HL in this scheme is straightforward. For a fixed shearing strength χT, as the interaction energy is proportional toĴ 2 z , an initial SCS with different particle number N will evolve into quite different kinds of spin states. For small N, the interaction is too small to twist the SCS. The evolved states are slightly squeezed and the corresponding uncertainty regions in the generalized Bloch sphere are still maintain approximately a circle. With small amount of entanglement generated, the optimal QFI is only a bit better than the SQL. While as N increases, the interaction becomes strong to twist the SCS and the evolved state gradually changes to the over-squeezed states or even the fragmented states. For the strongly fragmented states, the Husimi distributions are located around both north and south poles, leading to the optimal QFI close to the HL.
In order to find a well-defined precision scaling close to the HL, we consider the interaction strength χ can be tunable for different N. For example with Bose condensed atoms, the twisting strength χ can be adjusted by tuning the s-wave scattering lengths via Feshbach resonance [15,35,60] or adjusting the spatial overlap via spin-dependent forces [11,61]. It is worthy to mention that, for a given N, χ becomes a variational parameter but it does not change with time, which is different from the variation of Ω x (t) and Ω y (t). Thus, for a given evolution time T, the variational Hamiltonian becomeŝ With the same technique of machine optimization, we can get the results numerically, see the blue dots in figure 6. The blue solid line is the linear fitting with the blue dots. With the twisting strength χ tunable with respect to N, a well-defined precision scaling can be attained by machine optimization. The slope of fitting is about −0.964 and the correlation coefficient is greater than 0.9999, which indicates that the ultimate precision bound nearly achieves the HL scaling.

Imperfect detection
In realistic experiments, there are lots of imperfections that limit the final estimation precision. The detection noise is one of the most significant imperfections, which is related to the detector performance and leads to inaccurate counting of particles [36,38,62,63]. Generally, the half-population difference measurement is chosen as the observable onto the final state owing to its feasibility in implement. Measurement result can be written as Ĵ z = m=−N/2 N/2 P m (ω)m, where P m (ω) is the probability distribution of the evolved state projecting onto the Dicke basis |J, m . Assuming that the detection noise obeys Gaussian distribution N (m, σ 2 ), the probability distribution becomes σ is the standard deviation of the detection noise and A m is the normalization factor. Since the optimized states are highly entangled, in the presence of detection noise, they may become fragile and the measurement precision may be dramatically reduced. In order to find a robust input state under imperfect detection, we select the CFI with detection noise as the objective function. It can be written as Figure 7 shows the optimal CFI against detection noise σ. It is worth noting that the special case of σ = 0 returns to the ideal detection described by equation (3). Although the optimal CFI decreases as detection noise increases, the estimation precision can still be better than SQL when σ becomes relatively large (above the shaded areas). Moreover, the CFI still shows its characteristics of ω-independent. In figure 7, we results of ω = 0, π/4, π/2 are shown. They show the same trend, which indicates our scheme can still work for a wide range of ω in the presence of imperfect detection.

Summary and discussion
In summary, we have presented a proposal for optimized concurrent entanglement generation and sensing in an atomic ensemble by using the technique of machine optimization. For concurrent entanglement generation and sensing via TNT dynamics, we find that the measurement precision bounds are sensitive to the value of the estimated parameter ω if rotation along only one direction is modulated. Generally, the measurement precision would be degraded if ω becomes large. To compensate the influence of ω and achieve better measurement precision when |ω| > 0, we introduce another modulated rotation along different direction and find out the optimal control protocols by maximizing the QFI as well as the CFI. As a result, the optimal measurement precision bounds become almost independent on the value of the estimated parameter, which is more experimentally applaudable. Further, we also analyze the optimal measurement precision scaling versus particle number. We find that, by optimizing the interaction strength for different particle number and combining the bi-modulation of Ω x (t) and Ω y (t), along with machine optimization, a Heisenberg-limited scaling can be achieved. Finally, we consider the influence of imperfection in our optimization scheme and we show our scheme can still work in the presence of detection noises. Our scheme may be realized in state-of-the-art experiments. In particular, for an ensemble of Bose condensed atoms occupying two hyperfine levels, it is feasible to perform concurrent entanglement generation and sensing via the general TNT Hamiltonian (11) [16,37]. In this system, the modulated controls of Ω x (t) and Ω y (t) could be implemented by precise controls of the radio-frequency or microwave pulses in experiment. The twisting strength χ can be adjusted by tuning the s-wave scattering lengths via a narrow magnetic Feshbach resonance [15,35,60] or adjusting the spatial overlap via spin-dependent forces [11,61]. The unknown parameter ω can be the detuning from the atomic resonance frequency between two hyperfine levels or related to an estimated magnetic filed. In experiment, one can first prepare an SCS by performing a π/2 pulse on the state of all atoms in lower level. In the pre-estimation process [58,59], one can get a rough value of the unknown parameter under the free evolution of Hamiltonian H pre = ωĴ z . Then, one can put this rough estimated value into the optimization procedure and obtain the optimal control sequences of Ω x (t) and Ω y (t). Finally, one can turn on the twisting strength χ and let the system evolve under the optimized control.
From figure 3 it is illustrated that, under different values of ω, all initial SCSs evolve to the optimal ones, the 'ribbon-shaped' states with different rotated angle around the z axis on the generalized Bloch sphere. These states obtained by machine optimization of both Ω x (t) and Ω y (t) have the same structure so that have the same QFI. Thus the measurement precision bounds become independent on the estimated parameter, which improves the dynamic range. It is shown that machine optimization provides a powerful tool for designing optimal quantum metrology protocols. It can also be widely used for improving the performances of the key stages of parameter estimation, such as efficient preparation of highly entangled states [64][65][66][67] and optimal control processes for parameter extraction [63,68,69].

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

Appendix A. Procedure of our scheme
In the main text, we choose several objective functions, including QFI, CFI and CFI with detection noise. We need to traverse all control strategies to maximize the objective function. It is equivalent to finding the extreme value of a multivariate function. Specifically, we divide the total evolution time into N equal segments, and the control parameter Λ i in each segment is an independent variable of the objective function. The value of objective function is a multivariate function with respect to control parameters vector Λ. Therefore, finding the optimal control parameters is equivalent to solving this unconstrained nonlinear optimization problem.
The above analysis works on the premise that the parameter to be estimated is known, but the actual reality is the opposite. So we need to implement a pre-estimation on the unknown parameter first. The general procedure of our optimization scheme is shown in figure A1, which consists of two parts. The first part is pre-estimation. The input SCS evolves according to H pre = ωĴ z and let it undergo a conventional parameter estimation process. We obtain a pre-estimated value ω est . The second part is machine optimization and we replace ω with ω est in the Hamiltonian. It is worth noting that now the input SCS undergo a parameter estimation process with concurrent entanglement generation and sensing. After choosing the objective function, we initialize the control parameters and set a suitable tolerance. Specifically, we choose BFGS algorithm as the machine optimization algorithm. Finally, a set of optimal control parameters can be automatically obtained.

Appendix B. Optimization algorithms
Many numerical methods have been proposed to solve unconstrained nonlinear optimization problems, such as gradient descent method, Newton's method and quasi-Newton method. Gradient descent method is a first-order optimization algorithm, where the objective function updates the parameter values along the opposite direction of the gradient vector at each iteration. Newton's method is a second-order optimization algorithm that makes use of the Hessian matrix and it has faster convergence speed. However, Newton's method requires the calculation of the inverse of the Hessian matrix, which costs abundant computationally resources and may not be stable depending on the properties of the objective function. Quasi-Newton method is an improvement of Newton's method, by approximate the inverse of the Hessian matrix using the gradient.
We will introduce basic principles of Newton's method and quasi-Newton method in detail in the following. For iteration k, the Taylor expansion of a continuous multivariate function f(x) is It may not guarantee that the function value declines steadily or even converge, and eventually leads to the algorithm invalid. This problem can be solved by introducing the step factor λ k . The search direction of each iteration is still d k , but an one-dimensional search along this direction is implemented to find the optimal step factor, that is λ k = arg min λ∈R f (Λ k + λd k ).
Newton's method needs to calculate the inverse of the Hessian matrix in each iteration and it is not always feasible. Quasi-Newton method is proposed to overcome these problems, by replacing the Hessian matrix approximately with a suitable matrix B. This approximation greatly simplifies the calculation of Hessian matrix. Quasi-Newton method has been one of the most widely used methods for nonlinear optimization. Firstly, we take the gradient on both sides of equation (B1) and get Choosing x = x k , equation (B3) becomes We define s k = x k+1 − x k , y k = g k+1 − g k , it can be further written as This formula is called secant condition of quasi-Newton method and it restricted the Hessian matrix (or its inverse matrix) in the iterative process. So next, we just need to find a way to construct the approximate matrix B k of H −1 k . Different construction of B k lead to different kinds of quasi-Newton method, such as DFP algorithm, BFGS algorithm and Broyden's algorithm. In our paper, we select the BFGS algorithm to search for the local extremum values of the objective function. The update of the approximate matrix in BFGS algorithm is defined as The corresponding program block diagram of BFGS algorithm is shown in figure B1. The initial approximation matrixes B 0 is often taken to be just the identity matrix.