Time-adaptive phase estimation

Phase estimation has found an increasing number of applications in metrology and quantum computing in recent years. Although resources are considered differently in these two settings [3, 4, 5], the methods of each have led to new applications in the other. In metrology, classical averaging limits to estimates with phase uncertainty $\Delta\hat{\phi}$ scaling at best according to the standard quantum limit (SQL), $\Delta\hat{\phi}>1/\sqrt{N}$, with the number of resources $N$, while general strategies are fundamentally limited only by Heisenberg’s uncertainty principle and can obtain a $\sqrt{N}$ improvement over the SQL reaching the so-called Heisenberg limit (HL), $\Delta\hat{\phi}>\pi/N$ [6]. Entanglement was initially thought to be a key ingredient in schemes to reach Heisenberg scaling, $\Delta\hat{\phi}\propto 1/N$ [7], but when it was also found that the same scaling could be reached using sequentially prepared unentangled systems [8, 9, 10, 11, 12, 13, 14]¹¹1This is not necessarily true when noise is considered [88]., ideas from quantum computing [16, 17] soon led to new experimentally accessible metrology procedures with this scaling [18]. On the other hand, there have been several more recent proposals to use metrology methods to calibrate operations for quantum computing [19, 20, 21, 22, 1].

Many metrology proposals make use of a Bayesian approach to estimation which provides a natural framework to describe adaptive procedures, where the settings for future experiments are modified based on the results of previous measurements [23, 24, 25, 26, 27, 18, 28, 29, 30, 31]. While adaptive methods can lead to better performance, they are generally more complex than non-adaptive strategies and can also be more difficult to implement in some experiments. Remarkably, Higgins et al. found that they could reach near-optimal Heisenberg scaling using a non-adaptive procedure by optimising the number of measurements performed with different coherent applications of the unknown phase [32, 33].

While initial proposals like that of Higgins et al. [32] considered ideal settings with pure states and unitary operations, the role of noise and experimental imperfections has been increasingly studied over time [34]. This development has been both in the understanding of the fundamental limits to precision [35, 36, 37, 38, 39, 40, 41, 42, 43], and in devising better strategies to cope with non-ideal conditions, where the Bayesian framework with adaptive measurements has proved useful [44, 45, 46, 47, 48, 49, 50]. In some cases, adaptive methods have been used to design strategies better suited to the physics of particular experiments. In recent developments for the sensing of magnetic fields with NV centres, adaptive procedures were used to account for reduced visibility measurements [51, 52, 53, 54, 55, 56, 57, 58]²²2Readout can also be improved by adaptive methods [89, 90]. In many proposals the adaptive control acts on a phase that can be seen as an adjustment of the measurement basis, but the time of coherent interaction is chosen non-adaptively, e.g. as detailed in [32]. Since the interaction times proposed in [32] are optimised without noise, and thus may not be suited to experimental conditions, some experimenters perform optimisations using numerical simulations that include relevant noise in order to find interaction times better suited to their experiments[60, 61, 62, 53, 63]. Recently, Belliardo and Giovannetti [64] have analytically shown how to modify coherent interaction times in the non-adaptive procedure of Higgins et al. [32] to account for certain types of noise. Others have investigated the possibility of also choosing the interaction time, or equivalently, in the case of some optical measurements, the size of a so-called N00N-state [33]³³3This is an example where entanglement can be used to convert temporal resources into spacial resources [11]., adaptively: initially in proposals without noise [27, 33], or with mixed-state quantum computation [28], and later in proposals using numerical algorithms that can also account for noise and imperfections [66, 67, 68, 57, 69].

Commonly the performance of a strategy is only considered in the asymptotic regime, when the number of resources (i.e. number of physical systems and total estimation time) approaches infinity. Many experiments, however, could benefit from strategies that are optimised for finite resources. Phase estimation in metrology is typically studied in two settings. In one setting, often referred to as local, the goal is to achieve optimal sensitivity to phase fluctuations from an a priori known phase $\phi=\phi_{0}$. In the other setting, referred to as global, it is assumed that there is initially no a priori knowledge of the phase [47]. While the amount of a priori knowledge of the phase is not relevant in the asymptotic regime, it can significantly change the optimal estimation strategies for finite estimation times [70]. Therefore, it is useful to have methods that can account for arbitrary a priori phase knowledge.

In addition, the coherent interaction time in many proposals is optimised under the assumption that this time dominates the experiment time. While this is true in the asymptotic regime, many experiments may not reach this regime for all or at least a significant portion of the estimation time [54]. An extreme example is found in GaAs quantum dots where the measurement time is several orders of magnitude longer than the coherent evolution [71]; here the authors show an exponential improvement in the mean-square error (MSE) of parameter estimation by adaptively choosing the time of coherent evolution.

In this paper we present time-adaptive phase estimation (TAPE), a method that allows for the adaptive optimisation of the coherent interaction time and a control phase to be adjusted to the resources of the experiment. TAPE can provide strategies optimised both when the experiment time is proportional to the time of coherent interaction, and when the experiment time is proportional to the number of measurements; the method also allows for any resource allocation in between these two extreme cases. In addition TAPE can provide optimal strategies for arbitrary a priori knowledge of the phase, because the choice of measurement settings depends only on the prior knowledge of the phase after the last measurement, and not on any record of previous measurements. In contrast to earlier works investigating similar adaptive procedures [27, 33], we use a more general form for the measurement probabilities so that we can account for many types of noise or imperfections.

We propose and analyse several different objective functions for the adaptive parameter selection that lead to near-optimal performance with respect to known theoretical bounds, with and without noise. When the experiment time is proportional to the time of coherent interaction we reach uncertainties in the phase estimates that are a factor of $1.43$ larger than the HL. In addition, we find that TAPE is quite robust to errors that are not accounted for in the model of the estimator. This demonstrates that, similarly to proposals like robust phase estimation [1], our methods are well-suited to calibration of single-qubit operations in the context of quantum computing.

TAPE uses a numerical representation of the phase knowledge as a Fourier series that can easily describe arbitrary prior knowledge, at the possible expense of increased computation. Since the complexity of the Fourier representation (i.e. number of coefficients to track in memory) increases with phase knowledge, the computation time required for adaptive control of experimental parameters may become too large for some practical applications. To tackle this limitation we propose a method that reduces the interval over which the Fourier representation of the phase is used. This can significantly reduce the required memory and time of computation.

The evolution and measurement for a single step $s$ of sequential phase estimation can be described by the quantum circuit in Fig. 1 [72]

where $H$ is the Hadamard gate, $R_{z}(-\alpha)=\mathrm{e}^{i\alpha Z/2}$ with Pauli operator $Z$ is a z-rotation, and $\alpha\,\in[0,\pi]$ is a control phase to be optimised. $U|\phi\rangle=\mathrm{e}^{i\phi}|\phi\rangle$ with $|\phi\rangle\in\mathbb{C}^{2^{n}}$ an eigenstate of the unitary $U\in\mathbb{C}^{2^{n}\times 2^{n}}$ with unknown phase $\phi\,\in[0,2\pi]$. $k\in\mathbb{N}^{+}$ determines the number of applications of the unitary evolution. The parameters $\alpha$ and $k$ will be chosen adaptively for each measurement as detailed in section III. We denote the possible measurement outcomes at step $s$ by $\xi\in\{\pm 1\}$, and in order to describe a range of noise or measurement errors [54], we assume the outcome $\xi$ to occur with probability

where $\lambda_{k},\zeta_{k}\in[0,1]$ allow respectively for the description of asymmetry and reduced contrast in the measurement⁴⁴4While $\lambda_{k}<1$ describes a reduced probability for the outcome $\xi=-1$, the reverse situation can always be described by relabelling the outcomes.. The subscript $k$ of $\lambda_{k},\zeta_{k}$ indicates that these values can generally depend on $k$. Since the outcome $\xi$ and the control parameters $\alpha$ and $k$ depend on $s$ we should really denote them $\xi_{s},\alpha_{s},k_{s}$ but we omit the subscripts for simplicity of notation.

As in [27, 33, 51] we use a Fourier series to write the probability density describing our knowledge of the phase $\phi$ at step $s$

Given prior knowledge $p_{s-1}(\phi)$ and measurement outcome $\xi$, we update our state of knowledge using Bayes’ theorem

where

is the posterior probability of outcome $\xi$. Equation (3) specifies how to modify the coefficients $c_{n}^{(s-1)}\rightarrow c_{n}^{(s)}$ given the measurement outcome.

We choose to use the estimator

A nice feature of the Fourier series representation for $p_{s}(\phi)$ is that the estimator (4) depends only on one coefficient. Given an estimate $\hat{\phi}$ of $\phi$ (in general now, not necessarily given by (4)), we quantify the uncertainty in $\hat{\phi}$ by the square root of the Holevo variance [74]:

where $S(\hat{\phi})=|\langle\mathrm{e}^{i\hat{\phi}}\rangle|$ is the sharpness and the angular brackets indicate an average over the estimates $\hat{\phi}$. If the estimate is biased, one can use $S=\langle\cos(\hat{\phi}-\phi)\rangle$ instead, where $\phi$ is the value of the correct system phase [33]. In [26] the authors show that the estimator (4) minimises the Holevo variance (5).

In the previous section we have described how in general to use Bayes’ theorem to update our knowledge of the phase $\phi$ given the measurement outcomes, and how to obtain an estimate $\hat{\phi}$ from the density $p_{s}(\phi)$ at step $s$. The task of achieving minimal uncertainty in the estimates obtained from (4) is then dependent on the choices of the values of the control phase $\alpha$ and the number $k$ of applications of $U$. In a sequential procedure, the state coherently evolves according to $U^{k}$ for a time proportional to $k$, so that the choice of $k$ corresponds to choosing the time of coherent evolution. In the following we will assume such a sequential procedure, but we note that the results can also be applied in the case of some parallel procedures by using entanglement.

Given prior knowledge $p_{s-1}(\phi)$ at step $s-1$ we look to choose the optimal control phase for a given value of $k$ by maximising a function of $p_{s-1}(\phi)$, $k$, and $\alpha$ that quantifies the expected knowledge gain from the next measurement. One might then expect that the optimal choice of $k$ can be determined by maximising the expected knowledge gain over all possible values of $k$. This can be a good choice when the evolution time of $U^{k}$ is negligible compared to other times in the experiment, as in the example of GaAs quantum dots mentioned above [71]. But in general experiments with different values of $k$ require different resources in terms of execution time, so that the expected knowledge gains for different values of $k$ are not directly comparable⁵⁵5Even more generally, which resources are valuable depends on the setting and what the experimenter wants to optimise. e.g. they may have a restricted number of qubits to measure, so that number of measurements becomes the relevant resource. In this sense it is useful to have a method where the optimisation can be adjusted by the experimenter to describe best how they value their resources..

In [27, 33], the authors studied a method where the value of $k$ is chosen adaptively. In [27], $k$ corresponds to the number of photons in a N00N state, however as shown for the noiseless case in [33] the sequential procedure we discuss here is mathematically equivalent. This equivalence also holds under certain types of noise [13, 14]. Physically these two procedures are different but they are connected by the use of entanglement to convert between temporal and spatial resources [11].

In the case that the experiment time is proportional to $k$, the relevant resource is the number of applications of $U$. This way of considering resources, which is typically chosen in quantum metrology, is considered in [27, 33]. In these works the authors calculate the expected (differential) entropy ⁶⁶6The differential entropy is the entropy of a continuous random variable, but it lacks some important properties of the Shannon entropy for discrete random variables. See e.g. [91], chapter 8. In this manuscript we will usually write simply “entropy” when referring to the differential entropy. of the posterior after the next measurement using a Gaussian approximation. Motivated by the resource dependence on $k$ they choose the value of $k$ that minimises the expected entropy divided by $\ln(N)$. Here $N$ refers to the total number of resources used since the beginning of the estimation sequence, where an initially uniform prior is assumed (no prior knowledge of the phase). Since this method requires knowing $N$ assuming an initially uniform prior and because the priors in the later stages of the estimation must be approximately Gaussian, it is not well suited to incorporating arbitrary prior knowledge. Moreover, the approach of dividing the entropy by $\ln(N)$ is based on scaling of the uncertainty in the phase that is only valid without noise.

Here we study two different methods for adaptively choosing $k$ that have the following features:

1. 1.

   The resource dependence on $k$ is adjustable so that it can be chosen to match a given experiment. Since we consider sequential procedures, we take the resource requirement for performing an experiment with a particular value of $k$ to be the time required to perform that experiment.

2. 2.

   The choice of $k$ for a given measurement is determined only by the prior probability density at that step of the estimation sequence.

3. 3.

   The methods can be applied to different functions quantifying the expected knowledge gain associated with a particular choice of $k$. We study two such possible functions possessing different benefits: the expected sharpness gain and the expected (differential) entropy gain.

In the following we show how to calculate these two functions quantifying expected knowledge gain exactly; with the exact expressions in hand these functions can be computed for any prior, not only Gaussian. The expected knowledge gains are calculated from expressions that can describe experiments with noise. This allows the methods to determine good phase estimation procedures for noisy experiments, as well as in the noise-free case.

The representation of phase knowledge using the Fourier series (2) means that the number of non-zero coefficients $c_{n}$ to keep track of increases with our knowledge of the phase. If the number of non-zero coefficients at step $s$ is $\Gamma$, the computation for the Bayesian update (3) requires calculating $\Gamma+k$ new coefficients. Considering the case that we start with a uniform probability density, we initially have $c_{0}=1,\,c_{n}=0,\,\forall n\neq 0$ in (2). After performing some measurements with $N$ applications of $U$ in total, the last update will require calculating $\Gamma=N$ non-zero coefficients.

If we are in a situation which results in Heisenberg scaling, then $N\propto 1/\Delta\hat{\phi}$. If the scaling of $\Delta\hat{\phi}$ is slower than Heisenberg scaling, then it follows from (3) that the number of coefficients $\Gamma$ after reaching a given value of $\Delta\hat{\phi}$ will be greater than for Heisenberg scaling. However, if we assume a probability density $p(\phi)$ which on average preserves its functional form (i.e. the overall shape stays the same and only the width changes) as estimation proceeds and $\Delta\hat{\phi}$ decreases, then from standard Fourier analysis $\Gamma$ should scale as $1/\Delta\hat{\phi}$. This follows from the fact that we expect the square root of the Holevo variance of the Bayesian probability density, $\sqrt{V[p(\phi)]}$, to scale in the same way as $\Delta\hat{\phi}$ with $N$. This suggests that when we converge more slowly than Heisenberg scaling, the description of $p(\phi)$ obtained from (3) uses disproportionately many non-zero coefficients and becomes inefficient. For this case, we should seek a good representation of $p(\phi)$ with fewer coefficients to minimise computational overhead. In this sense we can expect that it should be possible for the number of coefficients to be inversely proportional to the uncertainty, $\Gamma\propto 1/\Delta\hat{\phi}$, even if we do not have Heisenberg scaling.

The computation for adaptively choosing $\alpha$ and $k$ scales more favourably with $\Gamma$ than the Bayesian update. Computing the expected sharpness gain (7) does not scale with the total number of coefficients $\Gamma$, while the expected entropy gain (9) scales as $\Gamma/k$. Since we expect the optimal value of $k$ to increase with $N$ (in the noise-free case), there may exist algorithms, offering a complexity independent of $N$, for choosing $k$ such that the expected entropy gain typically is maximised.

The practicality of using TAPE depends on the possibility of computing the Bayesian update and choice of control parameters ($\alpha$ and $k$) in a time shorter or equal to the time it takes to run the next measurement. Since the update computation time scales linearly with $\Gamma$, it will become longer than the measurement time for sufficiently large $\Gamma$. In the asymptotic limit when $N\rightarrow\infty$, the average measurement time, proportional to $k$ for sequential procedures, will also scale linearly with $N$ if we have Heisenberg scaling. If the slope of the linear increase in computation time is less than that of the measurement time, the computation will remain shorter than the measurement time on average for all $N$. However, this may not be the case for many experiments where TAPE could otherwise be useful, for example if decoherence prevents large $k$-value measurements from giving phase information, or if TAPE is used for parallel procedures. In these cases it is beneficial to find methods that reduce the time of computation. One approach is to change representation from a truncated Fourier series to another function such as a Gaussian which is described by a number of parameters that is independent of $N$ [81]. Here we propose an alternative method that allows the truncated Fourier series representation to be used for arbitrarily large $N$ while ensuring that the number of coefficients needed in the series scales much slower than $1/\Delta\hat{\phi}$. This way all the procedures of sections II and III can still be used for large $N$ with minimal changes, and the user is also free to adjust the generality of the description of $p(\phi)$ versus computational overhead as suits their purposes.

Given a probability density $p(\phi)$, we define a contraction of $p(\phi)$ as $(M,\phi_{0},q(\theta))$, such that $p(\phi_{0}+\theta/M)\approx q(\theta)$ for $0\leq\theta<2\pi$ and $p(\phi)\approx 0$ otherwise. We refer to $M\in\mathbb{N}^{+}$ as the magnification, to $\phi_{0}\in[0,2\pi)$ as the offset, and we represent $q(\theta)$ as a truncated Fourier series. Thus, by representing the density $p(\phi)$ as a contraction we use the truncated Fourier series representation for $\phi_{0}\leq\phi<\phi_{0}+2\pi/M$ and assume that $p(\phi)=0$ for values of $\phi$ outside this interval. We are effectively “zooming in” on a region around the expected value of $\phi$ and assuming that the probability of $\phi$ outside this region is negligible. Given a contraction $(M,\phi_{0},q(\theta))$ where $q(\theta)$ has coefficients $c_{n},\,\forall n=-\Gamma,\dots,\Gamma$ we can calculate a new contraction $(M^{\prime},\phi_{0}^{\prime},q^{\prime}(\theta^{\prime}))$, $M^{\prime}$ a multiple of $M$ with $m=M^{\prime}/M\in\mathbb{N}^{+}$ as

where $\hat{\theta}=\arg(c_{-1})$.

By using the truncated Fourier series on a reduced interval of size $2\pi/M$, a contraction $(M,\phi_{0},q(\theta))$ of $p(\phi)$ (with $\Gamma$ coefficients) uses only $\Gamma/M$ coefficients. The assumption that $p(\phi)=0$ outside the interval ${[\phi_{0},\phi_{0}+2\pi/M)}$ is often false with finite probability. However by performing more measurements this probability can always be reduced arbitrarily while increasing the number of coefficients $\Gamma$. Once this probability is small enough for the given application we can make contractions to ensure that the number of coefficients no longer exceeds $\Gamma$. Thus, with this method we must balance the probability that a contraction fails to represent the estimated phase against the maximum computation time per measurement.

In the Supplemental Material, section S.IV, we show that under the assumption of a Gaussian probability density $p(\phi)$, using contractions reduces the number of Fourier coefficients required by $\mathcal{O}(\sqrt{\log(\log(\Gamma))})$, while keeping the probability that $p(\phi)\neq 0$ outside the interval $[\phi_{0},\phi_{0}+2\pi/M)$ below a constant value (which can be chosen to be arbitrarily small). In general, the computational complexity will depend on the particular functional form of $p(\phi)$, but the Gaussian case demonstrates that the method can sometimes greatly reduce computation. If a much poorer complexity is found when faster computation is needed, it could be helpful to modify the adaptive strategy for choosing $\alpha$ and $k$ such that $p(\phi)$ becomes closer to a Gaussian before each contraction is performed.

The expected sharpness or entropy gain relations (7) and (9) can be used for the contracted Fourier series representation $q(\theta)$ with the replacement $\phi\rightarrow\theta$, although the values of the expected sharpness gain in that case will not correspond to those for the full density $p(\phi)$ ⁷⁷7Although the differential entropy is not invariant when changing the scale, the differential entropy gain is.. However the control parameters that achieve the maximum gain for the contraction can still be used to determine the values to use for the next measurement. Denoting the optimal control phase thus obtained as $\beta$ and the optimal number of applications of $U$ as $j$, then the corresponding values to use for the un-contracted density $p(\phi)$ can be calculated as $k=Mj$, $\alpha=\beta+k\phi_{0}$. In this way we can use contractions to reduce the time not only for the update computation but also for the search of optimal parameters $\alpha$ and $k$ using the expected entropy gain ⁸⁸8this also works when using sharpness gain, but there is no speedup.

In order to characterise the performance of the methods considered above we study the results of numerical simulations. We focus on estimation starting from a uniform prior (initially no information about the phase), though we note that our methods are also well-suited to any prior that can be well-represented with a Fourier series and the contractions described in section IV, given available computational resources and particular timing requirements of an experiment. We first consider the case of noise-free quantum metrology. In this context the time of the measurement is assumed to be proportional to the number of applications $k$ of the unitary evolution $U$; in the following we set the time of $U$ to 1 so that the time of a single measurement is equal to $k$. In this case, the number of resources $N$ in the SQL, $\Delta\hat{\phi}>1/\sqrt{N}$, and HL, $\Delta\hat{\phi}>\pi/N$, is the total number of applications of $U$ in the estimation procedure, and is equal to the total estimation time.

In figure 4 we plot the uncertainty in the phase estimates resulting from repeating the estimation procedure with randomly chosen values for the system phase $\phi$ ⁹⁹9The performance of estimation may be dependent on the system phase when starting from a uniform prior if the control phase for the first measurement is fixed (e.g. $\alpha=0$). By choosing a random control phase $\alpha\in[0,\pi]$ for the first measurement, the performance is the same for any value of the system phase. Thus, our results can also be interpreted as the performance for any particular value of the system phase assuming a random control phase is used in the first measurement.. More specifically, we plot the uncertainty multiplied by total estimation time for which any phase estimation reaching HS results in a constant value independent of estimation time; for the HL, the constant value is $\pi$. The uncertainty in the phase estimates $\Delta\hat{\phi}$ is calculated according to equation (5), where the sharpness is calculated as $S(\hat{\phi})=\langle\cos(\hat{\phi}-\phi)\rangle$ for a true system phase $\phi$, and the expectation value is approximated as an average over simulations. This choice can only lead to larger uncertainty values than the usual sharpness, $|\langle\mathrm{e}^{i(\hat{\phi}-\phi)}\rangle|$, since $\langle\cos(\hat{\phi}-\phi)\rangle=\mathfrak{Re}\{\langle\mathrm{e}^{i(\hat{\phi}-\phi)}\rangle\}$. For the results in figure 4 we find $\mathfrak{Im}\{\langle\mathrm{e}^{i(\hat{\phi}-\phi)}\rangle\}\ll\mathfrak{Re}\{\langle\mathrm{e}^{i(\hat{\phi}-\phi)}\rangle\}$, which shows that the estimator is unbiased and either choice of sharpness measure leads to essentially the same values for the uncertainty. Nevertheless, in all simulation results presented in this manuscript we have calculated the phase uncertainty $\Delta\hat{\phi}$ using $S(\hat{\phi})=\langle\cos(\hat{\phi}-\phi)\rangle$.

For the multi-step (gain rate) method, each plotted point is $(\Delta\hat{\phi})\times(\text{total estimation time})$ calculated from $6\times 10^{4}$ ($3\times 10^{6}$) realisations of estimation. We compare four adaptive methods for the choice of control parameters $\alpha$ and $k$. For the multi-step method we plot separately the results of maximising either the expected sharpness gain (7) or the expected entropy gain (9) for total estimation times $2^{m},\,m=1,\dots,12$. For the gain rate method we plot the results of maximising the expected sharpness gain (7) to determine the choices of $\alpha$ and $k$, but we find that maximising the expected entropy gain (9) performs worse (not shown), which is due to the properties of entropy mentioned in section III.1. We also note that when maximising the expected entropy gain only, we tend to obtain very accurate estimates with a small number of outliers that significantly increase the phase uncertainty.

Instead of studying further the gain rate method maximising only expected entropy gain, we plot the result of a hybrid strategy in which the choices of $\alpha$ and $k$ are determined using the expected entropy gain for at most the first half of the total estimation time, and the expected sharpness gain is used for the remaining time. This approach is motivated by the fact that although maximising ones information about the phase is desired, obtaining a precise estimate for the phase additionally requires a narrow probability density $p_{s}(\phi)$. Since multi-peak density functions can have the same entropy as single-peak functions with a broader peak, maximising the entropy doesn’t necessarily lead to a narrow density. We choose the hybrid method to study the performance of strategies that maximise the information gain, while still ensuring a narrow density $p_{s}(\phi)$ that is required for accurate phase estimation.

As shown by Gutiérrez-Rubio et al. [50], maximising the entropy gain leads to the maximum likelihood of estimating the correct system phase. Another way to understand the drawback of maximising the entropy gain only, is that maximising the likelihood of estimating the system phase $\phi$ does not generally minimise $\Delta\hat{\phi}$, which depends on the entire distribution of phase estimates, call it $q(\hat{\phi})$, rather than the value $q(\hat{\phi}=\phi)$ at a single point. Note that if one could maximise expected gain over an entire measurement sequence used for estimation, maximising the expected sharpness gain should minimise $\Delta\hat{\phi}$ [33]. Since we perform the optimisation only for the next measurement, it becomes interesting to study the hybrid method.

The two plotted methods based on gain rate are for total estimation times $2^{m},\,m=1,\dots,18$, and all methods make use of contractions (section IV) to reduce computation times; if the square root of the Holevo variance of the phase density, $\sqrt{V[p(\phi)]}$ is less than $\frac{\pi}{2^{13}M}$, where $M$ is the current magnification, then we perform at least one measurement using parameters $\alpha,k$ that maximise the expected sharpness gain, and then we perform a contraction with $m=2$. The additional measurement to maximise the expected sharpness gain is used in the hybrid method, and multi-step method maximising entropy to reduce the probability of the system phase being outside the reduced interval over which the Fourier series representation of the phase knowledge is used in the contraction. Since $k$-values are chosen adaptively, a sequence of optimal choices does not generally lead to the total times we choose to plot in figure 4. For the simulations plotted we therefore always restrict the optimisation over $k$ to values less than or equal to the total remaining time. For this reason we have performed simulations for each plotted point independently of the other points for different total times.

The results in figure 4 show that all methods perform better than the SQL. The gain rate method using sharpness gain or the hybrid method perform best, and we fit the results for total estimation times from $2^{12}$ to $2^{18}$ with a functional form $\Delta\hat{\phi}=A\pi N^{-\gamma}$. The values of $A$ and $\gamma$ obtained from the fits are summarised in table 1. The hybrid method clearly reaches HS ($\gamma=1$), while the method maximising sharpness only is extremely close to HS, and both methods are close to the theoretically lowest achievable uncertainty for the values simulated. In particular, the uncertainty for a total time of $2^{18}$ (i.e. the longest time simulated) is a factor of $1.422\pm 0.015$ from the HL when maximising sharpness gain only. For the same total time the hybrid method reaches an uncertainty within a factor of $1.429\pm 0.011$ from the HL.

We attribute the poorer performance of the multi-step method to the fact that the compared gains require only approximately equal time and to the fact that only a local optimum is used (see Supplemental Material, section S.III). Although we expect the performance could be improved by starting the search at a $k$-value dependent on the Holevo variance of the prior (cheap to compute) and additionally using a “search down” method analogous to algorithm S.4 (see Supplemental Material, section S.III), we choose rather to focus on the gain rate method since it demonstrates near-optimal performance in addition to being simpler and computationally more efficient.

In the Supplemental Material, section S.I, we have also compared the multi-step and gain rate methods in the case where we allow only the subset of $k$-values: $k\in\{2^{n}\},\,n=1,2,3,\dots$. In this case the gains compared by the multi-step method require exactly the same time, and the method performs similarly to the gain rate method. However even in this situation the multi-step method at best performs similarly to the hybrid method using the gain rate. For completeness, we have also simulated the gain rate methods studied here when a brute force search for $k$ is used rather than a Fibonacci search. In that case the hybrid method performs best reaching an uncertainty within $1.394\pm 0.010$ of the HL. A summary of simulation results for the best performing methods we have studied is given in table S.1 in the Supplemental Material, section S.I.

In figure 5 we consider the performance of the methods using the gain rate when the measurement time is independent of $k$. This situation describes the limit where the time required to perform the unitary evolution $U$ is negligible compared to other times in the experiment such as state preparation and readout. In this case the total time is equal to the number of measurements. We plot methods based on maximising either sharpness gain or entropy gain as well as a hybrid method, as in figure 4. In addition we compare the case where the gain maximisation is based on the correct rate in this context, i.e. $\mathbf{R}=\mathbf{G}/\mathbf{T}$ where $\mathbf{T}=\{1,1,\dots,1\}$ (so $\mathbf{R}=\mathbf{G}$), with the usual case in metrology: $\mathbf{T}\propto\mathbf{K}$. We plot the mean error resulting from $12\times 10^{6}$ realisations of estimation from 1 to 40 measurements.

The results in figure 5 show that the strategies using the correct rate ($\mathbf{R}=\mathbf{G}$) outperform those optimised for the usual metrology setting with $\mathbf{T}\propto\mathbf{K}$. In this case we find that the method that maximises the expected entropy gain performs best while that maximising expected sharpness gives the lowest accuracy in phase estimation; the hybrid method has a performance in between the two. We fit the results from $t=10$ to $40$ measurements with an exponential decay $\Delta\hat{\phi}=A\exp{(-\kappa t)}$. The values of the fitted parameters $A$ and $\kappa$ for each of the plotted curves in figure 5 are listed in table 2.

Although the performance is better when the correct rate is used, the improvement is not dramatic. This is perhaps not surprising since the quantum phase estimation algorithm (QPEA) [85, 72] is known to be optimal even when the number of applications of $U$ is not the relevant resource (indeed this is also the case in the setting of quantum computation) [3, 4]. Methods based on QPEA also lead to Heisenberg scaling in the context of quantum metrology [33]; that is, estimation procedures using similar allocation of $U$ per measurement have been shown to be optimal in both settings. Nevertheless, we see that in the case of the single-step optimisation we perform for sequential strategies we can obtain a slight improvement by using the appropriate rate in the optimisation. In this setting we attribute the better performance when maximising the information gain rather than the sharpness gain to the fact that higher $k$-values are cheap compared to the metrology setting, allowing strategies that are better quantified by the richer nature of the entropy.

For the remaining simulations we return to the metrology setting where $\mathbf{T}\propto\mathbf{K}$ since in this context we can compare performance with known theoretical bounds. In particular, we study the performance in the presence of noise.

In figure 6 we model a system with dephasing. We assume perfectly prepared $|+\rangle\equiv(|0\rangle+|1\rangle)/\sqrt{2}$ states, and $U=\mathrm{e}^{-iZ\phi/2}$, $Z$ the Pauli-z operator. After each application of $U$ we additionally apply a dephasing channel described by the Kraus operators $K_{0}=\sqrt{\frac{1+\eta}{2}}\openone$, $K_{1}=\sqrt{\frac{1-\eta}{2}}Z$. At the end we perform the measurement in the x-basis. In this case it is possible to have Heisenberg scaling only initially while increasing $k$ is beneficial, but as $k$ is increased dephasing eventually reduces the information available by the measurement and we are restricted to $1/\sqrt{N}$ scaling. The ultimate bound on precision can then be expressed as $c/\sqrt{N}$, where $c$ is a prefactor depending on the noise channel. In the case of dephasing this prefactor has been shown to be equal to at least $\sqrt{1-\eta^{2}}/\eta$ (but it is not known if this bound is tight) [39]. In figure 6 we consider the case where $\eta=0.995$; the bound $c/\sqrt{N}$ is also plotted. As before we simulate the gain rate method for sharpness gain, entropy gain, and the hybrid case described above. We compare these three methods when no decoherence is accounted for in the estimator model ($\zeta_{k}=1$ in equation (1)) and when the decoherence rate is already known ($\zeta_{k}=\mathrm{e}^{-(1-\eta)k}$). We plot the phase uncertainties calculated from $32\times 10^{4}$ realisations of the estimation starting from a uniform prior and system phase chosen uniformly at random.

When decoherence is not accounted for in the estimator model, maximising the entropy leads to the poorest performance since this chooses larger $k$-values for which the system decoheres, leading to reduced information in the measurement. For all methods that do not account for decoherence in the model, estimation works initially, but as $k$ is increased decoherence eventually introduces large random errors in the estimates. When decoherence is accounted for in the model we find that all methods perform well. When sufficiently many resources are used the method that maximises the rate of entropy gain performs similarly to the other methods. This can be explained by the fact that decoherence eventually limits the $k$-values chosen by the strategy and the density $p_{s}(\phi)$ is far less likely to have multiple peaks. We find that when the decoherence is accounted for, the performance after a total estimation time of $2^{15}$ is within $1.6715\pm 0.0021$ of the theoretical bound when maximising the expected entropy gain, $1.6807\pm 0.0063$ when maximising the expected sharpness gain, and $1.6734\pm 0.0021$ for the hybrid method. Our methods could also be combined with an optimised estimation of decoherence timescales as proposed by Arshad et al. [79].

The results in figure 6 demonstrate some robustness of the estimation methods to errors not accounted for in the model. In figure 7 we present the results of a last set of simulations to examine further the robustness of the methods. In this case we assume perfectly prepared $|0\rangle$ states and $U=\mathrm{e}^{iY\phi/2}$. After the $k$ applications of $U$ we apply a bit-flip channel with Kraus operators $K_{0}=\sqrt{1-p_{b}}\openone$, $K_{1}=\sqrt{p_{b}}X$, $X,Y$ Pauli-x,y, followed by spontaneous emission with Kraus operators

We choose to set $p_{b}=p/2$, $p_{s}=p$ and simulate the hybrid gain rate method with $p=0.1,0.2,0.3$. For each value of $p$ we simulate estimation using an estimator model with ($\lambda_{k}=\zeta_{k}=(1-p),\,\forall k$) and without ($\lambda_{k}=\zeta_{k}=1,\,\forall k$) errors included.

We see that when errors are included in the model the phase uncertainty is significantly larger than in the noise-free case, but we still retain better than classical scaling for longer estimation times. We fit the phase uncertainty with a functional form $\Delta\hat{\phi}=A\pi N^{-\gamma}$, with fit parameters $A$ and $\gamma$. By fitting subsets of points we see that $\gamma$ tends to increase for longer times, but we would need to simulate further to see if HS is eventually reached. When errors are not included in the model, only the $p=0.1$ case performs similarly to the SQL initially, while for larger errors the performance drops significantly. We see also that for longer estimation times the lack of errors in the model can prevent further reduction of the phase uncertainty. This suggests that it may be better to restrict to lower $k$-values in the presence of large unknown errors. Overall, these results demonstrate robustness of the method to large errors indicating that TAPE can also be useful for calibration of quantum systems, e.g. for quantum computing, as considered in ref. [1].

In the noise-free case we find average computation times are $\sim 1\,\mathrm{ms}$ for the computation required for a single shot using a standard PC (CPU model: Intel(R) Core(TM) i7-8565U CPU @ 1.80GHz) without any parallel computation. When noise is included in the model we find that the maximisation of expected knowledge gain can take a time up to $\sim 3\,\mathrm{ms}$ for the entropy gain, and $\sim 1\,\mathrm{ms}$ when maximising the expected sharpness gain. This makes the method suitable for trapped ion or neutral atom based quantum computing systems where the time of a single shot is typically a few milliseconds. Since the computation of expected knowledge gain for each $k$-value could be performed in parallel it may be possible to obtain sufficient speedups for some applications in superconducting circuits or NV centres. Further speedups may also be possible when maximising the expected sharpness gain if the optimal value of $\alpha$ where determined analytically. An alternative approach to achieve faster computation with the implementation we use here is to restrict measurements to use only a subset of possible $k$-values; a further discussion can be found in the Supplemental Material, section S.I.

For experimental systems with sufficiently low measurement rates, TAPE thus provides a near-optimal and very flexible method for phase estimation. In the metrology setting without noise, the best performance we know of for phase estimation using a classical sequential strategy is the adaptive method of Higgins et al. that is inspired by QPEA from quantum computing [18, 4]. Their adaptive method reaches a phase uncertainty (as quantified by the square root of the Holevo variance of the estimates) a factor of $1.56$ larger than the HL, while they have later devised a non-adaptive method, also inspired by QPEA, that performs similarly, demonstrating uncertainty less than a factor of $2.03$ larger the HL [32, 33]. The latter non-adaptive method has also been shown to be a robust method for calibrating single-qubit gates for quantum computation [1]. Non-classical strategies for phase estimation reaching lower uncertainty estimates are known, as e.g. the method of Pezzè and Smerzi [86, 87] which provides phase estimates with uncertainty 1.27 larger than the HL. However, classical sequential strategies are best suited to single-qubit gate calibration. In the metrology setting we have found that TAPE reaches an uncertainty that is a factor of $1.43$ larger than the HL using the hybrid method, demonstrating similar or better performance to known methods. The results presented in figures 6 and 7 also show significant robustness to errors, demonstrating that TAPE can be a good method for calibrating single-qubit operations.

In contrast to algorithms inspired by QPEA where a predetermined number of measurements is performed for values of $k\in 2^{n},\,n=1,2,3,\dots$, TAPE chooses the time of phase evolution (i.e. the $k$-value) adaptively. For the non-adaptive method that is used for robust phase estimation [1, 32], the number of measurements for each $k$-value is optimised beforehand for a given total number of phase applications $N=\sum_{s}k_{s}$ assuming no prior knowledge of the phase. If noise is present the optimisation needs to be modified as shown in [64]. But as noted therein, the noise models they consider “are to be thought more as toy models” that “capture some of the key features of those scenarios”. Depending on the type of noise present in a given experiment a further analysis and optimisation of the number of measurements to perform with each $k$-value will be required to minimise the phase uncertainty. By using the general form for the measurement probabilities (1), TAPE allows for the description of a wide range of noise to be included in the model of the estimator and directly provides near-optimal phase estimation procedures by accounting for the modelled noise in the optimisation for the control parameters $\alpha$ and $k$.

In addition, TAPE allows the exact experiment times to be easily included in the optimisation for the adaptive choice of $k$-value. This is very convenient for experiments where state preparation and readout times cannot be neglected in comparison to the time required to apply the unknown phase. Using the QPEA inspired methods above, a further optimisation of the number of measurements with each $k$-value would otherwise be needed to minimise the phase uncertainty.

In the QPEA inspired methods above, the value of $k$ to use for a particular step in the estimation sequence requires knowing how many measurements have been performed so far with each $k$-value. In TAPE, once the model of the noise is specified (by setting the values $\lambda_{k},\,\zeta_{k}\,\forall k$), and the experiment times required for each $k$-value are set, the values of the control parameters $\alpha$ and $k$ are determined using only the current prior knowledge density $p(\phi)$. This makes it easy to apply TAPE in situations where some prior knowledge may be available. As an example suppose we would like to use phase estimation for calibrating single-qubit operations on a quantum computer where internal parameters of the device can drift over time. One could use either TAPE or a QPEA inspired method to initially estimate parameters. However, it would be easy to include the drift rate in the model of the phase knowledge $p(\phi)$ as e.g. a broadening of the probability density function over time. Then one could use TAPE to perform a minimal number of measurements to keep track of the parameters needed for single-qubit operations over time.

Other proposals for phase estimation such as [27, 67, 31], and some discussed in [25] are more similar to TAPE in that they choose the value of $k$ adaptively. However, they do not demonstrate better performance in terms of uncertainty of phase estimates or flexibility in terms of using potentially available prior information or accounting for experimental resources and imperfections.

Between the different forms of TAPE compared we find that choosing the control parameters for the phase $\alpha$ and number of unknown phase applications $k$ based on the rate of knowledge gain gives near-optimal performance in several different settings, while requiring computation times that make it accessible to many experiments. In the context of noise-free quantum metrology we reach uncertainties in phase estimates within $1.43$ of the HL using a hybrid method maximising the rate of expected entropy and sharpness gains. In addition, we have found uncertainties within $1.39$ of the HL using the hybrid gain rate method performing a brute-force search over measurement setting $k$, rather than a Fibonacci search (Supplemental Material, section S.I, table S.1). Performing the computations for each $k$-value in parallel would allow this to be done in times comparable to those we find for performing a computation in series with the Fibonacci search method, or even faster.

In a setting where experiment times are proportional to the number of measurements rather than to the number of unknown phase applications $k$, we find that maximising the information gain only, leads to the best performance, while the hybrid strategy performs only slightly worse. The method is also able to find optimal strategies in the presence of different types of noise, and demonstrates significant robustness to errors. Combined with the fact that the optimisation can be easily tuned to the real times of experiments as a function of $k$ and can be used with arbitrary prior information, TAPE thus provides an extremely versatile phase estimation method that can directly give optimal performance in a wide range of experimental settings.

The core implementation of TAPE used for all simulations in this work is available at:

https://bitbucket.org/brennann/qtape/src/main/

Initial theory for phase estimation in the noise-free case for $k=1$ including derivation of an expression for the expected entropy gain was done by A.V.L., and the resulting adaptive method was implemented by V.N. for a trapped-ion experiment in the group of J.P.H. to perform adaptive Ramsey measurements. B.N. had the idea to choose different $k$-values adaptively using the same formalism, and worked out the theory in the general case allowing for a range of noise. B.N. devised the gain rate and multi-step methods and chose to study the sharpness as a measure of knowledge in addition to the entropy. The manuscript was written by B.N. with input from all the authors.