# Fragmented imaginary-time evolution for early-stage quantum signal
processors

Thais L. Silva  Federal University of Rio de Janeiro, Caixa Postal 68528, Rio
de Janeiro, RJ 21941-972, Brazil  Márcio M. Taddei  Federal University of Rio
de Janeiro, Caixa Postal 68528, Rio de Janeiro, RJ 21941-972, Brazil  ICFO -
Institut de Ciencies Fotòniques, The Barcelona Institute of Science and
Technology, 08860, Castelldefels, Barcelona, Spain  Stefano Carrazza  TIF Lab,
Dipartimento di Fisica, Università degli Studi di Milano and INFN Sezione di
Milano, Milan, Italy  Quantum Research Centre, Technology Innovation
Institute, Abu Dhabi, UAE  Leandro Aolita  Federal University of Rio de
Janeiro, Caixa Postal 68528, Rio de Janeiro, RJ 21941-972, Brazil  Quantum
Research Centre, Technology Innovation Institute, Abu Dhabi, UAE

###### Abstract

Simulating quantum imaginary-time evolution (QITE) is a major promise of
quantum computation. However, the known algorithms are either probabilistic
(repeat until success) with unpractically small success probabilities or
coherent (quantum amplitude amplification) but with circuit depths and
ancillary-qubit numbers unrealistically large in the mid term. Our main
contribution is a new generation of deterministic, high-precision QITE
algorithms significantly more amenable experimentally. These are based on a
surprisingly simple idea: partitioning the evolution into several fragments
that are sequentially run probabilistically. This causes a huge reduction in
wasted circuit depth every time a run fails. Indeed, the resulting overall
runtime is asymptotically better than in coherent approaches and the hardware
requirements even milder than in probabilistic ones, remarkably. More
technically, we present two QITE-circuit sub-routines with excellent
complexity scalings. One of them is optimal in ancillary-qubit overhead (one
single ancillary qubit throughout) whereas the other one is optimal in runtime
for small inverse temperature or high precision. The latter is shown by noting
that the runtime saturates a cooling-speed limit that is the imaginary-time
counterpart of the no fast-forwarding theorem of real-time simulations, which
we prove. Moreover, we also make two technical contributions to the quantum
signal processing formalism for operator-function synthesis (on which our sub-
routines are based) that are useful beyond QITE. Our findings are specially
relevant for the early fault-tolerance stages of quantum hardware.

##  I Introduction

QITE is central not only to ground-state optimisations [1, 2, 3, 4, 5] but
also to partition-function estimation and quantum Gibbs-state sampling [6, 7,
8, 9, 10, 11, 12, 13, 12, 14, 15, 16, 17, 18], i.e. the task of preparing
thermal quantum states at tunable inverse temperature $\beta$. This is both
fundamentally relevant and useful for notable algorithmic applications. For
instance, even though approximating ground states of generic Hamiltonians is
not expected to be efficient even on a quantum computer – as it can solve QMA-
complete problems [19] –, significant speed-ups over classical simulations are
possible. This has motivated several ground-state cooling algorithms (with and
without QITE), specially for combinatorial optimisations [20, 21, 22, 2, 23]
or molecular electronic structures [24, 25, 26, 1]. On the other hand, Gibbs-
state samplers are used as main sub-routines for quantum semi-definite
programme solvers [13, 12, 14] or for training [27, 28, 29] quantum machine-
learning models [30, 31], e.g. Moreover, QITE also enables quantizations [2]
of the METTS or Lanczos algorithms, which directly simulate certain thermal
properties without Gibbs-state sampling.

Quantum Gibbs states can be approximated by quantum Metropolis Markov-chains
[8, 9] or by variational circuits trained to minimise the free energy [16],
e.g. However, the former involve deep and complex circuits, whereas the latter
are highly limited by the variational Ansatz. In turn, heuristic QITE
algorithms for ground-state optimisations exist [1, 2, 3, 4, 5, 32, 33, 34].
There, one simulates pure-state QITE with a unitary circuit that depends on
the input state, the Hamitonian, and $\beta$. For small-$\beta$ steps, one can
determine the circuits by measurements on the input state at each step and
classical post-processing. One possibility is to optimise a variational
circuit on the measured data [1], but this is again limited by the
expressivity of the Ansatz. Another possibility is to invert a linear system
generated from the measurements [2, 3, 4, 5], but the size of such system (as
well as the number of measurements required) is exponential in the number of
qubits, unless restrictive locality assumptions are made.

Figure 1: High-level schematics of our algorithms. a) QITE primitives: A
system register $\mathcal{S}$ carries the input state $|\Psi\rangle$, whereas
an ancillary register $\mathcal{A}$ is initialised in a computational-basis
state $|0\rangle$. A unitary transformation $U_{F_{\beta}(H)}$, composed of a
sequence $\\{V_{k}\\}_{\in[q]}$ of $q$ gates, with $[q]\coloneqq\\{1,\ldots
q\\}$, is applied and then the ancillas are measured. Each gate makes one
query to the Hamiltonian oracle (not shown). The specific choice of gates in
the sequence is determined by a powerful formalism for operator-function
design, known as quantum signal processing, which allows us to synthesise the
QITE propagator $F_{\beta}(H)\coloneqq e^{-\beta(H-\lambda_{\rm min})}$
efficiently. More precisely, conditioned on detecting $|0\rangle$ on the
ancillas, the desired state
$\frac{F_{\beta}(H)|\Psi\rangle}{\left\|F_{\beta}(H)|\Psi\rangle\right\|}$ is
output up to controllable error. We refer to the circuit generating
$U_{F_{\beta}(H)}$ as a QITE primitive. We develop two QITE primitives, one
for each oracle model. Both primitives feature excellent scalings of $q$ with
$\beta$ and the error. b) Master QITE algorithms: The post-selection
probability – given approximately by
$\left\|F_{\beta}(H)|\Psi\rangle\right\|^{2}$ – can decrease with $\beta$ very
fast. Hence, for high $\beta$, probabilistic approaches based on repeat-until-
success fail for the vast majority of trials, thereby producing an enormous
waste in runtime. In turn, coherent approaches based on quantum amplitude
amplification provide a close-to-quadratic runtime speed-up, but at the
expense of enormous circuit depths, among other disadvantages. In contrast, we
introduce a master algorithm that concatenates $r$ QITE fragments of inverse
temperatures $\\{\Delta\beta_{l}\\}_{l\in[r]}$, with
$\sum_{l\in[r]}\Delta\beta_{l}=\beta$ and $\beta_{l}<\beta$ for all $l\in[r]$.
Each fragment is successively run probabilistically and has both a success
probability significantly higher and a query complexity significantly lower
than that of the entire evolution run at once. Indeed, the ancillary
measurements in b) detect potential errors in the operator synthesis after
each fragment, instead of just once at the end in a). This ends up yielding an
enormous saving in overall runtime (even beating coherent approaches for high
$\beta$) while at the same time preserving all the practical advantages of
probabilistic approaches for experimental implementations.

The most general, guaranteed-precision QITE algorithms are based on unitary
circuits followed by ancillary-qubit post-selection [6, 11, 14, 15, 18]. These
circuits – to which we refer as _QITE primitives_ – are efficient in $\beta$
as well as in the target precision. However, due to the intrinsically
probabilistic post-selection, they must be applied multiple times – by what we
refer to as _master QITE algorithms_ – to obtain a deterministic output.
Repeat-until-success master algorithms apply the primitive in parallel (i.e.
in independent probabilistic runs), thereby not inducing any increase in
circuit depth. However their overall complexity is inversely proportional to
the post-selection probability. Instead, coherent master algorithms [6, 11,
14, 15], based on amplitude amplification [35], have close-to-quadratically
smaller overall complexity. However, they require enormous circuit depths and
significantly more ancillas. In addition, no fundamental efficiency limit for
generic QITE algorithms is known.

Here, we introduce two efficient QITE primitives as well as a practical master
QITE algorithm (see Fig. 1); and prove a universal lower bound for the
complexity of QITE primitives. The first primitive is designed for
Hamiltonians given in the well-known block-encoding oracle model, whereas the
second one for a simplified model of real-time evolution oracles involving a
single time. Both primitives feature excellent query complexity (number of
oracle calls) and ancillary-qubit overhead. In fact, for the first primitive
the complexity is sub-additive in $\beta$ and $\log(\varepsilon^{-1})$, with
$\varepsilon$ the tolerated error. This scaling saturates our universal bound
when $\beta\ll\log(\varepsilon^{-1})$. Hence Primitive 1 is optimal in that
regime, which, interestingly, turns out crucial for our master algorithm. In
contrast, Primitive 2’s complexity is multiplicative in $\beta$ and
$\log(\varepsilon^{-1})$, but it requires a single ancilla throughout and its
oracle significantly fewer gates. This is appealing for intermediate-scale
quantum hardware. In turn, our master QITE algorithm breaks the evolution into
small-$\beta$ fragments and runs each fragment’s primitive probabilistically.
Surprisingly, this yields an overall runtime competitive with – and, in the
relevant regime of high $\beta$, even better than – that of coherent
approaches while, at the same time, preserving all the advantages of
probabilistic ones for experimental feasibility.

More technically, our primitives are based on the quantum signal processing
(QSP) framework [36, 37, 15, 14], to which we also make two contributions
relevant on their own. The first one is a generic bound on the approximation
error of Hermitian-operator functions by their truncated Chebyshev series, for
any analytical real function. We apply this to the Jacobi-Anger expansion of
the exponential function [38], on which we base our first primitive. This
expansion has been used to obtain optimal-scaling algorithms for real-time
evolution (RTE) [36, 37, 15]. Moreover, it has also been used for partition-
function estimation [18]. However, there, the different terms in the expansion
are synthesized separately and their coefficients are then used to classically
combine the measurement results. Our first primitive directly synthesizes the
whole Jacobi-Anger expansion of the QITE propagator. Our second contribution
(see also [39]) is a Fourier approximation of analytical real functions of
Hermitian operators for RTE oracles, on which we base our second primitive.
This is a novel QSP variant superior to previous Fourier-based approaches [14]
in that it requires a single ancilla, instead of multiple ones. This variant
does not require qubitization [37]. In turn, the universal runtime bound we
prove can be seen as an imaginary-time counterpart of the no fast-forwarding
theorem for RTE [40, 41, 42].

Finally, the complexity of our master algorithm depends on the _fragmentation
schedule_ , i.e. number $r$ of fragments and their relative sizes. On one
hand, for Primitive 1, we rigorously prove that, from a critical inverse
temperature $\beta_{\rm c}=\mathcal{O}(2^{N/2}\,N)$ on, the runtime is lower
than that with coherent QITE. This is shown by explicitly constructing
schedules with only $r=2$ fragments that do the job, remarkably. On the other
hand, that fragmented QITE outperforms coherent QITE is also observed for both
primitives through extensive numerical evidence. More precisely, we study the
overall runtime as a function of $\beta$ and $\varepsilon$, up to $N=15$
qubits, and for numerically-optimized schedules. These experiments involve
random instances of Hamiltonians encoding four computationally hard classes of
problems: Ising models associated to the $i$) MaxCut and $ii$) weighted MaxCut
problems [20, 21, 22]; $iii$) restricted quantum Boltzmann machines
(transverse-field Ising models) [30, 31]; and $iv$) a quantum generalization
(fully-connected Heisenberg models) of the Sherrington-Kirkpatrick model [43,
44] for spin glasses. We see a clear trend whereby, from $\beta_{\rm
c}=\mathcal{O}\big{(}2^{N/2}\big{)}$ on, fragmentation outperforms coherent
QITE for both primitives, for an optimal number of fragments $r\lesssim 6$.
The obtained values for $\beta_{\rm c}$ imply that our algorithm outperforms
coherent QITE in the computationally hardest range of $\beta$, particularly
relevant for Hamiltonians with an exponentially small spectral gap [45, 46].
Moreover, impressively, such advantages are attained at no cost in circuit
depth or number of ancillas, which are identical to those of probabilistic
QITE.

The paper is organised as follows. In Sec. II we introduce basic definitions
and notation. In Sec. III we present our main results: In Sec. III.1, we state
two theorems giving the complexities of the primitives. In Sec. III.2, we
prove a universal lower-bound on the complexity of generic QITE primitives. In
Sec. III.3, we present our master algorithm and prove its correctness and
complexity. In Sec. III.4, we present exhaustive numerical results showing the
advantages of fragmentation over quantum amplitude amplification for quantum
Gibbs-state sampling. Finally, in Sec. IV, we give our concluding remarks and
outlook. The technical details of our primitives and theorem proofs are left
for the Methods, in Secs. V and VI, where we also review the basics of QSP.

##  II Preliminaries

We consider an $N$-qubit system $\mathcal{S}$, of Hilbert space
$\mathbb{H}_{\mathcal{S}}$. QITE with respect to a Hamiltonian $H$ on
$\mathbb{H}_{\mathcal{S}}$ and over an imaginary time $-i\,\beta$ is
represented by the non-unitary operator $e^{-\beta H}$. This can be simulated
via post-selection with a unitary operator $U$ that encodes $e^{-\beta H}$ in
one of its matrix blocks [6, 11, 14, 15, 18]. We denote by
$\mathbb{H}_{\mathcal{A}}$ the Hilbert space of an ancillary register
$\mathcal{A}$, by
$\mathbb{H}_{\mathcal{SA}}\coloneqq\mathbb{H}_{\mathcal{S}}\otimes\mathbb{H}_{\mathcal{A}}$
the joint Hilbert space of $\mathcal{S}$ and $\mathcal{A}$, and by
$\left\|A\right\|$ the spectral norm of an operator $A$. The following
formalizes the encoding.

######  Definition 1.

(Block encodings). For sub-normalization $0\leq\alpha\leq 1$ and tolerated
error $\varepsilon>0$, a unitary operator $U_{A}$ on
$\mathbb{H}_{\mathcal{SA}}$ is an $(\varepsilon,\alpha)$-block-encoding of a
linear operator $A$ on $\mathbb{H}_{\mathcal{S}}$ if $\left\|\alpha\,A-\langle
0|\,U_{A}\,|0\rangle\right\|\leq\,\varepsilon$, for some
$|0\rangle\in\mathbb{H}_{\mathcal{A}}$. For $\varepsilon=0$ and
$(\varepsilon,\alpha)=(0,1)$ we use the short-hand terms perfect
$\alpha$-block-encoding and perfect block-encoding, respectively.

E.g., if $U_{A}$ is a perfect $\alpha$-block-encoding of $A$, measuring
$|0\rangle\in\mathbb{H}_{\mathcal{A}}$ on
$U_{A}|\Psi\rangle|0\rangle\in\mathbb{H}_{\mathcal{SA}}$, for any
$|\Psi\rangle\in\mathbb{H}_{\mathcal{S}}$, leaves $\mathcal{S}$ in the state
$\frac{A|\Psi\rangle}{\left\|A|\Psi\rangle\right\|}$. The probability of that
outcome is $\alpha^{2}\left\|A|\Psi\rangle\right\|^{2}$. Note that, since
$\left\|U_{A}\right\|=1$, a perfect $\alpha$-block-encoding is possible only
if $\alpha\left\|A\right\|\leq 1$. Hence, $\alpha$ allows one to encode
matrices even if their norm is greater than $1$. Typically, however, one
wishes $\alpha$ as high as possible, to avoid unnecessary reductions in post-
selection probability.

Our algorithms admit two types of oracle as input. The first one is based on
perfect block-encodings of $H$ and therefore requires $\left\|H\right\|\leq
1$. If $\left\|H\right\|>1$, however, the required normalisation can be
enforced by a simple spectrum rescaling. More precisely, for $\lambda_{-}$ and
$\lambda_{+}$ arbitrary lower and upper bounds, respectively, to the minimal
and maximal eigenvalues of $H$, $\lambda_{\rm min}$ and $\lambda_{\rm max}$,
the rescaled Hamiltonian
$H^{\prime}\coloneqq\frac{H-\bar{\lambda}\mathds{1}}{\Delta\lambda}$ fulfils
$\left\|H^{\prime}\right\|\leq 1$ by construction, with the short-hand
notation $\bar{\lambda}\coloneqq\frac{\lambda_{+}+\lambda_{-}}{2}$ and
$\Delta\lambda\coloneqq\frac{\lambda_{+}-\lambda_{-}}{2}$. Then, by
correspondingly rescaling the inverse temperature as
$\beta^{\prime}\coloneqq\Delta\lambda\,\beta$, one obtains the propagator
$e^{-\beta^{\prime}H^{\prime}}$, which induces the same physical
transformation as $e^{-\beta H}$. Hence, from now on, without loss of
generality we assume throughout that $\left\|H\right\|\leq 1$, i.e. that
$-1\leq\lambda_{\rm min}\leq\lambda_{\rm max}\leq 1$.

We are now in a good position to define our first oracle, $O_{1}$, which is
the basis of our first primitive, $P_{1}$. We denote by $\mathcal{A}_{1}$ the
entire ancillary register needed for $P_{1}$ and by
$\mathcal{A}_{O_{1}}\subset\mathcal{A}_{1}$ the specific ancillary qubits
required to implement $O_{1}$.

######  Definition 2.

(Block-encoding Hamiltonian oracles). We refer as a block-encoding oracle for
a Hamiltonian $H$ on $\mathbb{H}_{\mathcal{S}}$ to a controlled unitary
operator $O_{1}$ on $\mathbb{H}_{\mathcal{SA}_{O_{1}}}$ of the form
${O_{1}=U_{H}\otimes|0\rangle\langle 0|+\mathds{1}\otimes|1\rangle\langle
1|}$, where $\mathds{1}$ is the identity operator on
$\mathbb{H}_{\mathcal{S}}$, $\\{|0\rangle,|1\rangle\\}$ a computational basis
for the control qubit, and $U_{H}$ a perfect block encoding of $H$.

This is a powerful oracle paradigm used both in QITE [11, 14, 15, 18] and
real-time evolution [47, 36, 37, 15]. It encompasses, e.g., Hamiltonians given
by linear combinations of unitaries, $d$-sparse Hamiltonians (i.e. with at
most $d$ non-null matrix entries per row), and Hamiltonians given by states
[37]. Its complexity depends on $H$, but highly efficient implementations are
known. E.g., for $H$ a linear combination of $m$ unitaries, each one requiring
at most $c$ two-qubit gates, $O_{1}$ can be implemented with
$|\mathcal{A}_{O_{1}}|=\mathcal{O}(\log_{2}m)$ ancillary qubits and gate
complexity (i.e. total number of two-qubit gates)
$g_{O_{1}}=\mathcal{O}\big{(}m(c+\log_{2}m)\big{)}$ [47, 37].

The second oracle model that we consider encodes $H$ through the real-time
unitary evolution it generates.

######  Definition 3.

(Real-time evolution Hamiltonian oracle). We refer as a real-time evolution
oracle for a Hamiltonian $H$ on $\mathbb{H}_{\mathcal{S}}$ at a time
$t\in\mathbb{R}$ to a controlled-$e^{-itH}$ gate
$O_{2}=\mathds{1}\otimes|0\rangle\langle 0|+e^{-itH}\otimes|1\rangle\langle
1|$.

This is a simplified version of the models of [6, 14], e.g. There, controlled
real-time evolutions at multiple times are required, thus involving multiple
ancillas. In contrast, $O_{2}$ involves a single real time, so the ancillary
register $\mathcal{A}_{O_{2}}$ consists of $|\mathcal{A}_{O_{2}}|=1$ single
qubit (the control). In fact, we show below that no other ancilla is needed
for our second primitive, $P_{2}$, i.e. $\mathcal{A}_{2}=\mathcal{A}_{O_{2}}$.
This is advantageous for near-term implementations. There, one may for
instance apply product formulae [48, 49] to implement $O_{2}$ with gate
complexities $g_{O_{2}}$ that, for intermediate-scale systems, can be
considerably smaller than for $O_{1}$ (see App. A). Furthermore, this oracle
is also relevant to hybrid analogue-digital platforms, for which QSP schemes
have already been studied [50].

QITE algorithms based on post-selection rely on a unitary quantum circuit to
simulate a block encoding of the QITE propagator. We refer to such circuits as
QITE primitives.

######  Definition 4.

(QITE primitives). Let $\beta\geq 0$, $\varepsilon^{\prime}\geq 0$, and
$\alpha\leq 1$. A $(\beta,\varepsilon^{\prime},\alpha)$-QITE-primitive of
query complexity $q(\beta,\varepsilon^{\prime},\alpha)$ is a circuit $P$, with
$q(\beta,\varepsilon^{\prime},\alpha)$ calls to an oracle $O$ for $H$ or its
inverse $O^{\dagger}$, that generates an
$(\varepsilon^{\prime},\alpha)$-block-encoding $U_{F_{\beta}(H)}$ of
$F_{\beta}(H)\coloneqq e^{-\beta(H-\lambda_{\rm min})}$, for all $H$.

Note that $P$ is Hamiltonian agnostic, i.e. it admits any $H$ provided it is
properly encoded in the corresponding oracle. The factor
$e^{-\beta\lambda_{\rm min}}$ implies that $\|F_{\beta}(H)\|=1$, thus
maximizing the post-selection probability. However, if $\lambda_{\rm min}$ is
unknown, one can replace it by a suitable lower bound $\lambda_{-}\geq-1$.
This introduces only a constant sub-normalisation. In turn, the query
complexity is the gold-standard figure of merit for efficiency of oracle-based
algorithms. It quantifies the runtime of $P$ relative to that of an oracle
query. In fact, $P$ is time-efficient if its query complexity and gate
complexity per query $g_{P}$ are both in $\mathcal{O}\big{(}{\rm
poly}(N,\beta,1/\varepsilon^{\prime},\alpha)\big{)}$.

Importantly, normalisation causes the post-selection probability
$p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)$ of $P$ (on an input state
$|\Psi\rangle$) to propagate onto the error $\varepsilon$ in the output state,
making the latter in general greater than $\varepsilon^{\prime}$. The exact
dependance of $\varepsilon$ on $\varepsilon^{\prime}$ is dictated by
$p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)$. However, if
$\varepsilon^{\prime}\leq\,\varepsilon\,\sqrt{p_{\Psi}(\beta,\alpha)}/2$, with
$p_{\Psi}(\beta,\alpha)\coloneqq
p_{\Psi}(\beta,0,\alpha)=\alpha^{2}\,\left\|F_{\beta}(H)|\Psi\rangle\right\|^{2}$,
the output-state error is in $\mathcal{O}(\varepsilon)$ (see App. B). In turn,
the primitives must be incorporated into master procedures to attain
deterministic QITE. We call these master QITE algorithms.

######  Definition 5.

(Master QITE algorithms). Given $\varepsilon\geq 0$, $\beta\geq 0$,
$|\Psi\rangle\in\mathbb{H}_{\mathcal{S}}$, and
$(\beta^{\prime},\varepsilon^{\prime},\alpha^{\prime})$-QITE-primitives
$P_{\beta^{\prime},\varepsilon^{\prime},\alpha^{\prime}}$ querying oracles for
a Hamiltonian $H$, a $(\beta,\varepsilon)$-master-QITE-algorithm for $H$ on
$|\Psi\rangle$ is a procedure that outputs the state
$\frac{F_{\beta}(H)|\Psi\rangle}{\left\|F_{\beta}(H)|\Psi\rangle\right\|}$ up
to trace-distance error $\varepsilon$ with unit probability. Its overall query
complexity $Q(\beta,\varepsilon)$ is the sum over the query complexities of
each $P_{\beta^{\prime},\varepsilon^{\prime},\alpha^{\prime}}$ applied.

Until now, two variants of master QITE algorithms had been reported,
probabilistic and coherent. The former leverage repeat-until-success: perform
$\mathcal{O}(1/p_{\Psi}(\beta,\varepsilon^{\prime},\alpha))$ trials of
$P_{\beta,\varepsilon^{\prime},\alpha}$ (on independent systems) until getting
the desired output. In contrast, the latter are based on quantum amplitude
amplification [35]. There, $P_{\beta,\varepsilon^{\prime},\alpha}$ is
incorporated into a unitary amplification engine that is sequentially applied
(on the same system)
$\mathcal{O}\big{(}\sqrt{1/p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)}\big{)}$
times. Conveniently, $1/p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)$ is well-
approximated by $1/p_{\Psi}(\beta,\alpha)$ for $\varepsilon\ll 1$ (see App.
C). Hence, the overall complexity of both variants is given by the unified
expression

$Q_{\kappa}(\beta,\varepsilon,\alpha)=\mathcal{O}\bigg{(}\frac{1}{{\big{(}p_{\Psi}(\beta,\alpha)}\big{)}^{{\mu_{\kappa}}}}\,q(\beta,\varepsilon^{\prime},\alpha)\bigg{)},$
(1)

where $\kappa=$ prob or coh for probabilistic or coherent schemes,
respectively, $\mu_{\rm prob}=1$, $\mu_{\rm coh}=1/2$,
$p_{\Psi}(\beta,\alpha)\coloneqq\alpha^{2}\left\|F_{\beta}(H)|\Psi\rangle\right\|^{2}$,
and $\varepsilon^{\prime}\coloneqq\varepsilon\sqrt{p_{\Psi}(\beta,\alpha)}/2$.

Since $p_{\Psi}(\beta,\alpha)$ can decrease with $N$ exponentially, the
quadratic advantage in $1/p_{\Psi}(\beta,\alpha)$ of coherent approaches is
highly significant. However, while probabilistic schemes make no assumption on
the input-state source, coherent ones require a unitary oracle that prepares
$|\Psi\rangle$ from a reference state. Most importantly, coherent algorithms
have a circuit depth $\mathcal{O}(\sqrt{1/p_{\Psi}(\beta,\alpha)})$ times
greater than in probabilistic ones and require $\mathcal{O}(N)$ extra
ancillas. This makes coherent schemes impractical for intermediate-scale
quantum devices.

##  III Results

We first discuss the primitives, then the universal complexity lower bound,
and the master algorithm at last.

###  III.1 Quantum imaginary-time evolution primitives

We introduce two QITE primitives. Both of them possess the basic structure
shown in Fig. 1 a), where a sequence of gates $\\{V_{k}\\}_{\in[q]}$, with
$[q]\coloneqq\\{1,2,\ldots q\\}$, generates an approximate block-encoding of
the QITE propagator $F_{\beta}(H)$. The approximation consists of truncating
an expansion of the exponential function at finite order $q$. The gates are
determined by such expansion using quantum signal processing [36, 37, 15, 14]
(see Sec. V). Conceptually, the two primitives differ in the kind of expansion
and the type of oracle. Their circuit descriptions are given in the Methods.

Primitive 1 is the circuit output by Algorithm 2, in Sec. V.2, on inputs
$f(\lambda)=e^{-\beta(\lambda-\lambda_{\rm min})}$, with $\beta\geq 0$,
$\varepsilon^{\prime}\geq 0$, a block-encoding oracle $O_{1}$, and its inverse
$O_{1}^{\dagger}$. Recall that $|\mathcal{A}_{O_{1}}|$ is the ancillary-
register size and $g_{O_{1}}$ is the gate complexity of $O_{1}$. In Sec.
V.2.1, we prove the following.

######  Theorem 1.

(QITE primitive for block-encoding oracles). Circuit $P_{1}$ constructed by
Alg. 2 for $f(\lambda)=e^{-\beta(\lambda-\lambda_{\rm min})}$ is a
$(\beta,\varepsilon^{\prime},1)$-QITE-primitive with query complexity

$q_{{1}}(\beta,\varepsilon^{\prime})=\mathcal{O}\left(\frac{e\,\beta}{2}+\frac{\ln(1/\varepsilon^{\prime})}{\ln\left[e+2\ln(1/\varepsilon^{\prime})/(e\,\beta)\right]}\right),$
(2)

$|\mathcal{A}_{1}|=|\mathcal{A}_{O_{1}}|+1$ total ancillary qubits, and gate
complexity per query $g_{P_{1}}=\mathcal{O}(g_{O_{1}}+|\mathcal{A}_{O_{1}}|)$.
Moreover, Algorithm 2 has classical run-time $\mathcal{O}({\rm
poly}\big{(}q(\beta,\varepsilon^{\prime})\big{)}$.

The highlights of Eq. (2) are its sub-logarithmic dependence on $\varepsilon$
and its sub-additivity in $\beta$ and $\ln(1/\varepsilon^{\prime})$. We note
that a QITE primitive from [15] works for the same oracle model and has
complexity upper-bounded by
$\mathcal{O}\big{(}\sqrt{2\,\max[e^{2}\,\beta,\ln(2/\varepsilon^{\prime})]\,\ln(4/\varepsilon^{\prime})}\big{)}$.
Even though that is asymptotically better in $\beta$ than the bound in Eq.
(2), it underperforms it for all $\beta\lesssim 8\ln(4/\varepsilon^{\prime})$.
On the other hand, the complexity obtained in Ref. [18] lacks the denominator
in the second term of Eq. (2), also making the former worse than the latter
for low $\beta$ or high $\ln(1/\varepsilon^{\prime})$. While both the scalings
of Refs. [15, 18] tend to
$\mathcal{O}\big{(}\ln(1/\varepsilon^{\prime})\big{)}$ for $\beta\to 0$, Eq.
(2) tends to zero. In fact, in Sec. III.2 we show that the scaling in Eq. (2)
becomes optimal as $\beta$ decreases relative to
$\ln(1/\varepsilon^{\prime})$. We stress that the latter regime is crucial for
the master algorithm of Sec. III.3, whose first fragments require precisely
low inverse temperatures and high precisions. In turn, in the opposite regime
of high $\beta$, preliminary numerical observations [51] suggest that the
asymptotic scaling of the exact value of $q_{1}$ could actually be as good as
$q_{1}(\beta,\varepsilon^{\prime})=\mathcal{O}\big{(}\sqrt{\beta\,\ln(1/\varepsilon^{\prime})}\big{)}$,
i.e. similar to that from [15].

Primitive 1 is based on the Jacobi-Anger expansion [38]. This immediately
gives a truncated Chebyshev-polynomial series [52, 53] for $F_{\beta}(H)$,
which can be synthesized with quantum signal processing (see Sec. V). In Refs.
[42, 36, 37, 15, 18], the truncation error by this expansion was bounded [42]
using properties of the Bessel functions [54]. In contrast, here, we use a
generic upper bound (Lemma 9, in Sec. V.2) for truncated Chebyshev series of
arbitrary Hermitian-operator functions. This gives the same bound for
$F_{\beta}(H)$ as [42] but holds for any analytical real function, and is
hence useful for operator-function design in general.

Primitive 2 is the circuit output by Algorithm 3, in Sec. V.3, on inputs
$f(\lambda)=e^{-\beta(\lambda-\lambda_{\rm min})}$, with $\beta\geq 0$,
$\varepsilon^{\prime}\geq 0$, a real-time evolution oracle $O_{2}$ at time
$t=\frac{\pi}{2}\big{(}1+\frac{\gamma}{\beta}\big{)}^{-1}$, and its inverse
$O_{2}^{\dagger}$, with $\gamma\geq 0$ a parameter that will determine the
subnormalization. In Sec. V.3.1, we prove the following.

######  Theorem 2.

(QITE primitive for single real-time evolution oracles). Circuit $P_{2}$
constructed by Alg. 3 for $f(\lambda)=e^{-\beta(\lambda-\lambda_{\rm min})}$
is a $(\beta,\varepsilon^{\prime},\alpha)$-QITE-primitive with
$\alpha=e^{-\beta(1+\lambda_{\rm min})-\gamma}$ and query complexity

$q_{{2}}(\beta,\varepsilon^{\prime},\alpha)=\mathcal{O}\big{(}(\beta/\gamma+1)\ln(4/\varepsilon^{\prime})\big{)},$
(3)

$|\mathcal{A}_{2}|=1$, and $g_{P_{2}}=g_{O_{2}}+\mathcal{O}(1)$. Moreover,
Algorithm 3 has classical runtime $\mathcal{O}({\rm
poly}\big{(}q(\beta,\varepsilon^{\prime},\alpha)\big{)}$.

In contrast to Eq. (2), the relation between $\beta$ and
$\ln(1/\varepsilon^{\prime})$ in Eq. (3) is multiplicative. However, in
return, $P_{2}$ requires $|\mathcal{A}_{2}|=1$ ancillary qubit throughout,
remarkably. This is a drastic reduction with respect to not only block-encoded
oracles but also other oracles based on real-time evolution [6, 14]. In fact,
$|\mathcal{A}_{2}|=1$ is the minimum possible, because, since $F_{\beta}(H)$
is non-unitary, at least 1 ancilla is needed to block-encode it. Moreover, the
scaling of $g_{P_{2}}$ is optimal too, as it adds only a small, constant
number of gates per query to the intrinsic gate complexity $g_{O_{2}}$ of the
oracle. These features make $P_{2}$ specially well-suited for near-term
devices.

Our algorithms support any $\lambda_{\rm min}\in[-1,1]$. For $P_{2}$, this is
reflected by the sub-normalization factor $e^{-\beta(1+\lambda_{\rm min})}$,
which decreases as $\lambda_{\rm min}$ departs from $-1$. In turn, the other
factor, $e^{-\gamma}$, arises from the Gibbs phenomenon of Fourier series. The
theorem holds for all $\gamma\geq 0$, allowing one to trade success
probability for query complexity. However, in App. D, we obtain the optimal
$\gamma$ minimising Eq. (1), for both probabilistic ($\kappa=$prob) and
coherent ($\kappa=$coh) algorithms. Conveniently, for $\varepsilon^{\prime}\ll
1$, this is well-approximated by

$\gamma_{\kappa}(\beta)\coloneqq\frac{\beta}{2}\left(\sqrt{1+\frac{2}{\mu_{\kappa}\,\beta}}-1\right).$
(4)

Importantly, rather than a peculiarity of $P_{2}$, the favourable scalings of
$|\mathcal{A}_{2}|$ and $g_{P_{2}}$ are generic features of the type of
operator-function design behind it: An optimised Fourier-approximation
algorithm for arbitrary analytical real functions of Hermitian operators (see
Sec. V.3). This is a novel quantum-signal-processing variant interesting on
its own. Since it is based based on real-time evolution oracles, it requires
no qubitization [37]. This explains the scaling of $g_{P_{2}}$. In addition,
in contrast to other Fourier-based approaches, ours involves a single real
time instead of an error-dependent number of them [14]. Here, the different
frequency terms in the Fourier approximation are generated through the action
of $O_{2}$ multiple times, instead of using ancillas to make a linear
combination of several unitaries, one unitary for each frequency term. This
explains the scaling of $|\mathcal{A}_{2}|$.

Finally, Theos. 1 and 2 can be straightforwardly extended to the realistic
case of approximate oracles: In App. A, we show (for generic analytical
operator functions) that it suffices to take the oracle error (deviation from
an ideal oracle) as $\varepsilon^{\prime}_{O}=O(\varepsilon^{\prime}/q)$ to
keep the primitive’s error in $O(\varepsilon^{\prime})$.

###  III.2 Cooling-speed limits for oracle-based QITE algorithms

The most challenging applications of QITE involve small post-selection
probabilities, decreasing exponentially in $N$ in the worst cases. In an
effort to reduce the overall complexity [see Eq. (1)], this has fueled a long
race [6, 7, 11, 12, 14, 15] to improve $q(\beta,\varepsilon^{\prime})$, going
from the seminal $\mathcal{O}\big{(}\beta\,{\rm
poly}(1/\varepsilon^{\prime})\big{)}$ of [6] to the recent
$\mathcal{O}\big{(}\sqrt{2\max[e^{2}\,\beta,\ln(2/\varepsilon^{\prime})]\,\ln(4/\varepsilon^{\prime})}\big{)}$
of [15] or the additive scaling of Eq. (2). However, to our knowledge, no
runtime limit for QITE simulations has been established. This contrasts with
RTE, where fundamental runtime lower bounds are given by the “no-fast-
forwarding theorem” [40, 41, 42]. These are saturated by optimal RTE
algorithms [36, 37, 15]. Here we derive an analogous bound for imaginary time,
which we call cooling-speed limit in allusion to the use of QITE to cool
systems down to their ground state.

More precisely, in Sec. VI, we prove a universal efficiency limit for QITE
primitives based on block-encoded oracles. This is convenient as it directly
applies to our primitive with lowest query complexity, i.e. $P_{1}$.

######  Theorem 3 (Imaginary-time no-fast-forwarding theorem).

Let $\beta>0$ and $0<\varepsilon^{\prime}<\alpha/2$. Then, any
$(\beta,\varepsilon^{\prime},\alpha)$-QITE-primitive querying block-encoding
Hamiltonian oracles has query complexity at least
$q_{\min}(\beta,\varepsilon^{\prime},\alpha)\geq\tilde{q}$, where
$\tilde{q}\in\mathbb{R}_{>0}$ is the unique solution to the equation

$\left|\frac{1-e^{-\frac{\beta}{4\tilde{q}}}}{2}\right|^{2\,\tilde{q}}=\frac{2\,\varepsilon^{\prime}}{\alpha}\
.$ (5)

Even though the bound is only given implicitly, interesting conclusions can
readily be drawn. First, for any fixed $\beta$, the left-hand-side of Eq. (5)
decreases monotonically with $\tilde{q}$ (therefore the uniqueness of the
solution). Second, for any fixed $\varepsilon^{\prime}$ and $\alpha$,
$\tilde{q}$ grows monotonically with $\beta$. Third, and most important, Eq.
(5) is approximated by
$\big{(}\frac{\beta}{{8}\,\tilde{q}}\big{)}^{2\tilde{q}}=2\,\varepsilon^{\prime}/\alpha$
for $\beta\ll\tilde{q}$, as a Taylor expansion shows. The latter equation has
a known explicit solution [15], which, for $\alpha=1$, is given by Eq. (2).
Hence, for $\beta/\tilde{q}\to 0$, Eq. (2) tends to the optimal scaling. Note
that $\beta\ll\tilde{q}$ is equivalent to the first term Eq. (2) being much
smaller than its second term, which in turn implies that
$\varepsilon^{\prime}$ should be exponentially small in $\beta$. Thus, $P_{1}$
is close to optimal for small inverse temperatures or high precisions.
Interestingly, this is the regime at which the first fragments of our master
algorithm operate, as we see next.

###  III.3 Fragmented master QITE algorithm

Our master algorithm relies on the basic identity
$F_{\beta}(H)=\prod_{l=1}^{r}F_{\Delta\beta_{l}}(H)$ to partition the
evolution into $r\in\mathbb{N}$ fragments of inverse temperatures
$S_{r}\coloneqq\\{\Delta\beta_{l}>0\\}_{l\in[r]}$, such that
$\sum_{l\in[r]}\Delta\beta_{l}=\beta$. We refer to $S_{r}$ as the
_fragmentation schedule_. For each $l$, the algorithm repeats until success a
$(\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l})$-QITE-primitive
$P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}$ on the output state
$|\Psi_{l-1}\rangle$ of the $(l-1)$-th step, with $\varepsilon^{\prime}_{l}$
given in Eq (6). That is, if the post-selection succeeds, its output state
$|\Psi_{l}\rangle$ is input into the $(l+1)$-th fragment. Else, the algorithm
starts all over from the first fragment on
$|\Psi_{0}\rangle\coloneqq|\Psi\rangle$, until $|\Psi_{l-1}\rangle$ is
prepared and the $l$-th fragment can be run again. Alternatively, the
measurement on $\mathcal{A}$ after each fragment can be seen as monitoring
that the correct block of $U_{F_{\Delta\beta_{l}}(H)}$ is applied on each
$|\Psi_{l-1}\rangle$, in contrast to the single error detection after
$U_{F_{\beta}(H)}$ in the probabilistic master algorithm (see Fig. 1). Note
that the total number of trials (i.e. preparations of $|\Psi\rangle$)
coincides with the number of repetitions of the first fragment. The following
pseudocode summarizes all the algorithm:

1

input : $|\Psi\rangle$, $\beta\geq 0$, $\varepsilon\geq 0$, $S_{r}$, and QITE
primitives
$\\{P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}\\}_{l\in[r]}$
querying an oracle for $H$ (for $\Delta\beta_{l}\in S_{r}$,
$\varepsilon^{\prime}_{l}$ given in Eq (6), and $\alpha_{l}>0$).

output : $F_{\beta}(H)|\Psi\rangle/||F_{\beta}(H)|\Psi\rangle||$ up to trace-
distance error $\mathcal{O}(\varepsilon)$.

2 Set $l=1$, initialize $\mathcal{S}$ in the input state $|\Psi\rangle$ and
$\mathcal{A}$ in the reference state $|0\rangle$;

3 while  _$l\leq r$_ do

4 apply the circuit $P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}$
on $\mathcal{S}$ and $\mathcal{A}$;

5 measure $\mathcal{A}$ in computational basis;

6 if  _the outcome is $|0\rangle$_ then $l\to l+1$;

7 else break loop and go back to line 1;

8

9 end while

Algorithm 1 Fragmented QITE

The correctness and complexity of Algorithm 1 are established by the following
theorem, proven in App. E.

######  Theorem 4.

(Fragmented master QITE algorithm). If

$\varepsilon^{\prime}_{l}\leq\left\\{\begin{array}[]{lr}\frac{\varepsilon\,\prod_{k=1}^{r}\alpha_{k}}{2\times
4^{r-1}}\sqrt{p_{\Psi}(\beta)}&\text{ if }l=1,\\\
\frac{\varepsilon\,\prod_{k=l}^{r}\alpha_{k}}{4^{r-l+1}}\sqrt{\frac{p_{\Psi}(\beta)}{p_{\Psi}(\beta_{l-1})}}&\text{
if }l>1,\\\ \end{array}\right.$ (6)

for all $l\in[r]$, Algorithm 1 is a
$\big{(}\beta,\mathcal{O}(\varepsilon)\big{)}$-master-QITE-algorithm for $H$
on $|\Psi\rangle$ with average query complexity

${Q_{S_{r}}(\beta,\varepsilon)=\sum_{l=1}^{r}n_{l}\,q(\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l})},$
(7)

where
$n_{l}\coloneqq\frac{p_{\Psi}(\beta_{l-1})}{p_{\Psi}(\beta)\prod_{k=l}^{r}\alpha^{2}_{k}}$
is the average number of times that
$P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}$ is run, with
$\beta_{0}\coloneqq 0$, $\beta_{l}\coloneqq\sum_{k=1}^{l}\Delta\beta_{k}$ for
all $l\in[r]$, and
$p_{\Psi}(\beta_{l-1})\coloneqq\|F_{\beta_{l-1}}(H)|\Psi\rangle\|^{2}$ for all
$l\in[r+1]$.

Figure 2: Intuition behind the complexity reduction by fragmentation. The
overall complexity of the probabilistic master algorithm is dominated by the
area of the yellow rectangle. In contrast, the corresponding complexity of the
fragmented algorithm (here, for the exemplary case of $r=3$ fragments) is
dominated by the area of the blue-shaded rectangles. Up to logarithmic
corrections in the precision, the cumulative width of the blue-shaded
rectangles coincides with the width of the yellow one, of order $\beta$. In
contrast, while the height of the yellow rectangle is of order
$1/p_{\Psi}(\beta)$, the height of the blue-shaded ones decreases from order
$1/p_{\Psi}(\beta)$ till order $p_{\Psi}(\beta_{r-1})/p_{\Psi}(\beta)$, making
the blue-shaded area smaller than the yellow one. For high-enough $\beta$, the
reduction can be so strong that the complexity of the fragmented algorithm can
reach even that of the coherent algorithm (not shown). This intuition is
rigorously proven for Primitive 1 (in Theorem 5) and numerically verified to
exhaustion for both Primitives 1 and 2 (in Sec. III.4).

We note that, for Primitive 1, the average total number of trials coincides
with that of the probabilistic algorithm: $n_{1}=1/p_{\Psi}(\beta)=:n_{\rm
prob}$. This is important because the probabilistic algorithm consumes
$q_{1}(\beta,\varepsilon^{\prime})$ queries per trial, successful or not. In
contrast, the fragmented one consumes per trial
$q_{1}(\Delta\beta_{1},\varepsilon^{\prime}_{1})$ queries, plus
$q_{1}(\Delta\beta_{2},\varepsilon^{\prime}_{2})$ queries only if the first
post-selection succeeds, plus
$q_{1}(\Delta\beta_{3},\varepsilon^{\prime}_{3})$ queries only if the second
one succeeds too, and so on. Hence, the total waste in queries is lower with
fragmentation. See Fig. 2. The strength of the reduction depends on how fast
$p_{\Psi}(\beta_{l})$ (and so $n_{l}$) decreases with $l$; but, in any case,
it gets more drastic as $\beta$ increases. That is, the largest reductions are
expected at the hardest regime of $p_{\Psi}(\beta)\ll 1$. To maximize the
effect, one wishes $q_{1}(\Delta\beta_{l},\varepsilon^{\prime}_{l})$ to
decrease with $l$ as fast as possible. Note that Eq. (6) implies
$\varepsilon^{\prime}_{l}<\varepsilon^{\prime}_{l+1}$, which plays against the
latter wish. However, fortunately,
$q_{1}(\Delta\beta_{l},\varepsilon^{\prime}_{l})$ grows approximately linearly
in $\Delta\beta_{l}$ but sub-logarithmically in $1/\varepsilon^{\prime}_{l}$.
Hence, for sufficiently high $\beta$, one can make
$q_{1}(\Delta\beta_{l},\varepsilon^{\prime}_{l})$ arbitrarily smaller than
$q_{1}(\Delta\beta_{l+1},\varepsilon^{\prime}_{l+1})$ by choosing
$\Delta\beta_{l}$ sufficiently smaller than $\Delta\beta_{l+1}$.

Based on these heuristics, we next prove for Primitive 1 that Alg. 1 can not
only outperform the probabilistic algorithm but also – for sufficiently high
$\beta$ – even the coherent one, surprisingly. The proof is constructive: we
devise suitable schedules that give the desired advantage for fragmentation.
Remarkably, it is enough to consider only $r=2$ fragments. The result is valid
for any $|\Psi\rangle$ and $H$, under only mild assumptions on the success
probability $p_{\Psi}$ as a function of $\beta$. We denote the inverse
function of $p_{\Psi}$ by $p^{-1}_{\Psi}$. For simplicity, we state the
theorem explicitly for the restricted case of $H$ non-degenerate, with a
unique ground state $|\lambda_{\rm min}\rangle$ of overlap
${o}^{2}\coloneqq|\langle\lambda_{\rm min}|\Psi\rangle|^{2}$ with
$|\Psi\rangle$. However, it can be straightforwardly generalized to the
degenerate case by redefining ${o}^{2}$ as the overlap with the lowest-energy
subspace.

######  Theorem 5 (Fragmented QITE outperforms coherent QITE).

Let $|\lambda_{\rm min}\rangle\in\mathbb{H}_{\mathcal{S}}$ be the unique
ground state of $H$ and $|\Psi\rangle\in\mathbb{H}_{\mathcal{S}}$ such that
$0<{o}\leq 1/2.2$. Define the critical inverse temperature
$\beta_{c}=\frac{2}{{o}}\left[\frac{2}{e}\ln\left(\frac{8}{{o}\,\varepsilon}\right)+p^{-1}_{\Psi}(\frac{{o}}{2.2})\right]$.
Then, if ($H$ and $|\Psi\rangle$ are such that) $p_{\Psi}(\beta_{c}){\leq
1/4}$, there exists a two-fragment schedule $S_{2}$ for which, for $P_{1}$, it
holds that $Q_{S_{2}}(\beta,\varepsilon)<Q_{\text{coh}}(\beta,\varepsilon)$
for all $\beta\geq\beta_{c}$ and $0<\varepsilon<1$. In particular,
$S_{2}=\big{\\{}\Delta\beta_{1}={p^{-1}_{\Psi}\big{(}\frac{{o}}{2}\frac{1}{{\ln[e+2\ln(2/{o}\,\varepsilon)/e\beta]}}\big{)}},\Delta\beta_{2}=\beta-\Delta\beta_{1}\big{\\}}$
is a valid choice of such schedules.

The proof is given in App. G. The schedules constructed there have the sole
purpose of proving the existence of $\beta_{c}$ in general and are therefore
not necessarily optimal for each specific $H$ and $|\Psi\rangle$. For
instance, in App. H we study Gibbs-state sampling (i.e. for the maximally-
mixed state as input, with ${o}=2^{-N/2}$) for $H$ describing non-interacting
particles, where a closed-form expression for $p_{\Psi}(\beta)$ can be
obtained. For this simple case, the theorem yields
$\beta_{c}=\mathcal{O}\big{(}2^{N/2}\,N\big{)}$. However, in Sec. III.4 we
numerically optimize the schedules and obtain
$\beta_{c}=\mathcal{O}\big{(}2^{N/2}\big{)}$ for hard-to-simulate, interacting
systems. The proof in App. G exploits the additive dependence of $q_{{1}}$ on
$\beta$ and the logarithmic term in Eq. (2). Its extension to the
multiplicative case of $q_{2}$ is left for future work. Nevertheless, here, we
do consistently observe an advantage of fragmented QITE over coherent one for
$P_{2}$. More precisely, in Sec. III.4, we numerically find that also for
$P_{2}$ does fragmentation outperform coherent-QITE at Gibbs-state sampling,
with $\beta_{c}$ scaling with $N$ as in $P_{1}$ but with a somewhat larger
pre-factor (which is expectable, as $\alpha_{l}<1$ gives an exponential
dependance of $n_{l}$ on $r$ that worsens the performance). Either way, that
fragmentation can outperform quantum amplitude amplification at all is
remarkable, since the latter requires circuits
$\mathcal{O}(\sqrt{1/p_{\Psi}(\beta)})$ times deeper and $\mathcal{O}(N)$ more
ancillas than the former.

Importantly, our findings would have little practical relevance if $\beta_{\rm
c}$ was unphysically high. Fortunately, $\beta_{\rm
c}=\mathcal{O}\big{(}2^{N/2}\big{)}$ is in an intermediate regime useful for
important applications: E.g., Ground-state cooling (or, more generally, Gibbs-
state sampling at low temperatures) requires $\beta$ scaling inversely
proportionally to the spectral gap, which can be exponentially small in $N$
even for relatively simple Hamiltonians such as transverse-field Ising models
[45, 46]. Finally, as already mentioned, $P_{1}$ is particularly well-suited
for fragmentation. On the one hand, it displays $\alpha_{l}=1$ for all
$l\in[r]$. On the other hand, and most importantly, $q_{1}$ becomes optimal as
$\beta$ decreases relative to $\ln(1/\varepsilon^{\prime})$. This is
convenient to minimize Eq. (7), because the first fragments (specially the
first one) operate precisely at low $\Delta\beta_{l}$ and
$\varepsilon^{\prime}_{l}$, close to that optimality regime. The latter is
verified both analytically for the non-interacting case of App. H and
numerically for the examples of Sec. III.4 in App. I, where we consistently
observe that $\beta_{1}$ is typically only a tinny fraction of
$\ln(1/\varepsilon_{1}^{\prime})$. Colloquially speaking, the widths of the
first blue-shaded rectangles in Fig. 2 can be reduced more with $P_{1}$ than
with other primitives.

Figure 3: Runtimes and circuit depths of quantum Gibbs-state samplers running
on Primitive 1 versus inverse temperature. Red corresponds to the
probabilistic master QITE algorithm, green to the coherent one, blue to the
fragmented one with uniform schedule $S_{r}$ for the best $r$, and orange to
the fragmented one with a schedule $S_{r,a}$ as in Eq. (8) for the best $r$
and $a$ (see also Fig. 5). The three classes of Hamiltonians (expressions in
upper panels and lattice geometries in lower ones) correspond to the weighted
MaxCut problem, the quantum RBM, and a quantum version of the Sherrington-
Kirkpatrick spin glass model. Solid curves represent the means over 1000
random instances from each class, whereas (the thicknesses of) shaded curves
are the corresponding standard deviations. Also, note the logarithmic vertical
scale. The examples shown correspond to $N=12$ qubits and a tolerated error of
$\varepsilon=10^{-3}$, but qualitatively identical behaviors are observed for
all $N$ between 2 and 15 as well as for $\varepsilon=10^{-2}$ and
$\varepsilon=10^{-1}$. Besides, qualitatively similar behaviors are observed
for Primitive 2 (see App. J). Upper panels: average overall query complexity.
Both fragmented algorithms comfortably outperform the probabilistic one
already at small $\beta$. In addition, fragmentation with non-uniform schedule
outperforms even coherent QITE at a critical inverse temperature $\beta_{c}$.
Moreover, fragmentation with uniform schedule is in the same order of
magnitude as coherent QITE at $\beta=\beta_{c}$ and eventually also
outperforms it (not shown) at larger $\beta$ (between 5500 and 10000 depending
on the Hamiltonian class). The dependance of $\beta_{c}$ with $N$ is shown in
Fig. 4. Lower panels: average query depth. Defined as the maximum number of
queries per circuit run (i.e., not taking into account independent trials),
the query depth quantifies the circuit depth (relative to the depth per query)
required by the last run (the successful one). As expected, coherent QITE lies
orders of magnitude above probabilistic QITE; but, surprisingly, fragmented
QITE with non-uniform and uniform schedules are respectively almost identical
to and of the same order of the latter.

###  III.4 Fragmented quantum Gibbs-state samplers

We benchmark the performance of Alg. 1 at quantum Gibbs-state sampling by
comparing Eqs. (7) and (1) for four classes of spin-$1/2$ systems: Ising
models associated to the $i$) MaxCut and $ii$) weighted MaxCut problems [20,
21, 22]; $iii$) transverse-field Ising interactions on the restricted-
Boltzmann-machine (RBM) geometry [30, 31]; and $iv$) Heisenberg all-to-all
interactions, corresponding to a quantum generalization of the Sherrington-
Kirkpatrick model [43, 44] for spin glasses. All four classes feature long-
range frustation; and classically simulating their Gibbs states (for random
instances) is a computationally-hard task [55, 56, 57, 58, 59].

The Gibbs state $\varrho_{\beta}\coloneqq\frac{e^{-\beta(H-\lambda_{\rm
min})}}{Z_{\beta}}$ of $H$ at $\beta$, with
$Z_{\beta}\coloneqq\textrm{Tr}\left[e^{-\beta(H-\lambda_{\rm min})}\right]$
its partition function, can be prepared by QITE at $\beta/2$ on the maximally-
mixed state $\varrho_{0}=\frac{\mathds{1}}{Z_{0}}$, where $Z_{0}=2^{N}$.
Hence, the post-selection probability is
$p_{\Psi}(\beta/2,\alpha)=\alpha^{2}\,\frac{Z_{\beta}}{Z_{0}}$, where
$\alpha=1$ for $P_{1}$ and $e^{-\beta(1+\lambda_{\rm min})-\gamma}$ for
$P_{2}$. This, together with Eqs. (2) and (3), determine the overall query
complexities, with respect to $P_{1}$ and $P_{2}$, respectively, for the three
master algorithms: probabilistic [Eq. (1) for $\kappa=$ prob], coherent [Eq.
(1) for $\kappa=$ coh], and fragmented [Eq. (7)]. More technically, rather
than Eqs. (2) or (3) we use their ceiling functions, to guarantee that each
fragment’s query complexity is integer.

For $N$ up to 15 qubits, we draw 1000 random $H$’s within each class. For fair
comparison, we re-scale all $H$’s so that $\lambda_{\rm min}=-1$ and
$\lambda_{\rm max}=1$. For each of them, we calculate the complexities for
$\beta$ between 0 and 10000 and $\varepsilon=0.1$, $0.01$, or $0.001$.
Partition functions are evaluated by exact diagonalization of $H$. Evaluating
Eq. (7) requires in addition a choice of schedule. We propose

$S_{r,a}\coloneqq\left\\{\left[\left(\frac{l}{r}\right)^{a}-\left(\frac{l-1}{r}\right)^{a}\right]{\beta/2}\right\\}_{l\in[r]},$
(8)

for $a>1$, so that $\beta_{l}=\big{(}\frac{l}{r}\big{)}^{a}{\beta/2}$ for all
$l\in[r]$. This guarantees that
$\Delta\beta_{1}<\Delta\beta_{2}\ldots<\Delta\beta_{r}$ and allows us to
control the strength of the inequalities by varying $a$. For each problem
instance ($N$, $H$, and $\beta$), we sweep $r$ and, for each value of $r$, we
find the optimal $a$ through the Broyden–Fletcher–Goldfarb–Shanno (BFGS)
algorithm until minimizing $Q_{S_{r,a}}({\beta/2},\varepsilon)$ [60].

Figure 4: Critical inverse temperatures for $P_{1}$ versus number of qubits.
The error and Hamiltonian classes are the same as in Fig. 3, except for
MaxCut, defined as weighted MaxCut but with random $J_{i,j}\in\\{0,1\\}$ for
all $(i,j)$. Blue dots represent the means over 1000 instances from each
class, whereas dashed orange curves their fits over the Ansatz $\beta_{\rm
c}(N)=A\,2^{\eta\,N}+B$, with $A,B,\eta\in\mathbb{R}$. The fit results,
together with their root-mean-square deviations (RMSDs), are shown in the
insets. Similar scalings with $N$ are observed for $P_{2}$ (see App. J). In
all cases, $\beta_{\rm c}=\mathcal{O}\big{(}2^{N/2}\big{)}$ is satisfied.
Figure 5: Optimal fragmentation schedules for Primitive 1 versus inverse
temperature. System sizes are $N=5$ (blue), $N=10$ (orange), and $N=15$
(green). Solid curves represent the means over 1000 random weighted-MaxCut
Hamiltonians, whereas (the thicknesses of) shaded curves are the standard
deviations (see App. J). The tolerated error is $\varepsilon=10^{-3}$.
Qualitatively identical behaviors are observed for all $N$ between 2 and 15 as
well as for $\varepsilon=10^{-2}$ and $\varepsilon=10^{-1}$; and the same
holds for the other Hamiltonian classes in Fig. 4. The upper panel shows the
optimal number of fragments $r$ for uniform schedules $S_{r,1}$. The central
and lower panels respectively show the optimal $r$ and $a$ for non-uniform
schedules $S_{r,a}$. The dashed and dotted curves in the upper and lower
panels respectively represent fits over the Ansätze $r(\beta)=A\,\beta^{\eta}$
and $a(\beta)=A\,\beta^{\eta}$, with $A$ and $\eta\in\mathbb{R}$. The fit
results are shown in the insets. Remarkably, for non-uniform schedules, the
observed scaling for $r$ is constant not only with $\beta$ but also with $N$.

The overall complexities and circuit depths obtained (together with those for
uniform schedules, i.e. with fixed $a=1$) are shown in Fig. 3 for $P_{1}$; and
the scalings with $N$ of $\beta_{\rm c}$ in Fig. 4. Similar scalings for the
critical inverse temperature are obtained for $P_{2}$ but with somewhat higher
constant pre-factors (see App. J), which is expectable due to the non-unit
sub-normalization factors $\alpha_{k}$ in $n_{l}$. Summarizing, our numerical
experiments support the following observation.

######  Observation 6.

(Gibbs-state sampling with fragmented QITE). Let the primitives be of fixed
type, either $P_{1}$ or $P_{2}$. Then, for every $H$ and $\varepsilon>0$
studied, there exists $\beta_{\rm c}=\mathcal{O}\big{(}2^{N/2}\big{)}$ such
that, for all $\beta\geq\beta_{\rm c}$, there is a schedule
$S_{r}=S_{r}(\beta)$ that makes $Q_{S_{r}}({\beta/2},\varepsilon)\leq Q_{\rm
coh}({\beta/2},\varepsilon,\alpha)$. Moreover, the maximal circuit depth
required by fragmentation is asymptotically the same as that of probabilistic
QITE.

Apart from the notable fact that fragmentation outperforms coherent QITE for
both primitives, it is also remarkable that, long before $Q_{S_{r}}$ reaches
$Q_{\rm coh}$, at $\beta_{\rm c}$, $Q_{S_{r}}$ is already much smaller than
$Q_{\rm prob}$. Crucially, these advantages of fragmented QITE come at no cost
in circuit depth, since the query depth of fragmentation,
$\sum_{l=1}^{r}\,q(\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l})$, is
observed to almost coincide with that of repeat until success,
$q({\beta/2},\varepsilon^{\prime},\alpha)$, specially for high $\beta$. Note
that the latter needs not be the case: strictly speaking, neither $q_{1}$ nor
$q_{2}$ are additive in $\Delta\beta_{l}$ due to the non-linear dependance of
$\varepsilon^{\prime}_{l}$ on $\Delta\beta_{l}$.

Of course, the optimal schedules as functions of $\beta$ are a priori unknown.
Nevertheless, the trends we observe for the schedule proposals in Eq. (8) are
so compelling that they provide a sound basis for educated guesses in general:

######  Observation 7.

(Optimal schedules). For uniform and non-uniform schedules $S_{r,1}$ and
$S_{r,a}$, given by Eq. (8), the overall complexity for $P_{1}$ is
respectively minimised by $r=\mathcal{O}(\beta^{1/2})$ and
$(r,a)=\big{(}\mathcal{O}(1),\mathcal{O}(\beta^{1/3})\big{)}$ (see Fig. 5);
whereas for $P_{2}$ by $r=2$ and
$(r,a)=\big{(}\mathcal{O}(1),\mathcal{O}(\beta^{1/6})\big{)}$ (see App. J).

As expected from the exponential dependence on $r$ in Eq. (6), a slow growth
of $r$ with $\beta$ is observed for each $N$ to minimize
$Q_{S_{r}}({\beta/2},\varepsilon)$. This is indeed seen for $P_{1}$ with
uniform schedules (Fig. 5, upper panel). On the other hand, for $P_{2}$ with
uniform schedules, $r=2$ is observed (App. J) to minimize $Q_{S_{r}}$ but the
resulting complexity does not reach $Q_{\rm coh}$ over the scanned domain
($0\leq\beta\leq 10000$). However, for both $P_{1}$ (Fig. 5, central panel)
and $P_{2}$ (App. J) with non-uniform schedules (where fragmentation does
outperform amplitude amplification), the observed scaling of $r$ is constant
with both $\beta$ and $N$, remarkably. In turn, that $a$ grows with $\beta$
implies that each $\Delta\beta_{l}$ decreases relative to $\Delta\beta_{l+1}$
as $\beta$ grows. This is consistent with the intuition from Sec. III.3 that
each $\Delta\beta_{l}$ should be smaller than $\Delta\beta_{l+1}$. In
addition, see also App. I, where we consistently observe that, for the
obtained optimal schedules, $\Delta\beta_{1}$ is only a tinny fraction (around
0.1% to 2%) of $8\ln(4/\varepsilon_{1}^{\prime})$. In fact, for both
primitives, inserting the obtained $a(\beta)$ into Eq. (8), one sees that all
$\Delta\beta_{l}$’s (except the last one, $\Delta\beta_{r}$) also decrease in
absolute terms as $\beta$ grows. Yet, that $a$ grows slowly with $\beta$
guarantees that the $\Delta\beta_{l}$’s do not decrease too much. More
precisely, comparing with Eq. (6), we see that
$\Delta\beta_{l}>\varepsilon^{\prime}_{l}$ for all $l\in[r]$. This is an
important sanity check, because if $\Delta\beta_{l}<\varepsilon^{\prime}_{l}$,
the identity operator would readily provide an
$(\varepsilon^{\prime}_{l},1)$-block-encoding of $F_{\Delta\beta_{l}}(H)$,
hence rendering the obtained scaling for $a(\beta)$ meaningless.

##  IV Concluding discussion

We have presented two QITE primitives and a master QITE algorithm. The first
primitive is designed for block-encoding Hamiltonian oracles and has query
complexity (number of oracle calls) sub-additive in the inverse-temperature
$\beta$ and $\ln(\varepsilon^{-1})$, with $\varepsilon$ the error. This
scaling is better than all previously-known bounds [11, 15] for $\beta\lesssim
8\ln(4\,\varepsilon^{-1})$ and becomes provably optimal for
$\beta\ll\ln(\varepsilon^{-1})$. Optimality is proven by showing saturation of
a universal cooling-speed limit that is an imaginary-time counterpart of the
celebrated no fast-forwarding theorem for real-time simulations [40, 41, 42].
It is an open question what the optimal scaling is away from the saturation
regime. Coincidentally, the first steps of our master algorithm operate
precisely in that regime. On the other hand, the second primitive is designed
for a simplified model of real-time evolution oracles involving a single time.
Its complexity is multiplicative in $\beta$ and $\ln(\varepsilon^{-1})$, but
it requires a single ancillary qubit throughout and its oracle is
experimentally-friendlier than in previous QITE primitives. Interestingly,
preliminary numerical analysis [51] suggests that the asymptotic scaling with
$\beta$ of both primitives’ complexities could actually be significantly
better than in the analytical bounds above, for $P_{1}$ even reaching levels
as good as
$q_{1}(\beta,\varepsilon^{\prime})=\mathcal{O}\big{(}\sqrt{\beta\,\ln(1/\varepsilon^{\prime})}\big{)}$.

Our primitives are based on two technical contributions to quantum signal
processing (QSP) [36, 37, 15, 14] relevant on their own. The first one is a
bound on the approximation error of Hermitian-operator functions by their
truncated Chebyshev series, for any analytical real function. The second one
is a novel, Fourier-based QSP variant for real-time evolution oracles superior
to previous ones [14] in that it requires a single real time (and therefore a
single ancilla), instead of multiple ones. Moreover, it is also experimentally
friendly in that it requires no qubitization [37].

Primitive technicalities aside, the main conceptual contribution of this work
is the master QITE algorithm, which is conceptually simple, yet surprisingly
powerful. It is based on breaking the evolution into small-$\beta$ fragments.
This gives a large reduction in wasted queries and circuit depth, yielding an
overall runtime competitive with (and for high $\beta$ even better than) that
of coherent approaches based on quantum amplitude amplification (QAA). This is
remarkable since the latter requires in general $N$ extra ancillary qubits and
circuits $\mathcal{O}\big{(}1/\sqrt{p_{\Psi}(\beta,\alpha)}\big{)}$ times
deeper than the former. To put this in perspective, it is illustrative to
compare with quantum amplitude estimation (QAE). In its standard form, QAE has
similar hardware requirements as QAA [35]. However, recently, interesting
algorithms have appeared [61, 62] that perform partial QAE with circuit depths
that can interpolate between the probabilistic and coherent cases. In
contrast, here, we beat full QAA using circuit depths for most runs much lower
than in the bare probabilistic approach.

That fragmented QITE outperforms coherent QITE is proven rigorously for
Primitive 1 and also supported by exhaustive numerical evidence for both
primitives. Namely, our numerical experiments address random instances of
Ising, transverse-field Ising, and Heisenberg-like Hamiltonians encoding
computationally hard problems relevant for combinatorial optimisations,
generative machine learning, and statistical physics, e.g. We emphasize that
our analysis of is based on the analytical upper bounds on the query
complexity we obtained, instead of the complexities themselves. The
corresponding analysis for the actual (numerically obtained) query
complexities requires re-optimizing the fragmentation schedules. Preliminary
observations [51] in that direction are again promising, indicating that the
actual overall complexities may be orders of magnitude lower than in Fig. 3,
e.g. In any case, qualitatively similar interplays between fragmentation and
QAA are expected even for other types of primitives (beyond QITE) whose
complexity and post-selection probability have similar scalings. All these
exciting prospects are being explored for future work.

Our findings open a new research direction towards mid-term high-precision
quantum algorithms. In particular, the presented primitives, cooling-speed
limit, QSP methods, and master algorithm constitute a powerful toolbox for
quantum signal processors specially relevant for the transition from NISQ to
early prototypes of fault-tolerant hardware.

##  V Methods: quantum signal processing

Quantum signal processing (QSP) is a powerful method to obtain an
$\varepsilon^{\prime}$-approximate block encoding of an operator function
$f(H)\coloneqq\sum_{\lambda}f(\lambda)|\lambda\rangle\langle\lambda|$, where
$\\{|\lambda\rangle\in\mathbb{H}_{\mathcal{S}}\\}$ are the eigenvectors and
$\\{\lambda\\}$ the eigenvalues of a Hamiltonian $H$, from queries to an
oracle for $H$ [36]. We note that QSP can also be extended to non-Hermitian
operators [15], but here we restrict to the Hermitian case for simplicity. We
present two QSP methods for general functions one for each oracle model in
Defs. 2 and 3. Our QITE primitives are then obtained by particularizing these
methods to the case $f(H)=F_{\beta}(H)$, with
$F_{\beta}(H)=e^{-\beta(H-\lambda_{\rm min})}$.

###  V.1 Real-variable function design with single-qubit rotations

We start by reviewing how to approximate functions of one real variable with
single-qubit pulses.

####  V.1.1 Single-qubit QSP method 1

Consider the single qubit rotation $R_{1}(\theta,\phi)\coloneqq
e^{i\theta{X}}e^{i\phi{Z}}$, where ${X}$ and ${Z}$ are the first and third
Pauli matrices, respectively, and $\phi\in[0,2\pi]$. The angle
$\theta\in[-\pi,\pi]$ is the signal to be processed and the rotation
$e^{i\theta{X}}$ is called the iterate. One can show [63] that, given
$q\in\mathbb{N}_{\text{even}}$ and a sequence of angles
${\boldsymbol{\Phi}_{1}}=\big{(}\phi_{1},\cdots,\phi_{q+1}\big{)}\in\mathbb{R}^{q+1}$,
the sequence of rotations
$\mathcal{R}_{1}\left(\theta,{\boldsymbol{\Phi}_{1}}\right)\coloneqq
e^{i\phi_{q+1}{Z}}\prod_{k=1}^{q/2}R_{1}(-\theta,\phi_{2k})R_{1}(\theta,\phi_{2k-1})$
has matrix representation in the computational basis

$\mathcal{R}_{1}\left(\theta,{\boldsymbol{\Phi}_{1}}\right)=\left(\begin{array}[]{cc}B(\cos\theta)&i\,\sin\theta\,D(\cos\theta)\\\
i\,\sin\theta\,D^{*}(\cos\theta)&B^{*}(\cos\theta)\end{array}\right),$ (9)

where $B$ and $D$ are polynomials in $\cos\theta$ with complex coefficients
determined by $\boldsymbol{\Phi}_{1}$.

For target real polynomials $\mathscr{B}(\cos\theta)$ and
$\mathscr{D}(\cos\theta)$, we wish to find ${\boldsymbol{\Phi}_{1}}$ that
generates $B(\cos\theta)$ and $D(\cos\theta)$ with $\mathscr{B}(\cos\theta)$
and $\mathscr{D}(\cos\theta)$ as either their real or imaginary parts,
respectively. This can be done iff they satisfy (App. K)

$\mathscr{B}^{2}(\cos\theta)+\sin^{2}\theta\,\mathscr{D}^{2}(\cos\theta)\leq
1$ (10)

for all $\theta$, and have the form

$\begin{split}\mathscr{B}(\cos\theta)&=\sum_{k=0}^{q/2}b_{k}\cos(2k\theta)\\\
\sin\theta\,\mathscr{D}(\cos\theta)&=\sum_{k=1}^{q/2}d_{k}\sin\left(2k\theta\right),\end{split}$
(11)

with $b_{k}\in\mathbb{R}$ and $d_{k}\in\mathbb{R}$. Alternatively, Eq. (11)
can also be expressed in terms of Chebyshev polynomials of first
$T_{k}(\cos\theta)\coloneqq\cos(k\theta)$ and second
$U_{k}(\cos\theta)\coloneqq\sin\left((k+1)\theta\right)/\sin\theta$ kinds.
This can be used to obtain either Chebyshev or Fourier series of target
operator functions. If the target expansion satisfies Eqs. (10) and (11), the
angles ${\boldsymbol{\Phi}_{1}}$ can be computed classically in time
$\mathcal{O}\left(\text{poly}(q)\right)$ [63, 64, 65, 66].

####  V.1.2 Single-qubit QSP method 2

This method is inspired by a construction in Ref. [67] and shown in detail in
a companying paper [39]. The fundamental gate is
$R_{2}(x,\omega,\zeta,\eta,\varphi,\kappa)=e^{i\frac{\zeta+\eta}{2}Z}e^{-i\varphi
Y}e^{i\frac{\zeta-\eta}{2}Z}e^{i\omega xZ}e^{-i\kappa Y}$, which has five
adjustable parameters $\\{\omega,\boldsymbol{\xi}\\}\in\mathbb{R}^{5}$, where
$\boldsymbol{\xi}\coloneqq\\{\zeta,\eta,\varphi,\kappa\\}$. Here,
$x\in\mathbb{R}$ will play the role of the signal and $e^{i\omega xZ}$ that of
the iterate. In Ref. [67], it was observed that the gate sequence
$\mathcal{R}_{2}(x,\boldsymbol{\omega},\boldsymbol{\Phi}_{2})\coloneqq\prod_{k=0}^{q}R_{2}(x,\omega_{k},\boldsymbol{\xi}_{k})$,
with
$\boldsymbol{\omega}\coloneqq\\{\omega_{0},\cdots,\omega_{q}\\}\in\mathbb{R}^{q+1}$
and
$\boldsymbol{\Phi}_{2}=\\{\boldsymbol{\xi}_{0},\cdots,\boldsymbol{\xi}_{q}\\}\in\mathbb{R}^{4(q+1)}$,
can encode certain finite Fourier series into its matrix components. In [39],
not only it is formally proven that for any target series a unitary operator
can be built with it as one of its matrix elements but also we provide an
explicit, efficient recipe for finding the adequate choice of pulses
$\boldsymbol{\Phi}_{2}$. This is the content of the following lemma.

######  Lemma 8.

[Single-qubit Fourier series synthesis] Given
$\tilde{g}_{q}(x)=\sum_{m=-q/2}^{q/2}c_{m}\,e^{imx}$, with $q\in\mathbb{N}$
even, there exist $\boldsymbol{\omega}$ and $\boldsymbol{\Phi}_{2}$ such that
$\langle
0|\,\mathcal{R}_{2}(x,\boldsymbol{\omega},\boldsymbol{\Phi}_{2})\,|0\rangle=\tilde{g}_{q}(x)$
for all $|x|\leq\pi$ iff $|\tilde{g}_{q}(x)|\leq 1$ for all $|x|\leq\pi$.
Moreover, $\omega_{0}=0$ and $\omega_{k}=(-1)^{k}/2$, for all $1\leq k\leq q$,
and $\boldsymbol{\Phi}_{2}$ can be calculated classically from
$\\{c_{m}\\}_{m}$ in time $\mathcal{O}\left(\text{poly}(q)\right)$.

###  V.2 Operator-function design from block-encoded oracles

Here, we synthesize an ($\varepsilon^{\prime}$,1)-block-encoding of $f(H)$
from queries to an oracle for $H$ as in Def. 2. The algorithm can be seen as a
variant of the single-ancilla method from Ref. [37] with slightly different
pulses. The basic idea is to design a circuit, $P_{1}$, that generates a
perfect block-encoding $V_{\boldsymbol{\Phi}_{1}}$ of a target Chebyshev
expansion $\tilde{f}_{q}(H)\coloneqq\sum^{q/2}_{k=0}b_{k}\,T_{k}(H)$ that
$\varepsilon_{\text{tr}}$-approximates $f(H)$, for some
$0\leq\varepsilon_{\text{tr}}\leq\varepsilon^{\prime}$. This can be done by
adjusting $\boldsymbol{\Phi}_{1}$ as in Sec. V.1.1. Note that the
achievability condition (10) requires that $\|\tilde{f}_{q}(H)\|\leq 1$, but
we only guarantee $\|\tilde{f}_{q}(H)\|\leq 1+\varepsilon_{\text{tr}}$.
However, this can be easily accounted for introducing an inoffensive sub-
normalization $\alpha=(1+\varepsilon_{\text{tr}})^{-1}$ (see, e.g., Lemma 14
in Ref. [37]), which we neglect here throughout. Choosing $\tilde{f}_{q}$ as
the truncated Chebyshev series of $f$ with truncation error
$\varepsilon_{\text{tr}}\leq\varepsilon^{\prime}$, we obtain the desired
block-encoding of $f(H)$. For $f$ analytical, the error fulfills [53]

$\varepsilon_{\text{tr}}\leq\frac{\underset{\lambda\in[\lambda_{\text{min}},\lambda_{\text{max}}]}{\max}\left|f^{(q/2+1)}(\lambda)\right|}{2^{\frac{q}{2}}(q/2+1)!},$
(12)

with $f^{(q/2+1)}$ the $(q/2+1)$-th derivative of $f$. This allows one to
obtain the truncation order $q/2$ and Chebyshev coefficients
$\boldsymbol{b}\coloneqq\\{b_{k}\\}_{0\leq k\leq q/2}$ (see App. L). Then,
from $\boldsymbol{b}$, one can calculate the required $\boldsymbol{\Phi}_{1}$
(see App. M).

Next, we explicitly show how to generate $V_{\boldsymbol{\Phi}_{1}}$. Using
the short-hand notation
$|0_{\lambda}\rangle\coloneqq|\lambda\rangle|0\rangle\in\mathbb{H}_{\mathcal{SA}_{O_{1}}}$
and Def. 2, one writes
$O_{1}|0_{\lambda}\rangle=\lambda|0_{\lambda}\rangle+\sqrt{1-\lambda^{2}}|0_{\lambda}^{\perp}\rangle$
with $\langle 0_{\lambda}|0_{\lambda}^{\perp}\rangle=0$. This defines the
2-dimensional subspace
$\mathbb{H}_{\lambda}\coloneqq\text{span}\\{|0_{\lambda}\rangle,|0_{\lambda}^{\perp}\rangle\\}$.
To exploit the single-qubit formalism from Sec. V.1.1, one needs an iterate
that acts as an $SU(2)$ rotation within each $\mathbb{H}_{\lambda}$. In
general, $O_{1}$ itself is not appropriate for this due to leakage out of
$\mathbb{H}_{\lambda}$ by repeated applications of $O_{1}$. However, there is
a simple oracle transformation – qubitization – that maps $O_{1}$ into another
block-encoding $O^{\prime}_{1}$ of the same $H$ but with the desired property
[37]. The transformed oracle reads (see App. N)

$O^{\prime}_{1}=\bigoplus_{\lambda}e^{-i\theta_{\lambda}{Y}_{\lambda}},$ (13)

with $\theta_{\lambda}{:}=\cos^{-1}(\lambda)$ and $Y_{\lambda}\coloneqq
i(|0_{\lambda}^{\perp}\rangle\langle 0_{\lambda}|-|0_{\lambda}\rangle\langle
0_{\lambda}^{\perp}|)$.

Although the qubit resemblance could be considered in a direct analogy to QSP
for a single qubit, it leads to a more strict class of achievable functions
than if we resort to one additional qubit (single-ancilla QSP) [37]. This
extra ancilla controls the action of the oracle $O^{\prime}_{1}$ through the
iterate

$V_{0}=\mathds{1}\otimes|+\rangle\langle+|+O^{\prime}_{1}\otimes|-\rangle\langle-|$
(14)

on $\mathbb{H}_{\mathcal{SA}}$, where $|\pm\rangle$ are the eigenstates of the
Pauli operator ${X}$ for the QSP qubit ancilla. Throughout this section,
$\mathds{1}$ is the identity operator on $\mathbb{H}_{\mathcal{SA}_{O_{1}}}$
and $M$ denotes the single-qubit Hadamard gate. Let us define the operators
$V_{\phi}=V_{0}\left(\mathds{1}\otimes e^{i\phi{Z}}\right)$ and
$\bar{V}_{\phi}=V_{0}^{\dagger}\left(\mathds{1}\otimes e^{i\phi{Z}}\right)$
for a given phase $\phi\in[0,2\pi]$, which play the role of $R(\theta,\phi)$
of the previous sub-section with $\theta_{\lambda}$ playing the role of
$\theta$ for each $\lambda$. These operators can be phase iterated to generate

$V_{{\boldsymbol{\Phi}_{1}}}\coloneqq
W_{\text{out}}\left(\bar{V}_{\phi_{q}}{V}_{\phi_{q-1}}\cdots\bar{V}_{\phi_{2}}{V}_{\phi_{1}}\right)W_{\text{in}}$
(15)

on $\mathbb{H}_{\mathcal{SA}}$, with ancilla pre- and post-processing
unitaries $W_{\text{in}}\coloneqq\mathds{1}\otimes M$ and
$W_{\text{out}}:=\mathds{1}\otimes\big{(}M\,e^{i\phi_{q+1}{Z}}\big{)}$,
respectively, with $M$ the single-qubit Hadamard matrix. The resulting
circuit, $P_{1}$, is depicted in Figs. 6.a) and 6.b).

The following pseudocode gives the entire procedure.

input : analytical function $f:\text{Dom}\to$ Img, with
$[\lambda_{\text{min}},\lambda_{\text{max}}]\subseteq\text{Dom}$ and
$\text{Img}\subseteq[-1,1]$, error $\varepsilon^{\prime}>0$, oracle $O_{1}$
for $H$, and its inverse $O^{\dagger}_{1}$.

output : a unitary quantum circuit $P_{1}$.

1 Obtain a truncation order $q/2$ s.t.
$\varepsilon_{\text{tr}}\leq\varepsilon^{\prime}$;

2 calculate the Chebyshev coefficients $\boldsymbol{b}$;

3 calculate the rotation angles $\boldsymbol{\Phi}_{1}$ (App. M);

4 begin construction of $P_{1}$:

5 apply $W_{\text{in}}=\mathds{1}\otimes M$ on $\mathcal{A}$;

6 for  _$k=1$ to $k=q$_ do

7 apply $\mathds{1}\otimes e^{i\phi_{k}Z}$;

8 if  _$k$ is odd_ then apply $V_{0}$ from Eq. (14);

9 else apply $V_{0}^{\dagger}$;

10

11 end for

12 apply $W_{\text{out}}=\mathds{1}\otimes\big{(}M\,e^{i\phi_{q+1}{Z}}\big{)}$
on $\mathcal{A}$;

13 end

14

Algorithm 2 Operator-function design from block-encoded Hamiltonian oracles

Figure 6:  QSP primitives for generic operator function design. a) Both
circuits $P_{1}$ from Alg. 2 and $P_{2}$ from Alg. 3 have the same structure.
If the ancillas are initialised and post-selected in
$|0\rangle_{\mathcal{A}}$, the circuit prepares the system state
$\frac{\tilde{f}_{q}(H)|\Psi\rangle}{||\tilde{f}_{q}(H)|\Psi\rangle||}$, which
$\varepsilon$-approximates the target output
$\frac{f(H)|\Psi\rangle}{\|f(H)|\Psi\rangle\|}$. The details specific to
$P_{1}$ and $P_{2}$ are respectively shown in panels b) and c).
$W_{\text{in}}$ and $W_{\text{out}}$ are fixed ancillary unitaries, and $M$ is
a single-qubit Hadamard gate. The basic blocks $V_{k}$ in panel a) represent
the gates $V_{\phi_{k}}$ in b) and $V_{\boldsymbol{\xi}_{k}}$ in c). Each
$V_{\phi_{k}}$ involves one query to the qubitized oracle $O^{\prime}_{1}$,
which in turn requires one query to $O_{1}$ and one to its inverse
$O^{\dagger}_{1}$ (see Fig. 13). Whereas each $V_{\boldsymbol{\xi}_{k}}$
involves one query to the oracle $O_{2}$. $\bar{V}_{k}$ is defined as $V_{k}$
but with $O_{1}^{\prime\dagger}$ substituting $O^{\prime}_{1}$ or
$O_{2}^{\dagger}$ substituting $O_{2}$. Hence, the query complexities of
$P_{1}$ and $P_{2}$ are respectively $2q$ and $q$. The approximating function
$\tilde{f}_{q}$ is determined by the angles
${\boldsymbol{\Phi}_{1}}=\big{(}\phi_{1},\cdots,\phi_{q+1}\big{)}$ or
$\boldsymbol{\Phi}_{2}=\\{\boldsymbol{\xi}_{0},\cdots,\boldsymbol{\xi}_{q}\\}$
in the rotations $R_{z}(2\phi_{k})=e^{i\phi_{k}{Z}}$ or
$R_{y}(2\kappa_{k})=e^{i\kappa_{k}Y}$ and
$R_{zyz}(\boldsymbol{\xi}_{k})=R_{z}(\zeta_{k}+\eta_{k})\,R_{y}(2\varphi_{k})\,R_{z}(\zeta_{k}-\eta_{k})$,
with $\boldsymbol{\xi}_{k}=\\{\zeta_{k},\eta_{k},\varphi_{k},\kappa_{k}\\}$.
For $P_{1}$ and $P_{2}$, these angles are chosen such that $\tilde{f}_{q}$ is
a high-precision Chebyshev and Fourier approximation of $f$, respectively.

The correctness and complexity of Algorithm 2 are addressed by the following
lemma, proven in App. O .

######  Lemma 9.

Let $f$, $\varepsilon^{\prime}$, $O_{1}$, and $O^{\dagger}_{1}$ be as in the
input of Alg. 2. Then, for any $q$ s.t. $\varepsilon_{\text{tr}}$ in Eq. (12)
is no greater than $\varepsilon^{\prime}$, there exists
$\boldsymbol{\Phi}_{1}\in\mathbb{R}^{q+1}$ such that
$V_{\boldsymbol{\Phi}_{1}}$ in Eq. (15) is a ($\varepsilon^{\prime}$,1)-block-
encoding of $f(H)$. The circuit $P_{1}$ generating $V_{\boldsymbol{\Phi}_{1}}$
requires a single-qubit ancilla, $q$ queries to $O_{1}$ and $O^{\dagger}_{1}$
each, and $g_{P_{1}}=\mathcal{O}\left(g_{O_{1}}+|\mathcal{A}_{O_{1}}|\right)$
gates per query, with $g_{O_{1}}$ the gate complexity of $O_{1}$. Furthermore,
the classical runtime (calculations of $\boldsymbol{b}$ and
$\boldsymbol{\Phi}_{1}$) is within complexity $\mathcal{O}\big{(}{\rm
poly}(q/2)\big{)}$.

Some final comments about the input function are in place. The restriction of
$f$ being analytical is needed to determine the truncation order through Eq.
(12). In fact, to evaluate the RHS of the equation exactly, one needs in
general closed-form expression for $f$. However, if the required truncation
order is given in advance, the corresponding Chebyshev coefficients can be
obtained from $q/2+1$ evaluations of $f$ in specific points (the nodes of the
Chebyshev polynomials). In that case, a closed-form expression for $f$ is not
required and a classical oracle for evaluating it suffices. Moreover, it is
important to note that a satisfactory Chebyshev approximation is guaranteed to
exist for all bounded and continuous functions [52]. If the Chebyshev
expansion is given, then step 1 of Algorithm 2 can obviously be skipped and
$f$ is not required at all. We further note that Alg. 2 can also be applied
even to non-continuous functions over restricted domains without the
discontinuities. This is for instance the case of the inverse function, which
can be well-approximated over the sub-domain $[-1,-\delta]\cup[\delta,1]$ by a
pseudo-inverse polynomial of $\delta$-dependent degree [68].

####  V.2.1 QITE-primitive from a block-encoding oracle

QITE primitive 1 corresponds to the output of Alg. 2 for
$f=F_{\beta}=e^{-\beta(\lambda-\lambda_{\rm min})}$. The Chebyshev
coefficients can be readily obtained from the Jacobi-Anger expansion [38]

$e^{-\beta\lambda}=I_{0}(\beta)+2\sum_{k=1}^{\infty}(-1)^{k}I_{k}(\beta)\,T_{k}({\lambda}),$
(16)

where $I_{k}(\beta)$ is a modified Bessel function. The proof of Theorem 1
thus follows straightforwardly from Lemma 9.

###### Proof of Theorem 1..

The function $F_{\beta}:[\lambda_{\rm min},\lambda_{\rm
max}]\rightarrow(0,1]$, with
$F_{\beta}(\lambda)=e^{-\beta(\lambda-\lambda_{\rm min})}$ for all
$\lambda\in[\lambda_{\rm min},\lambda_{\rm max}]$ satisfies all the
assumptions of Lemma 9. Hence, on input $f=F_{\beta}$, Alg. 2 outputs an
$(\beta,\varepsilon^{\prime},1)$-QITE-primitive. By Eq. (12), the
corresponding truncation error is ($q^{\prime}\coloneqq q/2$)

$\displaystyle\varepsilon_{\text{tr}}\leq\frac{\beta^{q^{\prime}+1}}{2^{q^{\prime}}(q^{\prime}+1)!}$
$\displaystyle\leq\sqrt{\frac{2}{\pi(q^{\prime}+1)}}\left(\frac{e\beta}{2(q^{\prime}+1)}\right)^{q^{\prime}+1}$
$\displaystyle=\mathcal{O}\left(\left(\frac{\beta
e}{2q^{\prime}}\right)^{q^{\prime}}\right),$ (17)

where Stirling inequality has been invoked. (We note also that the first
inequality can also be obtained from explicit summation using Eq. (16) and the
properties of the Bessel functions [54].) Then, imposing
$\mathcal{O}\Big{(}\big{(}\frac{\beta
e}{2q^{\prime}}\big{)}^{q^{\prime}}\Big{)}\leq\varepsilon^{\prime}$ and
solving for $q^{\prime}$ [15] gives the query complexity of Eq. (2). ∎

###  V.3 Operator function design from real-time evolution oracles

Here, we synthesize an ($\varepsilon^{\prime},\alpha$)-block-encoding of
$f(H)$ from an oracle for $H$ as in Def. 3. We proceed as in Sec. V.2, but
with a circuit $P_{2}$ generating a perfect block-encoding
$V_{\boldsymbol{\Phi}_{2}}$ of a target Fourier expansion
$\tilde{g}_{q}(H)\coloneqq\sum_{m=-q/2}^{q/2}c_{m}e^{imHt}$ that
$\varepsilon_{\text{tr}}$-approximates $\alpha\,f(H)$, for some
$\varepsilon_{\text{tr}}\leq\varepsilon^{\prime}$, $\alpha\leq 1$, and a
suitable $t>0$. This is done by adjusting $\boldsymbol{\Phi}_{2}$ according to
Lemma 8. The function $\tilde{g}_{q}$ is a Fourier approximation of an
intermediary function $g$ such that $g(\lambda,t)\coloneqq
g(x_{\lambda})=\alpha\,f(\lambda)$, for $t$ chosen so that
$x_{\lambda}=\lambda\,t$ is in the interval of convergence of $\tilde{g}_{q}$
to $g$ for all $\lambda\in[\lambda_{\text{min}},\lambda_{\text{max}}]$. The
reason for this intermediary step here is to circumvent the well-known Gibbs
phenomenon, by virtue of which convergence of a Fourier expansion cannot in
general be guaranteed at the boundaries. In turn, the sub-normalization factor
$\alpha$ arises because our $\tilde{g}_{q}$ converges to $g$ only for
$|x_{\lambda}|<{\pi/2}$, whereas Lemma 8 requires that
$|\tilde{g}_{q}(x_{\lambda})|\leq 1$ for all $|x_{\lambda}|\leq\pi$. This
forces one to sub-normalize the expansion so as to guarantee normalization
over the entire domain. (As in Sec. V.2, the inoffensive sub-normalization
factor $(1+\varepsilon_{\text{tr}})^{-1}$ is neglected.)

More precisely, we employ (see Ref. [39]) a construction from Ref. [14] that,
given $0<\delta\leq\pi/2$ and a power series that
$\frac{\varepsilon_{\text{tr}}}{4}$-approximates $g$, gives
$\boldsymbol{c}\coloneqq\\{c_{m}\\}_{|m|\leq q/2}$ such that $\tilde{g}_{q}$
$\varepsilon_{\text{tr}}$-approximates $g$ for all
$x_{\lambda}\in[-\pi/2+\delta,\pi/2-\delta]$, if

$q\geq\left\lceil\frac{2\pi}{\delta}\ln\left(\frac{4}{\varepsilon_{\text{tr}}}\right)\right\rceil.$
(18)

For $f$ analytical, one can obtain the power series of $g$ from a truncated
Taylor series of $f$ using that $g(x_{\lambda})=\alpha f(\lambda)$. The
truncation order $L$ can be obtained from the remainder:

$\frac{\varepsilon_{\text{tr}}}{4}\leq\frac{\underset{\lambda\in[\lambda_{\text{min}},\lambda_{\text{max}}]}{\max}\left|\alpha\,f^{(L+1)}(\lambda)\right|}{(L+1)!}.$
(19)

In turn, the conditions
$[\lambda_{\text{min}},\lambda_{\text{max}}]\subseteq[-1,1]$ and
$x_{\lambda}\in[-\pi/2+\delta,\pi/2-\delta]$ lead to the natural choice
$t=\pi/2-\delta$. (This renders $\tilde{g}_{q}$ periodic in $x_{\lambda}$ with
period $2\pi$.) In addition, in Ref. [39] the sub-normalization constant
$\alpha$ is bounded in terms of the obtained $t$ and Taylor coefficients
$\boldsymbol{a}\coloneqq\\{a_{l}\\}_{0\leq l\leq L}$ of $f$. It suffices to
take $\alpha$ such that

$\sum_{l=0}^{L}\big{|}a_{l}/(1-2\,\delta/\pi)^{l}\big{|}\leq\alpha^{-1}.$ (20)

Note that $L$ and $\alpha$ are inter-dependent. One way to determine them is
to increase $L$ and iteratively adapt $\alpha$ until Eqs. (19) and (20) are
both satisfied. Alternatively, if the expansion converges sufficiently fast
(e.g., if
$\lim_{l\rightarrow\infty}|\frac{a_{l+1}}{a_{l}}|<1-\frac{2\delta}{\pi}$), one
can simply substitute $L$ in Eq. (20) by $\infty$. This is indeed the case
with QITE primitives (see Sec. V.3.1). There, the substitution introduces a
slight increase of unnecessary sub-normalization but makes the resulting
$\alpha$ independent of $L$, thus simplifying the analysis. Then, from the
obtained $\boldsymbol{c}$, one can finally calculate the required
$\boldsymbol{\Phi}_{2}$ [39].

Next, we explicitly show how to generate $V_{\boldsymbol{\Phi}_{2}}$. The
iterate is now taken simply as the oracle itself:
$O_{2}=\mathds{1}\otimes|0\rangle\langle 0|+e^{-iHt}\otimes|1\rangle\langle
1|$. Notice that, in contrast to $O_{1}$, $O_{2}$ readily acts as an $SU(2)$
rotation on each 2-dimensional subspace
$\text{span}\\{|\lambda\rangle|0\rangle,|\lambda\rangle|1\rangle\\}$. This
relaxes the need for qubitization. The basic QSP blocks for the unitary
operator $V_{{\boldsymbol{\Phi}_{2}}}\coloneqq
W_{\text{out}}\left(\bar{V}_{\boldsymbol{\xi}_{q}}{V}_{\boldsymbol{\xi}_{q-1}}\cdots\bar{V}_{\boldsymbol{\xi}_{2}}{V}_{\boldsymbol{\xi}_{1}}\right)W_{\text{in}}$
are

${V}_{\boldsymbol{\xi}_{k}}\coloneqq\left[\mathds{1}\otimes\left(e^{i\frac{\zeta_{k}+\eta_{k}}{2}Z}e^{-i\varphi_{k}Y}e^{i\frac{\zeta_{k}-\eta_{k}}{2}Z}\right)\right]O_{2}\left[\mathds{1}\otimes
e^{-i\kappa_{k}Y}\right],$ (21a) and
$\bar{V}_{\boldsymbol{\xi}_{k}}\coloneqq\left[\mathds{1}\otimes\left(e^{i\frac{\zeta_{k}+\eta_{k}}{2}Z}e^{-i\varphi_{k}Y}e^{i\frac{\zeta_{k}-\eta_{k}}{2}Z}\right)\right]O_{2}^{\dagger}\left[\mathds{1}\otimes
e^{-i\kappa_{k}Y}\right]$ (21b)

with
${\boldsymbol{\xi}_{k}}\coloneqq\\{\zeta_{k},\eta_{k},\varphi_{k},\kappa_{k}\\}$.
$V_{\boldsymbol{\xi}_{k}}$ and $\bar{V}_{\boldsymbol{\xi}_{k}}$ play a similar
role to $R_{2}(x,\omega_{k},\boldsymbol{\xi}_{k})$ in Sec. V.1.2 (with
$x_{\lambda}$ inside $O_{2}$ playing the role of $x$ there for each
$\lambda$). Here we take

$W_{\text{in}}=\mathds{1}\otimes\big{[}e^{i\frac{\zeta_{0}+\eta_{0}}{2}Z}e^{-i\varphi_{0}Y}e^{i\frac{\zeta_{0}-\eta_{0}}{2}Z}e^{-i\kappa_{0}Y}\big{]}$
(22)

and $W_{\text{out}}=\mathds{1}$. The circuit is depicted in Figs. 6.a) and
6.c).

The following pseudocode presents the entire procedure.

1

input : analytical function $f:\text{Dom}\to$ Img, with
$[\lambda_{\text{min}},\lambda_{\text{max}}]\subseteq\text{Dom}$ and
$\text{Img}\subseteq[-1,1]$, error $\varepsilon^{\prime}>0$,
$\delta\in(0,\pi/2)$, oracle $O_{2}$ for $H$ at a time $t={\pi/2-\delta}$, and
its inverse $O^{\dagger}_{2}$.

output : unitary quantum circuit $P_{2}$.

2 Obtain $\alpha$, $L$, and $\boldsymbol{a}$ s.t.
$\varepsilon_{\text{tr}}\leq\varepsilon^{\prime}$;

3 obtain $q$ s.t. Eq. (18) holds;

4 calculate the Fourier coefficients $\boldsymbol{c}$;

5 calculate the rotation angles $\boldsymbol{\Phi}_{2}$ [39];

6 begin construction of $P_{2}$:

7 apply $W_{\text{in}}$ from Eq. (22) on $\mathcal{A}$;

8 for  _$k=1$ to $k=q$_ do

9 if  _$k$ is odd_ then apply $V_{\boldsymbol{\xi}_{k}}$ from Eq. (21a);

10 else appy $\bar{V}_{\boldsymbol{\xi}_{k}}$ from Eq. (21b);

11

12 end for

13

14 end

15

Algorithm 3 Operator-function design from real-time evolution Hamiltonian
oracles

The correctness and complexity of Alg. 3 are addressed by the following lemma,
proved in Ref. [39].

######  Lemma 10.

Let $f$, $\varepsilon^{\prime}$, $\delta$, $O_{2}$, and $O_{2}^{\dagger}$ be
as in the input of Alg. 3. Then, for any $q$ satisfying Eq. (18) and $\alpha$
satisfying Eq. (20), there exists
$\boldsymbol{\Phi}_{2}\in\mathbb{R}^{4(q+1)}$ such that
$V_{\boldsymbol{\Phi}_{2}}$ is an ($\varepsilon^{\prime},\alpha$)-block-
encoding of $f(H)$. The circuit $P_{2}$ generating $V_{\boldsymbol{\Phi}_{2}}$
requires $q/2$ queries to $O_{2}$ and $O_{2}^{\dagger}$ each and $4+g_{O_{2}}$
gates per query, with $g_{O_{2}}$ the gate complexity of $O_{2}$. Moreover,
the classical runtime is within complexity
$\mathcal{O}\big{(}\text{poly}(L,q/2)\big{)}$.

Clearly, if a suitable power series for $f$ is a-priori available, analyticity
of $f$ is not required and step 1 in Alg. 3 can be skipped. Finally, we note
that it is always possible to avoid sub-normalization by introducing a
periodic extension of $f$ that is readily normalized over the entire domain of
the Fourier expansion. However, Eq. (18) is then no longer valid and one must
assess the query complexity on a case-by-case basis. This can for instance be
tackled numerically by variationally optimising the gate sequence
$\mathcal{R}_{2}(x,\boldsymbol{\omega},\boldsymbol{\Phi}_{2})$ to block-encode
the periodic extension of $f$ [67]. Either way, clearly, if a normalized
Fourier expansion is a-priori available, one can skip steps 1 to 3 in Alg. 3.

####  V.3.1 QITE-primitive from a real-time evolution oracle

QITE primitive 2 is the output of Alg. 3 for $f=F_{\beta}$. The proof of Theo.
2 thus follows straight from Lemma 10.

###### Proof of Theorem 2.

The function
$F_{\beta}:[\lambda_{\text{min}},\lambda_{\text{max}}]\rightarrow(0,1]$, with
$F_{\beta}(\lambda)=e^{-\beta(\lambda-\lambda_{\rm min})}$, satisfies the
assumptions of Lemma 10 with $L$ given by Eq. (19) for
$\max_{\lambda\in[\lambda_{\text{min}},\lambda_{\text{max}}]}\big{|}\alpha\,f^{(L+1)}(\lambda)\big{|}=\alpha\,\beta^{L+1}$.
Its Taylor coefficients are $a_{l}=\frac{e^{\beta\lambda_{\rm
min}}(-\beta)^{l}}{l!}$, for all $l\in\mathbb{N}$. To obtain $\alpha$, we note
that
$\sum_{l=0}^{L}|a_{l}/\big{(}1-\frac{2\delta}{\pi}\big{)}^{l}|\leq\sum_{l=0}^{\infty}|a_{l}/\big{(}1-\frac{2\delta}{\pi}\big{)}^{l}|=e^{\beta\lambda_{\rm
min}}e^{\frac{\beta}{1-2\delta/\pi}}$. Hence, by Eq. (20), we can take
$\alpha=e^{-\beta(\lambda_{\rm min}+\frac{1}{1-2\delta/\pi})}$. Introducing
$\gamma\coloneqq-\beta+\frac{\beta}{1-2\delta/\pi}=\frac{\beta\delta}{\pi/2-\delta}$,
we re-write $\alpha=e^{-\beta(1+\lambda_{\rm min})-\gamma}$. This allows us to
specify $e^{-\gamma}$ instead of the Fourier convergence interval, i.e. to
subordinate $\delta$ to the desired $\gamma$. This is done by fixing
$\delta=\frac{\pi}{2}\frac{1}{1+\frac{\beta}{\gamma}}$. This, together with
Eq. (18), leads to Eq. (3). ∎

##  VI Methods: Minimum query complexity of QITE primitives based on block-
encoding oracles

Our proof strategy for Theorem 5 is analogous to that of the no-fast-
forwarding theorem for real-time evolutions [40, 41, 42]. That is, it is based
on a reduction to QITE of the task of determining the parity ${\rm
par}(\boldsymbol{x})\coloneqq x_{0}\oplus x_{1}\oplus\cdots\oplus x_{N-1}$,
with $\oplus$ the bit-wise sum, of an unknown $N$-bit string
$\boldsymbol{x}\coloneqq x_{0}\,x_{1}\cdots x_{N-1}$ from a parity oracle
$U_{\boldsymbol{x}}$ for $\boldsymbol{x}$; together with known fundamental
complexity bounds for the latter task [69, 70]. More precisely, our proof
relies on three facts: $i$) No algorithm can find ${\rm par}(\boldsymbol{x})$
from $U_{\boldsymbol{x}}$ with fewer than a known number of queries to it [69,
70]; $ii$) a QITE primitive querying an oracle for an appropriate Hamiltonian
$H_{\boldsymbol{x}}$ gives an algorithm to find ${\rm par}(\boldsymbol{x})$;
and $iii$) a block-encoding oracle for $H_{\boldsymbol{x}}$ can be synthesized
from one call to $U_{\boldsymbol{x}}$. The three facts are established in the
following lemmas.

The first lemma, proven in [69], lower-bounds the complexity of any quantum
circuit able to obtain ${\rm par}(\boldsymbol{x})$ from queries to
$U_{\boldsymbol{x}}$. For our purposes, it can be stated as follows.

######  Lemma 11.

Let $C$ be a quantum circuit composed of $\boldsymbol{x}$-independent gates
and $q$ times the $\boldsymbol{x}$-dependent unitary

$U_{\boldsymbol{x}}\coloneqq\sum_{j=0}^{N}|j\rangle\langle j|\otimes
X^{x_{j}}\ ,$ (23)

with $\\{|j\rangle\\}_{j\in[N+1]}$ an orthogonal basis, such that, acting on
an $\boldsymbol{x}$-independent input state and upon measurement on an
$\boldsymbol{x}$-independent basis, outputs ${\rm par}(\boldsymbol{x})$ with
probability greater than $1/2$ for all $\boldsymbol{x}\in\\{0,1\\}^{N}$. Then,
$q\geq\lceil N/2\rceil$.

The second lemma, proven in App. P, reduces parity finding to QITE and is the
key technical contribution of this section.

######  Lemma 12.

Let $\beta>0$, $\varepsilon^{\prime}>0$, $\alpha\in(0,1]$, and
$\boldsymbol{x}$ an $N$-bit string such that

$\left|\frac{1-e^{-\frac{\beta}{2N}}}{2}\right|^{N}>\frac{2\,\varepsilon^{\prime}}{\alpha}\
.$ (24)

Then, there exists $H_{\boldsymbol{x}}$, with $\|H_{\boldsymbol{x}}\|\leq 1$,
such that a $(\beta,\varepsilon^{\prime},\alpha)$-QITE-primitive with calls to
a block-encoding oracle for $H_{\boldsymbol{x}}$, acting on an
$\boldsymbol{x}$-independent input state and upon measurement on an
$\boldsymbol{x}$-independent basis, outputs ${\rm par}(\boldsymbol{x})$ with
probability greater than $1/2$ for all $\boldsymbol{x}\in\\{0,1\\}^{N}$.

Finally, the missing link between Lemmas 11 and 12 is a sub-routine to query
$H_{\boldsymbol{x}}$ given queries to $U_{\boldsymbol{x}}$ in Eq. (23). This
is provided by the following lemma, proven in App. Q.

######  Lemma 13.

A block-encoding oracle for $H_{\boldsymbol{x}}$ can be generated from a
single query to $U_{\boldsymbol{x}}$ and $\mathcal{O}(N)$
$\boldsymbol{x}$-independent gates, for all $\boldsymbol{x}\in\\{0,1\\}^{N}$.
(See Fig. 15 for circuit.)

###### Proof of Theorem 5.

With these three Lemmas, the proof of Theo. 5 is straightforward. First, note
that the left-hand side of Eq. (24) decreases monotonically in $N$. Hence, for
any fixed $(\beta,\varepsilon^{\prime},\alpha)$, the largest $N\in\mathbb{N}$
that satisfies Eq. (24) is $N=\lfloor 2\tilde{q}\rfloor$, with
$\tilde{q}\in\mathbb{R}$ defined by Eq. (5). This, together with Lemmas 12 and
13, implies that any $(\beta,\varepsilon^{\prime},\alpha)$-QITE-primitive
synthesized from queries to the parity oracle $U_{\boldsymbol{x}}$ provides a
quantum circuit to determine the parity of any string $\boldsymbol{x}$ of
length $\lfloor 2\tilde{q}\rfloor$. Then, by virtue of Lemma 11, the query
complexity $q_{\min}(\beta,\varepsilon^{\prime},\alpha)\in\mathbb{N}$ of the
primitive cannot be smaller than $q\geq\lceil
N/2\rceil=\left\lceil\frac{\lfloor 2\tilde{q}\rfloor}{2}\right\rceil$. Note
that this number is the nearest integer to $\tilde{q}$. Ergo,
$q_{\min}(\beta,\varepsilon^{\prime},\alpha)\geq\tilde{q}$. ∎

Finally, a comment on why Theo. 5 does not hold for QITE primitives based on
real-time evolution (RTE) oracles is useful at this point. The reason is that,
by virtue of the RTE no-fast-forwarding theorem [40, 41, 42], a single call to
an RTE oracle suffices to find the parity with probability greater than $1/2$.
Hence, it is the query complexity RTE oracles themselves what is lower-bounded
by the parity considerations above, but not that of RTE-based QITE primitives.
It is an open question whether similar bounds can be obtained for RTE-based
QITE primitives by other arguments.

###### Acknowledgements.

We thank Martin Kliesch for for helpful comments that led us to write App. H.
We acknowledge financial support from the Serrapilheira Institute (grant
number Serra-1709-17173), and the Brazilian agencies CNPq (PQ grant No.
305420/2018-6) and FAPERJ (PDR10 E-26/202.802/2016 and JCN E-26/202.701/2018).
MMT acknowledges also the Government of Spain (FIS2020-TRANQI and Severo Ochoa
CEX2019-000910-S), Fundació Cellex, Fundació Mir-Puig, Generalitat de
Catalunya (CERCA, AGAUR SGR 1381) and ERC AdG CERQUTE.

## References

  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. Benjamin, and X. Yuan, Variational ansatz-based quantum simulation of imaginary time evolution, NPJ Quantum Information 5, 75 (2019).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nat. Phys. 16, 205 (2020).
  * Gomes  _et al._ [2020] N. Gomes, F. Zhang, N. F. Berthusen, C.-Z. Wang, K.-M. Ho, P. P. Orth, and Y. Yao, Efficient step-merged quantum imaginary time evolution algorithm for quantum chemistry, J. Chem. Theory Comput. 10, 6256 (2020).
  * Sun  _et al._ [2021] S.-N. Sun, M. Motta, R. N. Tazhigulov, A. T. Tan, G. K.-L. Chan, and A. J. Minnich, Quantum computation of finite-temperature static and dynamical properties of spin systems using quantum imaginary time evolution, PRX Quantum 2, 010317 (2021).
  * Nishi  _et al._ [2021] H. Nishi, T. Kosugi, and Y. Matsushita, Implementation of quantum imaginary-time evolution method on nisq devices by introducing nonlocal approximation, NPJ Quantum Information 7, 85 (2021).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Sampling from the Thermal Quantum Gibbs State and Evaluating Partition Functions with a Quantum Computer, Physical Review Letters 103, 220502 (2009), arXiv:0905.2199 .
  * Bilgin and Boixo [2010] E. Bilgin and S. Boixo, Preparing thermal states of quantum systems by dimension reduction, Phys. Rev. Lett. 105, 170405 (2010).
  * Temme  _et al._ [2011] K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, Quantum metropolis sampling, Nature 471, 87 (2011).
  * Yung and Aspuru-Guzik [2012] M.-H. Yung and A. Aspuru-Guzik, A quantum–quantum metropolis algorithm, Proceedings of the National Academy of Sciences 109, 754 (2012).
  * Kastoryano and Brandão [2016] M. J. Kastoryano and F. G. S. L. Brandão, Quantum gibbs samplers: the commuting case, Communications in Mathematical Physics 344, 915 (2016).
  * Chowdhury and Somma [2017] A. N. Chowdhury and R. D. Somma, Quantum algorithms for gibbs sampling and hitting-time estimation, Quantum Inf. Comput. 17, 41 (2017).
  * Brandão  _et al._ [2017] F. G. S. L. Brandão, A. Kalev, T. Li, C. Y.-Y. Lin, K. M. Svore, and X. Wu, Quantum sdp solvers: Large speed-ups, optimality, and applications to quantum learning, Eprint: 1710.02581 (2017).
  * Brandão and Svore [2017] F. G. S. L. Brandão and K. M. Svore, Quantum speed-ups for semidefiniteprogramming, Proceedings of 58th IEEEFOCS , 415 (2017).
  * van Apeldoorn  _et al._ [2020] J. van Apeldoorn, A. Gilyén, S. Gribling, and R. de Wolf, Quantum SDP-Solvers: Better upper and lower bounds, Quantum 4, 230 (2020).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) p. 193.
  * [16] Y. Wang, G. Li, and X. Wang, Variational quantum gibbs state preparation with a truncated taylor series, Eprint: 2005.08797 .
  * Kuwahara  _et al._ [2020] T. Kuwahara, K. Kato, and F. G. S. L. Brandão, Clustering of conditional mutual information for quantum gibbs states above a threshold temperature, Phys. Rev. Lett. 124, 220601 (2020).
  * Chowdhury  _et al._ [2021] A. N. Chowdhury, R. D. Somma, and Y. Subasi, Computing partition functions in the one clean qubit model, Phys. Rev. A 103, 032422 (2021).
  * Kempe  _et al._ [2006] J. Kempe, A. Kitaev, and O. Regev, The complexity of the local hamiltonian problem, SIAM Journal of Computing 35(5), 1070 (2006).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum computation by adiabatic evolution, MIT-CTP-2936 (2000).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, A quantum approximate optimization algorithm, MIT-CTP-4610 (2014).
  * Montanaro [2015] A. Montanaro, Quantum walk speedup of backtracking algorithms, Eprint: 1509.02374 (2015).
  * et al [2021] M. P. H. et al, Quantum approximate optimization of non-planar graph problems on a planar superconducting processor, Nat. Phys. (2021).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 4213 (2014).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Quantum chemistry calculations on a trapped-ion quantum simulator, Phys. Rev. X 8, 031022 (2018).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 3007 (2019).
  * Kieferova and Wiebe [2017] M. Kieferova and N. Wiebe, Tomography and generative data modeling via quantum boltzmann training, Phys. Rev. A 96, 062327 (2017).
  * Wiebe and Wossnig [2019] N. Wiebe and L. Wossnig, Generative training of quantum boltzmann machines with hidden units, Eprint: 1905.09902 (2019).
  * Zoufal  _et al._ [2020] C. Zoufal, A. Lucchi, and S. Woerner, Variational quantum boltzmann machines, Eprint: 2006.06004 (2020).
  * nd Peter Wittek  _et al._ [2017] J. B. nd Peter Wittek, N. Pancotti, N. W. Patrick Rebentrost, and S. Lloyd, Quantum machine learning, Nature 549, 195 (2017).
  * Amin  _et al._ [2018] M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, Quantum boltzmann machine, Phys. Rev. X 8, 021050 (2018).
  * Benedetti  _et al._ [2021] M. Benedetti, M. Fiorentini, and M. Lubasch, Hardware-efficient variational quantum algorithms for time evolution, Phys. Rev. Research 3, 033083 (2021).
  * Lin  _et al._ [2021] S.-H. Lin, R. Dilip, A. G. Green, A. Smith, and F. Pollmann, Real- and imaginary-time evolution with compressed quantum circuits, PRX Quantum 2, 010342 (2021).
  * Cao  _et al._ [2021] C. Cao, Z. An, S.-Y. Hou, D. L. Zhou, and B. Zeng, Quantum imaginary time evolution steered by reinforcement learning (2021), arXiv:2105.08696 [quant-ph] .
  * Brassard  _et al._ [2002] G. Brassard, P. Høyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, in  _Quantum Computation and Quantum Information_, Contemporary Mathematics, Vol. 305, edited by S. J. Lomonaco and H. E. Brandt (American Mathematical Society, Providence, Rhode Island, 2002) pp. 53–74, arXiv:quant-ph/0005055 [quant-ph] .
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Optimal hamiltonian simulation by quantum signal processing, Phys. Rev. Lett. 118, 010501 (2017).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian Simulation by Qubitization, Quantum 3, 163 (2019).
  * [38] M. Abramowitz, I. A. Stegun, and Others, “Handbook of mathematical functions,” Applied mathematics series 55, 62 (1966).
  * Silva  _et al._ [2022] T. L. Silva, L. Borges, and L. Aolita, Quantum signal processing for fourier approximations (2022), paper in preparation.
  * Berry  _et al._ [2007] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Efficient Quantum Algorithms for Simulating Sparse Hamiltonians, Communications in Mathematical Physics 270, 359 (2007).
  * Berry  _et al._ [2014] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Exponential improvement in precision for simulating sparse Hamiltonians, in  _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ (ACM, New York, NY, USA, 2014) pp. 283–292, arXiv:1312.1414 .
  * Berry  _et al._ [2015a] D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters, 2015 IEEE 56th Annual Symposium on Foundations of Computer Science , 792 (2015a), arXiv:1501.01715 .
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Solvable model of a spin glass, Phys. Rev. Lett. 35, 1792 (1975).
  * Panchenko [2012] D. Panchenko, The Sherrington-Kirkpatrick model: an overview, J. Stat. Phys. 149, 362 (2012).
  * Altshuler  _et al._ [2009] B. Altshuler, H. Krovi, and J. Roland, Adiabatic quantum optimization fails for random instances of np-complete problems, Arxiv: 0908.2782 (2009).
  * Altshuler  _et al._ [2010] B. Altshuler, H. Krovi, and J. Roland, Anderson localization makes adiabatic quantum optimization fail, Proceedings of the National Academy of Sciences 107(28), 12446 (2010).
  * Berry  _et al._ [2015b] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating Hamiltonian Dynamics with a Truncated Taylor Series, Physical Review Letters 114, 090502 (2015b), arXiv:arXiv:1412.4687 .
  * Campbell [2019] E. Campbell, Random compiler for fast hamiltonian simulation, Phys. Rev. Lett. 123, 070503 (2019).
  * Childs  _et al._ [2019] A. M. Childs, A. Ostrander, and Y. Su, Faster quantum simulation by randomization, Quantum 3 (2019).
  * Lloyd  _et al._ [2021] S. Lloyd, B. T. Kiani, D. R. M. Arvidsson-Shukur, S. Bosch, G. D. Palma, W. M. Kaminsky, Z.-W. Liu, and M. Marvian, Hamiltonian singular value transformation and inverse block encoding, ArXiv: 2104.01410 (2021).
  * [51] L. Borges, T. L. Silva, and L. Aolita, in preparation (2022) .
  * Fraser [1965] W. Fraser, A survey of methods of computing minimax and near-minimax polynomial approximations for functions of a single independent variable, Journal of the Association for Computing machinery 12, 295 (1965).
  * Elliott  _et al._ [1987] D. Elliott, D. Paget, G. Phillips, and P. Taylor, Error of truncated Chebyshev series and other near minimax polynomial approximations, Journal of Approximation Theory 50, 49 (1987).
  * Luke [1972] Y. L. Luke, Inequalities for generalized hypergeometric functions, Journal of Approximation Theory 5, 41 (1972).
  * Karp [1972] R. M. Karp, Reducibility among combinatorial problems, in Miller, R.; Thacher, J. W. (eds.), Complexity of Computer Computation, Plenum Press , 85 (1972).
  * Trevisan  _et al._ [2000] L. Trevisan, G. Sorkin, M. Sudan, and D. Williamson, Gadgets, approximation, and linear programming, Proceedings of the 37th IEEE Symposium on Foundations of Computer Science , 617 (2000).
  * Long and Servedio [2010] P. M. Long and R. A. Servedio, Restricted boltzmann machines are hard to approximately evaluate or simulate, ICML’10: Proceedings of the 27th International Conference on International Conference on Machine Learning , 703 (2010).
  * Montanari [2019] A. Montanari, Optimization of the Sherrington-Kirkpatrick Hamiltonian, Proceedings of the 60th Annual Symposium on Foundations of Computer Science (FOCS ’19) , 1417 (2019).
  * Farhi  _et al._ [2019] E. Farhi, J. Goldstone, S. Gutmann, and L. Zhou, The quantum approximate optimization algorithm and the Sherrington-Kirkpatrick model at infinite size, Arxiv: 1910.08187 (2019).
  * Silva  _et al._ [2021] T. L. Silva, M. M. Taddei, S. Carrazza, and L. Aolita, scarrazza/qite: v1.0.0 (2021).
  * Giurgica-Tiron  _et al._ [a] T. Giurgica-Tiron, I. Kerenidis, F. Labib, A. Prakash, and W. Zeng, Low depth algorithms for quantum amplitude estimation, arXiv:2012.03348 (a).
  * Giurgica-Tiron  _et al._ [b] T. Giurgica-Tiron, S. Johri, I. Kerenidis, J. Nguyen, N. Pisenti, A. Prakash, K. Sosnova, K. Wright, and W. Zeng, Low depth amplitude estimation on a trapped ion quantum computer, ArXiv:2109.09685 (b).
  * Low  _et al._ [2016] G. H. Low, T. J. Yoder, and I. L. Chuang, Methodology of resonant equiangular composite quantum gates, Phys. Rev. X 6, 041067 (2016).
  * Haah [2019] J. Haah, Product Decomposition of Periodic Functions in Quantum Signal Processing, Quantum 3, 190 (2019).
  * Chao  _et al._ [2020] R. Chao, D. Ding, A. Gilyen, C. Huang, and M. Szegedy, Finding angles for quantum signal processing with machine precision (2020), arXiv:2003.02831 [quant-ph] .
  * Dong  _et al._ [2020] Y. Dong, X. Meng, K. B. Whaley, and L. Lin, Efficient phase factor evaluation in quantum signal processing (2020), arXiv:2002.11649 [quant-ph] .
  * Pérez-Salinas  _et al._ [2021] A. Pérez-Salinas, D. López-Núñez, A. García-Sáez, P. Forn-Díaz, and J. I. Latorre, One qubit as a universal approximant, Phys. Rev. A 104, 012405 (2021).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017).
  * Farhi  _et al._ [1998] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Limit on the Speed of Quantum Computation in Determining Parity, Physical Review Letters 81, 5442 (1998), arXiv:quant-ph/9802045 .
  * Beals  _et al._ [1998] R. Beals, H. Buhrman, R. Cleve, M. Mosca, and R. de Wolf, Quantum lower bounds by polynomials, in  _Proceedings 39th Annual Symposium on Foundations of Computer Science (Cat. No.98CB36280)_, Vol. 48 (IEEE Comput. Soc, 1998) pp. 352–361, arXiv:9802049 [quant-ph] .

##  Appendix A Operator-function design from imperfect oracles

Both algorithms for operator-function design we present assume access to ideal
oracles as given by Def. 2 and Def. 3. This is not the case in an experimental
scenario where only approximate oracles are available. This leads to a total
error in the algorithm that is proportional to the number of approximate
oracle calls. The following lemma deals with that error for generic operator
functions.

######  Lemma 14.

(Primitives from imperfect oracles) Let a circuit with $q$ queries to an
oracle $O$ for $H$ generate an $(\varepsilon,\alpha)$-block-encoding
$V_{f(H)}$ of $f(H)$, for arbitrary $f$. If $O$ is substituted by $\tilde{O}$,
with $\|O-\tilde{O}\|\leq\varepsilon_{O}$, the circuit generates an
$(\tilde{\varepsilon},\alpha)$-block-encoding $\tilde{V}_{f(H)}$ of $f(H)$
with $\tilde{\varepsilon}\leq\varepsilon+q\,\varepsilon_{O}$.

###### Proof.

For $q\,\varepsilon_{O}\ll 1$, straightforward matrix multiplication shows
that the total error is
$||V_{f(H)}-\tilde{V}_{f(H)}||=\mathcal{O}(q\,\varepsilon_{O})$. This implies
$||f(H)-\langle
0|\tilde{V}_{f(H)}|0\rangle||\leq\varepsilon_{P}+q\,\varepsilon_{O}$ by virtue
of the triangle inequality. This concludes the proof. ∎

As an exemplary application of Lemma 14, we calculate the gate complexity
$g_{\tilde{O}_{2}}$ required to implement an approximate real-time evolution
oracle $\tilde{O}_{2}$ for QITE primitive $P_{2}$ by Trotter-Suzuki-like
simulation methods. To attain a total error $\varepsilon^{\prime}$, we choose
$\varepsilon=\varepsilon^{\prime}/2$ and
$\varepsilon_{O}=\varepsilon^{\prime}/(2q)$. For real time $t$, the gate
complexity of state-of-the-art algorithms [48, 49] based on first-order
product formulae is $\mathcal{O}(t^{2}/\varepsilon_{O})$. Substituting for
$\varepsilon_{O}$, and then for $t$ and $q$ with the real time and query
complexity given in Theorem 2, we obtain the oracle gate complexity
$g_{\tilde{O}_{2}}=\mathcal{O}\left(\frac{\pi^{2}\beta}{\gamma}\frac{\ln(8/\varepsilon^{\prime})}{2\,\varepsilon^{\prime}(1+\gamma/\beta)}\right)$.
This scaling is exponentially worse in $1/\varepsilon^{\prime}$ than what we
would get if we simulated $O_{2}$ with more sophisticated real-time evolution
algorithms [42, 47, 36, 37, 15]. However, those algorithms assume more
powerful oracles, which require significantly more ancillary qubits. The fact
that $P_{2}$ can be synthesised with only 1 ancillary qubit throughout and
with no qubitization makes it appealing for near-term implementations on
intermediate-sized quantum hardware.

##  Appendix B Post-selection probability propagation onto output state error

Since $\alpha\,F_{\beta}(H)$ is not unitary, the (operator-norm) error
$\varepsilon^{\prime}$ in its approximation by an
$(\alpha,\varepsilon^{\prime})$-block-encoding gets amplified at the post-
selection due to state renormalization. The following allows us to control the
error in the output state.

######  Lemma 15.

(Error propagation from block encoding to output state) Let
$p_{\Psi}(\beta,\alpha)\leq 1$ be the post-selection probability of an
$(\alpha,\varepsilon^{\prime})$-block-encoding $U_{F_{\beta}(H)}$ of
$F_{\beta}(H)$ on an input state $|\Psi\rangle\in\mathbb{H}_{\mathcal{S}}$.
Then, if
$\varepsilon^{\prime}\leq\,\varepsilon\,\sqrt{p_{\Psi}(\beta,\alpha)}/2$ and
$\varepsilon\ll 1$, the output-state trace-distance error is
$\mathcal{O}(\varepsilon)$.

###### Proof.

By definition, $\|\langle 0|\
U_{F_{\beta}(H)}|0\rangle-\alpha\,F_{\beta}(H)\|\leq\varepsilon^{\prime}$.
Then, for some $|\Xi\rangle\in\mathbb{H}_{S}$, with
$\||\Xi\rangle\|\leq\varepsilon^{\prime}$, it is

$\langle 0|\ U_{F_{\beta}(H)}|0\rangle\
|\Psi\rangle=\alpha\,F_{\beta}(H)|\Psi\rangle+|\Xi\rangle$ (25a) and, so, for
some $\epsilon\in\mathbb{R}$, with
$|\epsilon|\leq\||\Xi\rangle\|\leq\varepsilon^{\prime}$, we get $\|\langle 0|\
U_{F_{\beta}(H)}|0\rangle\
|\Psi\rangle\|=\alpha\,\|F_{\beta}(H)|\Psi\rangle\|+\epsilon.$ (25b)

Next, we Taylor-expand the output state in terms of $\varepsilon$. To that
end, note first that

$\frac{\langle 0|\ U_{F_{\beta}(H)}|0\rangle\ |\Psi\rangle}{\|\langle 0|\
U_{F_{\beta}(H)}|0\rangle\
|\Psi\rangle\|}=\frac{\alpha\,F_{\beta}(H)|\Psi\rangle+|\Xi\rangle}{\alpha\,\|F_{\beta}(H)|\Psi\rangle\|+\epsilon}\\\
=\frac{F_{\beta}(H)|\Psi\rangle}{\|F_{\beta}(H)|\Psi\rangle\|}\left(1-\frac{\epsilon}{\alpha\,\|F_{\beta}(H)|\Psi\rangle\|}\right)\\\
+\frac{|\Xi\rangle}{\alpha\,\|F_{\beta}(H)|\Psi\rangle\|}+\mathcal{O}\left(\frac{\varepsilon^{\prime}\,|\Xi\rangle}{\alpha^{2}\,\|F_{\beta}(H)|\Psi\rangle\|^{2}}\right)\
.$ (26)

So, the output state-vector error is upper-bounded as

$\displaystyle\left\|\frac{\langle 0|\ U_{F_{\beta}(H)}|0\rangle\
|\Psi\rangle}{\|\langle 0|\ U_{F_{\beta}(H)}|0\rangle\
|\Psi\rangle\|}\right.-$
$\displaystyle\left.\frac{F_{\beta}(H)|\Psi\rangle}{\|F_{\beta}(H)|\Psi\rangle\|}\right\|$
(27)
$\displaystyle\leq\frac{2\,\varepsilon^{\prime}}{\alpha\,\|F_{\beta}(H)|\Psi\rangle\|}$
$\displaystyle+\mathcal{O}\left(\frac{\varepsilon^{\prime
2}}{\alpha^{2}\,\|F_{\beta}(H)|\Psi\rangle\|^{2}}\right).$

Hence, if
$\varepsilon^{\prime}\leq\varepsilon\,\alpha\,\|F_{\beta}(H)|\Psi\rangle\|/2=\varepsilon\,\sqrt{p_{\Psi}(\beta,\alpha)}/2$,
the error in $l_{2}$-norm distance of the output state vector is at most
$\varepsilon+\mathcal{O}(\varepsilon^{2})$, which equals
$\mathcal{O}(\varepsilon)$ for $\varepsilon\ll 1$.

Figure 7:  Probabilistic versus coherent master QITE algorithms. a) The
probabilistic approach repeatedly applies the unitary $U_{F_{\beta}(H)}$
generated by the primitive (on independent preparations of
$|\Psi\rangle|0\rangle\in\mathbb{H}_{\mathcal{SA}}$) until the post-selection
is successful, i.e. the measurement on the ancillas returns $|0\rangle$ as
outcome. This takes on average
$\mathcal{O}\big{(}1/p_{\Psi}(\beta,\alpha)\big{)}$ repetitions. No a priori
knowledge of the input state $|\Psi\rangle$ is required (it can be fully
generic, even mixed). b) The coherent approach is based on quantum amplitude
amplification. It operates only on pure input states. So, if the input state
is mixed, an extra ancillary register $\mathcal{A}_{\mathrm{pur}}$ of
$|\mathcal{A}_{\mathrm{pur}}|=|\mathcal{S}|=N$ qubits is required to purify
it. This is the case in quantum Gibbs state sampling, where $|\Psi\rangle$ is
a purification of the maximally mixed state on $\mathbb{H}_{\mathcal{S}}$. The
primitive is repeated applied sequentially (on the same preparation of
$|\Psi\rangle|0\rangle$), interleaved with reflection operators
$R_{|\Phi_{0}\rangle}$ and $R_{|0\rangle}$ around the states
$|\Phi_{0}\rangle\coloneqq
U_{F_{\beta}(H)}|\Psi\rangle|0\rangle\in\mathbb{H}_{\mathcal{SA}}$ and
$|0\rangle\in\mathbb{H}_{\mathcal{A}}$, respectively. In practice, this
requires full a priori knowledge of $|\Psi\rangle$. The total number of
repetitions of the primitive is
$\mathcal{O}\big{(}1/\sqrt{p_{\Psi}(\beta,\alpha)}\big{)}$, after which the
desired output is obtained with probability close to 1. Hence, the coherent
master algorithm displays a significantly lower overall query complexity than
the probabilistic one. However, in return, the former requires much larger
circuit depth than the latter.

With the $l_{2}$-norm distance between the two state vectors, we can upper-
bound the trace distance between their corresponding rank-1 density operators
using well-known inequalities. First, note that the $l_{2}$-norm distance
between two arbitrary normalised vectors $|\phi\rangle$ and $|\psi\rangle$ can
be expressed as

$\||\phi\rangle-|\psi\rangle\|=\sqrt{2}\sqrt{1-{\rm
Re}\langle\psi|\phi\rangle}.$ (28)

Second, recall that the trace distance

$\||\phi\rangle\langle\phi|-|\psi\rangle\langle\psi|\|_{{\rm
tr}}\coloneqq\frac{1}{2}{\rm
tr}\sqrt{\left(|\phi\rangle\langle\phi|-|\psi\rangle\langle\psi|\right)^{2}}\
.$ (29)

between rank-1 density operators $|\phi\rangle\langle\phi|$ and
$|\psi\rangle\langle\psi|$ can be written in terms of the overlap as

$\left\||\phi\rangle\langle\phi|-|\psi\rangle\langle\psi|\right\|_{{\rm
tr}}=\sqrt{1-|\langle\psi|\phi\rangle|^{2}}\ .$ (30)

Then, note that, for all $|\phi\rangle$ and $|\psi\rangle$, the RHS of Eq.
(28) upper-bounds the RHS of Eq. (30). Hence,
$\||\phi\rangle\langle\phi|-|\psi\rangle\langle\psi|\|_{{\rm
tr}}\leq\||\phi\rangle-|\psi\rangle\|$. This, together with Eq. (27) gives the
promised trace-norm error of the output state. ∎

##  Appendix C Traditional master QITE algorithms

The average number of times a QITE primitive is applied in probabilistic and
coherent master QITE algorithms is
$1/p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)$ and
$1/\sqrt{p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)}$, respectively; see Fig.
7. Conveniently, for $\varepsilon\ll 1$,
$1/p_{\Psi}(\beta,\varepsilon^{\prime},\alpha)$ can be approximated by the
more practical expression $1/p_{\Psi}(\beta,\alpha)$ up to error
$\mathcal{O}(\varepsilon)$. This follows from a similar Taylor expansion to
that in Eq. (26). The probabilistic algorithm applies the primitive on
independent input state preparations [Fig. 7-a)]. In contrast, the coherent
one (Fig. 7-b) leverages quantum amplitude amplification [35], which we
briefly discuss next.

The coherent amplification process is realized by repeatedly applying to
$|\Phi_{0}\rangle\coloneqq U_{F_{\beta}(H)}|\Psi\rangle|0\rangle$ the unitary
operator

$U_{\mathrm{eng}}\coloneqq R_{|\Phi_{0}\rangle}\,R_{|0\rangle},$ (31)

where $R_{|\Phi_{0}\rangle}$ and $R_{|0\rangle}$ are respectively the
reflection operators around
$|\Phi_{0}\rangle\in\mathbb{H}_{\mathcal{S}}\otimes\mathbb{H}_{\mathcal{A}_{\rm
pur}}\otimes\mathbb{H}_{\mathcal{A}}$ and
$|0\rangle\in\mathbb{H}_{\mathcal{A}}$. The former reflection can in turn be
decomposed as
$R_{|\Phi_{0}\rangle}=U_{F_{\beta}(H)}\,R_{|\Psi\rangle|0\rangle}\,U^{\dagger}_{F_{\beta}(H)}$,
where $U^{\dagger}_{F_{\beta}(H)}$ is the inverse of the block-encoding
$U_{F_{\beta}(H)}$ of the QITE propagator $F_{\beta}(H)$ and
$R_{|\Psi\rangle|0\rangle}$ the reflection around $|\Psi\rangle|0\rangle$.
Unitary $U_{\mathrm{eng}}$ is sometimes referred to as the amplification
engine. Importantly, $U_{\mathrm{eng}}$ acts on the 2-dimensional subspace
spanned by $|\Phi_{0}\rangle$ and $|\Phi_{\rm
target}\rangle\propto|0\rangle\langle 0|\,|\Phi_{0}\rangle$ as an SU$(2)$
rotation. For $k\in\mathbb{N}$, it gives

$\displaystyle U_{\mathrm{eng}}^{k}\,|\Phi_{0}\rangle$ $\displaystyle=$
$\displaystyle\sin[(2k+1)\theta]\,|\Phi_{\rm target}\rangle$ (32)
$\displaystyle+$ $\displaystyle\cos[(2k+1)\theta]\,|\Phi_{\perp}\rangle,$

where $|\Phi_{\perp}\rangle\propto|1\rangle\langle 1|\,|\Phi_{0}\rangle$ and
$\sin(\theta)\coloneqq\alpha\|F_{\beta}(H)|\Psi\rangle\|$. Hence, taking
$k=k_{\rm opt}=\mathcal{O}(1/\theta)$ yields $\sin[(2k+1)\theta]\approx 1$ and
therefore probability close to 1 for desired output. For $\theta\ll 1$, this
entails
$k=\mathcal{O}(1/\alpha\|F_{\beta}(H)|\Psi\rangle\|)=\mathcal{O}(1/\sqrt{p_{\Psi}(\beta,\alpha)})$
repetitions of the primitive, as in Eq. (1).

Finally, since $\alpha^{2}\|F_{\beta}(H)|\Psi\rangle\|$ is in general unknown,
one has no a priori knowledge of $k_{\rm opt}$. However, fortunately, this can
be accounted for with successive attempts with $k$ randomly chosen within a
range of values that grows exponentially with the number of attempts (see [35,
Theorem 3]). Remarkably, the resulting average number of applications of the
primitive remains within $\mathcal{O}(1/\sqrt{p_{\Psi}(\beta,\alpha)})$.

##  Appendix D Optimal sub-normalisation for QITE Primitive 2

In this appendix, we prove that Eq. (4) gives approximately the optimal value
for the subnormalization of $P_{2}$, such that the overall query complexity
given by Eq. (1) is minimized. Notice that, by diminishing the value of
$\gamma$ we increase the success probability of $P_{2}$, decresing the number
of times it needs to be realized. At the same time, it leads to an increment
on the query complexity of $P_{2}$ given by Eq. (3). Therefore, the optimal
subnormalization is a tradeoff between these two contributions.

From Eqs. (1) and (3), given $\beta$, $\varepsilon$, and the master algorithm
type $\kappa$, we get that the optimal $\gamma$ minimizes

$\begin{split}Q_{\kappa}(\beta,\varepsilon,\gamma)=\frac{e^{2\mu_{\kappa}\gamma}}{\|F_{\beta}(H)|\Psi\rangle\|^{2\mu_{\kappa}}}&\left(\frac{\beta}{\gamma}+1\right)\\\
&\bigg{[}\ln\bigg{(}\frac{8}{\|F_{\beta}(H)|\Psi\rangle\|\varepsilon}\bigg{)}+g\gamma\bigg{]},\end{split}$
(33)

with $g=1$ introduced for convenience. Here we have used
$p_{\Psi}(\beta,\gamma)=e^{-2\gamma}\|F_{\beta}(H)|\Psi\rangle\|$ and
$\varepsilon^{\prime}=p_{\Psi}(\beta,\gamma)\varepsilon/2$.

Instead, let $\gamma=\gamma_{\kappa}(\beta)$ be such that
$Q_{\kappa}(\beta,\varepsilon,\gamma)|_{g=0}$ is minimized. By solving
$\frac{\partial}{\partial\gamma}Q_{\kappa}(\beta,\varepsilon,\gamma)\big{|}_{g=0}=0$
we obtain $\gamma_{\kappa}(\beta)$ given by Eq. (4). In order to prove that
that expression is a good approximation for the actual optimal value
$\gamma^{(\text{opt})}_{\kappa}(\beta)$, first notice that, with $g=1$, it is
easy to verify that
$Q_{\kappa}(\beta,\varepsilon,\gamma)|_{g=1}>Q_{\kappa}(\beta,\varepsilon,\gamma_{\kappa}(\beta))|_{g=1}$
if $\gamma>\gamma_{\kappa}(\beta)$. Therefore,
$\gamma^{(\text{opt})}_{\kappa}(\beta)<\gamma_{\kappa}(\beta)$. Defining
$\Delta Q\coloneqq
Q_{\kappa}(\beta,\varepsilon,\gamma^{(\text{opt})}_{\kappa}(\beta))|_{g=1}-Q_{\kappa}(\beta,\varepsilon,\gamma_{\kappa}(\beta))|_{g=1}$,
a Taylor expansion of $Q_{\kappa}(\beta,\varepsilon,\gamma)|_{g=1}$ at
$\gamma=\gamma_{\kappa}(\beta)$ shows that, up to first order in
$\Delta\gamma=\gamma^{(\text{opt})}_{\kappa}(\beta)-\gamma_{\kappa}(\beta)$,
we have

$\frac{\Delta
Q}{Q_{\kappa}(\beta,\varepsilon,\gamma_{\kappa}(\beta))|_{g=1}}\simeq\frac{\Delta\gamma}{\ln\bigg{(}\frac{8e^{2\gamma_{\kappa}(\beta)}}{\|F_{\beta}(H)|\Psi\rangle\|\varepsilon}\bigg{)}},$
(34)

which is much smaller than $1$ provided that $\varepsilon^{\prime}\ll 1$,
since $\gamma_{\kappa}(\beta)<1$ and, therefore, $\Delta\gamma<1$.

##  Appendix E Proof of Theorem 4

Let us first prove a convenient auxiliary lemma. Ideally, the input and output
states of the $l$-th fragment of should be
$|\Psi_{l-1}\rangle=\frac{F_{\beta_{l-1}}(H)|\Psi\rangle}{\|F_{\beta_{l-1}}(H)|\Psi\rangle\|}$
and
$|\Psi_{l}\rangle=\frac{F_{\Delta\beta_{l}}(H)|\Psi_{l-1}\rangle}{\|F_{\Delta\beta_{l}}(H)|\Psi_{l-1}\rangle\|}=\frac{F_{\beta_{l}}(H)|\Psi\rangle}{\|F_{\beta_{l}}(H)|\Psi\rangle\|}$,
respectively. However, the actual operator that each
$P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}$ perfectly block-
encodes is

$\alpha_{l}\tilde{F}_{\Delta\beta_{l}}(H)\coloneqq\langle 0|\
U_{F_{\Delta\beta_{l}}(H)}|0\rangle=\alpha\,F_{\Delta\beta_{l}}(H)+E_{l},$
(35)

for some $E_{l}$ on $\mathbb{H}_{\mathcal{S}}$ with
$||E_{l}||=:\varepsilon_{l}^{\prime}$. Hence, in analogy to Eq. (25b), the
actual input and output states are
$|\tilde{\Psi}_{l-1}\rangle=|{\Psi}_{l-1}\rangle+|\Xi_{l-1}\rangle$ and
$|\tilde{\Psi}_{l}\rangle=|{\Psi}_{l}\rangle+|\Xi_{l}\rangle$, for some “error
states” $|\Xi_{l-1}\rangle$ and $|\Xi_{l}\rangle$, respectively. (Note that
$|\Xi_{0}\rangle=0$.) The following lemma controls the output-state error in
terms of the errors in the input state and block encoding.

######  Lemma 16.

(Error propagation from input state and block encoding to output state) For
all $l\in[r]$, let $\varepsilon_{l-1}\coloneqq\||\Xi_{l-1}\rangle\|$ and
$\varepsilon_{l}^{\prime}$ be the input-state and block-encoding errors,
respectively; and $\varepsilon_{l}$ the tolerated output-state error. If
$\varepsilon_{l-1}+\varepsilon_{l}^{\prime}\leq\varepsilon_{l}\,\|\alpha_{l}\,F_{\Delta\beta_{l}}(H)|\Psi_{l-1}\rangle\|/2$
and $\varepsilon_{l}\ll 1$, then
$\||{\Xi}_{l}\rangle\|=\mathcal{O}(\varepsilon_{l})$.

###### Proof.

Similar to the proof of Lemma 15 but where also the input state is
approximate. ∎

Note that the proof of Lemma 15 (and therefore also that of Lemma 16) makes no
use of the specific form of $F_{\beta}$ (or $F_{\Delta\beta_{l}}$).
Consequently, both lemmas hold for $(\alpha,\varepsilon^{\prime})$-block-
encodings of any operator function, not just the QITE propagator. We are now
in a good position to prove Theorem 4.

###### Proof of Theorem 4.

We begin by proving soundness. The proof consists of showing that Eqs. (6)
imply that Lemma 16 holds for all $l\in[r]$ and gives
$\varepsilon_{r}\coloneqq\||{\Xi}_{r}\rangle\|=\mathcal{O}(\varepsilon)$. To
see this, apply the lemma from $l=1$ till $l=r$, with $\varepsilon_{0}=0$, and
iteratively use the property

$\displaystyle p_{\Psi_{l-1}}(\Delta\beta_{l})$ $\displaystyle=$
$\displaystyle\|F_{\Delta\beta_{l}}(H)|\Psi_{l-1}\rangle\|$ (36)
$\displaystyle=$
$\displaystyle\frac{\|F_{\beta_{l}}(H)|\Psi\rangle\|}{\|F_{\beta_{l-1}}(H)|\Psi\rangle\|}$
$\displaystyle=$
$\displaystyle\frac{p_{\Psi}(\beta_{l})}{p_{\Psi}(\beta_{l-1})}.$

We next prove complexity. To this end, we must show the validity of Eq. (7).
By Def. 5, the overall query complexity of Algorithm 1 is the sum over the
query complexities of each primitive applied. Each primitive
$P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}$ has complexity
$q(\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l})$. Hence, we must show
that the average number of times $n_{l}$ that each
$P_{\Delta\beta_{l},\varepsilon^{\prime}_{l},\alpha_{l}}$ is run is given by
$\frac{p_{\Psi}(\beta_{l-1})}{p_{\Psi}(\beta)\prod_{k=l}^{r}\alpha^{2}_{k}}$.
To see this, note that, by definition, it is
$n_{r}\coloneqq\big{(}{\alpha_{r}}^{2}\,p_{\Psi_{r-1}}(\Delta\beta_{r})\big{)}^{-1}$,
$n_{r-1}\coloneqq
n_{r}\times\big{(}{\alpha_{r-1}}^{2}\,p_{\Psi_{r-2}}(\Delta\beta_{r-1})\big{)}^{-1}$,
$\ldots$, and $n_{1}\coloneqq
n_{2}\times\big{(}{\alpha_{1}}^{2}\,p_{\Psi}(\Delta\beta_{1})\big{)}^{-1}$.
Then, use the latter together with Eq. (36) to get
$n_{l}=\frac{p_{\Psi}(\beta_{l-1})}{p_{\Psi}(\beta)\prod_{k=l}^{r}\alpha^{2}_{k}}$
for all $l\in[r]$. ∎

Interestingly, we note that the only specific detail about $F_{\beta}$ that
the proof of Theorem 4 uses is the fact that
$F_{\beta}(H)=\prod_{l=1}^{r}F_{\Delta\beta_{l}}(H)$ (this is required for Eq.
(36) to hold). Consequently, Theorem 4 applies not only to fragmented QITE but
actually also to the fragmentation of any operator function in terms of
suitable factors.

##  Appendix F Validity of the query complexity bounds of QITE primitives for
low beta

As discussed in Secs. III.3 and III.4, the first steps of fragmented QITE
involve small inverse temperatures, often with $\Delta\beta_{1}<1$. However,
Eqs. (2) and (3) are in principle only asymptotic upper bounds for the actual
query complexities. Hence, it is licit to question how valid they are to
access the complexity of the first steps of Alg. 1. In this appendix, we
discuss the validity and tightness for all $\beta>0$ of

$\tilde{q}_{1}(\beta,\varepsilon^{\prime})\coloneqq
2\left(\frac{e\beta}{2}+\frac{\ln\left(1/\varepsilon^{\prime}\right)}{\ln\left(e+2\ln(1/\varepsilon^{\prime})/e\beta\right)}\right)$
(37)

and

$\tilde{q}_{2}(\beta,\varepsilon^{\prime})\coloneqq
4\left(\beta/\gamma+1\right)\ln(4/\varepsilon^{\prime})$ (38)

as exact expressions for the query complexities of primitives 1 and 2,
respectively. These formulas are the ones used in our numeric experiments.
First of all, notice that although these expressions are continuous functions,
the actual query complexities used by Algs. 2 and 3 take only even integer
values. Thus, we take the value of the query for each round of fragmentation
as $2\lceil\tilde{q}_{k}(\beta,\varepsilon^{\prime})/2\rceil$, $k=1,2$. In
particular, for any value of $\beta$ such that
$0<\tilde{q}_{1}(\beta,{\varepsilon^{\prime})}\leq 2$ or
$0<\tilde{q}_{2}(\beta,{\varepsilon^{\prime})}\leq 2$ at least two queries to
the corresponding oracle is used.

For $P_{2}$ there is actually no concern because the reasoning of Sec. V.3
that led to Eq. (3) is valid for any beta. This is because
$q_{2}=2\lceil\tilde{q}_{2}(\beta,\varepsilon^{\prime})\rceil$ is the exact
expression furnished by Lemma 37 in Ref. [14], used to prove Lemma 10, for
$\alpha F_{\beta}(\lambda)$. Therefore, although this formula is not strictly
tight, because it comes from successive approximations [14], it is exact and
valid for any $\beta>0$ in the sense that making exactly
$\tilde{q}_{2}(\beta,\varepsilon^{\prime})$ queries to the oracle ensures that
the error is below the tolerated.

$P_{1}$, on the other hand, requires a further analysis since, according to
Eq. (2), ${q}_{1}(\beta,\varepsilon^{\prime})$ is known to be equal to
$\tilde{q}_{1}(\beta,\varepsilon^{\prime})$ only in big-$\mathcal{O}$ notation
and, particularly, for large $\beta$. Next, we show that
$\tilde{q}_{1}(\beta,{\varepsilon^{\prime})}$ is, in fact, an over-estimation
of the actual $q_{1}(\beta,\varepsilon^{\prime})$ needed to guarantee error
$\varepsilon^{\prime}$.