# Heisenberg-limited quantum phase estimation of multiple eigenvalues with up
to two control qubits

A. Dutkiewicz  Instituut-Lorentz, Universiteit Leiden, 2300 RA Leiden, The
Netherlands  B.M. Terhal  QuTech, Delft University of Technology, P.O. Box
5046, 2600 GA Delft, The Netherlands and JARA Institute for Quantum
Information, Forschungszentrum Juelich, D-52425 Juelich, Germany  T. O’Brien
Google Quantum AI, 80636 Munich, Germany  Instituut-Lorentz, Universiteit
Leiden, 2300 RA Leiden, The Netherlands

###### Abstract

Quantum phase estimation is a cornerstone in quantum algorithm design,
allowing for the inference of eigenvalues of exponentially-large sparse
matrices. The maximum rate at which these eigenvalues may be learned, –known
as the Heisenberg limit–, is constrained by bounds on the circuit depth
required to simulate an arbitrary Hamiltonian. Single-control qubit variants
of quantum phase estimation have garnered interest in recent years due to
lower circuit depth and minimal qubit overhead. In this work we show that
these methods can achieve the Heisenberg limit, also when one is unable to
prepare eigenstates of the system. Given a quantum subroutine which provides
samples of a ‘phase function’ $g(k)=\sum_{j}A_{j}e^{i\phi_{j}k}$ with unknown
eigenvalue phases $\phi_{j}$ and probabilities $A_{j}$ at quantum cost $O(k)$,
we show how to estimate the phases $\\{\phi_{j}\\}$ with accuracy (root-mean-
square) error $\delta$ for total quantum cost $T=O(\delta^{-1})$. Our scheme
combines the idea of Heisenberg-limited multi-order quantum phase estimation
for a single eigenvalue phase [1, 2] with subroutines with so-called dense
quantum phase estimation which uses classical processing via time-series
analysis for the QEEP problem [3] or the matrix pencil method. For our
algorithm which adaptively fixes the choice for $k$ in $g(k)$ we prove
Heisenberg-limited scaling when we use the time-series/QEEP subroutine. We
present numerical evidence that using the matrix pencil technique the
algorithm can achieve Heisenberg-limited scaling as well.

###### Contents

  1. 1 Introduction
    1. 1.1 Outline
  2. 2 The classical and quantum tasks of phase estimation
    1. 2.1 Phase distance on the circle
    2. 2.2 Limits for single-eigenvalue phase estimation
  3. 3 Defining the task of multiple-eigenvalue phase estimation
    1. 3.1 Methods of dense signal phase estimation
  4. 4 Multiple eigenvalues: multi-order estimation and the phase matching problem
    1. 4.1 Heisenberg-limited algorithm for multiple phases
    2. 4.2 Bounding the error with and without failures
    3. 4.3 Algorithm 4.1 achieves the Heisenberg limit
  5. 5 Numerical implementation
  6. 6 Conclusion
  7. A Proof of Lemma 4.2
  8. B Range of shifted phase estimates

##  1 Introduction

For quantum computers to overcome the $50^{+}$ year head start in research and
development enjoyed by their classical competition, quantum algorithms must
eke out every inch of quantum speedup over their classical counterparts.
Quantum phase estimation (QPE) of a unitary operator $U$, –a BQP-complete
problem [4]–, plays a central or support role in many promising quantum
applications [5, 6, 7].

However, not all flavours of quantum phase estimation are equally powerful. To
perform QPE with accuracy (root-mean-square) error $\delta$, initial
implementations of quantum phase estimation [8] used a
$O(\log(\delta^{-1}))$-qubit control register, and required computation time
scaling as $T=O(\delta^{-2})$, known as the Sampling Noise Limit.

In fact, this ‘textbook’ QPE [8, 9] estimates the first $n$ bits of the binary
expansion of the phase with quantum cost $T=O(2^{n})$, with high probability
and thus one might naively think that $T=O(1/\delta)$, i.e. Heisenberg-
limited, due to the error induced by missing out on the remaining bits of the
phase. However, for textbook QPE, as noted in [1], one gets sampling noise
limit scaling, $T=O(\delta^{-2})$, if we do the proper analysis for the mean-
square-error $\delta^{2}$ which includes contributions from outlying
(unlikely) data.

Separate from this, the requirement to perform applications of $U$ conditional
on a large entangled control register is technically challenging and has
strong coherence requirements. It has been known for a long time [10] that the
control register in QPE can be replaced by a single qubit using classical
feedback and re-preparation of the control qubit, also known as iterative QPE
[11]. For estimating the phase of a single eigenstate, –assuming the
preparation of this eigenstate–, the sampling noise limit can thus simply be
achieved using a single control qubit. In [1] iterative-QPE was extended to
achieve the Heisenberg limit $T=O(\delta^{-1})$. This Heisenberg limit can be
shown, via Cramer-Rao bounds [1], to be a lower bound on the cost of phase
estimation, assuming one cannot fast-forward the unitary $U$ [12]. This type
of estimation have additionally demonstrated a relative robustness to error
[2, 13], which is of interest to NISQ applications. Other analyses of QPE use
maximum-likelihood [14], or Bayesian [13, 15] inference.

The requirement to prepare eigenstates of the unitary $U$ is however not
possible for most applications. In [16] it was demonstrated numerically that
this requirement is unnecessary; one may infer single eigenvalues from mixed
or superposed initial eigenstates using standard classical signal processing
techniques [17]. Under some additional constraints on the system, it was
recently found that these techniques could be performed in the absence of any
control qubits or the need to apply controlled unitary operations [18, 19], a
further significant saving. Due to the need to ‘densely sample’ the phase
function $g(k)=\sum_{j}A_{j}e^{i\phi_{j}k}$, and a lack of optimization of the
classical post-processing techniques, Ref. [16] only achieved sampling-noise-
limited scaling, but not Heisenberg-limited scaling. By dense sampling we mean
that we draw samples from $g(k)$ with $k$ a sequence of integers,
$k=0,1,\ldots,K$ (as opposed to, say, choosing $k=2^{d}$, i.e. exponentially
increasing which is used in textbook QPE and iterative QPE). By a clever
adjustment of the quantum phase estimation problem to target estimation of the
spectral function, Eq. (3), of the input state, Ref. [3] was able to prove
rigorous results, with bounds that were subsequently improved in Ref. [20].
Still, these results fall short of reaching “the Heisenberg limit for the
problem of estimating multiple phases", however that should be defined.

In this work, we demonstrate single-control qubit quantum phase estimation at
a so-called Heisenberg limit. We extend the methods used in Refs. [1, 2] that
obtain Heisenberg-limited scaling for single eigenphases to the multiple-phase
setting by the use of a multi-order scheme and phase matching subroutines
between different orders. We show that to make this phase matching unambiguous
requires the sampling scheme to be adaptive, i.e. the next choice for $k$ of
$g(k)$ depends on the current phase estimates. At each order the multi-order
algorithm requires input from a dense phase estimation method: for a given
order $k$ we use samples from $g(k{\sf k})$ with ${\sf k}=0,1,\ldots,K$. As we
require the freedom to choose $k$ a real number, to be applicable to a
completely general $U$ we must invoke the fractional query algorithm of Ref.
[21], which requires one additional control bit. Using the time-series or QEEP
analysis of Ref. [3] as classical processing subroutine, we are able to obtain
a rigorous proof of Heisenberg-limited scaling of our multi-order scheme.
Using the matrix pencil method analysed in Ref. [16], as such a dense
subroutine, we are able to show numerical results consistent with the
Heisenberg limit, with a performance improvement over the time-series analysis
results.

In essence, our paper is concerned with what choices of $k$ in $g(k)$ and what
classical processing are needed to enable Heisenberg-limited scaling, i.e.
scaling which minimizes the total number of applications $T$ of (controlled)
$U$ given a targeted error $\delta$ with which to estimate multiple eigenvalue
phases of $U$ present in some input state $\left\lvert{\Psi}\right\rangle$. It
can thus be viewed as purely solving a problem of classical signal processing.
This does not mean that such questions are trivial: for example, the question
of how to estimates phases if one is allowed to only get single samples from
$g(k)$ for a set of randomly chosen $k$ relates to the dihedral hidden
subgroup problem in quantum information theory [22]. We note that other work,
based on a Monte Carlo extension of [3], achieving Heisenberg scaling (up to
polylog factors) was recently presented in [23]. Another recent work also
demonstrated numerical evidence for Heisenberg-limited phase estimation using
Bayesian methods [24]. We also note that the information-theoretically optimal
method in [14] which picks random $k$ in $g(k)$ has a classical processing
cost which is linear in the quantum cost $T$, and Theorem 1 in [14] can be
converted to bound the mean-squared-error in estimating a single phase, see
comments below Theorem 2.6 in Section 2.2. In principle this information-
theoretic method can be extended to the case of $n_{\phi}$ eigenvalue phases,
but the classical processing cost will be exponential in $n_{\phi}$ as one
iterates over the possible values of the phases, while our Algorithm 4.1 has a
polynomial quantum and classical cost in terms of the number of phases
$n_{\phi}$.

###  1.1 Outline

We begin this paper in Sec. 2 by defining a few mathematical objects. We
separate the quantum part of a phase estimation problem, namely sampling of
the phase (or signal) function $g(k)$ in Eq. (2) given an input unitary $U$
and input state $|\Psi\rangle$ in Definition 2.3 through running some quantum
circuits, and the classical processing of samples from $g(k)$ to extract the
eigenvalue data of $U$. In Sec. 2.1 we state and prove some needed properties
of the distance between periodic phases. In Sec. 2.2 we prove several Cramer-
Rao bounds on the scaling of the error versus the total quantum cost for the
estimation of a single eigenvalue phase, Theorem 2.6. We state the previous
result on getting Heisenberg-limited scaling for a single eigenvalue phase in
Algorithm 2.7.

In Sec. 3 we properly define a multi-eigenvalue phase estimation problem (with
assumptions) and state properties of the so-called QEEP algorithm in Theorem
3.2. Here we prove some error bounds when one uses this algorithm in the form
of Alg. 3.4 in Lemma 3.5. Alg. 3.4 is the fixed-order subroutine in our final
Algorithm 4.1 which achieves Heisenberg scaling.

In Section 4 we present our Heisenberg-limited algorithm. We discuss a
critical aliasing problem to be solved which occurs when estimating multiple
eigenvalues. We show that an adaptive choice for $k$ in $g(k)$ can solve this
issue and we prove that such adaptive choice always exists in Lemma 4.2. In
Lemma 4.3 and Lemma 4.4 we prove some properties about Algorithm 4.1 which
will be helpful in proving the final Theorem 4.5.

Thus in Theorem 4.5 we prove Algorithm 4.1 achieves Heisenberg-limited
scaling, which is the main result of this work. In Sec. 5 we numerically
compare this rigorous implementation to an implementation using the matrix
pencil method, used in Ref. [16], for which we are unable to find a rigorous
proof of Heisenberg scaling. We finish the paper with a discussion, Section 6.

##  2 The classical and quantum tasks of phase estimation

One may separate quantum phase estimation into the extraction of a signal
which consists of oscillations at eigenvalue frequencies $\phi_{j}$ at a
chosen time $k$, and the processing of this signal to resolve the frequencies.
Let us first define the following:

######  Definition 2.1 (Signal or Phase Function).

Let $U\in U(2^{N})$ be an $N$-qubit unitary operator, and
$|\Psi\rangle\in\mathbb{C}^{2^{N}}$ an $N$-qubit state. We label the
eigenstates $|\phi_{j}\rangle$ of $U$ by their phase —
$U|\phi_{j}\rangle=e^{i\phi_{j}}|\phi_{j}\rangle$. We can decompose
$|\Psi\rangle$ in terms of these eigenstates,

$|\Psi\rangle=\sum_{j}a_{j}|\phi_{j}\rangle,$ (1)

and write the probability $A_{j}:=|a_{j}|^{2},\sum_{j}A_{j}=1$. We define the
phase function, – also called the signal–, $g(k)$ for $k\in\mathbb{R}$ of a
state $|\Psi\rangle$ under $U$ as

$g(k)=\sum_{j}A_{j}e^{ik\phi_{j}}.$ (2)

The spectral function $A(\phi)$ is defined as

$A(\phi)=\sum_{j}A_{j}\delta(\phi-\phi_{j}).$ (3)

Note that $\int_{0}^{2\pi}d\phi A(\phi)=1$, and $g(k)=\int_{0}^{2\pi}d\phi
e^{ik\phi}A(\phi)$; i.e. the phase function sets the Fourier coefficients of
the spectral function.

Note that one may change seamlessly between the description of a unitary $U$
and its eigenvalues and a Hermitian operator $H$ and its eigenvalues using the
transform $U=e^{iHt}$ for an appropriate choice of $t$. Note that since
$g(-k)=g^{*}(k)$ we can restrict ourselves to $k\geq 0$.

One may consider algorithms estimating $g(k)$ at integer $k\in\mathbb{Z}^{+}$,
which require the quantum circuits using controlled-$U^{k}$ in Fig. 1 with
$k\in\mathbb{Z}^{+}$. In our final Alg. 4.1 we will however use
$k\in\mathbb{R}^{+}$ (in practice $k\in\mathbb{Q}^{+}$). This generalization
can be pulled back to $k\in\mathbb{Z}^{+}$ without additional overhead in a
query-based model via a result from [21]:

######  Lemma 2.2.

[21] Fix $k\in\mathbb{R}^{+}$. Given a quantum circuit that implements a
unitary $U$ with circuit depth $t_{U}$, one may construct a quantum circuit to
implement $U^{k}$ with error $\varepsilon$ with circuit depth $\lfloor
k\rfloor
t_{U}+O\left(\frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)}t_{U}\right)$.

Crucially, the additional term is additive, and does not affect the scaling in
$k$. However, the gadget required to implement $U^{k}$ does require an
additional control bit. (In the case of performing Hamiltonian simulation
$U=e^{iHt}$ this additional control is not needed.)

The following task summarizes the quantum subroutine for phase estimation
which is to be executed with the quantum circuits in Fig. 1:

######  Definition 2.3 (Phase Function Estimation, PFE).

Let $U$ be an $N$-qubit unitary operator and $|\Psi\rangle$ an $N$-qubit
quantum state. Assume

  1. 1.

A polynomial-size (in $N$) quantum circuit implementation of $U$ (conditional
on a control qubit) and,

  2. 2.

A polynomial-size (in $N$) quantum circuit that prepares $|\Psi\rangle$.

Given a $k\in\mathbb{Z}^{+}$, error $\epsilon>0$, and confidence $0<p\leq
1\in\mathbb{R}$, PFE outputs an estimate $\tilde{g}(k)$ of the phase function
$g(k)$ of $|\Psi\rangle$ under $U$, with
$\mathbb{P}(|\tilde{g}(k)-g(k)|\leq\epsilon)\geq 1-p$ with quantum cost
$T={\rm poly}(N)M|k|$ where $M$ is the number of repetitions of both
experiments in Fig. 1 and $M=\Theta(|\ln(1-p)|\epsilon^{-2})$ via a Chernoff
bound. Lemma 2.2 shows that an identical statement holds when
$k\in\mathbb{R}^{+}$ with cost $T$ scaling as ${\rm poly}(N)M(\lfloor
k\rfloor+O(1))$.

Note that estimating the quantum cost of the subroutine in Def. 2.3 as linear
in $k$ is consistent with general no-fast forwarding statements [12] which
state that for general Hamiltonians one cannot implement $U^{t}=\exp(itH)$ in
time sub-linear in $t$. It is expected that phase function estimation is hard
to do efficiently on a classical computer as it allows one, via classical
post-processing, to sample from the eigenvalue distribution from the input
state which can be reformulated as a BQP-complete problem [4].

Figure 1: For an input state
$\left\lvert{\Psi}\right\rangle=\sum_{j}a_{j}\left\lvert{\phi_{j}}\right\rangle$
the probability for the ancilla measurement outcome to be $\pm 1$ is
$\mathbb{P}(\pm 1)=\frac{1}{2}\sum_{j}A_{j}[1\pm\cos(k\phi_{j})]$ (left
circuit measuring in the $X$-basis) and $\mathbb{P}(\pm
1)=\frac{1}{2}\sum_{j}A_{j}[1\mp\sin(k\phi_{j})]$ (right circuit measuring in
the $Y$-basis).

The quantum subroutine for PFE proceeds by executing the circuits in Fig. 1
for the given $k$. The control qubit is prepared in the
$\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ state, and is used to control $k$
applications of the unitary $U$ on the system register prepared in
$|\Psi\rangle$. The reduced density matrix of the control qubit then takes the
form

$\rho=\frac{1}{2}\left(\begin{array}[]{cc}1&g(k)\\\
g^{*}(k)&1\end{array}\right).$ (4)

The phase function $g(k)$ is extracted by state tomography of the control
qubit: one estimates the real and imaginary parts of
$g(k)=g^{(r)}(k)+ig^{(i)}(k)$ by $M$ repetitions of the circuit in Fig. 1 and
measurements of the control qubit in the X- or Y-basis respectively. We will
ignore any dependence of phase function estimation on $N$ in the rest of the
paper.

###  2.1 Phase distance on the circle

Quantum phase estimation describes a series of protocols to estimate the
eigenphases $\phi_{j}$. As these eigenvalues are defined on the circle
$[0,2\pi)$, we need a notion of distance which respects this periodicity:

######  Definition 2.4.

For $x\in\mathbb{R}$ we define the distance $|\cdot|_{T}\in[0,\pi]$ as

$\displaystyle|x|_{T}:=\min_{m\in\mathbb{Z}}(|\Delta|),$ $\displaystyle\mbox{
with }\;x=\Delta+2\pi m,$ $\displaystyle\mbox{ with }\Delta\in[-\pi,\pi).$ (5)

Clearly, the distance obeys the triangle inequality: for
$x_{1},x_{2}\in\mathbb{R}$

$|x_{1}+x_{2}|_{T}\leq|x_{1}|_{T}+|x_{2}|_{T}.$ (6)

The following Lemma addresses a technical issue in the proof of performance of
Algorithm 4.1. For integer $k\in\mathbb{Z}^{+}$, we have
$|kx|_{T}=\min_{m\in\mathbb{Z}}|k\Delta+2\pi
m|=k\min_{m\in\mathbb{Z}}|\Delta+\frac{2\pi m}{k}|_{T}$, implying Eq. (8)
below directly. However, for $k\in\mathbb{R}^{+}$ (or rational numbers
$k\in\mathbb{Q}^{+}$) we need to specify a range of $x$ for which such a
statement holds, that is:

######  Lemma 2.5.

Suppose $k\in\mathbb{R}^{+}>1$ and $\theta,\phi\in[0,2\pi)$. If

$\frac{\pi}{k}\leq\phi\leq\frac{\pi(2\lfloor k\rfloor-1)}{k},$ (7)

we have for any $\theta$

$\min_{n\in\\{0,\ldots,\lfloor
k\rfloor-1\\}}\left|\phi-\frac{\theta}{k}-\frac{2\pi
n}{k}\right|_{T}=\frac{1}{k}\left|k\phi-\theta\right|_{T}.$ (8)

###### Proof.

Let $x=\phi-\frac{\theta}{k}$, then Eq. (7) and $\theta\in[0,2\pi)$ imply that

$-\pi\leq kx\leq 2\pi\lfloor k\rfloor-\pi.$ (9)

and thus $kx=\Delta+2\pi m$ with $m\in\\{0,\ldots\lfloor k\rfloor-1\\}$ and
$\Delta\in[-\pi,\pi)$ and $|k\phi-\theta|_{T}=|\Delta|$. Hence
$x=\frac{\Delta}{k}+\frac{2\pi m}{k}$ with $m\in\\{0,\ldots\lfloor
k\rfloor-1\\}$, implying Eq. (8) where the minimum can be achieved by $m=n$. ∎

###  2.2 Limits for single-eigenvalue phase estimation

For the special case of estimating a single eigenvalue phase, the Cramér-Rao
theorem can be used to lower bound the quantum cost $T$ to learn the phase
with accuracy $\delta$, known as the Heisenberg limit. The problem of
estimating multiple phases $\phi_{j}$, in the presence of unknown
probabilities $A_{j}$, is not easily amenable to such Fisher information
analysis as it is a multi-parameter estimation problem. However, it can be
expected that the cost $T$ of this task is at least as high as that of single
phase estimation, hence it is of interest to review these bounds here. The
following theorem also proves a dense signal limit which sits in between
Heisenberg and sampling noise scaling:

######  Theorem 2.6 (The Heisenberg, Dense Signal and Sampling Limits).

The (root-mean-square) error $\delta$ of an estimator $\tilde{\phi}$ of the
eigenvalue phase $\phi$ employing the circuits in Fig. 1 on a eigenstate of U
is always lower bounded as

$\mbox{ {\rm Heisenberg Limit}}\colon\delta\geq cT^{-1}$ (10)

where $T$ is the quantum cost of implementing the circuits. If we choose to
use only quantum circuits with $k=1$, the sampling noise limit holds:

$\mbox{ {\rm Sampling Noise Limit}}\colon\delta\geq T^{-1/2}$ (11)

If we choose circuits with $k=0,1,\ldots,K$ with a fixed number of repetitions
$M$ for each circuit we are bound by a so-called ‘dense signal’ limit:

$\mbox{ {\rm Dense Signal Limit}}\colon\delta\geq cT^{-3/4}$ (12)

In these statements $c$ is some constant.

###### Proof.

Let $\tilde{\phi}$ be an estimator of $\phi$ which is inferred from the data
${\bf x}$. Here the data ${\bf x}$ is the string of outcomes of the ancilla
qubit measurements for all the experiments using the left and right circuits
in Fig. 1. We have

$\delta^{2}=\sum_{\bf x}\mathbb{P}({\bf x}|\phi)(\phi-\tilde{\phi}({\bf
x}))^{2}\geq I^{-1}(\phi),$ (13)

by the Cramér-Rao theorem [25, 26] where the Fisher information is defined as

$I(\phi)=\sum_{{\bf x}}\mathbb{P}({\bf
x}|\phi)\left\\{\frac{\partial}{\partial\phi}\ln\left[\mathbb{P}({\bf
x}|\phi)\right]\right\\}^{2}.$ (14)

Thus $I(\phi)$ limits the information we may learn about $\phi$ given a
dataset ${\bf x}$ drawn from $\mathbb{P}({\bf x}|\phi)$ and we can calculate
$I(\phi)$. Let $M_{k}^{r}$ be the number of experiments, using the circuit
with the $X$ measurement, and $M_{k}^{i}$ be the number of experiments using
the circuit with the $Y$ measurement with a certain chosen $k$. The Fisher
information for all independent experiments together is additive, i.e.
$I(\phi)=\sum_{k}[M_{k}^{r}I(\phi|k,r)+M_{i}^{r}I(\phi|k,i)]$ with
$I(\phi|k,r)$ and $I(\phi|k,i)$ the Fisher information of a single experiment
and $\sum_{k}$ is the sum over the chosen set of $k$s. For a single experiment
we can calculate, using Eq. (14) and the probability for the output bit given
in Fig. 1, that $I(\phi|k,r)=I(\phi|k,i)=k^{2}$ and thus

$I(\phi)=\sum_{k}k^{2}(M_{k}^{r}+M_{k}^{i}).$ (15)

At the same time the total quantum cost of all experiments is

$T=\sum_{k}k(M_{k}^{r}+M_{k}^{i}).$ (16)

The key insight here is that the relative dependence on $k$ is different
between $T$ and $I$ and this implies that the trend of the number of
experimental runs $M_{k}^{r/i}$ as a function of $k$ will affect the maximum
rate of estimation. If we choose only $k=1$ we see that
$\delta^{2}\geq\frac{1}{T}$ which is the sampling noise limit.

The biggest value for $I(\phi)$ for a given $T$ is obtained when we choose a
single largest possible $k=K$ so that $T=KM_{K}$ and
$I(\phi)=K^{2}M_{K}=T^{2}/M_{K}$ with $M_{K}=M_{K}^{r}+M_{K}^{i}$ . This
implies a Heisenberg limit, i.e. $\delta\geq cT^{-1}$ where $c$ is some
constant depending on $M_{K}$. If we however make the ‘dense signal’ choice,
that is, $M^{r}_{k}=M^{i}_{k}=M$ for $k=1,2\ldots K$, then

$I(\phi)=\frac{M}{3}K(K+1)(1+2K),$ (17)

while the total quantum cost is $T=MK(K+1)$. To leading order in $K$ we then
have $\delta\geq
I(\phi)^{-\frac{1}{2}}=\sqrt{\frac{3}{2}}M^{\frac{1}{4}}T^{-\frac{3}{4}}$,
which is the dense signal limit.

∎

Remark: We note that a randomized version of the dense signal choice can
potentially scale in near-Heisenberg-limited fashion. In this method, one
would draw $k$ at random from $1,\ldots,K$ and repeat this $S$ times to
generate random variables $k_{1},\ldots,k_{S}$ and repeat experiments with
fixed $M$ for each such $k_{i}$. With the right choice of $S\times M={\rm
polylog}(K)$, one can argue, using the Cramer-Rao lower bound analysis above,
that the expected error $\mathbb{E}(\delta)\geq\frac{{\rm
polylog}(\mathbb{E}(T))}{\mathbb{E}(T)}$ where $\mathbb{E}(T)$ is the expected
quantum cost. The algorithm in [14] uses such strategy with $M=1$. Clearly,
the sampling noise limit can be achieved by choosing $k=1$ in the circuits of
Fig. 1. However, one can ask whether the dense signal limit or the Heisenberg
limit can also be achieved, in particular when we demand that the classical
post-processing is computationally efficient, meaning that this processing is
polynomial in the quantum cost $T$. For the dense signal limit one needs a
classical method to process the estimates of $g(k)$ at $k=1,\ldots,K$ to
estimate $\tilde{\phi}$. Using perturbation theory in the noise, the matrix
pencil method has been claimed to achieve this for a single eigenvalue [27].

Achieving the Heisenberg limit is non-trivial due to phase aliasing: the phase
function $g(k)$ obtained by the experiments using $U^{k}$ remains invariant if
the phase $\phi$ is shifted by $\frac{2\pi}{k}$. This implies that a strategy
of estimating $\phi$ from a single point $g(k)$ at large $k$ will fail unless
$\phi$ is already known to sit within a window of width $\frac{2\pi}{k}$. This
issue is circumvented by sampling $g(k)$ at multiple orders $k=2^{d}$ to
‘gradually zoom in’ on $\phi$. To get Heisenberg scaling, one lets the number
of samples $M$ and thus the confidence, to depend on the order, so that the
most significant bits of $\phi$ are determined with the highest confidence.
Methods for doing this were first introduced in Ref. [1], and improved in Ref.
[2] for the purpose of gate calibration. Here we state the result:

######  Algorithm 2.7 (Heisenberg Algorithm For Single Eigenvalue Phase [1,
2]).

Given an targeted error $\delta>0$, and numbers
$\alpha,\gamma\in\mathbb{Z}^{+}$. The Heisenberg algorithm which outputs an
estimate $\tilde{\phi}$ for $\phi$ proceeds as follows:

  1. 1.

Fix $d_{f}=\lceil\log_{2}(1/\delta)\rceil$.

  2. 2.

For $d=0,1\ldots,d_{f}$:

    1. (a)

Use the PFE subroutine, Def. 2.3, circuits in Fig. 1, to obtain an estimate
$\tilde{g}(k)$ of $g(k)$ for $k=2^{d}$ using $M_{d}=\alpha+\gamma(d_{f}+1-d)$
repetitions of both experiments.

    2. (b)

Compute $\tilde{\theta}^{(d)}=\mathrm{Arg}[\tilde{g}(2^{d})]\in[0,2\pi)$

    3. (c)

If $d=0$, set $\tilde{\phi}^{(0)}=\tilde{\theta}^{(0)}$.

    4. (d)

Else, set $\tilde{\phi}^{(d)}$ to be the unique value in the interval
$[\tilde{\phi}^{(d-1)}-\frac{\pi}{2^{d}},\tilde{\phi}^{(d-1)}+\frac{\pi}{2^{d}})$
(with periodic boundaries) such that

$2^{d}\tilde{\phi}^{(d)}=\tilde{\theta}^{(d)}\mod 2\pi.$ (18)

  3. 3.

Return $\tilde{\phi}=\tilde{\phi}^{(d_{f})}$ as an estimate for $\phi$.

It was proved in [2] that for some choices of $\alpha$ and $\gamma$ the root-
mean-square error $\delta$ on the final estimate $\phi^{(d_{f})}$ is at most
$cT^{-1}$ for a constant $c$ and total cost $T=2\sum_{d=0}^{d_{f}}2^{d}M_{d}$,
thus reaching the Heisenberg limit.

##  3 Defining the task of multiple-eigenvalue phase estimation

In this section we define the goal of estimating multiple eigenvalue phases of
some unitary $U$. When the input state $\left\lvert{\Psi}\right\rangle$ is
supported on multiple eigenstates, choosing a single $k$ does not suffice,
simply since knowing Eq. (2) at a single point $k$ does not give a unique
solution $\\{A_{j},\phi_{j}\\}$ [16]. A simple way to circumvent this problem
is thus to estimate ‘densely’: estimating $g(k)$ for all integers $0\leq k\leq
K$ would allow us to fit up to $O(K)$ $(\phi_{j},A_{j})$ pairs. However, this
does not saturate the Heisenberg limit as shown in Theorem 2.6, hence we need
to come up with a different method.

Separate from this, the full eigenspectrum of an arbitrary $N$-qubit unitary
$U$ has up to $2^{N}$ unique values, making it impossible to describe in
polynomial time in $N$. In addition, the spectral content of the input state
$\left\lvert{\Psi}\right\rangle$ could be very dense, with many eigenvalues
clustered together instead of separated by gaps, and the probability for these
eigenvalues, $A_{j}$, could be sharply concentrated or uniformly spread. To
deal with general input states, Ref. [3] thus formulated the quantum
eigenvalue estimation problem (QEEP): instead of estimating individual phases,
the focus is on estimating the spectral function $A(\phi)$ in Eq. (3) with
some resolution.

######  Definition 3.1 (Quantum Eigenvalue Estimation Problem, QEEP).

Let $A(\phi)$ be the spectral function defined in Eq. (3) for a unitary $U$
and $\left\lvert{\Psi}\right\rangle$. Given is an error parameter
$\epsilon>0$, a confidence parameter $1\geq p>0$, and a set of non-negative
(approximate indicator) functions $f^{l}(\phi)$ for $\phi\in[0,2\pi)$ for
$l=0,\ldots,L-1$, $L=\lceil\frac{2\pi}{\epsilon}\rceil$, where $f^{l}(.)$ has
support on only the interval bin

$\mathcal{V}_{l}=[(l-1)\epsilon,(l+1)\epsilon]_{T},$ (19)

and $f^{l}(\phi)+f^{(l-1)}(\phi)=1$ for all
$\phi\in\mathcal{V}_{l}\cap\mathcal{V}_{l-1}$. Assuming access to the PFE
subroutine, Def. 2.3, the goal is to output an approximation $\tilde{b}_{l}$
for $l=0,\ldots,L-1$ to the integral

$b_{l}=\int_{0}^{2\pi}d\phi A(\phi)f^{l}(\phi),$ (20)

which satisfies

$\sum_{l=0}^{L-1}|\tilde{b}_{l}-b_{l}|\leq\epsilon,$ (21)

with probability at least $p$.

Note that the bins $\mathcal{V}_{l}$ have width $2\epsilon$ and overlap on a
region of width $\epsilon$ and $\sum_{l=0}^{L-1}b_{l}=1$.

###  3.1 Methods of dense signal phase estimation

For estimating multiple phases there are two classical processing methods to
infer phases using estimates of $g(k)$ with $k=0,1,\ldots,K$. These are the
matrix pencil method first studied for QPE in Ref. [16], see Algorithm 5.1,
and the ‘time series analysis’ proposed in Ref. [3] to solve the QEEP problem
in Definition 3.1. We can examine the scaling of the total quantum cost in
light of the dense signal limit of Theorem 2.6. Ref. [3] proved the following

######  Theorem 3.2 (QEEP Algorithm [3]).

One can solve the QEEP problem in Definition 3.1 with $\\{f^{l}(.)\\}$ a set
of ‘bump’ functions

$\displaystyle f^{l}(\phi)=$
$\displaystyle\frac{2a}{\epsilon}\int_{l\epsilon-\frac{\epsilon}{2}}^{l\epsilon+\frac{\epsilon}{2}}d\phi^{\prime}$
$\displaystyle\times\exp\left\\{-\left[1-\frac{4}{\epsilon^{2}}(\phi-\phi^{\prime})^{2}\right]^{-1}\right\\},$
(22)

with normalization constant $a\approx 2.252$, using PFE in Def. 2.3 with
$k=0,1,\ldots,K$ with $K=O\left(\epsilon^{-1}\ln^{2}(\epsilon^{-1})\right)$
and $M=\tilde{O}(|\ln(1-p)|\epsilon^{-4})$ for each $k=1,2\ldots,K$. The total
quantum cost for $U$ is then bounded as
$T=O(MK^{2})=\tilde{O}(|\ln(1-p)|\epsilon^{-6})$.

We note that the approximate indicator functions $f^{l}(.)$ in Eq. (22) are
designed to have a quickly decaying Fourier series, which is required to
achieve polynomial-time scaling. We also refer the reader to Ref. [20], which
has extended the result by relaxing the requirement that
$f^{l}(\phi)+f^{(l-1)}(\phi)=1$ on the interval
$\mathcal{V}_{l}\cap\mathcal{V}_{l-1}$. The idea of this QEEP algorithm is as
follows. Since $b_{l}=\sum_{j}A_{j}f^{l}(\phi_{j})$, using Eq. (3),
periodically extending $f^{l}(\phi)$ beyond $[0,2\pi)$ and Fourier decomposing
gives
$b_{l}=\sum_{j}A_{j}\sum_{k\in\mathbb{Z}}e^{ik\phi_{j}}\tilde{f}^{l}(k)=\sum_{k\in\mathbb{Z}}g(k)\tilde{f}^{l}(k)$.
At the same time, the fact that $f^{l}(.)$ is an indicator function ensures
that $b_{l}\approx\sum_{\phi_{j}\in\mathcal{V}_{l}}A_{l}$. Thus knowledge of
$g(k)$ and the Fourier coefficients $\tilde{f}^{l}(k)$ for a range of $k$
allows one to estimate the weights $b_{l}$. The requirement to estimate the
spectral function to within a $1$-norm $\epsilon$, Eq. (21), is very
stringent, hence the scaling of $T$ with error $\epsilon$ is quite costly,
$T=O(\epsilon^{-6})$. It is possible that one can improve the scaling by re-
examining the analysis in [3].

Given estimates of $g(k)$ for $k=0,\ldots,K$ (dense signal), one can ask for
error bounds of processing this data to estimate phases $\phi_{j}$ and their
associated probabilities $A_{j}$ using the matrix pencil technique. Bounding
the error of the matrix pencil method naturally requires a minimal gap
$\Delta$ between these phases, $\Delta=\min_{i\neq j}|\phi_{i}-\phi_{j}|_{T}$.
Namely, a proven necessary condition to get convergent estimates for
$\phi_{j}$ and $A_{j}$ is that $K\geq c/\Delta$ [28] for some constant $c$. We
do however not want to impose the presence of a gap in this work as this is
unnecessarily restrictive.

Here is our precise definition of the problem to be solved and achieving the
Heisenberg limit for multiple eigenvalues requires solving this problem with a
total quantum cost $T=O(\delta^{-1})$:

######  Definition 3.3 (Multiple eigenvalue estimation problem).

Fix an error bound $\delta>0$, a probability bound $A$. For a unitary $U$ and
state $|\Psi\rangle$, we assume that $A_{j}>A$ for exactly $n_{\phi}$ phases
$\phi_{j}$ and $A_{j}=0$ for all other phases so that $n_{\phi}\leq A^{-1}$.
Let $g(k)=\sum_{j}A_{j}e^{ik\phi_{j}}$ be the phase function in Def. 2.1 and
assume access to the PFE quantum subroutine in Def. 2.3 for any
$k\in\mathbb{R}^{+}$, generating data ${\bf x}$. The task is to output a set
$\\{\tilde{\phi}_{l}\\}$ of $n_{\phi}$ or fewer estimates of the phases
$\\{\phi_{j}\\}$ such that, if we take the closest estimate
$\tilde{\phi}_{j}^{(\mathrm{closest})}({\bf x})$ of each phase $\phi_{j}$
given the data ${\bf x}$,

$\tilde{\phi}_{j}^{(\mathrm{closest})}({\bf
x})=\operatorname*{arg\,min}_{\tilde{\phi}_{l}}(|\tilde{\phi}_{l}({\bf
x})-\phi_{j}|_{T}),$ (23)

the accuracy error

$\delta_{j}=\sqrt{\sum_{\mathbf{x}}\mathbb{P}\left(\mathbf{x}|\\{\phi_{l},A_{l}\\}\right)\left|\tilde{\phi}_{j}^{(\mathrm{closest})}(\mathbf{x})-\phi_{j}\right|_{T}^{2}},$
(24)

is bounded by $\delta_{j}\leq\delta$ for all $j=1,\ldots,n_{\phi}$.

Remarks: Def. 3.3 allows us the freedom to assign a single estimate to
multiple phases when calculating the final mean-squared-error. The classical
computational efficiency of the algorithm which solves this problem should be
polynomial in $n_{\phi}$ and $T$.

We intend to solve the problem considered in Def. 3.3 by identifying phases
from the output of the QEEP problem in Def. 3.1 that have probability above
the probability bound $A$. To extend this to the Heisenberg limit we will use
this QEEP algorithm in Theorem 3.2 for a unitary $V=U^{k}$ with exponentially
growing $k$, thus $V$ has eigenphases

$\theta_{j}=k\phi_{j}\mod 2\pi.$ (25)

Hence in Algorithm 3.4 and Lemma 3.5 below we refer to estimating eigenvalue
phases $\theta$ of such $V=U^{k}$.

Algorithm 3.4 is an extension to the QEEP algorithm allowing us to output
maximally $n_{\phi}$ phases by solving the QEEP problem. Lemma 3.5 then states
what error is obtained on the estimates for $\theta$ when the QEEP algorithm
3.2 succeeds with probability $p$.

Let us first motivate Algorithm 3.4. One may identify phases with sufficient
probability in the output of the QEEP algorithm as the bins $l$ with
$b_{l}>b_{\mathrm{cutoff}}$ with $b_{l}$ in Eq. (20). Then to convert this
into an estimate of a phase $\theta_{j}$, for each such bin above cut-off we
could output the estimate $\epsilon l$, i.e. in the middle of the
corresponding bin $\mathcal{V}_{l}$.

We calculate appropriate values of $b_{\mathrm{cutoff}}$ and the QEEP error
$\epsilon$ to guarantee that we output an estimate for each $\theta_{j}$ with
$A_{j}>A$ with a provable confidence, and to guarantee no estimate in the
absence of any $\theta_{j}$. Def. 3.1 states that when there exists such a
$\theta_{j}\in\mathcal{V}_{l}\cap\mathcal{V}_{l-1}$, we are guaranteed with
confidence $p$ that $\tilde{b}_{l}+\tilde{b}_{l-1}+\epsilon>A$ (as
$b_{l}+b_{l-1}>A$). To guarantee that at least one of $\tilde{b}_{l}$ or
$\tilde{b}_{l-1}$ is larger than $b_{\mathrm{cutoff}}$ with the same
confidence, we thus require $b_{\mathrm{cutoff}}\leq(A-\epsilon)/2$.

Similarly, Def. 3.1 states that when there exists no
$\theta_{j}\in\mathcal{V}_{l}$ with $A_{j}>0$, we are guaranteed with
confidence $p$ that $b_{l}<\epsilon$. To prevent a spurious estimate in this
case, we require $b_{\mathrm{cutoff}}\geq\epsilon$. Solving this to maximise
$\epsilon$ (which minimizes the cost of the QEEP routine) yields
$b_{\mathrm{cutoff}}=\epsilon=A/3$ which is what we use in Alg. 3.4.

A further small complication exists in solving the problem posed in Def. 3.3:
we require that one outputs at most $n_{\phi}$ phases. This will be satisfied
if we can guarantee at most one phase estimate per $\theta_{j}$ with
$A_{j}>A$, as we know there are $n_{\theta}=n_{\phi}$ such estimates. As the
bins $\mathcal{V}_{l}$ (Def. (19)) overlap, a phase $\theta_{j}$ may
participate in up to two neighbouring bins — corresponding to amplitudes
$b_{l},b_{l+1}>A/3$. To ensure that this does not result in two estimates
being generated, one could take all contiguous sets
$b_{l_{1}},b_{l_{1}+1},\ldots,b_{l_{2}}>A/3$ and prune away every second index
$l$. In doing this we need to respect the periodicity of the
$\mathcal{V}_{l}$: $\mathcal{V}_{L-1}\cap\mathcal{V}_{0}\neq\emptyset$, so
these contiguous sets may wrap around the circle. Pruning every second index
when starting from the middle of one of these contiguous sets may result in
two neighbouring $l,l+1$ being removed, which is not what we desire. Instead,
in the following pseudocode, after generating the set of all $l$ with
sufficient $b_{l}$, we find the first gap (in $l$) between these regions
(corresponding to the first $b_{l}<A/3$). We then iterate (from this point
$l_{\min}$ to $L-1$ and then from $0$ to $l_{\min}$) over the $b_{l}$, and
remove each $l$ from our set if $b_{l-1}>A/3$ and $l-1$ was not itself
removed.

######  Algorithm 3.4 (Conservative QEEP Eigenvalue Extraction).

Fix a probability bound $A$, an error bound $0<\epsilon<\frac{A}{3}$, and a
confidence bound $0<p<1$. Assume access to a QEEP Algorithm 3.2 for a unitary
$V$. Define $+_{L}$ and $-_{L}$ as addition and subtraction modulo $L$ (so
$0-_{L}1=L-1$). The algorithm proceeds as follows:

  1. 1.

Use the QEEP subroutine with error $\epsilon$, Alg. 3.2, and confidence $p$ to
generate an estimate $\tilde{b}_{l}$ for $b_{l}$ as defined in Eq. (20)

  2. 2.

Construct the set

$S=\\{l\in\\{0,\ldots,L-1\\}|\tilde{b}_{l}\geq A/3\\}.$ (26)

  3. 3.

Find the smallest $l\in\\{0,\ldots,L\\}$ with $l\not\in S$ and call it
$l_{\min}$.

  4. 4.

For $l^{\prime}=l_{min},\ldots,l_{min}+_{L}(L-1)$:

    1. (a)

if $l^{\prime}\in S$ and $l^{\prime}-_{L}1\in S$, remove $l^{\prime}$ from
$S$.

  5. 5.

Return the set $\\{\tilde{\theta}_{l}=l\epsilon\\}_{l\in S}$ as a set of
estimates of eigenphases of $V$.

The following Lemma then relates the performance of Alg. 3.4 to the
performance of the QEEP subroutine. Here $\\{\theta_{j}\\}$ is an ordered list
of $n_{\phi}$ phases, while $\\{\tilde{\theta}_{l}\\}$ is an ordered list
(that we will show contains at most $n_{\phi}$ phases).

######  Lemma 3.5.

Fix a confidence bound $0<p<1$, a probability bound $A$, a number of phases
$n_{\theta}\geq A^{-1}$, and an error bound $0<\epsilon<\frac{A}{3}$. Let
$g(k)=\sum_{j}A_{j}e^{ik\theta_{j}}$ be the phase function for a unitary $V$.
Let $\\{\tilde{\theta}_{l}\\}$ be a set of estimates of $\\{\theta_{j}\\}$
generated by Alg. 3.4 with error $\epsilon$ and confidence bound $p$. With
probability at least $p$, the following statements are true:

  1. 1.

For each phase $\theta_{j}$ with $A_{j}>0$, there exists at least one estimate
$\tilde{\theta}_{l}$ such that

$|\theta_{j}-\tilde{\theta}_{l}|_{T}\leq 2\epsilon.$

  2. 2.

For each estimate $\tilde{\theta}_{l}$ there exists at least one phase
$\theta_{j}$ with $A_{j}>0$ such that

$|\theta_{j}-\tilde{\theta}_{l}|_{T}\leq 2\epsilon.$

  3. 3.

The number of estimates $|\\{\tilde{\theta}_{l}\\}|\leq n_{\theta}$.

###### Proof.

Our proof follows by showing that the output from Alg. 3.4 satisfies these
statements whenever Eq. (21) holds. This yields our confidence bound as Eq.
(21) holds with probability $p$ by Def. 3.1.

To see that statement 2 holds when Eq. (21) is satisfied, note that if
$\mathcal{V}_{l}$ contains no phases $\theta_{j}$ with $A_{j}>A$ and Eq. (21)
is satisfied, $\tilde{b}_{l}<\epsilon<A/3$, and $l$ will not be added to the
set $S$ in Alg. 3.4. This implies that when Eq. (21) holds, if $l\in S$ there
exists some $\theta_{j}\in\mathcal{V}_{l}$ with $A_{j}>A$, in which case
$|\theta_{j}-\tilde{\theta}_{l}|_{T}=|\theta_{j}-\epsilon l|_{T}\leq
2\epsilon$.

To see that statement 3 holds when Eq. (21) is satisfied, note that
$\theta_{j}\in\mathcal{V}_{l}\cap\mathcal{V}_{l+_{L}1}$ for exactly one $l$
(and $\theta_{j}\notin\mathcal{V}_{m}$ for $m\neq l,l+_{L}1$). Then, when Eq.
(21) holds, $\tilde{b}_{l}+\tilde{b}_{l+_{L}1}>2A/3$, so
$\max(\tilde{b}_{l},\tilde{b}_{l+_{L}1})>A/3$, and either $l$, $l+_{L}1$ or
$l$ and $l+_{L}1$ will be added to the set $S$ during step 2 of Alg. 3.4 for
each phase $\theta_{j}$. Then, step 4 of Alg. 3.4 will remove $l+_{L}1$ if $l$
remains in $S$, so each phase $\theta_{j}$ can contribute to only one final
estimate $\tilde{\theta}_{l}=l\epsilon$, and the number of estimates is
bounded from above by the number of phases.

To see that statement 1 holds, we use the point in the previous paragraph
that, when Eq. (21) is satisfied, each phase $\theta_{j}$ with $A_{j}>A$ adds
either $l$, $l+_{L}1$, or $l$ and $l+_{L}1$ to the set $S$ during step 2. of
Alg. 3.4. Then in step 4. of Alg. 3.4, $l$ is removed from $S$ only if
$l-_{L}1$ remains in $S$, and $l+_{L}1$ is removed from $S$ only if $l$
remains in $S$. This implies that the distance from $\theta_{j}$ to an
estimate is bounded by

$\max_{\theta_{j}\in\mathcal{V}_{l}\cap\mathcal{V}_{l+_{L}1}}\hskip
5.69046pt\max_{l^{\prime}=l-_{L}1,l,l+_{L}1}|\theta_{j}-\epsilon
l^{\prime}|_{T}=\\\ =\epsilon(l+_{L}1)-\epsilon(l-_{L}1)=2\epsilon,$ (27)

as required. ∎

##  4 Multiple eigenvalues: multi-order estimation and the phase matching
problem

To achieve Heisenberg-limited scaling for multiple phases, we combine the
dense signal algorithms which can resolve multiple phases in the previous
section with the single-phase Heisenberg limited algorithm, Algorithm 2.7,
which achieves the correct scaling.

A natural way to achieve such combination is to estimate phases
$\\{\theta^{(d)}_{j}\\}$ of $V=U^{2^{d}}$ for multiple orders
$d=0,1\ldots,d_{f}$ via a dense signal method (e.g. Alg. 3.4), and then
combine them in the same manner as in Algorithm 2.7.

If we would manage to get an estimate of $\theta_{j}^{(d)}=\phi_{j}2^{d}$ at
each order $d$ with error $\epsilon$ and be able to combine these estimates in
an unequivocal manner, then reaching the Heisenberg limit for multiple phases
may be feasible. Note that the error in the final $d_{f}^{\rm th}$ estimate in
this case would be $\delta\sim\epsilon/2^{d_{f}}$ with $\epsilon$ in Algorithm
3.4. One could thus achieve arbitrarily small $\delta$ for fixed $\epsilon$ by
making $d_{f}$ arbitrarily large. This allows us to use (possibly non-optimal)
routines such as the QEEP algorithm since the scaling with $\epsilon$ does not
propagate to a scaling in $\delta$ for a sufficiently small $\epsilon$.

However, the combination of phase information at different orders in the case
of multiple phases may not be feasible when we use $V=U^{2^{d}}$ for
increasing $d$ as in Algorithm 2.7. The reason is that if we have multiple
phases, the previous order estimates provide sets of ‘ballpark’ intervals and
it may not clear or unambiguous which interval to choose in order to convert a
new estimate $\tilde{\theta}^{(d)}_{j}$ to an updated $\tilde{\phi}^{(d)}_{j}$
(as in step 2d of Algorithm 2.7). We would like the choice of the next order
to be such that this ‘matching with a previous estimate’ can be done
unambiguously.

For this, we will estimate the multiple eigenphases of $V=U^{k_{d}}$ for
$k_{d}=\prod_{d^{\prime}=1}^{d}\kappa_{d}$ with $\kappa_{d}\geq 2$ a, possibly
non-integer, multiplier. We can make use of non-integer $\kappa_{d}$
efficiently by virtue of Lem. 2.2. For this algorithm to have a means of
associating each $d$th order estimate $\tilde{\theta}^{(d)}_{j}$ with a single
previous-order estimates $\tilde{\phi}^{(d-1)}_{j}$, we use an adaptive
strategy for choosing the next multiplier $\kappa_{d}$ in Alg. 4.1. That is,
the algorithm will determine a $\kappa_{d}$ based on the estimates
$\tilde{\phi}_{j}^{(d-1)}$ from the previous round. Although this scheme
requires some classical processing of the experimental data before the
experiment is finished, it is not an adaptive scheme in the same sense as
iterative QPE [11], as we do not require feedback within the coherent lifetime
of a single experiment.

The generalization from using $U^{2^{d}}$ to $U^{k_{d}}$ for
$k_{d}\in\mathbb{R}^{+}$ presents one small additional complication. In order
to prove bounds on the estimation at each order we will require invoking Lemma
2.5. However, this requires that our phases $\phi_{j}$ satisfy Eq. (7) (unless
$k_{d}\in\mathbb{Z}^{+}$ ). If a phase $\phi_{j}$ does not satisfy Eq. (7),
one can construct a situation where two corresponding estimates
$\tilde{\phi}_{j}^{(d)}$ are found on either side of the branch cut at $2\pi$,
and where we cannot guarantee that our algorithm would choose the ‘correct’
one (without knowledge of the hidden $\phi_{j}$). To solve this issue, we note
that one may shift the phases of $U$ by a constant $\chi$ by performing phase
estimation on $Ue^{-i\chi}$ instead of $U$. This need not even be done on the
quantum device, as one simply multiplies estimates of $g(k)$ by $e^{-ik\chi}$.
As we assume the existence of only $n_{\phi}$ phases, we can always find some
$Ue^{i\chi}$ with phases in some window $[\phi_{\min},\phi_{\max})$ with
$\phi_{\min}\geq\frac{\pi}{k}$, $\phi_{\max}\leq\frac{\pi(2\lfloor
k\rfloor-1)}{k}$ when $k\geq 3n_{\phi}$. This will allow us to invoke Lem. 2.5
to match estimates of eigenphases of $Ue^{i\chi}$ and estimates of eigenphases
of $(Ue^{i\chi})^{k}$ as we require. We also note that the above issue can be
circumvented when $U=e^{iHt}$ by a suitable choice of $t$.

###  4.1 Heisenberg-limited algorithm for multiple phases

We now describe our Heisenberg-limited phase estimation algorithm. This
algorithm targets a final error $\delta=O(\delta_{c})$, where $\delta_{c}$ is
a fixed input to the algorithm itself (We will calculate the constant of
proportionality in the proof of Theorem 4.5). The Heisenberg limit will be
achieved by making this $\delta_{c}$ smaller while keeping the error
$\epsilon$ of the phase extraction subroutine, Alg. 3.4, constant.

######  Algorithm 4.1.

[Adaptive multi-order phase estimation algorithm] We assume access to the
conservative QEEP eigenvalue extraction algorithm, Alg. 3.4 for a unitary
$V=U^{k}$ (for arbitrary $k\in\mathbb{R}^{+}$), and an initial state
$|\Psi\rangle$. Fix a final error $\delta_{c}$, a probability bound $A$, a
number of phases $n_{\phi}\leq A^{-1}$, and error parameters $\epsilon_{0}$
and $\epsilon$ bounded as

$\epsilon_{0}\leq\epsilon_{\mathrm{crit},0}\equiv\frac{2\pi}{300n_{\phi}^{4}}$
(28)

and

$\epsilon\leq\epsilon_{{\rm crit}}\equiv\frac{2\pi}{300n_{\phi}^{2}}.$ (29)

Let the confidence parameter $p_{d}$ be

$p_{d}=1-e^{-\alpha}\left(\frac{k_{d}\delta_{c}}{\pi}\right)^{\gamma},$ (30)

given some real numbers $\alpha>0$, $\gamma>2$ and $k_{d}$ to be chosen below.
The algorithm proceeds as follows:

  1. 1.

Let $d=0$ and $k_{d=0}=1$. Use Alg. 3.4 to find a set of first estimates
$\\{\tilde{\phi}^{(0)}_{j}\\}$ of eigenvalues of $U$ with error parameter
$\epsilon_{0}$ in Eq. (28), probability bound $A$, and confidence $p_{d=0}$ in
Eq. (30). If this set is empty or has more than $n_{\phi}$ elements, return
$\\{0\\}$ (this is a failure mode).

  2. 2.

Find the point $\zeta\in[0,2\pi)$ defined by

$\zeta=\operatorname*{arg\,max}_{\zeta^{\prime}\in[0,2\pi]}\;\min_{j}\left|\tilde{\phi}^{(0)}_{j}-\zeta^{\prime}\right|_{T},$
(31)

i.e. $\zeta$ is the midway point in the largest gap between the phase
estimates $\tilde{\phi}^{(0)}_{j}$. Let

$d_{\zeta}=\min_{j}\left|\tilde{\phi}^{(0)}_{j}-\zeta\right|_{T},$ (32)

i.e. $d_{\zeta}$ is half the size of the largest gap. Shift the unitary
$U\rightarrow Ue^{-i(\zeta+\frac{1}{2}d_{\zeta}-8\epsilon_{0})}$,
$\tilde{\phi}_{j}^{(0)}\rightarrow\tilde{\phi}^{(0)}_{j}-\zeta-
d_{\zeta}/2+8\epsilon_{0}\mod 2\pi$.

  3. 3.

Choose $\kappa_{1}=k_{1}$ with $k_{1}\in[3n_{\phi},3n_{\phi}+1]$ such that for
all $\tilde{\phi}^{(0)}_{j}\neq\tilde{\phi}^{(0)}_{l}$, either

$|\tilde{\phi}_{j}^{(0)}k_{1}-\tilde{\phi}_{l}^{(0)}k_{1}|_{T}>4\epsilon_{0}(1+k_{1}).\leavevmode\nobreak\
$ (33)

or

$|\tilde{\phi}_{j}^{(0)}-\tilde{\phi}_{l}^{(0)}|_{T}<\frac{\pi}{k_{1}}.\leavevmode\nobreak\
$ (34)

  4. 4.

While $k_{d}<\frac{2\epsilon}{\delta_{c}}$ :

    1. (a)

Set $d\rightarrow d+1$.

    2. (b)

Use Alg. 3.4 to find a set of estimates $\\{\tilde{\theta}^{(d)}_{l}\\}$ of
eigenvalues of $V=U^{k_{d}}$ with error parameter $\epsilon$ in Eq. (29),
probability bound $A$, and confidence $p_{d}$ in Eq. (30).

    3. (c)

If there exists some $\tilde{\phi}_{j}^{(d-1)}$ such that

$\min_{l}|k_{d}\tilde{\phi}_{j}^{(d-1)}-\tilde{\theta}_{l}^{(d)}|_{T}>2\epsilon(1+\kappa_{d}),$
(35)

or there exists some $\tilde{\theta}_{l}^{(d)}$ such that

$\min_{j}|k_{d}\tilde{\phi}_{j}^{(d-1)}-\tilde{\theta}_{l}^{(d)}|_{T}>2\epsilon(1+\kappa_{d}),$
(36)

or the number of estimates $|\\{\tilde{\theta}_{l}^{(d)}\\}|>n_{\phi}$, return
$\\{\tilde{\phi}_{j}^{(d-1)}+\zeta+d_{\zeta}/2-8\epsilon_{0}\mod 2\pi\\}$.
This is a failure mode.

    4. (d)

If not, for each $\tilde{\theta}_{l}^{(d)}$, find the estimate
$\tilde{\phi}_{j}^{(d-1)}$ and an integer $n\in[0,k_{d})$ which minimizes

$|\tilde{\phi}_{j}^{(d-1)}-(\tilde{\theta}_{l}^{(d)}+2\pi n)/k_{d}|_{T},$ (37)

and set $\\{\tilde{\phi}_{l}^{(d)}\\}=\\{(\tilde{\theta}^{(d)}_{l}+2\pi
n)/k_{d}\\}$.

    5. (e)

If any $\tilde{\phi}_{j}^{(d)}\in[0,\frac{\pi}{k_{d}})\cup(\frac{\pi(2\lfloor
k_{d}\rfloor-1)}{k_{d}},2\pi]$, return
$\\{\tilde{\phi}_{j}^{(d-1)}+\zeta+d_{\zeta}/2-8\epsilon_{0}\mod 2\pi\\}$.
This is a failure mode.

    6. (f)

Choose the multiplier $\kappa_{d+1}\in[2,3]$ such that for all
$\tilde{\phi}^{(d)}_{j}\neq\tilde{\phi}^{(d)}_{l}$, either

$\displaystyle|\tilde{\phi}_{j}^{(d)}k_{d}\kappa_{d+1}-\tilde{\phi}_{l}^{(d)}k_{d}\kappa_{d+1}|_{T}$
$\displaystyle>4\epsilon(1+\kappa_{d+1}).\leavevmode\nobreak\ $ (38)

or

$|\tilde{\phi}_{j}^{(d)}-\tilde{\phi}_{l}^{(d)}|_{T}<\frac{\pi-2\epsilon(1+\kappa_{d+1})}{k_{d}\kappa_{d+1}},\leavevmode\nobreak\
$ (39)

and set $k_{d+1}=k_{d}\kappa_{d+1}$.

  5. 5.

Return $\\{\tilde{\phi}_{j}^{(d)}+\zeta+d_{\zeta}/2-8\epsilon_{0}\mod
2\pi\\}$.

In the rest of this section, we prove that Alg. 4.1 can achieve the Heisenberg
limit. We note that Alg. 3.4 which is invoked in this algorithm assumes that
the error parameter $\epsilon$ or $\epsilon_{0}$ is at most $A/3\leq
1/(3n_{\phi})$ which is fulfilled by Eq. (28) and Eq. (29).

We note that the first use of the QEEP subroutine possibly needs to use a
smaller error parameter, called $\epsilon_{0}$, bounded by Eq. (28) than
subsequent uses where $\epsilon$ needs to be only bounded by Eq. (29). This
relates to a technical issue, namely $k_{1}\geq 3n_{\phi}$ in order for Lemma
B.1 in Appendix B and thus Lemma 2.5 to apply. For later rounds $k_{d}\geq
3n_{\phi}$ automatically, not putting a constraint on the multiplier
$\kappa_{d}$, and indirectly allowing to relax the region of valid choices for
$\epsilon$. The first step in proving the performance of Algorithm 4.1 is to
show that the multipliers can be chosen in the first (step 3 in Alg. 4.1), and
subsequent rounds (step 4f in Alg. 4.1), which obey the desired conditions.
This is accomplished by the following Lemma which is proved in Appendix A.

######  Lemma 4.2.

Let $\\{\tilde{\phi}_{j}^{(0)}\\}\in[0,2\pi)$ be a set of at most $n_{\phi}$
phases. Assuming Eq. (28), for a randomly chosen
$k_{1}\in[3n_{\phi},3n_{\phi}+1]$ with probability at least $1/2$, either Eq.
(33) or Eq. (34) holds for all
$\tilde{\phi}_{j}^{(0)}\neq\tilde{\phi}_{l}^{(0)}$. Fix a $k_{d}$. Let
$\\{\tilde{\phi}^{(d)}_{j}\\}\in[0,2\pi)$ be a set of at most $n_{\phi}$
phases. Assuming Eq. (29), for a randomly chosen $\kappa_{d+1}\in[2,3]$ with
probability at least $3/4$, either Eq. (38) or Eq. (39) holds for all
$\tilde{\phi}_{l}^{(d)}\neq\tilde{\phi}_{j}^{(d)}$.

Remarks: The probability with which a multiplier can be found which obeys the
desired property is rather arbitrary in this Lemma and can be increased by
choosing a smaller $\epsilon$. Note that it is easy to verify whether for a
randomly chosen multiplier the desired conditions hold or not. The validity of
this Lemma importantly does not depend on whether the phase estimates are
actually accurate, it only depends on the number of phases $n_{\phi}$. In
practice, we do not generate a random multiplier $\kappa_{d+1}$ through this
Lemma, but simply exhaustively search for a valid $\kappa_{d+1}$ starting at
the maximal value, see Section 5.

The reason to adaptively choose the multiplier $\kappa_{d+1}$ for
$d=0,1,\ldots$ is that two (estimated) phases in principle need to lead to
separate estimates at the next order: this is expressed in Eq. (38). An
exception to this occurs when the (estimated) phases are still close enough,
as in Eq. (39), so that their next-order refined estimates could merge at the
next order, see Fig. 2. Phase estimates can thus split and merge over the
multiple orders. They split when sufficient accuracy is available at the next
order to distinguish them, they can stay or are allowed to merge when such
accuracy is not yet needed at the given order.

In what follows below we will assume, just for simplicity of the proof, that
the error parameter $\epsilon$ is bounded by $\epsilon_{{\rm crit},0}$ in Eq.
(28) for all rounds, and $\epsilon$ is the same for all rounds, including the
first one.

Figure 2: Schematic of the execution of Alg. 4.1 to estimate three hidden
phases, $\phi_{0}$, $\phi_{1}$ and $\phi_{2}$ (dashed lines). The algorithm
progresses from bottom to top as the estimation order $d$ increases. At each
order $d$, the phase extraction subroutine (Alg. 3.4) promises to return an
estimate $\tilde{\theta}^{(d)}_{\rightarrow j}$ of each
$\theta^{(d)}_{j}=k_{d}\phi_{j}\mod 2\pi$ that corresponds to an estimate
$\tilde{\phi}^{(d)}_{\rightarrow j}$ (solid lines) lying within the promised
estimate region about $\phi_{j}$ (coloured boxes). By matching phases at
subsequent orders (arrows), the algorithm is able to converge to an ever-more
accurate estimate of each $\phi_{j}$. The phase extraction subroutine only
promises that each region will contain at least one phase (and that the total
number of estimates at each order is bounded by $n_{\phi}$) - when two regions
overlap, the subroutine may merge the phases to give a single estimate (green
and purple lines). Estimates at subsequent orders may continue to separate and
even re-merge until the regions separate, at which point the algorithm
promises with high confidence a precise estimate for each hidden phase. The
estimate $\tilde{\phi}_{\rightarrow j}^{(d)}$ at each order is only known
$\mod(2\pi)/k_{d}$, leading to a set of potential aliases (dotted lines at
$d=4$) for each phase. We do not know a priori which alias is correct, and
must rely on the fact that the true estimate $\tilde{\phi}_{\rightarrow
j}^{(d)}$ needs to be close to a previous estimate $\tilde{\phi}_{\rightarrow
j}^{(d-1)}$. By carefully choosing each $k_{d}$, we can guarantee that no
alias will satisfy this condition (so long as Alg. 3.4 succeeds), and our
phase matching will be unambiguous.

###  4.2 Bounding the error with and without failures

In this section we state and prove the two key intermediate lemmas, Lemma 4.3
and Lemma 4.4 on our way towards proving that Alg. 4.1 reaches the Heisenberg-
limit. Together, these lemmas allow us to bound the error in Alg. 4.1,
–assuming that the phase extraction subroutine succeeds for the first $d$
rounds–, to within $O(\epsilon/k_{d})$.

These Lemmas deal with the issue of ‘aliasing’ or the correct matching of new
estimates with older estimates which is solved by the specific choice of
$\kappa_{d+1}$ in step 4f of Alg. 4.1, see also Fig. 2. It is important to
note that there is no 1-1 relation between these estimates and the actual
phases since the number of estimates is at most the number of phases.

Let $d_{f}$ be the last order executed in Alg. 4.1, i.e. the last order for
which we go through step 4b, construct the estimates
$\\{\tilde{\theta}_{l}^{(d_{f})}\\}$ and pass the tests at step 4c and 4e and
output $\\{\tilde{\phi}_{l}^{(d_{f})}\\}$. When none of the failure modes is
encountered, $d_{f}$ is set by the first $k_{d}$ such that
$k_{d}\geq\frac{\epsilon}{\delta_{c}}$ (since the next $k_{d+1}\geq
2\epsilon/\delta_{c}$ as $\kappa_{d+1}\geq 2$). Since $\kappa_{d}\geq 2$, we
observe that

$d_{f}\leq\log_{2}\left(\frac{2\epsilon}{\delta_{c}}\right).$ (40)

In Corollary 4.3.1 we argue that when the QEEP subroutines, Alg. 3.4, succeed
up to order $d_{f}$, we indeed never exit via these failure modes.

######  Lemma 4.3.

If each invocation of the QEEP subroutine, Alg. 3.4, succeeds in Alg. 4.1 up
to order $d_{f}$, then in this last round $d_{f}$ in step 4d it holds that

  * •

(Property 1a) For every phase $\phi_{j}$ there exists an estimate
$\tilde{\phi}_{l}^{(d_{f})}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d_{f})}|_{T}\leq 2\epsilon/k_{d_{f}}.$

  * •

(Property 1b) For every estimate $\tilde{\phi}_{l}^{(d_{f})}$ there exists a
phase $\phi_{j}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d_{f})}|_{T}\leq 2\epsilon/k_{d_{f}}.$

###### Proof.

We prove this Lemma by induction. Consider the first round $d=0$
($k_{d=0}=1$), i.e. step 1 of Alg. 4.1. If the QEEP subroutine, Alg. 3.4,
succeeds (with probability $p_{0}$) then Lemma 3.5 holds, namely for each
$\phi_{j}$ there exists an estimate $\tilde{\phi}_{l}^{(0)}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(0)}|_{T}\leq 2\epsilon.$ (41)

and for each estimate $\tilde{\phi}_{l}^{(0)}$ there exists at least one
$\phi_{j}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(0)}|_{T}\leq 2\epsilon.$ (42)

Hence the statement to be proven holds at $d=0$. Now consider step 2 of Alg.
4.1 and invoke Lemma B.1 for which the assumptions are fulfilled by Eqs.
(41),(42). This implies that with the choice of $k_{1}\geq 3n_{\phi}$ in step
3 in Alg. 4.1 the shifted phases and their 0th-order estimates obey the
technical condition in Lemma 2.5 and we can use Eq. (8). In the next steps we
work with these shifted phases but for simplicity we don’t use any new
notation and refer to them as $\phi_{j}$ and estimates
$\tilde{\phi}_{j}^{(d)}$ etc.

Now assume the statement to be proven holds at order $d$ i.e. let
$\\{\tilde{\phi}_{l}^{(d)}\\}$ be a set of at most $n_{\phi}$ estimates of the
phases $\\{\phi_{j}\\}$ with

  * •

(Assumption 1a) For every phase $\phi_{j}$ there exists an estimate
$\tilde{\phi}_{l}^{(d)}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d)}|_{T}\leq 2\epsilon/k_{d}.$

  * •

(Assumption 1b) For every estimate $\tilde{\phi}_{l}^{(d)}$ there exists a
phase $\phi_{j}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d)}|_{T}\leq 2\epsilon/k_{d}.$

Note that these assumptions certainly imply that one can apply Lemma 2.5 to
the estimates $\tilde{\phi}_{l}^{(d)}$. That is, given that the real phases
$\phi_{j}$ are $2\epsilon/k_{d}$ close to these estimates and that the
(shifted) $\phi_{j}$ obey Eq. (106), it implies that Eq. (8) can be used with
$k\geq 3n_{\phi}$ (which is the case for all rounds $d\geq 1$).

We consider the QEEP subroutine, Alg. 3.4, with a given choice of $\kappa_{d}$
obeying the conditions in step 3 (for $d=1$) and step 4f (for higher $d$),
executed in step 4b with confidence $p_{d}$. In the math below we refer to
conditions on $\kappa_{d>1}$, namely Eq. (38) and Eq. (39), but the conditions
on $\kappa_{1}$ in Eq. (33) and Eq. (34) are of identical form (so we don’t
make a separate argument for the $d=0\rightarrow d=1$ induction step).

Let thus $\\{\tilde{\theta}^{(d+1)}_{l}\\}$ be a set of estimates of the
eigenphases $\\{\theta^{(d+1)}_{j}\\}$ of $U^{k_{d+1}}$ corresponding to the
set $\\{\phi_{j}\\}$, that is,

$\theta^{(d+1)}_{j}=k_{d+1}\phi_{j}\mod 2\pi.$ (43)

and $k_{d+1}=k_{d}\kappa_{d+1}$. By assuming that Alg. 3.4 succeeds we can
invoke Lemma 3.5, namely

  * •

(Assumption 2a) For every phase $\theta_{j}^{(d+1)}$ there exists an estimate
$\tilde{\theta}_{l}^{(d+1)}$ such that

$|\theta_{j}^{(d+1)}-\tilde{\theta}_{l}^{(d+1)}|_{T}\leq 2\epsilon.$

  * •

(Assumption 2b) For every estimate $\tilde{\theta}_{l}^{(d+1)}$ there exists a
phase $\theta_{j}^{(d+1)}$ such that

$|\theta_{j}^{(d+1)}-\tilde{\theta}_{l}^{(d+1)}|_{T}\leq 2\epsilon.$

To prove the induction step, we thus need to show that the set
$\tilde{\phi}_{j}^{(d+1)}$ generated by step 4d of Alg. 4.1 satisfies the
following two properties

  * •

(Property 1a) For every phase $\phi_{j}$ there exists an estimate
$\tilde{\phi}_{l}^{(d+1)}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d+1)}|_{T}\leq 2\epsilon/k_{d+1}.$

  * •

(Property 1b) For every estimate $\tilde{\phi}_{l}^{(d+1)}$ there exists a
phase $\phi_{j}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d+1)}|_{T}\leq 2\epsilon/k_{d+1}.$

First consider Assumption 2a. Assumption 2a implies that for every phase
$\phi_{j}$ there exists a $\tilde{\theta}_{l}^{(d+1)}$ such that

$|k_{d+1}\phi_{j}-\tilde{\theta}_{l}^{(d+1)}|_{T}\leq 2\epsilon.$ (44)

In this proof we will use the label $l=\rightarrow j$ for this
$\tilde{\theta}_{l}^{(d+1)}$ associated with $\phi_{j}$. Thus, also using Eq.
(8), Eq. (44) is equivalent to

$\displaystyle\min_{n\in\\{0,\ldots,\lfloor
k_{d+1}\rfloor-1\\}}|\phi_{j}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+$
$\displaystyle 2\pi n)/k_{d+1}|_{T}$
$\displaystyle\leq\frac{2\epsilon}{k_{d+1}},$ (45)

with $n^{\rm ideal}_{j,\rightarrow j}$ the integer which achieves the minimum,
i.e.

$\displaystyle n^{\rm ideal}_{j,\rightarrow j}=$
$\displaystyle\operatorname*{arg\,min}_{n\in\\{0,\ldots,\lfloor
k_{d+1}\rfloor-1\\}}\Big{\\{}$
$\displaystyle|\phi_{j}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi
n)/k_{d+1}|_{T}\Big{\\}}.$ (46)

Similarly, by Assumption 1a, there is some $\tilde{\phi}_{l=\rightarrow
j}^{(d)}$ which is $2\epsilon/k_{d}$-close to $\phi_{j}$, again using a label
which shows this association.

Consider the optimization in Eq. (37) at step 4d in Alg. 4.1. We define

$\displaystyle\xi_{l}$ $\displaystyle=\min_{j}\xi_{j,l},$ (47)
$\displaystyle\xi_{j,l}$
$\displaystyle\equiv\left|\tilde{\phi}^{(d)}_{j}-(\tilde{\theta}_{l}^{(d+1)}+2\pi
n_{j,l})/k_{d+1}\right|_{T}$ (48)

with

$\displaystyle n_{j,l}=$
$\displaystyle\operatorname*{arg\,min}_{n\in\\{0,\ldots,\lfloor
k_{d+1}\rfloor-1\\}}\Big{\\{}$
$\displaystyle\left|\tilde{\phi}^{(d)}_{j}-(\tilde{\theta}_{l}^{(d+1)}+2\pi
n)/k_{d+1}\right|_{T}\Big{\\}},$ (49) $\displaystyle n_{l}=$
$\displaystyle\operatorname*{arg\,min}_{n_{j,l}}\xi_{j,l}.$ (50)

The goal is thus to prove that for each $\phi_{j}$, using the corresponding
$\tilde{\theta}_{\rightarrow j}^{(d+1)}$, we have $n_{\rightarrow j}=n^{\rm
ideal}_{j,\rightarrow j}$ which directly implies Property 1a.

We can bound using Eq. (6) and then Eq. (8), Assumptions 1a and 2a and the
optimality of $n_{\rightarrow j,\rightarrow j}$,

$\displaystyle\xi_{\rightarrow j,\rightarrow j}$ $\displaystyle\hskip
5.69046pt=\left|\tilde{\phi}^{(d)}_{\rightarrow
j}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi n_{\rightarrow j,\rightarrow
j})/k_{d+1}\right|_{T}$ $\displaystyle\hskip
5.69046pt\leq\left|\tilde{\phi}^{(d)}_{\rightarrow j}-\phi_{j}\right|_{T}$
$\displaystyle\hskip
5.69046pt\qquad+\left|\phi_{j}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi
n^{\rm ideal}_{j,\rightarrow j})/k_{d+1}\right|_{T}$ $\displaystyle\hskip
5.69046pt\leq\frac{2\epsilon}{k_{d}}+\frac{1}{k_{d+1}}\left|k_{d+1}\phi_{j}-\tilde{\theta}^{(d+1)}_{\rightarrow
j}\right|_{T}$ $\displaystyle\hskip
5.69046pt=\frac{2\epsilon}{k_{d}}+\frac{1}{k_{d+1}}\left|\theta_{j}^{(d+1)}-\tilde{\theta}_{\rightarrow
j}^{(d+1)}\right|_{T}$ $\displaystyle\hskip
5.69046pt\leq\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}.$ (51)

Now if Eq. (38) holds for some other $m\neq\rightarrow j$, we claim on the
other hand that

$\xi_{m,\rightarrow j}>\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}},$ (52)

hence matching $\theta_{\rightarrow j}^{(d+1)}$ with such
$\tilde{\phi}_{m}^{(d)}$, with $m\neq\rightarrow j$ is non-optimal and will
not be chosen in the Algorithm. To indeed see that Eq. (38) implies Eq. (52),
we can calculate

$\displaystyle\frac{4\epsilon(1+\kappa_{d+1})}{k_{d+1}}<\frac{1}{k_{d+1}}\left|k_{d+1}\tilde{\phi}_{\rightarrow
j}^{(d)}-k_{d+1}\tilde{\phi}_{m}^{(d)}\right|_{T}$ $\displaystyle\hskip
56.9055pt\leq\frac{1}{k_{d+1}}\left|k_{d+1}\tilde{\phi}_{\rightarrow
j}^{(d)}-\tilde{\theta}_{\rightarrow j}^{(d+1)}\right|_{T}$
$\displaystyle\hskip
56.9055pt+\frac{1}{k_{d+1}}\left|\tilde{\theta}_{\rightarrow
j}^{(d+1)}-k_{d+1}\tilde{\phi}_{m}^{(d)}\right|_{T}$
$\displaystyle\leq\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}+\min_{n\in\\{0,\ldots,\lfloor
k_{d+1}\rfloor-1\\}}\Big{\\{}$ $\displaystyle\hskip
56.9055pt\Big{|}\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow
j}^{(d+1)}+2\pi n)/k_{d+1}\Big{|}_{T}\Big{\\}}$ (53)
$\displaystyle\Rightarrow\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}<$
$\displaystyle\min_{n\in\\{0,\ldots,\lfloor
k_{d+1}\rfloor-1\\}}\left|\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow
j}^{(d+1)}+2\pi n)/k_{d+1}\right|_{T}$ $\displaystyle\hskip
113.81102pt=\xi_{m,\rightarrow j}.$ (54)

Here we have used that Eq. (8) holds for the estimate
$\tilde{\phi}_{m}^{(d)}$.

Alternatively, for those $m\neq\rightarrow j$ for which Eq. (39) holds, we
claim that

$n_{m,\rightarrow j}=n_{\rightarrow j,\rightarrow j},$ (55)

hence for those $\tilde{\phi}^{(d)}_{m}\neq\tilde{\phi}^{(d)}_{\rightarrow j}$
the algorithm produces a single new estimate equal to
$(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi n_{\rightarrow j,\rightarrow
j})/k_{d+1}$.

To see that Eq. (39) implies Eq. (55) indeed, note that it is sufficient to
prove that

$\left|\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi
n_{\rightarrow j,\rightarrow j})/k_{d+1}\right|_{T}<\frac{\pi}{k_{d+1}},$ (56)

as one can then show that for $n^{\prime}\neq n_{\rightarrow j,\rightarrow
j}\in\\{0,\ldots,\lfloor k_{d+1}\rfloor-1\\}$ that

$\displaystyle\frac{2\pi}{k_{d+1}}\leq\left|\frac{2\pi}{k_{d+1}}(n_{\rightarrow
j,\rightarrow j}-n^{\prime})\right|_{T}$ $\displaystyle\hskip
14.22636pt=|(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi n_{\rightarrow
j,\rightarrow j})/k_{d+1}$ $\displaystyle\hskip
14.22636pt\qquad\qquad-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi
n^{\prime})/k_{d+1}|_{T}$ $\displaystyle\hskip
14.22636pt\leq\left|(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi
n_{\rightarrow j,\rightarrow j})/k_{d+1}-\tilde{\phi}_{m}^{(d)}\right|_{T}$
$\displaystyle\hskip
14.22636pt\qquad+\left|\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow
j}^{(d+1)}+2\pi n^{\prime})/k_{d+1}\right|_{T}$ $\displaystyle\hskip
14.22636pt<\frac{\pi}{k_{d+1}}+\left|\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow
j}^{(d+1)}+2\pi n^{\prime})/k_{d+1}\right|_{T}$ (57)
$\displaystyle\Rightarrow\frac{\pi}{k_{d+1}}<\left|\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow
j}^{(d+1)}+2\pi n^{\prime})/k_{d+1}\right|_{T},$ (58)

so $n_{\rightarrow j,\rightarrow j}$ is optimal. Using Eq. (6), Eq. (39) and
Eq. (62), we can prove Eq. (56) since

$\displaystyle\left|\tilde{\phi}_{m}^{(d)}-(\tilde{\theta}_{\rightarrow
j}^{(d+1)}+2\pi n_{\rightarrow j,\rightarrow j})/k_{d+1}\right|_{T}$
$\displaystyle\hskip
5.69046pt\leq\left|\tilde{\phi}_{m}^{(d)}-\tilde{\phi}_{\rightarrow
j}^{(d)}\right|$ $\displaystyle\hskip
11.38092pt+\left|\tilde{\phi}_{\rightarrow
j}^{(d)}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi n_{\rightarrow
j,\rightarrow j})/k_{d+1}\right|_{T}$ $\displaystyle\hskip
5.69046pt<\frac{\pi-2\epsilon(1+\kappa_{d+1})}{k_{d+1}}+\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}$
$\displaystyle\hskip 113.81102pt=\frac{\pi}{k_{d+1}}.$ (59)

We have thus shown that for each $\phi_{j}$, there is a
$\tilde{\theta}_{\rightarrow j}^{(d+1)}$, such that step 4d will output
$(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi n_{\rightarrow j,\rightarrow
j})/k_{d+1}$ with $n_{\rightarrow j,\rightarrow j}$ defined in Eq. (50),
related to the previous order estimate $\phi_{\rightarrow j}^{(d)}$ which was
already close to $\phi_{j}^{(d)}$. The last step is to show that
$n_{\rightarrow j,\rightarrow j}=n_{j,\rightarrow j}^{\rm ideal}$ using
Property 1a. It holds that

$\displaystyle|\phi_{j}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi
n_{\rightarrow j,\rightarrow j})/k_{d+1}|_{T}$ $\displaystyle\hskip
5.69046pt\leq|\phi_{j}-\tilde{\phi}^{(d)}_{\rightarrow j}|_{T}$
$\displaystyle\hskip 14.22636pt+|\tilde{\phi}^{(d)}_{\rightarrow
j}-(\tilde{\theta}_{\rightarrow j}^{(d+1)}+2\pi n_{\rightarrow j,\rightarrow
j})/k_{d+1}|_{T}$ $\displaystyle\hskip
5.69046pt\leq\frac{2\epsilon}{k_{d}}+\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}<\frac{\pi}{k_{d+1}},$
(60)

where we used that $4\epsilon(\kappa_{d+1}+1)<\pi$. Indeed for $d>1$,
$\epsilon<\frac{\pi}{16}$ ($\kappa_{d+1}\leq 3$) and for $d=0$,
$\epsilon<\frac{\pi}{4(3n_{\phi}+2)}$ ($\kappa_{1}\leq 3n_{\phi}+1$) given the
upper bounds on $\epsilon$ in Eqs. (29) and (28). This implies through the
same argument as in Eq. (58) that $n_{j,\rightarrow j}^{\rm ideal}$ achieving
the minimum in Eq. (45) equals $n_{\rightarrow j,\rightarrow j}$ and hence we
obtain Property 1a.

Now let’s prove Property 1b. Given a $\tilde{\theta}_{l}^{(d+1)}$, let
$\phi_{j=\rightarrow l}$ be the real phase for which, by Assumption 2b, it
holds that

$|k_{d+1}\phi_{\rightarrow l}-\tilde{\theta}_{l}^{(d+1)}|_{T}\\\
=k_{d+1}|\phi_{\rightarrow l}-(\tilde{\theta}_{l}^{(d+1)}+2\pi n^{\rm
ideal}_{\rightarrow l,l})/k_{d+1}|_{T}\\\ \leq 2\epsilon.$ (61)

To prove Property 1b, we need to show that $n_{l}=n^{\rm ideal}_{\rightarrow
l,l}$ with $n_{l}$ defined in Eq. (50). Let also $\tilde{\phi}_{\rightarrow
j}^{(d)}$ be the previous order $2\epsilon/k_{d}$-close estimate to
$\phi_{j=\rightarrow l}$ by Assumption 1a. The idea is that
$\tilde{\phi}_{\rightarrow j}^{(d)}=\tilde{\phi}_{\rightarrow(\rightarrow
l)}^{(d)}$ is matched to $\tilde{\theta}_{l}^{(d+1)}$ in the optimization step
of the algorithm, so that this leads to a better estimate for the phase
$\phi_{j=\rightarrow l}$.

Given a $\tilde{\theta}_{l}^{(d+1)}$ we can deduce, as before, that

$\displaystyle\xi_{\rightarrow j,l}$
$\displaystyle=\left|\tilde{\phi}^{(d)}_{\rightarrow
j}-(\tilde{\theta}_{l}^{(d+1)}+2\pi n_{\rightarrow j,l})/k_{d+1}\right|_{T}$
$\displaystyle\leq\left|\tilde{\phi}^{(d)}_{\rightarrow j}-\phi_{\rightarrow
l}\right|_{T}$ $\displaystyle\hskip 2.84544pt+\left|\phi_{\rightarrow
l}-(\tilde{\theta}_{l}^{(d+1)}+2\pi n^{\rm ideal}_{\rightarrow
l,l})/k_{d+1}\right|_{T}$
$\displaystyle\leq\frac{2\epsilon}{k_{d}}+\frac{1}{k_{d+1}}\left|k_{d+1}\phi_{\rightarrow
l}-\tilde{\theta}^{(d+1)}_{l}\right|_{T}$
$\displaystyle=\frac{2\epsilon}{k_{d}}+\frac{1}{k_{d+1}}\left|\theta_{\rightarrow
l}^{(d+1)}-\tilde{\theta}_{l}^{(d+1)}\right|_{T}$
$\displaystyle\leq\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}.$ (62)

Using previous arguments, all other $\xi_{m,l}$ are either larger or give the
same integer $n_{\rightarrow j,l}$ and thus $n_{l}=n_{\rightarrow j,l}$. In
addition, we can bound, using this equality and Assumption 1a

$|\phi_{\rightarrow l}-(\tilde{\theta}_{l}^{(d+1)}+2\pi n_{\rightarrow
j,l})/k_{d+1}|_{T}\\\ \leq|\phi_{\rightarrow
l}-\tilde{\phi}^{(d)}_{\rightarrow j}|_{T}+|\tilde{\phi}^{(d)}_{\rightarrow
j}-(\tilde{\theta}_{l}^{(d+1)}+2\pi n_{\rightarrow j,l})/k_{d+1}|_{T}\\\
\leq\frac{2\epsilon}{k_{d}}+\frac{2\epsilon(1+\kappa_{d+1})}{k_{d+1}}\\\
<\frac{\pi}{k_{d+1}},$ (63)

implying that $n_{l}=n_{\rightarrow j,l}=n^{\rm ideal}_{\rightarrow l,l}$ as
desired. ∎

Algorithm 4.1 has a few failure modes, namely steps 1, 4c and 4e where we exit
and return an estimate of lower order. Arguments in Lemma 4.3 show that these
failure modes are only encountered when the QEEP subroutine, Alg. 3.4, fails
at some order. For clarity, we collect these observations in a single
Corollary:

######  Corollary 4.3.1.

If each invocation of the QEEP subroutine, Alg. 3.4, succeeds in Alg. 4.1, we
never exit at step 1, 4c or 4e.

###### Proof.

Consider step 1 of Alg. 4.1 applying Alg. 3.4 which obeys Lemma 3.5, showing
that success of Alg. 3.4 implies that the number of phases is at most
$n_{\phi}$. By assumption there is at least one phase with $A_{j}>0$, and
hence success means that Alg. 3.4 cannot return the empty set due to statement
1. of Lemma 3.5. Hence if Alg. 3.4 succeeds we do not exit at step 1 of Alg.
4.1. Now consider step 4c of Alg. 4.1: again success of Alg. 3.4 implies that
the number of estimates does not exceed $n_{\phi}$. Consider Eq. (35) and Eq.
(36); we wish to show that these will not hold if Alg. 3.4 succeeds up to
order $d-1$. If Alg. 3.4 succeeds up to order $d-1$, the phase estimates
$\tilde{\phi}_{j}^{(d-1)}$ obey Eq. (8) for $k_{d}$, hence the condition in
Eq. (36) equals, for each $\tilde{\theta}_{l}^{(d)}$

$\displaystyle\min_{j}\min_{n}|\tilde{\phi}_{j}^{(d-1)}-(\tilde{\theta}_{l}^{(d)}+2\pi
n)/k_{d}|_{T}$ $\displaystyle\hskip
113.81102pt>\frac{2\epsilon(1+\kappa_{d})}{k_{d}},$ (64)

and we argued previously, via induction, that this does not happen when Alg.
3.4 succeeds up to order $d$, as $\min_{j}\xi_{j,l}$ is upper-bounded as in
Eq. (62) for all $d^{\prime}\leq d$. Similarly, Eq. (35) implies the existence
of a $\tilde{\phi}_{j}^{(d-1)}$ with

$\displaystyle\min_{l}\min_{n}|\tilde{\phi}_{j}^{(d-1)}-(\tilde{\theta}_{l}^{(d)}+2\pi
n)/k_{d}|_{T}$ $\displaystyle\hskip
113.81102pt>\frac{2\epsilon(1+\kappa_{d})}{k_{d}}.$ (65)

which can not happen due to the success of Alg. 3.4 which implies the bound in
Eq. (51). Consider lastly step 4e which exits if the current $d$th order
estimates do not lie in the region for which Eq. (7) holds with given $k_{d}$.
We have argued in Lemma 4.3 that, assuming success of the subroutines
implementing Alg. 3.4, that Eq. (7) holds for the phase estimates at all
orders. ∎

Now let us consider failures of the QEEP subroutine, Alg. 3.4, which do not
lead to exiting. Let’s imagine that the first failure occurs at some order
$d_{0}$. Now we want to make sure that continuing with higher orders after
such failure still leads to an error of order $\sim\epsilon/k_{d_{0}-1}$, even
though the failure (or any subsequent failure) is not detected.

To show this, we check that if Alg. 4.1 exits at some later round, namely
during $d=d_{f}+1$ and outputs estimates $\tilde{\phi}_{j}^{(d_{f})}$ that
these will be sufficiently close to the estimates right before failure, that
is, the set of phases $\tilde{\phi}_{j}^{(d_{0}-1)}$.

Then, by Lem. 4.3, these estimates will also be sufficiently close to the true
phases $\phi_{j}$.

######  Lemma 4.4.

Let Alg. 4.1 exit at order $d=d_{f}+1$ and let the QEEP subroutine, Alg. 3.4,
of step 4b first fail at $d=d_{0}\leq d_{f}+1$. For each $\phi_{j}$, there
will be an estimate $\tilde{\phi}_{l}^{(d_{f})}$, produced at step 4d in Alg.
4.1 which satisfies

$\left|\phi_{j}-\tilde{\phi}_{l}^{(d_{f})}\right|_{T}\leq\frac{14\epsilon}{k_{d_{0}-1}}.$
(66)

Vice-versa, for each estimate $\tilde{\phi}_{l}^{(d_{f})}$ there exists a
phase $\phi_{j}$ such that

$\left|\phi_{j}-\tilde{\phi}_{l}^{(d_{f})}\right|_{T}\leq\frac{14\epsilon}{k_{d_{0}-1}}.$
(67)

###### Proof.

Since each QEEP subroutine, Alg. 3.4, in Alg. 4.1 succeeds up to order
$d_{0}-1$, Lemma 4.3 guarantees that

  * •

(Property 1a) For every phase $\phi_{j}$ there exists an estimate
$\tilde{\phi}_{l}^{(d_{0}-1)}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d_{0}-1)}|_{T}\leq\frac{2\epsilon}{k_{d_{0}-1}}.$

  * •

(Property 1b) For every estimate $\tilde{\phi}_{l}^{(d_{0}-1)}$ there exists a
phase $\phi_{j}$ such that

$|\phi_{j}-\tilde{\phi}_{l}^{(d_{0}-1)}|_{T}\leq\frac{2\epsilon}{k_{d_{0}-1}}.$

Then, since the algorithm does not exit at step 4c through Eqs. (35) or step
4e for any order $d=d_{0},\ldots d_{f}$, for each estimate
$\tilde{\phi}_{l}^{(d-1)}$ we can associate some
$\tilde{\theta}_{m_{l}}^{(d)}$ that satisfies

$\frac{1}{k_{d}}\left|k_{d}\tilde{\phi}_{l}^{(d-1)}-\tilde{\theta}_{m_{l}}^{(d)}\right|_{T}=\\\
\min_{n\in\\{0,\ldots,\lfloor
k_{d+1}\rfloor-1\\}}\left|\tilde{\phi}_{l}^{(d-1)}-(\tilde{\theta}_{m_{l}}^{(d)}+2\pi
n)/k_{d}\right|_{T}\\\ \leq\frac{2\epsilon(1+\kappa_{d})}{k_{d}},$ (68)

where the second equality follows from being allowed to apply Eq. (8) (which
is validated by passing the test at step 4e). This implies that in step 4d of
Alg. 4.1 at round $d$, for a given $\tilde{\theta}_{m_{l}}^{(d)}$, the
optimization of $\xi_{n,m_{l}}$ over $n$ will pick the integer $n_{l,m_{l}}$,
i.e. the integer associated with matching $\tilde{\theta}_{m_{l}}^{(d)}$ with
$\tilde{\phi}_{l}^{(d-1)}$. Next, similar as in the proof of Lemma 4.3, we can
consider the possibility of matching to other estimates
$\tilde{\phi}_{k}^{(d-1)}\neq\tilde{\phi}_{l}^{(d-1)}$. Since $\kappa_{d}$ is
chosen in step 4f of Alg. 4.1, we claim that either Eq. (38) holds, in which
case

$\min_{n\in\\{0,\ldots,\lfloor
k_{d}\rfloor-1\\}}\left|\tilde{\phi}_{k}^{(d-1)}-\frac{2\pi
n+\tilde{\theta}_{m_{l}}^{(d)}}{k_{d}}\right|_{T}\\\
>\frac{2\epsilon}{k_{d}}(1+\kappa_{d}),$ (69)

hence this $\xi_{k,m_{l}}$ is not optimal, or that Eq. (39) holds, in which
case

$n_{l,m_{l}}=\\\ \operatorname*{arg\,min}_{n\in\\{0,\ldots,\lfloor
k_{d}\rfloor-1\\}}\left|\tilde{\phi}_{k}^{(d-1)}-\frac{2\pi
n+\tilde{\theta}_{m_{l}}^{(d)}}{k_{d}}\right|_{T}\\\ =n_{k,m_{l}}.$ (70)

The proofs of these claims are exactly the same as in the proof of Lemma 4.3,
i.e. using Eqs. (56), (58),(59).

Now let us prove Eq. (66). Given a phase $\phi_{j}$, we can use Property (1a)
to find an associated estimate $\tilde{\phi}_{\rightarrow j}^{(d_{0}-1)}$
within $2\epsilon/k_{d_{0}-1}$. Then for this estimate let
$\theta^{(d_{0})}_{m_{\rightarrow j}}$ be the matched estimate in the next
round for which Eq. (68) holds, so that the round produces a new estimate
$\tilde{\phi}^{(d_{0})}_{\rightarrow
j}=\frac{\tilde{\theta}^{(d_{0})}_{m_{\rightarrow j}}+2\pi n_{\rightarrow
j,m_{\rightarrow j}}}{k_{d_{0}}}$ (which we label with $\rightarrow j$ again)
for which

$|\tilde{\phi}_{\rightarrow j}^{(d_{0}-1)}-\tilde{\phi}^{(d_{0})}_{\rightarrow
j}|_{T}\leq\frac{2\epsilon}{k_{d_{0}}}(1+\kappa_{d_{0}}).$ (71)

Then again for $\tilde{\phi}^{(d_{0})}_{\rightarrow j}$ there exists some
matching $\theta^{(d_{0}+1)}_{m_{\rightarrow j}}$ etc. and this generates a
series of estimates $\tilde{\phi}^{(d)}_{\rightarrow j}$ up to order $d_{f}$.
For a given $\phi_{j}$ we can then bound, using this series of estimates and
$\kappa_{d}\geq 2$ for all $d$,

$\displaystyle\left|\tilde{\phi}_{j}^{(d_{f})}-\phi_{j}\right|_{T}\leq\left|\tilde{\phi}_{\rightarrow
j}^{(d_{0}-1)}-\phi_{j}\right|_{T}$ $\displaystyle\hskip
28.45274pt\qquad+\sum_{d=d_{0}}^{d_{f}}\left|\tilde{\phi}_{\rightarrow
j}^{(d-1)}-\tilde{\phi}_{\rightarrow j}^{(d)}\right|_{T}$ $\displaystyle\hskip
14.22636pt\leq\frac{2\epsilon}{k_{d_{0}-1}}+\sum_{d=d_{0}}^{d_{f}}\frac{2\epsilon(1+\kappa_{d})}{k_{d}}$
$\displaystyle\hskip
14.22636pt=\frac{2\epsilon}{k_{d_{0}-1}}\left(1+\sum_{d=d_{0}}^{d_{f}}\frac{1+\kappa_{d}}{\kappa_{d_{0}}\kappa_{d_{0}+1}\ldots\kappa_{d}}\right)$
$\displaystyle\hskip
14.22636pt\leq\frac{2\epsilon}{k_{d_{0}-1}}\left(1+\sum_{n=0}^{\infty}\frac{3}{2^{n}}\right)$
$\displaystyle\hskip 14.22636pt=\frac{14\epsilon}{k_{d_{0}-1}}.$ (72)

Now let’s prove Eq. (67) and start with an estimate
$\tilde{\phi}_{l}^{(d_{f})}$ which was obtained from some
$\tilde{\theta}_{l}^{(d_{f})}$ matched with a previous estimate
$\tilde{\phi}_{l}^{(d_{f}-1)}$ (just for convenience we again use the same
label) such that
$|\tilde{\phi}_{l}^{(d_{f})}-\tilde{\phi}_{l}^{(d_{f}-1)}|_{T}\leq\frac{2\epsilon}{k_{d_{f}}}(1+\kappa_{d_{f}})$,
using that we do not exit through Eq. (36). Then again for this previous
estimate $\tilde{\phi}_{l}^{(d_{f}-1)}$ we can repeat the argument and create
a sequence of estimates up to $\tilde{\phi}^{(d_{0}-1)}_{l}$. For the last
estimate, we invoke Property (1b), namely that there is a nearby $\phi_{j}$.
Then we can upperbound for this $\phi_{j}$:
$|\tilde{\phi}_{l}^{(d_{f})}-\phi_{j}|_{T}\leq|\tilde{\phi}_{l}^{(d_{0}-1)}-\phi_{j}|_{T}+\sum_{d=d_{0}}^{d_{f}}|\tilde{\phi}_{l}^{(d)}-\tilde{\phi}_{l}^{(d-1)}|_{T}$
etc., exactly as in Eq. (72), leading to Eq. (67). ∎

###  4.3 Algorithm 4.1 achieves the Heisenberg limit

We have seen that the success of the QEEP subroutines in Alg. 4.1 leads to an
error scaling as $\epsilon/k_{d_{f}}\sim\delta_{c}$. Now we must choose the
success probability $p_{d}$ of these subroutines in Eq. (30), depending on
$\alpha,\gamma$ so that the total mean-squared-error is bounded by some
$\delta^{2}=O(\delta_{c}^{2})$ while the quantum cost $T=O(\delta^{-1})$. We
note that the next theorem contains no logarithmic factors in $\delta^{-1}$,
as in [23], but achieves pure Heisenberg scaling.

######  Theorem 4.5.

Algorithm 4.1 solves the multiple eigenvalue estimation problem in Def. 3.3
with accuracy error $\delta$ and total quantum cost $T=O(\delta^{-1})$, given
$A,n_{\phi}$ and a fixed $\epsilon_{0}$ and $\epsilon$ obeying Eqs. (29) and
(28), and some choice for the constants $\alpha>0$ and $\gamma>2$.

Remarks: Note that the dependence on the number of phases $n_{\phi}$ is not
made explicit in the statement of this Theorem, but this dependence will be
polynomial in $n_{\phi}$, not necessarily a very low-order polynomial. This
dependence comes through the choice for $\epsilon_{0}$ (and $\epsilon$) via
Eq. (28) (resp. Eq. (29)) which sets the error and thus the running time of
the QEEP Algorithm 3.2.

###### Proof.

Our proof is motivated by the analysis in [2] for a single phase $\phi$
leading to Theorem 2.7. The idea is to bound the mean-squared-error in the
final estimation of $\phi$ by summing over error contributions at each order
$d$ at which the phase extraction subroutine may fail (with probability
$1-p_{d}$).

In our case the multipliers $\kappa_{d}$ (and $k_{d}$) at each order are not
fixed (as in Theorem 2.7) but depend on phase estimates at previous orders and
thus measurement data at previous orders. Our confidence parameter $p_{d}$ in
Eq. (30), which determines the number of repeats of experiments, and hence the
cost, in Alg. 3.4, depends on $k_{d}$ and is thus a random variable depending
on previous measurement data. All measurement data are denoted by ${\bf x}$
and thus we have random variables $k_{d}(\\{\kappa_{d^{\prime}}({\bf
x})\\}_{d^{\prime}=1}^{d})$ and $p_{d}(\\{\kappa_{d^{\prime}}({\bf
x})\\}_{d^{\prime}=1}^{d})$.

Consider the mean-squared-error $\delta_{j}^{2}$ for the $j$th phase
$\phi_{j}$ in Eq. (24) in Definition 3.3. We have three error contributions to
consider given a choice for the random variable $k_{d}$.

  1. 1.

With probability $1-p_{0}$ the subroutine Alg. 3.4 in Alg. 4.1 fails at step 1
($d=0$). In this case, as we always return some estimate, $\delta_{j}$ is
bounded for all $j$ by $\pi$.

  2. 2.

With probability at most $(1-p_{d_{0}})\prod_{d=0}^{d_{0}-1}p_{d}\leq
1-p_{d_{0}}=e^{-\alpha}\left(\frac{k_{d_{0}}\delta_{c}}{\pi}\right)^{\gamma}$,
the subroutine Alg. 3.4 in Alg. 4.1 fails for the first time at some order
$1\leq d_{0}\leq d_{f}$, and the algorithm proceeds in any way afterwards (by
possibly exiting or not). In this case, Lemma 4.4 bounds $\delta_{j}$ for all
$j$ by $\frac{14\epsilon}{k_{d_{0}-1}}$ or Lemma 4.3 bounds $\delta_{j}$ for
all $j$ by $\frac{2\epsilon}{k_{d_{0}-1}}\leq\frac{14\epsilon}{k_{d_{0}-1}}$.

  3. 3.

With probability less than $\prod_{d=0}^{d_{f}}p_{d}<1$ the subroutine Alg.
3.4 in Alg. 4.1 succeeds up to the final round $d_{f}$, and Lemma 4.3 implies
that $\delta_{j}\leq 2\epsilon/k_{d_{f}}\leq 2\delta_{c}$ for all $j$ as
$k_{d_{f}}\geq\epsilon/\delta_{c}$.

We can now bound the mean-squared-error as a sum over the above three
contributions weighted by their relevant probabilities:

$\displaystyle\delta_{j}^{2}$ $\displaystyle\leq(1-p_{0})\pi^{2}+\sum_{{\bf
x}}\mathbb{P}({\bf
x})\left[\sum_{d_{0}=1}^{d_{f}-1}\mathbb{P}(\kappa_{1},\ldots,\kappa_{d_{0}}|{\bf
x})(1-p_{d_{0}}(\\{\kappa_{d}^{\prime}\\}_{d^{\prime}=1}^{d_{0}}))\left[\frac{14\epsilon}{k_{d_{0}-1}(\\{\kappa_{d}^{\prime}\\}_{d^{\prime}=1}^{d_{0}-1})}\right]^{2}\right]+4\delta_{c}^{2}$
$\displaystyle=\pi^{2}e^{-\alpha}\left(\frac{\delta_{c}}{\pi}\right)^{\gamma}+\sum_{{\bf
x}}\mathbb{P}({\bf
x})\left[\sum_{d_{0}=1}^{d_{f}-1}\mathbb{P}(\kappa_{1},\ldots,\kappa_{d_{0}}|{\bf
x})e^{-\alpha}\left(\frac{k_{d_{0}}\delta_{c}}{\pi}\right)^{\gamma}\frac{196\epsilon^{2}}{k_{d_{0}-1}^{2}}\right]+4\delta_{c}^{2}$
$\displaystyle\leq\pi
e^{-\alpha}\left(\frac{\delta_{c}}{\pi}\right)^{\gamma}+0.15\times\delta_{c}^{2}\sum_{{\bf
x}}\mathbb{P}({\bf x})\mathbb{P}(\kappa_{1},\ldots,\kappa_{d_{0}}|{\bf
x})\sum_{d_{0}=1}^{d_{f}-1}\left(\frac{k_{d_{0}}\delta_{c}}{\pi}\right)^{\gamma-2}+4\delta_{c}^{2}.$
(73)

Here we have removed the dependency of $k_{d_{0}}$ and $k_{d_{0}-1}$ on the
previous multipliers for notational simplicity. For $d=1$ we have
$e^{-\alpha}\kappa_{d_{0}}^{2}196(\epsilon/\pi)^{2}=e^{-\alpha}k_{1}^{2}196(\epsilon_{0}/\pi)^{2}\leq\frac{196\times
16\times 4}{(300)^{2}}\leq 0.15$ due to Eq. (28). For $d>1$,
$e^{-\alpha}\kappa_{d_{0}}^{2}196(\epsilon/\pi)^{2}\leq\frac{9\times 196\times
4}{(300)^{2}}=0.08$, due to Eq. (29).

To evaluate the middle term, we write
$k_{d_{0}}=k_{d_{f}}\frac{k_{d_{0}}}{k_{d_{f}}}$, and note that as
$k_{d}=\prod_{d^{\prime}=1}^{d}\kappa_{d^{\prime}}$, we have
$\frac{k_{d_{0}}}{k_{d_{f}}}\leq 2^{d_{0}-d_{f}}$ as the multiplier
$\kappa_{d}\geq 2$. As
$k_{d_{f}}\leq\frac{2\epsilon}{\delta_{c}}<\frac{\pi}{2\delta_{c}}$, we have

$\displaystyle\sum_{d_{0}=1}^{d_{f}-1}\left(\frac{k_{d_{0}}\delta_{c}}{\pi}\right)^{\gamma-2}$
$\displaystyle\leq\frac{1}{2^{\gamma-2}}\sum_{d_{0}=1}^{d_{f}-1}(2^{\gamma-2})^{(d_{0}-d_{f})}$
$\displaystyle\leq\frac{2^{4-\gamma}}{2^{\gamma}-4},$ (74)

where the last inequality holds since $\gamma>2$. By letting the upper bound
be independent of the $\kappa_{d}$s, we can remove the dependence on ${\bf x}$
in Eq. (73), using that $\sum_{{\bf x}}\mathbb{P}({\bf
x})\mathbb{P}(\kappa_{1},\ldots,\kappa_{d_{0}}|{\bf x})=1$. This yields a
final bound on $\delta_{j}$ of

$\delta_{j}^{2}\leq\delta_{c}^{2}\left[\pi^{1-\gamma}e^{-\alpha}\delta_{c}^{\gamma-2}+\frac{0.15\times
2^{4-\gamma}}{2^{\gamma}-4}+4]\right].$ (75)

As $\gamma>2$, this scales as $\delta_{c}^{2}$ as $\delta_{c}\rightarrow 0$.

Let us now calculate the cost of executing Alg. 4.1 in terms of the number of
unitary applications. Again this depends on the choice of multiplier
$\kappa_{d}$ at each step. Let us fix a sequences of $k_{d}$s, and let $d_{f}$
be the final round of estimation in this algorithm, i.e. the final round for
which we invoked the quantum subroutine in Alg. 3.4. At each order $d$ we use
$V^{\sf k}=U^{k_{d}{\sf k}}$, where ${\sf k}=0,1,\ldots,K$, with $2M_{d}$
samples where $K$ is a function of $\epsilon$ as in Theorem 3.2. By Lemma 2.2,
the cost of each experiment is $\lfloor k_{d}{\sf k}\rfloor+O(M_{d})$.

We can calculate

$\displaystyle T=$ $\displaystyle\sum_{{\bf x}}\mathbb{P}({\bf
x})\mathbb{P}(\kappa_{1},\ldots,\kappa_{d_{0}}|{\bf x})$ $\displaystyle\hskip
28.45274pt\times\sum_{d=0}^{d_{f}}\sum_{{\sf k}=1}^{K}(2M_{d}\lfloor k_{d}{\sf
k}\rfloor+O(M_{d}^{2}))$ $\displaystyle\leq$
$\displaystyle\sum_{d=0}^{d_{f}}M_{d}k_{d}K(K+1)$ $\displaystyle\hskip
56.9055pt+(d_{f}+1)KO(M_{d=0}^{2}),$ (76)

where we have used $M_{d=0}\geq M_{d}$. The QEEP algorithm in Theorem 3.2
requires $M_{d}=\tilde{O}(|\ln(1-p_{d})|\epsilon^{-4})$ with $p_{d}$ in Eq.
(30) and $K=\tilde{O}(\epsilon^{-1})$.

The last term in Eq. (76) only has ${\rm poly}(\ln\delta_{c}^{-1})$
dependence, using Eq. (40) for $d_{f}$ and the expression of $M_{0}$ in terms
of $\delta_{c}$, and we ignore it in the rest of the proof. We may bound

$\displaystyle T$
$\displaystyle\leq\tilde{O}(\epsilon^{-6})\sum_{d=0}^{d_{f}}k_{d}\left|-\alpha+\gamma\ln\left(\frac{k_{d}\delta_{c}}{\pi}\right)\right|$
$\displaystyle=\tilde{O}(\epsilon^{-6})\sum_{d=0}^{d_{f}}k_{d}\left[\alpha-\gamma\ln\left(\frac{k_{d}\delta_{c}}{\pi}\right)\right].$
(77)

We again bound $k_{d}=\frac{k_{d}}{k_{d_{f}}}k_{d_{f}}\leq
2^{d-d_{f}}\frac{\pi}{2\delta_{c}}$, which yields

$\displaystyle T$
$\displaystyle\leq\delta_{c}^{-1}\tilde{O}(\epsilon^{-6})\Bigg{[}\frac{\pi}{2}$
$\displaystyle\times\sum_{d=0}^{d_{f}}2^{d-d_{f}}\big{[}\alpha-\gamma(d-d_{f}-1)\ln(2)\big{]}\Bigg{]}$
$\displaystyle\leq\delta_{c}^{-1}\tilde{O}(\epsilon^{-6})(\alpha+2\gamma\ln(2)).$
(78)

Combining our bounds then yields

$\displaystyle\delta$
$\displaystyle\leq\delta_{c}\left[e^{-\alpha}\pi^{1-\gamma}+4+\frac{0.15\times
2^{4-\gamma}}{2^{\gamma}-4}\right]^{\frac{1}{2}}$ (79) $\displaystyle\leq
T^{-1}\left[e^{-\alpha}\pi^{1-\gamma}+4+\frac{0.15\times
2^{4-\gamma}}{2^{\gamma}-4}\right]^{\frac{1}{2}}$
$\displaystyle\qquad\times\tilde{O}(\epsilon^{-6})[\alpha+2\gamma\ln(2)]$
$\displaystyle=O(T^{-1}),$ (80)

which is the Heisenberg limit. Note that a dependence on the number of phase
$n_{\phi}$ enters the scaling of $\delta$ via $\epsilon$ which needs to be
bounded as in Eqs. (28) and (29). ∎

##  5 Numerical implementation

Thm. 4.5 requires using the QEEP algorithm (Theorem 3.2 and Algorithm 3.4) in
order to obtain provable bounds. Instead of analytic bounds, we now turn to a
numerical demonstration, giving the opportunity to implement and test
Algorithm 4.1 with a few modifications. We test the algorithm using two
different sub-routines, one based on the matrix pencil method [27], and one
based on the QEEP time-series analysis of Theorem 3.2, as described in
Algorithm 3.4. Code to implement all simulations can be found at
https://github.com/alicjadut/qpe.

To improve the practical performance of Alg. 4.1, we make the following two
small changes. Firstly, instead of choosing $\kappa_{d+1}$ in step 4f in the
ranges declared in Lem. 4.2, we choose the largest $\kappa_{d+1}$ consistent
with Eq. (38) and Eq. (39) for all
$\tilde{\phi}_{j}^{(d)},\tilde{\phi}_{l}^{(d)}$. We note that the maximum such
$\kappa_{d+1}$ is bounded above by $\frac{\pi}{2\epsilon}-1$, as the left-hand
side of Eq. (38) is bounded above by $2\pi$ and the left-hand side of Eq. (39)
is bounded below by $0$. (In practice, tighter bounds can be found by checking
the boundaries of the regions $R_{jl}^{(n)}$ defined in Eq. (88), and we find
the largest possible $\kappa_{d+1}$ by iterating backwards through these
boundaries till a gap is found.) Secondly, as the bounds for $\epsilon$ and
$\epsilon_{0}$ in Lem. 4.2 are rather loose, and our performance scales rather
badly in both, we choose the largest $\epsilon=\epsilon_{0}$ that allows all
simulations to find a value of $\kappa_{d+1}>2$ at each order.

When using the matrix pencil processing subroutine, we follow the
implementation described in Ref. [16]:

######  Algorithm 5.1.

The matrix pencil method takes as input estimates of the phase function
$g({\sf k})=\sum_{j}A_{j}e^{i{\sf k}\theta_{j}}$ for a unitary $V$ at points
${\sf k}=0,1,\ldots K$ and a probability bound $A$, and proceeds as follows:

  1. 1.

Construct the $L_{K}\times(2K-L_{K}+1)$ Hankel matrices $G^{(0)}$, $G^{(1)}$,
where $G^{(a)}_{i,j}=g(i+j+a-K)$ for
$i\in\\{0,1,...L_{K}-1\\},j\in\\{0,1,...,2K-L_{K}\\}$, $a=0,1$, with
$L_{K}=\lfloor(K+1)/2\rfloor$, and using $g(-{\sf k})=g^{*}({\sf k})$.

  2. 2.

Construct the $L_{K}\times L_{K}$ shift matrix $T$ by least-squares
minimization of the matrix $2$-norm $\|TG^{(0)}-G^{(1)}\|$.

  3. 3.

Calculate the eigenvalues of $T$,
$\lambda_{j}=|\lambda_{j}|e^{i\tilde{\theta}_{j}}$ and from there the phase
estimates $\tilde{\theta}_{j}$.

  4. 4.

Calculate the probability estimates $\tilde{A}_{j}$ by least-squares
minimization of the vector $2$-norm $\|BA-g\|$, where $B$ is the $(K+1)\times
L_{K}$ matrix

$B_{k,j}=\lambda_{j}^{{\sf k}},$ (81)

and $g=[g(0),...g(K)]^{T}$.

  5. 5.

Return the phase estimates $\tilde{\theta}_{j}$ for which the corresponding
probability estimate $\tilde{A}_{j}\geq A$.

To use this algorithm as a subroutine in Alg. 4.1 (in place of Alg. 3.4), we
implement it on the matrix $V=U^{k_{d}}$, which requires implementing $V^{\sf
k}=U^{{\sf k}k_{d}}$ for a range of integer ${\sf k}$ on a quantum device.

To isolate the performance of the estimation routine from the generation of
the signal itself, we do not test our protocols on data generated from
simulating or approximating a particular unitary. Instead, we test the ability
of the protocols to estimate $n_{\phi}=2$ and $n_{\phi}=4$ randomly-chosen
phases $\phi_{j}\in[0,2\pi]$ when sampling from the true phase function
$g(k)$. We take all phases with equal weight — $A_{j}=1/n_{\phi}$. We simulate
the sampling from $g(k)$ in Algorithm 5.1 or Algorithm 3.2 for some
$V=U^{k_{d}}$ by simulating the readout of a control qubit with the reduced
density matrix of Eq. (4). (In practice this would be generated by the quantum
circuit in Fig. 1.) We first draw $M_{d}$ i.i.d. samples from the two
Bernoulli distributions

$\displaystyle\mathbb{P}_{k}^{r}(+1)=\frac{1}{2}\sum_{j=1}^{n_{\phi}}A_{j}(1+\cos(\theta_{j}k)),$
(82)
$\displaystyle\mathbb{P}_{k}^{i}(+1)=\frac{1}{2}\sum_{j=1}^{n_{\phi}}A_{j}(1-\sin(\theta_{j}k)),$
(83)

where $\theta_{j}=k_{d}\phi_{j}\mod 2\pi$ are the eigenvalues of $V$. Then, we
return the fraction of $+1$s drawn as estimates for the real and imaginary
parts of $g(k)$ respectively. This is then repeated at all points $k=k_{d}{\sf
k}$ for ${\sf k}=0,1,\ldots,K$. Following the discussion in Sec. 2.2 and using
the notation from Eq. (76), we sum the total quantum cost for the algorithm
over all requested $g(k)$ queries; $T=2\sum_{d}\sum_{{\sf k}=1}^{K}{\sf
k}k_{d}M_{d}$. (We ignore the sub-leading correction from the final term in
Eq. (76) as this will not affect the scaling of our result.) For the signal
length $K$ and number of points $M_{d}$ to sample each $g(k)$ at, we follow
the bounds given in Ref. [3] (both when using the QEEP and matrix pencil
subroutines):

  * •

signal length: $K=\lceil 0.1L\ln^{2}L\rceil$, with
$L=\lceil\frac{2\pi}{\epsilon}\rceil$ the number of bins used in the QEEP
subroutine (Def. 3.1).

  * •

number of measurements of each circuit:
$M_{d}=\left\lceil\left|\ln\left(1-p_{d}\right)\right|\epsilon^{-4}\right\rceil$.

Here, $p_{d}$ is given in Eq. (30) for a given $k_{d}$. This equation requires
fixing a choice of $\alpha$ and $\gamma$ — across all experiments we take
$\alpha=2$ and $\gamma=2.1$.

Figure 3:  Convergence of the algorithm with total quantum cost $T$. Phase
estimates were obtained with either QEEP (blue) or matrix pencil (green)
subroutines with parameters described in the text. Individual points show the
error $\delta=|\tilde{\phi}_{j}-\phi_{j}|_{T}$ on individual phases in each
simulation. This data is binned in the x-axis, and for each bin a root-mean-
square error and standard deviation (error bars) are plotted in the x- and
y-direction. Dotted lines show a fit of these means to $\delta\sim T^{-1}$.

To demonstrate that our methods achieve the Heisenberg limit, in Fig. 3 we
plot the error as a function of the total quantum cost for a set of
simulations using the methods described above. Each simulation draws a
different set of $n_{\phi}$ random phases, and a final error
$\delta_{c}\in[10^{-5},10^{-2}]$ (for each choice of $\delta_{c}$ we use the
same $50$ sets of phases). We plot the error for each phase estimate
separately in Fig. 3 (i.e. each simulation corresponds to $n_{\phi}$ points in
the plot). As both the total quantum cost and the error is different between
simulations, we bin all experiments within a range of $T$ values, and
calculate the root mean square error and root mean square total quantum cost.
This gives a good approximation to the accuracy error defined in Def. 3.3 for
the restricted data set used.

For the QEEP subroutine, we observe a clear fit of the data (blue points) to a
$\delta\sim T^{-1}$ trend, as expected from Thm. 4.5, but with a rather large
constant factor; we find $T\sim 10^{10}\delta^{-1}$ for estimating $2$ phases
and $T\sim 10^{15}\delta^{-1}$ for estimating $4$. Further optimization of the
QEEP algorithm for these purposes may yet improve this constant factor.
However, as the methods of Ref. [3] were not designed for estimating
individual phases, it may be expected that this method performs somewhat badly
for this purpose, so we have not pursued this further.

Simulations using the matrix pencil subroutine outperform simulations using
the QEEP subroutine by a factor of $10^{4}-10^{6}$, and clearly demonstrate
Heisenberg-limited scaling $\delta\sim T^{-1}$ as well. We take this result
instead of an analytic proof as strong numerical evidence for Heisenberg-
limited scaling when a version of Alg. 4.1 is constructed using the matrix
pencil method as a subroutine. We notice that the error in two phase estimates
in the first bin of the matrix pencil method for $n_{\phi}=2$ is significantly
above the remainder of the population (by about a factor $100\times$), which
blows up our error bars for this bin. Further investigation shows that the two
phases in question are from the same simulation, and separated by only
$1.5\times 10^{-4}$. By contrast, for the simulation in question (at $d=1$)
our algorithm sets $k_{1}K\sim 3\times 10^{3}<(1.5\times 10^{-4})^{-1}$ (where
$k_{1}K$ is the largest value of the phase function $g(k)$ sampled during this
simulation). This implies that our signal lies within the region where
improving our estimation accuracy by increasing the number of shots $M_{d}$ is
exponentially hard [28]. In latter simulations with these two phases we see
that our estimation error regresses to similar results as all other estimates.

##  6 Conclusion

In this work we studied Heisenberg-limited quantum phase estimation using a
single control qubit. In this form of phase estimation, we rely on classical
signal processing to extract eigenvalue data from the phase function $g(k)$ in
Eq. (2).

It has been an open question whether these methods can achieve the Heisenberg
limit in the case of multiple phases: Ref. [23] answered this question up to
log factors with a Heisenberg-limited Monte Carlo algorithm, providing a
sampling of the spectral function $A(\phi)$ in Eq. (3) from which to estimate
the phases. In this work we also answered this question in the affirmative
exactly with a new adaptive multi-order phase estimation algorithm, for which
we prove Heisenberg scaling if the algorithm uses a QEEP phase extraction
subroutine. We numerically show the performance of this algorithm, also when
instead of using a QEEP subroutine, one uses the matrix pencil method to
extract phase estimates from the phase function $g(k)$.

In obtaining our results we encountered at least two details of quantum phase
estimation that we have not seen discussed in the literature. The first is the
dense signal limit, Eq. (12) in Thm. 2.6: sampling $g(k)$ at all integer point
$k=1,\ldots,K$ is sub-optimal regardless of what method is used to process the
data. However, we also briefly argued that by picking points among
$k=1,\ldots,K$ at random one may go beyond this, and one could interpret this
as allowing the results obtained in [23] in which such randomized choices for
$k$ are taken.

The second is the need for adaptive choices of $k_{d}$ to solve the phase
matching problem. It is unclear to us how far this problem extends; although
Lemma 4.2 provides a practical solution, others may still exist. Another open
question with respect to Algorithm 4.1 is whether one can remove the need to
choose real-valued multipliers $\kappa_{d}$ and restrict to integer choices.
Restricting $\kappa_{d}\in\mathbb{N}$ would significantly simplify some
technical issues, i.e. the applicability of Lemma 2.5 and the need for
shifting phases in step 2 of Alg. 4.1, but we don’t know how to prove a
version of Lemma 4.2 for $\kappa_{d}\in\mathbb{N}$. In fact, we don’t know
whether there is a fundamental difference in performance between only using
data obtained with integer $k$ in $g(k)$ versus data obtain with real-valued,
–in practice rational–, $k$ in $g(k)$.

We have assumed in our problem description, Def. 3.3, that the spectrum in the
input state is discrete consisting of $n_{\phi}$ phases with probability above
some cut-off. In practice this condition may not be fulfilled and thus
studying the performance of the algorithms on more typical spectra induced by
many-body Hamiltonians and easy-to-prepare input states will be of interest.

A direction for future research is to make this algorithm efficient in
practice (i.e. improve the parameter dependence and the practical run time) or
devise yet-different Heisenberg-scaling algorithms and examine their
performance in the presence of experimentally-noisy signals $g(k)$.

## Acknowledgements

AD acknowledges financial support in the form of a scholarship from the Den
Adel Fund.

## References

  * Higgins et al. [2009] B. L. Higgins, D. W. Berry, S. D. Bartlett, M. W. Mitchell, H. M. Wiseman, and G. J. Pryde.  Demonstrating Heisenberg-limited unambiguous phase estimation without adaptive measurements.  _New J. Phys._ , 11(7):073023, 2009.  doi: https://doi.org/10.1088/1367-2630/11/7/073023.  URL https://arxiv.org/abs/0809.3308.
  * Kimmel et al. [2015] Shelby Kimmel, Guang Hao Low, and Theodore J. Yoder.  Robust calibration of a universal single-qubit gate-set via robust phase estimation.  _Phys. Rev. A_ , 92:062315, 2015.  doi: https://doi.org/10.1103/PhysRevA.92.062315.  URL https://arxiv.org/abs/1502.02677.
  * Somma [2019] Rolando D. Somma.  Quantum eigenvalue estimation via time series analysis.  _New J. Phys._ , 21:123025, 2019.  doi: https://doi.org/10.1088/1367-2630/ab5c60.  URL https://iopscience.iop.org/article/10.1088/1367-2630/ab5c60/pdf.
  * Wocjan and Zhang [2006] Pawel Wocjan and Shengyu Zhang.  Several natural BQP-complete problems.  _ArXiv:quant-ph/0606179_ , 2006.  URL https://arxiv.org/abs/quant-ph/0606179.
  * Shor [1997] Peter W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  _SIAM J. Sci. Stat. Comp._ , 26:1484, 1997.  doi: https://doi.org/10.1137/S0097539795293172.  URL https://arxiv.org/abs/quant-ph/9508027.
  * Harrow et al. [2009] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for solving linear systems of equations.  _Phys. Rev. Lett._ , 15(103):150502, 2009.  doi: https://doi.org/10.1103/PhysRevLett.103.150502.  URL https://arxiv.org/abs/0811.3171.
  * Whitfield et al. [2011] James D. Whitfield, Jacob Biamonte, and Alán Aspuru-Guzik.  Simulation of electronic structure Hamiltonians using quantum computers.  _Mol. Phys._ , 109:735–750, 2011.  doi: https://doi.org/10.1080/00268976.2011.552441.  URL https://arxiv.org/abs/1001.3855.
  * Nielsen and Chuang [2000] M.A. Nielsen and I.L. Chuang.  _Quantum Computation and Quantum Information_.  Cambridge Series on Information and the Natural Sciences. Cambridge University Press, 2000.  ISBN 9780521635035.  doi: https://doi.org/10.1017/CBO9780511976667.  URL https://books.google.de/books?id=65FqEKQOfP8C.
  * Cleve et al. [1998] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca.  Quantum algorithms revisited.  _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ , 454(1969):339–354, 1998.  doi: https://doi.org/10.1098/rspa.1998.0164.  URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1998.0164.
  * Griffiths and Niu [1996] Robert B. Griffiths and Chi-Sheng Niu.  Semiclassical Fourier transform for quantum computation.  _Physical Review Letters_ , 76(17):3228–3231, Apr 1996.  ISSN 1079-7114.  doi: https://doi.org/10.1103/physrevlett.76.3228.  URL http://dx.doi.org/10.1103/PhysRevLett.76.3228.
  * Kitaev [1995] A. Yu. Kitaev.  Quantum measurements and the Abelian stabilizer problem.  _ArXiv:quant-ph/9511026_ , 1995.  URL https://arxiv.org/abs/quant-ph/9511026.
  * Berry et al. [2007] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders.  Efficient quantum algorithms for simulating sparse Hamiltonians.  _Comm. Math. Phys._ , 270(359), 2007.  doi: https://doi.org/10.1007/s00220-006-0150-x.  URL https://arxiv.org/abs/quant-ph/0508139.
  * Wiebe and Granade [2016] Nathan Wiebe and Chris Granade.  Efficient Bayesian phase estimation.  _Phys. Rev. Lett._ , 117:010503, 2016.  doi: https://doi.org/10.1103/PhysRevLett.117.010503.  URL https://arxiv.org/abs/1508.00869.
  * Svore et al. [2013] Krysta M. Svore, Matthew B. Hastings, and Michael Freedman.  Faster phase estimation.  _Quant. Inf. Comp._ , 14(3-4):306–328, 2013.  URL https://arxiv.org/abs/1304.0741.
  * van den Berg [2020] Ewout van den Berg.  Efficient Bayesian phase estimation using mixed priors.  _ArXiv:2007.11629_ , 2020.  URL https://arxiv.org/abs/2007.11629.
  * O’Brien et al. [2019] Thomas E O’Brien, Brian Tarasinski, and Barbara M Terhal.  Quantum phase estimation of multiple eigenvalues for small-scale (noisy) experiments.  _New J. Phys._ , 21:023022, 2019.  doi: https://doi.org/10.1088/1367-2630/aafb8e.  URL https://iopscience.iop.org/article/10.1088/1367-2630/aafb8e.
  * Rife and Boorstyn [1974] David C. Rife and Robert R. Boorstyn.  Single-tone parameter estimation from discrete-time observations.  _IEEE Trans. Inf. Th._ , 20(5):591–598, 1974\.  doi: https://doi.org/10.1109/TIT.1974.1055282.  URL https://ieeexplore.ieee.org/document/1055282.
  * Lu et al. [2020] Sirui Lu, Mari Carmen Bañuls, and J. Ignacio Cirac.  Algorithms for quantum simulation at finite energies.  _PRX Quantum_ , 2:020321, 2020.  doi: 10.1103/PRXQuantum.2.020321.  URL https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.020321.
  * O’Brien et al. [2020] T.E. O’Brien, S. Polla, N.C. Rubin, W.J. Huggins, S. McArdle, S. Boixo, J.R. McClean, and R. Babbush.  Error mitigation via verified phase estimation.  _ArXiv:2010.02538_ , 2020.  doi: https://doi.org/10.1103/PRXQuantum.2.020317.  URL https://arxiv.org/abs/2010.02538.
  * Roggero [2020] Alessandro Roggero.  Spectral density estimation with the Gaussian integral transform.  _ArXiv:2004.04889_ , 2020.  doi: https://doi.org/10.1103/PhysRevA.102.022409.  URL https://arxiv.org/abs/2004.04889.
  * Berry et al. [2014] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma.  Exponential improvement in precision for simulating sparse Hamiltonians.  _Proc. ACM STOC_ , pages 283–282, 2014.  doi: https://doi.org/10.1145/2591796.2591854.  URL https://arxiv.org/abs/1312.1414.
  * Regev [2004] O. Regev.  A subexponential time algorithm for the dihedral hidden subgroup problem with polynomial space.  _ArXiv:quant-ph/0406151_ , 2004.  URL https://arxiv.org/abs/quant-ph/0406151.
  * Lin and Tong [2021] Lin Lin and Yu Tong.  Heisenberg-limited ground state energy estimation for early fault-tolerant quantum computers.  _ArXiv:2102.11340_ , 2021.  URL https://arxiv.org/abs/2102.11340.
  * Gebhart et al. [2020] Valentin Gebhart, Augusto Smerzi, and Luca Pezzè.  Heisenberg-limited bayesian multiphase estimation algorithm.  _ArXiv:2010.09075_ , 2020.  URL https://arxiv.org/abs/2010.09075.
  * Cramér [1946] Harald Cramér.  _Mathematical Methods of Statistics_.  Princeton University Press, 1946.  ISBN 0691080046.  doi: https://doi.org/10.1515/9781400883868.  URL https://archive.org/details/in.ernet.dli.2015.223699.
  * Rao [1945] Calyampudi Radakrishna Rao.  Information and the accuracy attainable in the estimation of statistical parameters.  _Bull. Calcutta Math. Soc._ , 37:81–89, 1945.  doi: https://doi.org/10.1007/978-1-4612-0919-5_16.  URL https://link.springer.com/chapter/10.1007/978-1-4612-0919-5_16.
  * Hua and Sarkar [1990] Y. Hua and T. K. Sarkar.  Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise.  _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , 38(5):814–824, 1990.  doi: https://doi.org/10.1109/29.56027.  URL https://ieeexplore.ieee.org/document/56027.
  * Moitra [2015] Ankur Moitra.  Super-resolution, extremal functions and the condition number of Vandermonde matrices.  In _Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing_ , STOC ’15, page 821–830, New York, NY, USA, 2015. Association for Computing Machinery.  ISBN 9781450335362.  doi: https://doi.org/10.1145/2746539.2746561.  URL https://doi.org/10.1145/2746539.2746561.

##  Appendix A Proof of Lemma 4.2

Let us first prove the existence of $\kappa_{d+1}\in[2,\kappa_{\rm max}]$ with
$\kappa_{\rm max}=3$ that satisfies our conditions (the proof for $k_{1}$ is
similar), for small enough $\epsilon$ in Eq. (29). Note that Eq. (29) implies

$\epsilon\leq\frac{\pi}{300}\approx 0.01.$ (84)

Given some pair $\tilde{\phi}^{(d)}_{j}\neq\tilde{\phi}^{(d)}_{l}$, $j<l$, let
$\Delta_{j,l}=|\tilde{\phi}^{(d)}_{j}-\tilde{\phi}^{(d)}_{l}|$ and let
$R_{j,l}$ be a set of $\kappa_{d+1}$ such that neither Eq. (39) nor Eq. (38)
holds for the chosen phases, that is,

$\displaystyle R_{j,l}=\bigg{\\{}$
$\displaystyle\kappa_{d+1}\in[2,\kappa_{\max}]:$
$\displaystyle|\Delta_{j,l}|_{T}\geq\frac{\pi-2\epsilon(1+\kappa_{d+1})}{k_{d}\kappa_{d+1}}$
$\displaystyle\land$ $\displaystyle|k_{d}\kappa_{d+1}\Delta_{j,l}|_{T}\leq
4\epsilon(1+\kappa_{d+1})\bigg{\\}}.$ (85)

We call $\cup_{j,l}R_{j,l}$ the forbidden region and want to show that we can
choose a value for $\kappa_{d+1}\in[2,3]$ outside this forbidden region if
$\epsilon$ is sufficiently small. We do this by bounding the size of the
forbidden region above and showing that this is smaller than the region
$[2,3]$, leaving room to choose $\kappa_{d+1}$.

Note that $R_{j,l}$ is nonempty only if

$\displaystyle k_{d}\Delta_{j,l}$ $\displaystyle\geq k_{d}|\Delta_{j,l}|_{T}$
$\displaystyle\geq\frac{\pi-2\epsilon(1+\kappa_{d+1})}{\kappa_{d+1}}$
$\displaystyle\geq\frac{\pi-2\epsilon(1+\kappa_{\max})}{\kappa_{\rm max}}$
$\displaystyle>\frac{73\pi}{225},$ (86)

using Eq. (84).

We may write the set $R_{j,l}$ as

$\displaystyle R_{j,l}=$
$\displaystyle\left[\max\left(2,\frac{\pi-2\epsilon}{k_{d}|\Delta_{j,l}|_{T}+2\epsilon}\right),\kappa_{\max}\right]$
$\displaystyle\cap\bigcup_{n\in\mathbb{N}}R_{j,l}^{(n)},$ (87)

where $R_{j,l}^{(n)}$ is the set of $\kappa_{d+1}$ for which

$\displaystyle|k_{d}\kappa_{d+1}\Delta_{j,l}|_{T}$
$\displaystyle=\bigg{|}k_{d}\kappa_{d+1}\Delta_{j,l}-2\pi n\bigg{|}$
$\displaystyle\leq 4\epsilon(1+\kappa_{d+1}),$ (88)

for some $n\in\mathbb{N}$.

Solving this equation for $\kappa_{d+1}$ yields

$\displaystyle R_{j,l}^{(n)}=$ $\displaystyle\left[\frac{2\pi
n-4\epsilon}{k_{d}\Delta_{j,l}+4\epsilon},\frac{2\pi
n+4\epsilon}{k_{d}\Delta_{j,l}-4\epsilon}\right].$ (89)

The size of each interval $R_{j,l}^{(n)}$ can then be calculated

$\left|R_{j,l}^{(n)}\right|=\frac{8\epsilon(2n\pi+k_{d}\Delta_{j,l})}{k_{d}^{2}\Delta_{j,l}^{2}-16\epsilon^{2}}.$
(90)

We can bound

$\displaystyle\left|R_{j,l}\right|$
$\displaystyle\leq\left|\bigcup_{n=1}^{n_{\max}}R^{(n)}_{j,l}\right|\leq\sum_{n=1}^{n_{\max}}\left|R^{(n)}_{j,l}\right|$
$\displaystyle=\frac{8\epsilon}{k_{d}^{2}\Delta_{j,l}^{2}-16\epsilon^{2}}\bigg{(}\pi
n_{\max}(n_{\max}+1)$
$\displaystyle\qquad\qquad\qquad+k_{d}\Delta_{j,l}n_{\max}\bigg{)}.$ (91)

Here, $n_{\max}=n_{\rm max}(j,l)$ is the largest index of a set
$R^{(n)}_{j,l}$ in Eq. (88) for which $\kappa_{d+1}\in[2,\kappa_{\max}]$.
Since $\kappa_{d+1}\leq 3$, Eq. (88) implies that

$n_{\max}\leq\frac{3(k_{d}\Delta_{j,l}+4\epsilon)+4\epsilon}{2\pi},$ (92)

Now using Eq. (86) and Eq. (84) gives

$\displaystyle\big{|}$ $\displaystyle
R_{j,l}\big{|}\leq\frac{1}{1-(4\epsilon/k_{d}\Delta_{j,l})^{2}}$
$\displaystyle\times\frac{30\epsilon(864000\epsilon^{2}+248160\epsilon\pi+11899\pi^{2})}{5329\pi^{3}}$
$\displaystyle\leq\frac{3\epsilon(864000\epsilon^{2}+248160\epsilon\pi+11899\pi^{2})}{532\pi^{3}}$
$\displaystyle\leq 0.24,$ (93)

where the last inequality used Eq. (84). As there are $n_{\phi}\geq 2$ phases
the length of the total forbidden region $\cup_{j,l}R_{j,l}$ is bounded from
above by $\frac{n_{\phi}^{2}}{2}\left|R_{j,l}\right|$. We want to this
interval to be, say, at most $1/4$, so that by choosing $\kappa_{d+1}$
randomly we have a 75$\%$ change of not landing in the forbidden interval. For
larger $n_{\phi}$ we thus should use

$\epsilon\leq\epsilon_{\mathrm{crit}}=\frac{2\pi}{300n_{\phi}^{2}},$ (94)

leading to Eq. (29).

We now repeat the above approach for the special case of finding the
multiplier in the first round $\kappa_{1}=k_{1}$. Consider thus
$k_{1}\in[3n_{\phi},\kappa_{\max}]$ with $\kappa_{\rm max}=3n_{\phi}+1$. The
key difference here is that there is a stricter lower bound on this multiplier
$\kappa_{\rm max}$ so that $\epsilon$ needs to be chosen smaller, depending on
$n_{\phi}$, namely we choose

$\epsilon_{0}\leq\epsilon_{\mathrm{crit},0}=\frac{2\pi}{300n_{\phi}^{4}},$
(95)

as expressed in Eq. (28).

Given some pair $\tilde{\phi}^{(0)}_{j}\neq\tilde{\phi}^{(0)}_{l}$, $j<l$, and
again let $\Delta_{j,l}=|\tilde{\phi}^{(0)}_{j}-\tilde{\phi}^{(0)}_{l}|$.
Then, let $R_{j,l}$ be the set of $k_{1}$ such that neither Eq. (39) nor Eq.
(38) holds for the chosen phases. That is,

$\displaystyle R_{j,l}=\bigg{\\{}$ $\displaystyle
k_{1}\in\left[3n_{\phi},\kappa_{\max}\right]:$
$\displaystyle|\Delta_{j,l}|_{T}\geq\frac{\pi-2\epsilon_{0}(1+k_{1})}{k_{1}}$
$\displaystyle\land$ $\displaystyle|k_{1}\Delta_{j,l}|_{T}\leq
4\epsilon_{0}(1+k_{1})\bigg{\\}}.$ (96)

We again call $\cup_{j,l}R_{j,l}$ the forbidden region and want to show that
we can choose a value for $k_{1}\in[3n_{\phi},3n_{\phi}+1]$ outside this
forbidden region, assuming that $\epsilon_{0}$ is small enough. Note that the
logic of the first few inequalities in Eq. (86) still holds in this new
calculation, leading to

$\displaystyle\Delta_{j,l}$
$\displaystyle\geq\frac{\pi-2\epsilon_{0}(1+\kappa_{\rm max})}{\kappa_{\rm
max}}$ $\displaystyle\hskip
28.45274pt\geq\frac{\pi(1-\frac{2+3n_{\phi}}{75n_{\phi}^{4}})}{1+3n_{\phi}},$
(97)

where the second inequality used Eq. (95) and the value for $\kappa_{\rm
max}$. Note that for large $n_{\phi}$ this allows $\Delta_{j,l}$ to decrease
like $\sim 1/n_{\phi}$, while previously $\Delta_{j,l}$ was lowerbounded by a
constant, Eq. (86).

This time, we may write the set $R_{j,l}$ as

$\displaystyle R_{j,l}=$
$\displaystyle\left[\max\left(3n_{\phi},\frac{\pi-2\epsilon_{0}}{|\Delta_{j,l}|_{T}+2\epsilon_{0}}\right),\kappa_{\max}\right]$
$\displaystyle\cap\bigcup_{n\in\mathbb{N}}R_{j,l}^{(n)},$ (98)

where $R_{j,l}^{(n)}$ is the set of $k_{1}$ satisfying

$\bigg{|}k_{1}\Delta_{j,l}-2\pi n\bigg{|}\leq 4\epsilon_{0}(1+k_{1}),$ (99)

for some $n\in\mathbb{N}$. Solving this equation for $k_{1}$ yields

$\displaystyle R_{j,l}^{(n)}=$ $\displaystyle\left[\frac{2\pi
n-4\epsilon_{0}}{\Delta_{j,l}+4\epsilon_{0}},\frac{2\pi
n+4\epsilon_{0}}{\Delta_{j,l}-4\epsilon_{0}}\right].$ (100)

with length

$\left|R_{j,l}^{(n)}\right|=\frac{8\epsilon_{0}(2n\pi+\Delta_{j,l})}{\Delta_{j,l}^{2}-16\epsilon_{0}^{2}}.$
(101)

We can then bound

$\displaystyle\left|R_{j,l}\right|$
$\displaystyle\leq\left|\bigcup_{n=n_{\min}}^{n_{\max}}R^{(n)}_{j,l}\right|\leq\sum_{n=n_{\min}}^{n_{\max}}\left|R^{(n)}_{j,l}\right|$
$\displaystyle=\frac{8\epsilon_{0}}{\Delta_{j,l}^{2}-16\epsilon_{0}^{2}}\bigg{[}\pi(n_{\max}^{2}-n_{\min}^{2})$
$\displaystyle\qquad+(\Delta_{j,l}+\pi)(n_{\max}-n_{\min})$
$\displaystyle\qquad+2\pi n_{\min}+\Delta_{j,l}\bigg{]}.$ (102)

Here, $n_{\max}=n_{\rm max}(j,l)$ and $n_{\min}=n_{\min}(j,l)$ are the largest
and smallest indices of sets $R^{(n)}_{j,l}$ in Eq. (99) for which
$k_{1}\in[3n_{\phi},3n_{\phi}+1]$. Finding the minimal and the maximal value
for $n$ in Eq. (99) given the bounds on $k_{1}$ gives

$\frac{3n_{\phi}(\Delta_{j,l}-4\epsilon_{0})-4\epsilon_{0}}{2\pi}\leq n_{\rm
min}\leq n_{\max}\\\
\leq\frac{(3n_{\phi}+1)(\Delta_{j,l}+4\epsilon_{0})+4\epsilon_{0}}{2\pi}.$
(103)

As there are $n_{\phi}$ phases, the length of the total forbidden region
$\cup_{j,l}R_{j,l}$ is bounded from above by
$\frac{n_{\phi}^{2}}{2}\left|R_{j,l}\right|$. By plugging the bound for
$\epsilon_{0}$ in Eq. (95) and the bound for $\Delta_{j,l}$ in Eq. (97)
together into Eq. (102), one can verify that

$|R_{j,l}|\frac{n_{\phi}^{2}}{2}\leq\frac{-16-168n_{\phi}-621n_{\phi}^{2}+552n_{\phi}^{4}+14400n_{\phi}^{5}+44550n_{\phi}^{6}+40500n_{\phi}^{7}+73125n_{\phi}^{8}+270000n_{\phi}^{9}+202500n_{\phi}^{10}}{450n_{\phi}^{3}(-4-9n_{\phi}-100n_{\phi}^{3}-150n_{\phi}^{4}+1875n_{\phi}^{7})},$
(104)

which can be verified to be less than 0.5 for all $n_{\phi}$. Hence a random
choice for $k_{1}$ in the interval $[3n_{\phi},3n_{\phi}+1]$ gives at least a
50% chance to not land in the forbidden region. ∎

##  Appendix B Range of shifted phase estimates

In this Appendix we prove that when the unitary $U$ is shifted as in step 2 of
Alg. 4.1, and as long as the output of the phase extraction subroutine (Alg.
3.4) meets the promises given in Lem. 3.5, all phase estimates of $U^{k_{d}}$
will lie in the region for which Lem. 2.5 holds. This allows us to invoke Lem.
2.5 as required during Lem. 4.3 and Lem. 4.4. Note that if we were to shift
the spectrum such that the middle of the largest gap would sit at 0, we would
do $U\rightarrow Ue^{-i\zeta}$. However, the ‘stay-away-from-the-boundary’
condition of Lem. 2.5 is not symmetric, hence we shift by a different amount
which also depends on the error $\epsilon_{0}$ in the phase estimates.

######  Lemma B.1.

Let $\\{\phi_{j}\\}$ be the list of eigenphases of unitary $U$, and let
$n_{\phi}$, $\\{\tilde{\phi}_{l}^{(0)}\\}$, $\zeta$, $d_{\zeta}$ be as defined
in steps 1 and 2 of Alg. 4.1. Assume:

  * •

(Assumption 1a) For every phase $\phi_{j}$, there exists an estimate
$\tilde{\phi}_{l}^{(0)}$ such that $|\phi_{j}-\tilde{\phi}_{l}^{(0)}|_{T}\leq
2\epsilon_{0}$.

  * •

(Assumption 1b) For every estimate $\tilde{\phi}_{l}^{(0)}$, there exists a
phase $\phi_{j}$ such that $|\phi_{j}-\tilde{\phi}_{l}^{(0)}|_{T}\leq
2\epsilon_{0}$.

Then for all $k\geq 3n_{\phi}$, for the estimated shifted eigenphases
$\\{\tilde{\varphi}_{l}^{(0)}\\}$ it holds that

$\frac{\pi}{k}+16\epsilon_{0}\leq\tilde{\varphi}_{l}^{(0)}\leq\frac{\pi(2\lfloor
k\rfloor-1)}{k}-16\epsilon_{0}.$ (105)

which implies Eq. (8) for $\tilde{\varphi}_{l}^{(0)}=\phi$. In addition, the
eigenphases $\\{\varphi_{j}\\}$ of the shifted unitary
$Ue^{-i(\zeta+d_{\zeta}/2-8\epsilon_{0})}$ satisfy

$\frac{\pi}{k}+14\epsilon_{0}\leq\varphi_{j}\leq\frac{\pi(2\lfloor
k\rfloor-1)}{k}-14\epsilon_{0}.$ (106)

which again implies Eq. (8) for $\varphi_{j}=\phi$.

###### Proof.

Let

$\tilde{\varphi}^{(0)}_{l}=\left(\tilde{\phi}_{l}^{(0)}-\zeta-\frac{d_{\zeta}}{2}+8\epsilon_{0}\right)\textrm{
mod }2\pi.$ (107)

Let us first show that

$\tilde{\varphi}^{(0)}_{l}\in\left[\frac{\pi}{2n_{\phi}}+8\epsilon_{0},2\pi-\frac{3\pi}{2n_{\phi}}+8\epsilon_{0}\right].$
(108)

By definition of $\zeta$ (as the midway point in the largest gap) and
$d_{\zeta}$ (as half the largest gap) we have

$\left(\tilde{\phi}_{l}^{(0)}-\zeta-\frac{d_{\zeta}}{2}\right)\textrm{ mod
}2\pi\in\left[\frac{d_{\zeta}}{2},2\pi-\frac{3d_{\zeta}}{2}\right].$ (109)

We have $d_{\zeta}\geq\pi/n_{\phi}$ (with equality corresponding to $n_{\phi}$
uniformly distributed estimates). By Eq. (28) it follows that

$\epsilon_{0}\leq\frac{\pi}{48n_{\phi}},$ (110)

and thus $\epsilon_{0}<\frac{d_{\zeta}}{16}$, leading to Eq. (108). By the
assumptions the shifted phases $\varphi_{j}$ lie within $2\epsilon_{0}$ from
the estimates $\tilde{\varphi}^{(0)}_{l}$. Thus for each $\varphi_{j}$ there
exists a $\tilde{\varphi}^{(0)}_{l}$ such that

$\displaystyle\varphi_{j}-14\epsilon_{0}$
$\displaystyle\geq\tilde{\varphi}^{(0)}_{l}-16\epsilon_{0}$
$\displaystyle\geq\frac{\pi}{2n_{\phi}}-8\epsilon_{0}$
$\displaystyle\geq\frac{\pi}{2n_{\phi}}-\frac{8\pi}{48n_{\phi}}$
$\displaystyle=\frac{\pi}{3n_{\phi}}\geq\frac{\pi}{k}.$ (111)

and

$\displaystyle\varphi_{j}+14\epsilon_{0}$
$\displaystyle\leq\tilde{\varphi}^{(0)}_{l}+16\epsilon_{0}$ $\displaystyle\leq
2\pi-3\left(\frac{\pi}{2n_{\phi}}-8\epsilon_{0}\right)$ $\displaystyle\leq
2\pi-\frac{3\pi}{k}$ $\displaystyle\leq\frac{\pi(2\lfloor k\rfloor-1)}{k}.$
(112)

where we have used Eq. (108), Eq. (110), $k\geq 3n_{\phi}$ and $\lfloor
k\rfloor>k-1$. This implies Eq. (106) and also Eq. (105). ∎