# Calculating response functions of coupled oscillators using quantum phase
estimation

Sven Danz<EMAIL_ADDRESS>Institute for Quantum Computing Analytics
(PGI-12), Forschungszentrum Jülich, 52425 Jülich, Germany  Theoretical
Physics, Saarland University, 66123 Saarbrücken, Germany  Mario Berta
Institute for Quantum Information, RWTH Aachen University, 52074 Aachen,
Germany  Department of Computing, Imperial College London, London SW7 2RH, UK
Stefan Schröder  Fraunhofer Institute for Production Technology IPT, 52074
Aachen, Germany  Pascal Kienast  Fraunhofer Institute for Production
Technology IPT, 52074 Aachen, Germany  Frank K. Wilhelm  Institute for Quantum
Computing Analytics (PGI-12), Forschungszentrum Jülich, 52425 Jülich, Germany
Theoretical Physics, Saarland University, 66123 Saarbrücken, Germany
Alessandro Ciani  Institute for Quantum Computing Analytics (PGI-12),
Forschungszentrum Jülich, 52425 Jülich, Germany

###### Abstract

We study the problem of estimating frequency response functions of systems of
coupled, classical harmonic oscillators using a quantum computer. The
functional form of these response functions can be mapped to a corresponding
eigenproblem of a Hermitian matrix $H$, thus suggesting the use of quantum
phase estimation. Our proposed quantum algorithm operates in the standard
$s$-sparse, oracle-based query access model. For a network of $N$ oscillators
with maximum norm $\lVert H\rVert_{\mathrm{max}}$, and when the eigenvalue
tolerance $\varepsilon$ is much smaller than the minimum eigenvalue gap, we
use $\mathcal{O}(\log(Ns\lVert H\rVert_{\mathrm{max}}/\varepsilon)$
algorithmic qubits and obtain a rigorous worst-case query complexity upper
bound $\mathcal{O}(s\lVert H\rVert_{\mathrm{max}}/(\delta^{2}\varepsilon))$ up
to logarithmic factors, where $\delta$ denotes the desired precision on the
coefficients appearing in the response functions. Crucially, our proposal does
not suffer from the infamous state preparation bottleneck and can as such
potentially achieve large quantum speedups compared to relevant classical
methods. As a proof-of-principle of exponential quantum speedup, we show that
a simple adaptation of our algorithm solves the random glued-trees problem in
polynomial time. We discuss practical limitations as well as potential
improvements for quantifying finite size, end-to-end complexities for
application to relevant instances.

##  I Introduction

It is part of the experience of scientists and engineers that many relevant
problems appearing in nature can be formulated as the problem of finding the
eigenvalues and eigenvectors of a matrix. Consequently, a natural question to
ask is whether quantum computers can provide some advantages in solving
eigenproblems, compared to standard classical algorithms. The answer to this
question is the subject of active research, but it is highly nontrivial and
depends on the fine details of the problem itself.

In this paper, we focus on a particular category of eigenproblems, namely
those associated with a system of $N$ coupled harmonic oscillators. Needless
to say, systems of coupled harmonic oscillators are ubiquitous in nature, as
they can be used to effectively model the low energy dynamics of apparently
very different physical systems, from electrical circuits to the vibration of
molecules. In essence, our goal is to perform a modal analysis of the system,
focusing on the calculation of response functions in frequency (or Laplace)
domain. Response functions describe how much an external perturbation applied
to a certain oscillator, influences the dynamics of another oscillator in the
system. These functions have an analytical functional form that is related to
an eigenproblem defined on the system of coupled oscillators. We thus focus on
the question of how these response functions can be extracted using a quantum
algorithm and in which cases we might expect significant quantum speedups.

###  I.1 Related work

The idea of using a quantum computer for eigenproblems dates back to the early
days of quantum algorithm development [1]. In the quantum chemistry setting,
where one is typically interested in characterizing the low-lying eigenvalues
and eigenstates of a molecular Hamiltonian, it is in fact considered one of
the most promising applications of quantum computers [2, 3, 4].

Broadly speaking, one can distinguish between two categories of quantum
algorithms for eigenproblems, namely algorithms based on the variational
quantum eigensolver [5, 6], which lack rigorous scaling guarantees, and those
based on the quantum phase estimation (QPE) subroutine [7, 8, 9]. In this
work, we focus on the latter, and in particular on a formulation of standard
QPE for eigenvalue estimation that relies on the block-encoding of the matrix
of interest in a (larger) unitary [10, 11, 12, 13, 14, 4, 15, 16]. A key
feature of the problem that we exploit is that to estimate the response
functions, we do not need to invoke a state preparation subroutine that
prepares a desired eigenvector as a quantum state, for instance via controlled
rotations [17], reject sampling method [18], eigenvalue transform [19] or
adiabatic state preparation (see discussion in [4]). These kinds of state
preparation subroutines hide additional complexity and can ruin potential
quantum speedups when properly taken into account. Instead, as we will show,
in our case the omnipresent state preparation problem in QPE reduces to the
preparation of a product state, from which we need to determine the relevant
eigenvalues and associated weights via importance sampling. The simultaneous
estimation of multiple eigenvalues using either QPE or other quantum
subroutines has been studied in the literature using different approaches [20,
21, 22, 23].

The works that are most related to ours on resolving coupled harmonic
oscillators on a quantum computers are Ref. [24] and Ref. [25]. Ref. [24]
analyzes the problem of obtaining the classical time-evolution of a system of
coupled harmonic oscillators. To do this, the authors employ a particular
encoding of the positions and momenta of the oscillators in the amplitudes of
an $n$-qubit quantum system. On this system, a Hamiltonian can be defined that
mimics the time-evolution of the system of coupled oscillators at any time.
With this approach, while one cannot efficiently access the full solution, it
is possible to estimate efficiently relevant quantities such as the kinetic
energy of a subset of oscillators. Moreover, Ref. [24] also shows that the
random glued-trees problem originally studied in Ref. [26] can be mapped to a
problem on coupled oscillators and solved efficiently using the quantum
approach. We complement this result by showing that the glued-trees problem
can also be solved efficiently with our QPE approach by a simple adaptation of
the algorithm for the estimation of response functions. An interesting open
question is whether response functions can also be estimated using the time-
evolution approach of Ref. [24].

Ref. [25] instead focuses on a normal mode analysis of systems of coupled
oscillators, in the same spirit as our work, and essentially using the same
kind of embedding of the problem in quantum systems. However, Ref. [25] does
not discuss the application to the estimation of response functions and
assumes that a good approximation of the eigenvectors can be prepared in a
quantum system. While this would also allow to estimate response functions in
principle, we show that it is not necessary for this task.

Finally, our work is also connected to quantum approaches to solve the wave
equation [27] or more generally for finite element simulations of solid
structures [28, 29]. In fact, as we detail in the next section, once
discretized, these problems can be mapped to a system of coupled harmonic
oscillators.

###  I.2 Background on applications

The solution of ordinary and partial differential or integral equations
describing technical and physical processes, such as the deformation of
components under given loads or heat conduction in solids, is rarely
successful analytically and must be solved numerically. The necessary transfer
to a finite-dimensional equivalent problem requires a discretization process
using methods such as the finite element or finite volume method [30]. In our
work, we are particularly interested in problems for which this discretization
step leads to an effective model of coupled harmonic oscillators. A prominent
example of these problems are those that emerge in the manufacturing industry
as we detail below.

Due to their time efficiency, finite element models are often used in practice
to simulate machining processes, such as milling. Milling is a metal cutting
manufacturing process that uses the circular cutting motion of a tool, usually
with multiple cutting edges, to produce a vast variety of surfaces on a
workpiece (see Fig. 1).

Figure 1:  In-process-workpiece (IPW) of a single blade demonstrator (left)
and the detailed representation of the finite element mesh. The cutting tool
is shown in blue. With the cutter location-dependent IPW, material properties
(density, Young’s modulus, Poisson’s ratio), and boundary conditions as input,
the finite element modal analysis can be conducted and the frequency response
function (FRF) calculated by solving the corresponding eigenproblem. Typical
solvers used for this task are based on the QR algorithm and have complexity
$\mathcal{O}(N^{3})$ [31, 32, 33].

In all milling processes, in contrast to other processes, such as turning and
drilling for instance, the cutting edges are not constantly engaged. Still, at
least one cutting interruption per cutting edge occurs with each revolution of
the tool [34]. Due to the constant cutting interruptions, depending on the
milling cutter speed, a dynamic excitation of the workpiece and the tool can
occur, which can have a negative effect on the surface quality in the form of
vibration marks. Dynamic process stability simulations are carried out to
analyze the vibrations and to improve the process design for the milling of
thin-walled components [35, 36, 37, 38, 39]. The ultimate goal of these
simulations is to obtain response functions over a certain frequency range,
i.e., the typical frequency range of the cutting tool, that provide
information about the effect that the cutting tool at a certain location has
on the workpiece at other, possibly different, locations. For this particular
case, the response function is a so-called compliance (with units
$\mathrm{m}/\mathrm{N}$) as a function of the frequency of the cutting edge,
as shown in Fig. 1, which effectively quantifies the vibration of the
workpiece caused by the cutting tool.

Ref. [40] analyzed a workflow for milling dynamics simulation, examining
computational time and identifying numerical problems that can be suited for
quantum algorithms. The milling dynamics workflow considers the dynamic
behavior of the workpiece, which continuously varies due to changes in the
stiffness and mass of the workpiece caused by the material removal and the
position-dependent force excitation by the milling tool [41]. This analysis
highlighted the fact that finite element modal analysis is typically a
computationally intensive simulation application in milling of thin-walled
aerospace components, and as such could benefit from possible advantages that
quantum algorithms could offer for this task.

###  I.3 Structure

The rest of the paper is structured as follows. In Sec. II, we discuss the
connection between the computation of the response function of coupled
oscillators and the solution of an eigenproblem. We further describe briefly
our quantum algorithm and its scaling behavior, which is summarized in 1. We
finish the section with a critical assessment of the limitations and possible
improvements of our quantum approach. In Sec. III, we describe our proposed
quantum algorithm in more detail, for which a pseudocode is provided in
Algorithm 1. Sec. IV discusses how the random-glued trees problem can be
solved using the QPE subroutine. We then summarize our results and provide an
outlook in Sec. V. The appendices contain a more in-depth derivations of the
analytical form of the response function (Appendix A), descriptions of the
quantum subroutines used in the algorithm (Appendix B), a detailed error
analysis (Appendix C), a modified version of the algorithm for non-local
response functions (Appendix D), and a scaling analysis of the algorithm
applied to the random glued-trees problem (Appendix E).

##  II Problem statement and main results

###  II.1 Linear response functions as an eigenproblem

Figure 2:  A representation of coupled oscillators. Point-like masses $m_{u}$
in blue are connected by springs with spring constant $\kappa_{wu}$. The force
$f_{v}$ applied to mass $v$ yields a displacement $x_{u}$ at mass $u$.

Let us consider a system of masses coupled by spring constants as shown in
Fig. 2 and let $G=(V,E)$ be a graph with $V$ the set of vertices and $E$ the
set of edges, that describes the network of oscillators. With each vertex
$u\in V$ we associate a harmonic oscillator with mass $m_{u}>0$ and a coupling
to a common wall with spring constant $\kappa_{u}>0$. We denote by $N=\lvert
V\rvert$ the number of oscillators. With each edge $(u,v)\in E$ in the graph,
instead, we associate a coupling spring constant $\kappa_{uv}>0$. In what
follows, we implicitly assume an ordering of the vertices, i.e., we associate
with each $u\in V$ a number between $1$ and $N$. For notational simplicity, we
also denote this number by $u$, since its meaning is clear from the context.
Additionally, we denote by $N_{G}[u]$ the closed neighborhood of $u\in V$,
i.e., the set of vertices connected to $u$ by an edge plus $u$ itself. We call
$s_{u}=|N_{G}[u]|$ and define the sparsity $s$ as $s=\mathrm{max}_{u\in
V}s_{u}$. For later convenience, we also define the set of vertices
$\tilde{N}_{G}[u]$ such that $N_{G}[u]\subset\tilde{N}_{G}[u]$ and
$|\tilde{N}_{G}[u]|=s$. Basically, $\tilde{N}_{G}[u]$ is the neighborhood of
$u$ enlarged with dummy vertices to ensure that the total number of elements
in the set is $s$.

####  II.1.1 Unperturbed dynamics

We can compactly write the dynamical equations that govern the motion of the
coupled oscillator systems, by defining two $N\times N$ matrices, namely a
diagonal mass matrix $\bm{M}$ with matrix elements

$M_{uv}=m_{u}\delta_{uv},$ (1)

and a stiffness matrix $\bm{K}$ with matrix elements

$K_{uv}=\begin{cases}\kappa_{u}+\sum_{w\in
N_{G}[u]\setminus\\{u\\}}\kappa_{wu}\quad&\mathrm{if}\quad u=v,\\\
-\kappa_{uv}\quad&\mathrm{if}\quad(u,v)\in E,\\\
0\quad&\mathrm{otherwise}.\end{cases}$ (2)

Additionally, we denote by $\bm{x}$ the $N$-dimensional column vector whose
elements are the oscillator positions $x_{u}$ with respect to the equilibrium
state.

The system evolves according to the dynamical equation

$\bm{M}\frac{d^{2}\bm{x}}{dt^{2}}=-\bm{K}\bm{x}.$ (3)

In order to transform this problem into an eigenproblem, let us first define
the matrix

$\bm{H}=\bm{M}^{-1/2}\bm{K}\bm{M}^{-1/2},$ (4)

which, as we will see, will play the role of a Hamiltonian in our quantum
algorithm. Note that since $\bm{M}$ is positive definite, also $\sqrt{\bm{M}}$
is positive definite and thus invertible. The matrix $\bm{H}$ is Hermitian and
positive-semidefinite and we denote its eigenvalues by $\lambda_{j}\geq 0$ for
$j\in\\{1,\dots,N\\}$. We now show that solving the dynamical equation Eq. 3
is equivalent to diagonalizing the matrix $\bm{H}$. We denote by $\bm{W}$ an
orthogonal matrix that diagonalizes $\bm{H}$,

$\bm{\Lambda}=\bm{W}^{T}\bm{H}\bm{W},$ (5)

where $\bm{\Lambda}$ is a diagonal matrix with the eigenvalues of $\bm{H}$ on
the diagonal, i.e., $\Lambda_{jj}=\lambda_{j}$. Defining the vector of normal
mode variables $\bm{y}$ as

$\bm{y}=\bm{W}^{T}\bm{M}^{1/2}\bm{x},$ (6)

we get an uncoupled system of harmonic oscillators satisfying

$\frac{d^{2}\bm{y}}{dt^{2}}=-\bm{\Lambda}\bm{y}\implies\frac{d^{2}y_{j}}{dt^{2}}=-\lambda_{j}y_{j}.$
(7)

Since $\lambda_{j}\geq 0$ they are usually denoted as
$\lambda_{j}=\omega_{j}^{2}$, with $\omega_{j}$ the resonance frequency of the
normal mode $j$. Note that Eq. 6 implies that

$x_{u}=\frac{1}{\sqrt{m_{u}}}\sum_{j=1}^{N}W_{uj}y_{j}.$ (8)

####  II.1.2 Adding an external force

Let us now add an external force to Eq. (3). In particular, we assume that a
time-dependent force $f_{v}(t)$ is applied to the oscillator $v\in V$ and we
gather all the forces in an $N$-dimensional column vector $\bm{f}(t)$.
Including the forcing term, the dynamical equation reads

$\bm{M}\frac{d^{2}\bm{x}}{dt^{2}}=-\bm{K}\bm{x}+\bm{f}.$ (9)

Eq. (9) gives rise to a convolution relation between $\bm{f}(t)$ and
$\bm{x}(t)$

$\bm{x}(t)=\int_{-\infty}^{+\infty}d\tau\bm{g}(t-\tau)\bm{f}(\tau),$ (10)

where $\bm{g}(t)$ is an $N\times N$ matrix that we call the matrix response
function. The matrix elements $g_{uv}(t)$ can simply be interpreted as the
response of the oscillator $u$ when a forcing term on the oscillator $v$,
$f_{v}(t)$, is a Dirac delta function, i.e., $f_{v}(t)=\delta(t)$, while no
force is applied to the other oscillators 111The response function is also
called the Green’s function for the system of coupled oscillators.. This
interpretation also suggests that causality implies $\bm{g}(t)=0$ if $t<0$.
Thus, we can simply take the upper limit of the integral in Eq. 10 to be $t$
rather than $+\infty$. Eq. 10 is the defining equation of a linear time-
invariant system [43].

Now, let us introduce the bilateral Laplace transform of a generic function
$\bm{z}(t)$ as

$\bm{Z}(s)=\int_{-\infty}^{+\infty}dte^{-st}\bm{z}(t),\quad
s=\sigma+i\omega\in\mathbb{C}.$ (11)

The Laplace transform is defined only in the region of convergence, that is
the complex domain where the integral in Eq. (11) converges. From the
convolution theorem for the Laplace transform Eq. (10) in Laplace domain
becomes

$\bm{X}(s)=\bm{G}(s)\bm{F}(s),$ (12)

with $\bm{G}(s)$ the Laplace transform of the response function in time domain
$\bm{g}(t)$, while $\bm{X}(s)$ and $\bm{F}(s)$ the Laplace transform of the
position vector $\bm{x}(t)$ and the vector of forces $\bm{f}(t)$,
respectively.

In the main text, we focus on a quantum algorithm to determine the diagonal
elements of $\bm{G}(s)$ that take the form

$G_{uu}(s)=\frac{1}{m_{u}}\sum_{j=1}^{N}\frac{W_{uj}^{2}}{\lambda_{j}+s^{2}},$
(13)

which we call the local response functions. However, it is also possible to
obtain the off-diagonal elements of $\bm{G}(s)$ or non-local response function
at the price of adding a Hadamard test, as we show in Appendix D. We provide a
derivation of Eq. 13, as well as the more general definition for local and
non-local response, in Appendix A.

Eq. 13 shows that the response function $G_{uu}(s)$ is completely determined
once we solve the eigenproblem associated with $\bm{H}$. In fact, all we need
to compute the response function is the eigenvalues $\lambda_{j}$ and the
coefficients $W_{uj}^{2}$, with $W_{uj}$ the matrix elements of the orthogonal
matrix $\bm{W}$ that diagonalizes $\bm{H}$. As discussed in Sec. I.2, the
response function is the relevant quantity to study in practical applications.
State-of-the-art classical algorithms that are used in practice for the
eigenproblem at hand, such as the QR algorithm, need at least
$\mathcal{O}(N^{3})$ arithmetic operations. In fact, most of the fastest
classical algorithms for eigenvalue computation rely on the Householder
transformation which transforms a matrix into the Hessenberg form (almost
triangular). This transformation needs $4N^{3}/3$ operations (see Chap. 11 p.
474 in [32]).

In the following quantum treatment, we refer to the $\bm{H}$ matrix in Eq. 4
as the Hamiltonian. From now on we avoid using the bold symbols, simply
denoting the Hamiltonian as $H$, and do the same for the other matrices.

###  II.2 Main idea

In this work, we describe a quantum algorithm that computes the local and non-
local response functions of a system of coupled oscillators. We start by
providing the basic intuition behind the algorithm for the case of the local
response function in Eq. 13, while we refer to Sec. III and the appendices for
more technical details. In order to compute $G_{uu}(s)$, we need to determine
the eigenvalues $\lambda_{j}$ and the coefficients $W_{uj}^{2}$. In our
approach, each vertex $u\in V$ is associated with a bitstring $\ket{u}$ of
length $n$ in the computational basis, where the number of qubits is
$n=\lceil\log_{2}N\rceil$. For simplicity, from now on we assume that $N$ is a
power of two, i.e., $N=2^{n}$. Since, $W$ is the orthogonal transformation
that diagonalizes $H$, we can express the bitstring $\ket{u}$ as a linear
combination of the eigenstates $\ket{\lambda_{j}}$ of $H$:

$\ket{u}=\sum_{j=1}^{N}W_{uj}\ket{\lambda_{j}}.$ (14)

Thus, we can interpret $W_{uj}$ as weights in the previous decomposition. Eqs.
(13) and (14) convey the core idea behind the algorithm. If we prepare the
bitstring $\ket{u}$ and perform measurement in the eigenbasis
$\ket{\lambda_{j}}$ of $H$ we will get the eigenvalues $\lambda_{j}$ with
probability $W_{uj}^{2}$. By repeating this procedure many times we can
estimate the coefficients $W_{uj}^{2}$ via statistical averages and obtain an
estimate for the response function $G_{uu}(s)$. Note that there is an implicit
importance sampling in this procedure: the eigenstates that contribute the
most to $G_{uu}(s)$ are also those whose eigenvalue will be measured more
often. In order to implement the measurement in the eigenbasis we use QPE. The
full algorithm is described in detail in Sec. III.

For the implementation of the quantum algorithm, we assume that we have oracle
access to the matrix elements of $H$ and that $H$ is $s$-sparse, i.e., that
each row has at most $s$ non-zero elements. In particular, we assume that we
have access to two matrix oracles. The first one is the position oracle
$O_{P}$, that acts on the $2n$-qubit state $\ket{u}\ket{0}^{\otimes n}$ as

$O_{P}\ket{u}\ket{0}^{\otimes
n}=\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}\ket{u}\ket{v}.$ (15)

$O_{P}$ stores in the second register a linear superposition of all the column
indices that correspond to non-zero matrix elements at row $u$ and the
additional dummy indices in $\tilde{N}_{G}[u]$ if necessary.

The second oracle is the angle oracle $O_{\vartheta}$ which contains
information about the matrix elements of $H$. It encodes the angle

$\vartheta_{uv}=\arccos\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}},$
(16)

and the sign of $H_{uv}$ in an $(r+1)$-qubit ancilla register

$O_{\vartheta}\ket{u,v}\ket{0}^{\otimes
r+1}=\ket{u,v}\ket{\vartheta_{uv},\Theta(-H_{uv})},$ (17)

with $\Theta(\cdot)$ the Heaviside step function and $r$ the number of bits
used to represent $\vartheta_{uv}$. We renormalize $H_{uv}$ in Eq. 16 with
$\norm{H}_{\mathrm{max}}=\max_{u,v}\absolutevalue{H_{uv}}$. Specifically,
since $0\leq\vartheta_{uv}\leq\pi/2$, we assume the following $r$-bit
representation of $\vartheta_{uv}$

$\vartheta_{uv}=\frac{\pi}{2}\sum_{l=1}^{r}\vartheta_{uv}^{(l)}2^{-l},$ (18)

with $\vartheta_{uv}^{(l)}\in\\{0,1\\}$. Thus, when we write
$\ket{\vartheta_{uv}}$ in an $r$-qubit register, this means

$\ket{\vartheta_{uv}}=\ket{\vartheta_{uv}^{(1)}}\ket{\vartheta_{uv}^{(2)}}\dots\ket{\vartheta_{uv}^{(r)}}.$
(19)

Finally, we assume that we have a quantum adder $\mathrm{ADD}$ available,
which adds the information in one quantum register to another, defined as

$\text{ADD}\ket{a}\ket{b}=\ket{a}\ket{b+a\mod 2^{n}},$ (20)

where $n$ is the number of qubits in the second register 222The number of
qubits in the first register does not need to be equal to the one in the
second register.. The inverse of this operation is a quantum subtractor
$\text{ADD}^{\dagger}=\text{SUB}$. Circuit realizations of the quantum adder
are well known in the literature [45, 46, 47, 48]. Together with the matrix
oracles $O_{P}$ and $O_{\vartheta}$, the quantum adder $\mathrm{ADD}$ is
needed for the implementation of the block-encoding described in detail in
Appendix B (see Fig. 5 in particular) and used in Sec. III.

The computational cost of the quantum adder will be neglected in the following
query complexity analysis. This is because the implementation of the oracle
$O_{\vartheta}$ usually requires multiple queries to the quantum adder. Hence,
the additional use of one ADD per oracle is negligible in terms of required
qubit number and computation time.

Table 1: Parameters used in the complexity analysis. Symbol | Meaning
---|---
$N$ | Number of rows or columns of $H$
$s$ | Sparsity of $H$
$\norm{H}_{\mathrm{max}}$ | Maximum absolute value of the entries of $H$
$\Delta_{\lambda}^{(u)}$ | Minimum eigenvalue gap at oscillator $u$
$\varepsilon$ | Additive eigenvalue tolerance
$\delta$ | Additive weight tolerance
$\zeta$ | Probability of failure
$N_{u}$ | Number of nonzero $W_{uj}$ at oscillator $u$

###  II.3 Complexity analysis

We now summarize the main results in terms of scaling with respect to relevant
parameters. We start by defining the additional parameters that enter in the
scaling analysis. First, we define the minimum, positive gap in the eigenvalue
spectrum as

$\Delta_{\lambda}=\mathrm{min}\\{|\lambda_{i}-\lambda_{j}|:i,j\in\\{1,\dots,N\\},\,|\lambda_{i}-\lambda_{j}|>0\\}.$
(21)

Additionally, we define the minimum, positive eigenvalue gap at oscillator $u$
as

$\Delta_{\lambda}^{(u)}=\mathrm{min}\\{|\lambda_{i}-\lambda_{j}|:i,j\in\\{1,\dots,N\\},\\\
\,|\lambda_{i}-\lambda_{j}|>0,W_{ui}^{2}>0,W_{uj}^{2}>0\\},$ (22)

with $W_{ui}$ and $W_{uj}$ being the coefficients in Eq. 14. This is the
parameter that matters in our algorithm since, in the presence of
degeneracies, weights associated with the same eigenvalue would just be
gathered together in Eq. 13 and eigenvectors that have no support at
oscillator $u$ give zero contribution. Obviously,
$\Delta_{\lambda}^{(u)}\geq\Delta_{\lambda}$. The scaling further depends on
the required, additive tolerances $\varepsilon$ and $\delta$ for the
eigenvalues $\lambda_{j}$ and weights $W_{uj}^{2}$, respectively, which we
treat as independent parameters of our choice 333The smaller the error on the
eigenvalues, the smaller the error on the weights, and vice versa. However,
for our analysis it is easier to think about them as independent tolerances..
Moreover, we define the parameter $N_{u}$ to be the number of eigenvectors
that have support at oscillator $u$, that is for which $W_{uj}^{2}>0$. Note
that in the worst case scenario $N_{u}=N$ from the true value. Finally, we
denote by $\zeta$ the probability that the estimate of $W_{uj}^{2}$ deviates
by more than $\delta$. All the parameters used in the complexity analysis are
summarized in Table 1.

Our main result is described by the following theorem.

######  Theorem 1.

Let $H$ be a $N\times N$ Hermitian matrix that describes a system of coupled
oscillators as in Eq. 4 and let $G=(V,E)$ be the graph associated with the
network of oscillators. Let $H$ be $s$-sparse with maximum norm
$\norm{H}_{\mathrm{max}}$. Given an oscillator $u\in V$, there exists a
quantum algorithm that approximates the eigenvalues $\lambda_{j}$ of $H$ and
the weights $W_{uj}^{2}$ in Eq. 13, with additive error $\varepsilon$ and
$\delta$, respectively, that requires

$n_{\mathrm{tot}}=\mathcal{O}\left(\log_{2}\left(Ns\norm{H}_{\text{max}}\max\left(\frac{1}{\varepsilon},\frac{1}{\delta\Delta_{\lambda}^{(u)}}\right)\right)\right).$
(23)

qubits and

$N_{\mathrm{queries}}^{(\mathrm{tot})}=\mathcal{O}\Biggl{(}\frac{s\norm{H}_{\mathrm{max}}}{\delta^{2}}\ln\left(\frac{N_{u}}{\zeta}\right)\max\left(\frac{1}{\varepsilon},\frac{1}{\delta\Delta_{\lambda}^{(u)}}\right)\Biggl{)},$
(24)

total queries of the matrix oracles $O_{P}$ and $O_{\vartheta}$ in Eq. 15, Eq.
17, respectively.

A detailed error analysis of the QPE-based algorithm behind 1 is presented in
presented in Sec. III (see Appendix C for more details). In what follows, we
provide further comments on the complexity of our result.

First of all, we note that Eq. 23 implies that $\varepsilon$ should be taken
at least $\varepsilon<\Delta_{\lambda}^{(u)}$, since $\delta<1$. Moreover, in
Eq. 23, we are not considering the ancilla qubits necessary to implement the
oracles $O_{P}$, $O_{\vartheta}$, since these are problem-dependent. If we use
$r$ bits to represent the matrix elements, the required number of bits in
classical algorithms is usually $\mathcal{O}(rN^{2})$, while for $s$-sparse
matrices it can be $\mathcal{O}(rsN)$ in the best case. The gate complexity of
our quantum algorithm depends on the complexity of implementing the oracles
$O_{P}$ and $O_{\vartheta}$, which are, as mentioned above, problem-dependent.
Generally, as long as the matrix elements of a matrix can be accessed
efficiently classically, this should also hold on a quantum computer, although
at the price of having several ancillas to make the computation reversible. As
an example, for the case of a linear chain of oscillators, we provide a full
compilation of the oracles into elementary gates in the code available in the
repository associated with this manuscript (see the “Data Availability”
section). As such, we choose to present the complexity of our algorithm in
terms of the number of oracle queries $N_{\mathrm{queries}}$. The query
complexity of a single run of QPE is (see Sec. C.1.3)

$N_{\mathrm{queries}}=\mathcal{O}\left(s\norm{H}_{\text{max}}\max\left(\frac{1}{\varepsilon},\frac{1}{\delta\Delta_{\lambda}^{(u)}}\right)\right).$
(25)

This is mostly due to the repetition of the controlled-$V$ operator in the
core of the phase estimation (cf. Fig. 3 and Algorithm 1).

QPE is a probabilistic algorithm from which we sample multiple times to
estimate the eigenvalues and the associated probabilities. As we show in Sec.
C.2, the number of samples $N_{S}$ needed to estimate each weight $W_{uj}^{2}$
via empirical averages with error $\delta$ and failure probability $\zeta$
scales as

$N_{S}=\mathcal{O}\left(\frac{1}{\delta^{2}}\ln\left(\frac{N_{u}}{\zeta}\right)\right).$
(26)

In the worst case scenario when $N_{u}=N=2^{n}$, the scaling would be
logarithmic with the size of the matrix $N$, i.e., linear with the number of
qubits $n$. We highlight the fact that this is not equivalent to requiring
that the probability distribution $W_{uj}^{2}$ is estimated with error
$\delta$ in total variation distance. Using the definition of total variation
distance and Hoeffding’s inequality similarly to Sec. C.2, one can show that
in this case the scaling would be linear, and not logarithmic, with $N_{u}$.
Moreover, we are not requiring to find all the $N_{u}$ eigenvalues that
contribute in the response function at oscillator $u$. This would require a
number of samples scaling at least as $N_{u}$, but what matters for us is the
evaluation of the weights $W_{uj}^{2}$ with error $\delta$, since the weights
appear at the numerator in Eq. 13. This leads to the inverse quadratic scaling
in $\delta$ in Eq. 26.

Combining Eq. 25 and Eq. 26, we obtain the total number of queries as in Eq.
24.

###  II.4 Critical assessment

####  II.4.1 Limitations of the algorithm

Focusing on the eigenvalue tolerance $\varepsilon$, the maximum number of
oracle queries required for the described algorithm is dominated by the
inverse linear scaling in the renormalized eigenvalue tolerance

$\frac{\varepsilon}{s\norm{H}_{\text{max}}}.$

Assuming the worst-case scenario in which,
$\Delta_{\lambda}^{(u)}=\Delta_{\lambda}$ the desired tolerance $\varepsilon$
is required to satisfy $\varepsilon<\Delta_{\lambda}$, which guarantees a
separation of all the eigenvalues. We further know that the $1$-norm
$\norm{H}_{1}=\max_{v}\sum_{u}\absolutevalue{H_{uv}}$ is smaller than
$s\norm{H}_{\mathrm{max}}$. Thus, we have

$\frac{\varepsilon}{s\norm{H}_{\mathrm{max}}}\leq\frac{\Delta_{\lambda}}{s\norm{H}_{\mathrm{max}}}\leq\frac{\Delta_{\lambda}}{\norm{H}_{\mathrm{1}}}\leq\frac{\Delta_{\lambda}}{\lambda_{\mathrm{max}}},$
(27)

where we used the upper limit for the maximum eigenvalue
$\lambda_{\mathrm{max}}\leq\norm{H}_{1}$. Hence, the number of oracle queries
scales with the ratio between the smallest eigenvalue gap and the maximum
eigenvalue. This, in the general case, increases polynomially with the number
of eigenvalues that have to fit between $0$ and $\lambda_{\mathrm{max}}$, and
accordingly it increases polynomially with the number of oscillators $N$.

The kinds of problems for which we can expect a polylogarithmic scaling in $N$
of the total query complexity necessarily need to satisfy
$\lambda_{\mathrm{max}}/\Delta_{\lambda}^{(u)}=\mathcal{O}(\mathrm{polylog}(N))$
in our rigorous worst case scaling analysis. On the one hand, the random
glued-trees problem we discuss in Sec. IV has this property, which ultimately
originates from its very large degeneracy. This shows that there exists well-
defined problems for which a large exponential speedup is possible. However,
the task there is not the determination of response functions. On the other
hand, if we simply consider a $\mathrm{1D}$ periodic chain with $N$
oscillators and unit spring constants and masses, the eigenvalues can be
obtained analytically and are given by [50]

$\lambda_{j}=\frac{2}{\sqrt{N}}\left(1-\cos\left(\frac{2\pi(j-1)}{N}\right)\right),$

for $j\in\\{1,\dots,N\\}$. For this problem $s\norm{H}_{\mathrm{max}}$ is
bounded, namely $s\norm{H}_{\mathrm{max}}=6$, but for large $N$ the gap for
this problem scales as $\Delta_{\lambda}=\mathcal{O}(N^{-5/2})$, yielding
eventually a polynomial scaling of the number of queries with $N$. Notice that
these considerations have no influence on the scaling with the tolerance on
the weights $\delta$, which is $1/\delta^{3}$ if we take it as an independent
parameter (see Eq. 24).

####  II.4.2 Possible improvements

Possible improvements in the scaling in the $\delta$ parameter could be
achieved combining our algorithm with the quantum amplitude estimation
subroutine [51]. Additionally, in the previous discussion, we have not
restricted the domain of interest of the eigenvalues that are in general
between $0$ and $\lambda_{\mathrm{max}}$. However, in practice we might be
interested in evaluating the response function only in a certain frequency
range for which only eigenvalues $\lambda_{j}\in[\lambda_{a},\lambda_{b}]$
matter. In fact, we might know in advance the typical frequency spectrum of
the forcing term, as it happens in the milling case discussed in Sec. I.2. In
this case, we can approximate the response function as

$G_{uu}(s)\approx\frac{1}{m_{u}}\sum_{j\in\mathcal{J}}\frac{W_{uj}^{2}}{\lambda_{j}+s^{2}},$
(28)

for $\mathrm{Im}(s)\in[\lambda_{a},\lambda_{b}]$ and
$\mathcal{J}=\\{j:\lambda_{j}\in[\lambda_{a},\lambda_{b}]\\}$. Thus, it might
be beneficial to apply a filter to the quantum state $\ket{u}$ to remove all
the contributions from eigenvalues that are outside the interval of interest
$[\lambda_{a},\lambda_{b}]$, and produce an initial, filtered state
$\ket{u_{f}}$ such that

$\ket{u_{f}}\sim\sum_{\lambda_{j}\in[\lambda_{a},\lambda_{b}]}W_{uj}\ket{\lambda_{j}}.$
(29)

Quantum eigenvalue filtering has been already studied in the literature [52,
13, 53] and could find an application in our problem. We note that also for
classical algorithms it is possible to restrict the eigenvalue search in a
specific range using filter diagonalization methods [54, 55, 56, 57].

Needless to say, our algorithm relies on the efficient implementation of the
matrix oracles. As discussed in Sec. II.3, this usually requires the
implementation of reversible, coherent arithmetic that comes at the cost of
additional ancilla qubits. Alternatively, one could adapt state preparation
techniques that avoid the need of coherent arithmetic [58, 19] to this task
444In fact, the oracles are ultimately used to prepare the state in Eq. B.5.
In case the matrix elements are already stored in a classical memory, these
data would need to be accessed by the oracle using a QRAM [60], which can lead
to severe overheads destroying any exponential speedup (see discussion in [4,
61]). For the milling use case discussed in Sec. I.2, this would mean that the
discretization of the problem would either need to happen in the quantum
computer or it must be given in a format that the quantum computer can access
efficiently. Rather than an intrinsic limitation, we believe that this is an
interesting scientific question and underexplored area of research.

Finally, we point out that our algorithm has some connection with the vector
fitting algorithm [62], where the goal is, given $K$ observations of
$G_{uu}(s)$ at $s=i\omega^{(k)}$, $k\in\\{1,\dots,K\\}$, to output the best
estimate of the weights $W_{uj}^{2}$, usually called residues, and the
eigenvalues $\lambda_{j}$. In our case, we want to achieve the same task, but
using the data obtained from running QPE multiple times to estimate the
weights and the eigenvalues.

##  III The quantum algorithm

In this section, we detail how the local response function $G_{uu}(s)$ in Eq.
13 can be obtained using QPE. In particular, we describe the algorithm, and
study its scaling, using the standard QPE [7, 9], that makes use of the
inverse quantum Fourier transform. A similar analysis can also be performed
using other versions of QPE that use a single ancilla qubit for the phase
register [63, 8, 64, 13] or coherent QPE with fewer qubits [65]. As we have
seen in Sec. II.1, in order to obtain the response function, we need to solve
the eigenproblem associated with the Hermitian matrix in Eq. 4.

[wire types= q,q,q,q,b,b,classical gap=0.07cm] [wires=4, label style=anchor=
east, xshift=-0cm]$\text{Phase:}\ket{0}^{\otimes m}$& H 5 … [4,disable auto
height]$\text{QFT}^{\dagger}$
H 4 …

H … 2
[label style=anchor= east, xshift = -0.2cm ]$\text{State:}\ket{0}^{\otimes n}$
∏_i=0^n-1X_i^u_i [2, disable auto height]V[2, disable auto height]V^2 … [2,
disable auto height]V^2^m-1
[label style=anchor= east, xshift=-0.2cm]$\text{Ancilla:}\ket{0}^{\otimes
a+l}$ …

Figure 3: QPE circuit for the computation of the local response function. The
state preparation is realized by a series of NOT-gates. QPE is performed using
the walk operator $V$ defined in Eq. 35 followed by the inverse quantum
Fourier transform ($\mathrm{QFT}^{\dagger}$) [7]. Sampling from the phase
register and repeating many times allows us to estimate the relevant
eigenvalues $\lambda_{j}$ and the weights $W_{uj}^{2}$ in Eq. 13.

###  III.1 Matrix encoding

The algorithm is based on the use of QPE as eigenvalue estimator [1]. In our
analysis, we focus on a version of QPE inspired by quantum walks [66, 67] and
similar to the one of Refs. [15, 16]. The basic idea of this method is that if
we want to estimate the eigenvalues $\lambda$ of a Hamiltonian $H$, we could
as well estimate a function $f(\lambda)$ as long as the function is
invertible. In particular, in our case instead of performing QPE with unitary
$e^{iH}$, we perform it using a unitary $V$ with eigenvalues $e^{\pm
i\arccos(\alpha\lambda)}$ with $\alpha$ a normalization factor. The appeal of
this method is that $V$ can be implemented exactly using the standard oracles,
which simplifies the error analysis. The full QPE is shown in Fig. 3. In what
follows, we briefly describe the method in the language of block-encoding and
qubitization [10, 12, 13, 14, 68, 4].

Let us assume that we have a unitary $U_{H}$ acting on $a+n$ qubits that
realizes an $(\alpha,a,0)$-block-encoding of our $n$-qubit Hamiltonian $H$.
Mathematically, this translates into the requirement

$\bra{0}^{\otimes a}U_{H}\ket{0}^{\otimes a}=\alpha H,$ (30)

where $a$ is the number of additional ancilla qubits used for the block-
encoding. The eigenstates of $H$ with eigenvalue $\lambda$ will be denoted by
$\ket{\lambda}$ in braket notation. Without loss of generality, we assume that
$U_{H}$ is Hermitian, i.e., $U_{H}=U_{H}^{\dagger}$, which implies
$U_{H}^{2}=I_{a+n}$ 555If $U_{H}$ is not Hermitian, one can add one ancilla
qubit and define the Hermitian unitary
$U_{H}^{(\mathrm{new})}=\outerproduct{+}{-}\otimes
U_{H}+\outerproduct{-}{+}\otimes U_{H}^{\dagger}$, which is a
$(\alpha,a+1,0)$-block-encoding of $H$.. We define the operator $\Pi$ as the
projector onto the subspace with all ancillas in the zero state

$\Pi=\outerproduct{0}{0}^{\otimes a}\otimes I_{n},$ (31)

where $I_{n}$ denotes the $n$-qubit identity.

The crucial property that qubitization uses is that the action of a Hermitian
block-encoding $U_{H}$ on the state $\ket{0}^{\otimes a}\ket{\lambda}$ can be
written as (see Section 10.4 in Ref. [4])

$U_{H}\ket{0}^{\otimes a}\ket{\lambda}=\alpha\lambda\ket{0}^{\otimes
a}\ket{\lambda}+\sqrt{1-\alpha^{2}\lambda^{2}}\ket{0\lambda^{\perp}},$ (32)

where $\ket{0\lambda^{\perp}}$ is a state such that $(\bra{0}^{\otimes
a}\bra{\psi})\ket{0\lambda^{\perp}}=0$ $\forall$ $n$-qubit state $\ket{\psi}$
and
$\innerproduct{0\lambda^{\prime\perp}}{0\lambda^{\perp}}=\delta_{\lambda\lambda^{\prime}}$.
Additionally,

$U_{H}\ket{0\lambda^{\perp}}=\sqrt{1-\alpha^{2}\lambda^{2}}\ket{0}^{\otimes
a}\ket{\lambda}-\alpha\lambda\ket{0\lambda^{\perp}}.$ (33)

Eq. 32 and Eq. 33 imply that $\forall\lambda$ we can identify a qubit-like
subspace $\mathcal{H}^{(\lambda)}=\mathrm{span}\\{\ket{0}^{\otimes
a}\ket{\lambda},\ket{0\lambda^{\perp}}\\}$, so that the action of $U_{H}$ onto
this subspace can be represented by a $2\times 2$ matrix $U_{H}^{(\lambda)}$
given by

$U_{H}^{(\lambda)}=\begin{pmatrix}\alpha\lambda&\sqrt{1-\alpha^{2}\lambda^{2}}\\\
\sqrt{1-\alpha^{2}\lambda^{2}}&-\alpha\lambda\end{pmatrix}.$ (34)

However, the eigenvalues of this matrix are $\pm 1$, since $U_{H}$ is a
reflection and thus, they do not carry information about the eigenvalues
$\lambda$ of $H$. To remedy this issue we define the “walk” operator $V$ that
we use in the QPE subroutine depicted in Fig. 3 as

$V=U_{H}(2\Pi-I_{a+n}).$ (35)

In fact, the operator $2\Pi-I_{a+n}$ also preserves the subspaces
$\mathcal{H}^{(\lambda)}$ and its action on $\mathcal{H}^{(\lambda)}$ is
simply represented by a Pauli $Z$ matrix:

$(2\Pi-I_{a+n})^{(\lambda)}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}.$ (36)

Consequently, the action of $V$ on $\mathcal{H}^{(\lambda)}$ is given by:

$V^{(\lambda)}=U_{H}^{(\lambda)}(2\Pi-I_{a+n})^{(\lambda)}=\\\
\begin{pmatrix}\alpha\lambda&-\sqrt{1-\alpha^{2}\lambda^{2}}\\\
\sqrt{1-\alpha^{2}\lambda^{2}}&\alpha\lambda\end{pmatrix}.$ (37)

The matrix $V^{(\lambda)}$ has eigenvalues

$\mu_{\pm}=\alpha\lambda\pm i\sqrt{1-\alpha^{2}\lambda^{2}}=e^{\pm
i\arccos(\alpha\lambda)},$ (38)

with corresponding eigenvectors

$\ket{\mu_{\pm}}=\frac{1}{\sqrt{2}}(\ket{0}^{\otimes a}\ket{\lambda}\mp
i\ket{0\lambda^{\perp}}).$ (39)

Thus, if we start from a state
$\ket{\psi}=\sum_{\lambda}c_{\lambda}\ket{0}^{\otimes a}\ket{\lambda}$ and
perform “ideal” phase estimation using the operator $V$, the system will
collapse in either $\ket{\mu_{+}}$ or $\ket{\mu_{-}}$ and the value
$+\arccos(\alpha\lambda)$ or $-\arccos(\alpha\lambda)$, respectively, can be
read out in the phase register. In either case, a subsequent measurement of
the ancilla register in the computational basis would give all zeros with
probability $1/2$ and thus, prepare the eigenstate $\ket{0}^{\otimes
a}\ket{\lambda}$.

The procedure we have just described works for an arbitrary Hermitian block-
encoding of $U_{H}$. In Appendix B, we show for completeness how such a block-
encoding can be realized if we have access to the oracles in Eq. 15 and Eq.
17. In this case, $a=n+2$ and the parameter $\alpha$ is given by

$\alpha=\frac{1}{s\lVert H\rVert_{\mathrm{max}}},$ (40)

where $s$ denotes the sparsity of the matrix $H$.

###  III.2 Obtaining the response function

In this section, we describe the full quantum algorithm to obtain the diagonal
response function in Eq. 13. As we see from Fig. 3, we have three main
registers:

  1. 1.

An $m$-qubit register that, at the end of the circuit, stores the information
about the eigenphase;

  2. 2.

An $n$-qubit register that represents the Hilbert space associated with the
Hamiltonian $H$;

  3. 3.

An $(a+l)$-qubit ancilla register, which includes the $a$ ancillas for the
block-encoding described in Sec. III.1 as well as $l$ ancilla qubits needed
for the implementation of the oracles (see Appendix B).

We start by preparing a state $\ket{u}$ associated with one of the vertices
$u\in V$ in the state register (cf. Fig. 3) by using a product of single-qubit
Pauli-$X$ gates

$\prod_{i=0}^{n-1}X^{u_{i}}_{i}\ket{0}^{\otimes
n}=\ket{u}=\sum_{j=1}^{N}W_{uj}\ket{\lambda_{j}},$ (41)

where $X_{i}$ denotes the Pauli-$X$ gate on qubit $i$. This is equivalent to
an amplitude encoding of the weights $W_{uj}$ in the basis of the eigenstates
$\ket{\lambda_{j}}$ of $H$. The choice of $u$ affects the set of weights
$W_{uj}$ and with it the probability to collapse into the corresponding
eigenstate $\ket{\lambda_{j}}$ upon measurement in this basis.

Including the $a$ ancilla qubits needed for the block-encoding, we can rewrite
the state in Eq. 41 in the eigenbasis of the unitary $V$ given in Eq. 39 as

$\ket{0}^{\otimes
a}\ket{u}=\frac{1}{\sqrt{2}}\sum_{j=1}^{N}W_{uj}\left(\ket{\mu_{+}^{(j)}}+\ket{\mu_{-}^{(j)}}\right).$
(42)

After the execution of QPE with $V$ as unitary operator controlled by the $m$
qubits of the phase register [7], we get:

$\xrightarrow{\text{QPE}}\frac{1}{\sqrt{2}}\sum_{j=1}^{N}W_{uj}\sum_{x=0}^{M-1}\left(a_{x+}^{(j)}\ket{\mu_{+}^{(j)}}\right.\\\
\left.+a_{x-}^{(j)}\ket{\mu_{-}^{(j)}}\right)\ket{x},$ (43)

where $M=2^{m}$ and $a_{x\pm}^{(j)}$ is the amplitude applied by the QPE
subroutine

$a_{x\pm}^{(j)}=\frac{1}{M}\sum_{z=0}^{M-1}e^{\frac{2\pi iz}{M}\left(\pm
y_{j}-x\right)},$ (44)

with

$y_{j}=\frac{M}{2\pi}\arccos\left(\frac{\lambda_{j}}{s\norm{H}_{\text{max}}}\right).$
(45)

The coefficients $a_{x\pm}^{(j)}$ are peaked around $x=\pm y_{j}\mod M$,
respectively. This tells us that the probability $P(x)$ to measure $x$ in the
phase register will peak around $\pm y_{j}\mod M$ for all
$j\in\\{1,\dots,N\\}$, and it is given by

$P(x)=\sum_{j=1}^{N}\frac{W_{uj}^{2}}{2}\left(\absolutevalue{a_{x+}^{(j)}}^{2}+\absolutevalue{a_{x-}^{(j)}}^{2}\right),$
(46)

where we used the orthonormality of $|\mu_{\pm}^{(j)}\rangle$.

After sampling multiple times from the phase register, we can obtain an
estimate of the probabilities Eq. 46, from which we need to extract the
eigenvalues $\lambda_{j}$ and the weights $W_{uj}^{2}$. Here, we assume that
we have sufficient resolution to determine the peaks and their centers.
Possible practical methods to determine the eigenvalues are discussed in Refs.
[20, 64, 70], and could be tested in practice for our problem. The center
$\tilde{y}_{j}$ of the peak can be translated into the eigenvalues by
inverting Eq. 45. The probability to sample from a certain eigenstate is
proportional to the weight $W_{uj}^{2}$ of interest. The last can be estimated
by adding the probabilities of the $2Q$ values closest to $\tilde{y}_{j}$

$W_{uj}^{2}\approx 2\sum_{q=-Q}^{Q-1}P(\lceil y_{j}\rceil+q),$ (47)

where, we assumed that $\tilde{y}_{j}\in\left[\lfloor y_{j}\rfloor,\lceil
y_{j}\rceil\right]$. The additional factor of 2 is due to the $1/2$ in Eq. 46.
The necessary size of $Q$ is further discussed in Sec. C.1.2. The full routine
is gathered in Algorithm 1.

Data:

  * •

Oscillator $u$

  * •

Mass $m_{u}$

  * •

Renormalization $s\norm{H}_{\mathrm{max}}$

  * •

Minimum eigenvalue gap $\Delta_{\lambda}^{(u)}$

  * •

Tolerances $\varepsilon$, $\delta$, $\zeta$

  * •

Controlled quantum walk operator $cV$ which contains the oracles
$O_{\vartheta}$ and $O_{P}$

  * •

Number of oscillators $N$

  * •

Number of distinguishable eigenvalues $N_{u}$

Result: Local response function $G_{uu}(s)$

begin

$m$, $n$, $a+l\leftarrow$ compute the required register sizes

$N_{S}\leftarrow$ compute the required sample size

for  _$t\in\\{1,\dots,N_{S}\\}$_ do

$\texttt{reg}_{\mathrm{phase}}\leftarrow\sum_{i=0}^{2^{m}-1}\ket{i}$

$\texttt{reg}_{\mathrm{state}}\leftarrow\ket{u}$

$\texttt{reg}_{\mathrm{ancilla}}\leftarrow\ket{0}^{\otimes a+l}$

for  _$i\in\\{0,\dots,m-1\\}$_ do

$cV^{2^{i}}$ to
$\texttt{[reg}_{\mathrm{phase}}^{(i)},\texttt{reg}_{\mathrm{state}},\texttt{reg}_{\mathrm{ancilla}}\texttt{]}$

$\mathrm{QFT}^{\dagger}$ to $\texttt{reg}_{\mathrm{phase}}$

$y_{t}\leftarrow$ measure $\texttt{reg}_{\mathrm{phase}}$

define $P(y)=\sum_{t=1}^{N_{S}}\delta(y-y_{t})/N_{S}$

$\bm{y}\leftarrow$ identify $N_{u}$ peaks in $P(y)$

for  _$j\in\\{j:y_{j}\in\bm{y}\\}$_ do

$\lambda_{j}\leftarrow s\norm{H}_{\mathrm{max}}\cos(2^{1-m}\pi y_{j})$

$W_{uj}^{2}\leftarrow 2\sum_{q=-Q}^{Q-1}P(\lceil y_{j}\rceil+q)$

$G_{uu}(s)\leftarrow
m_{u}^{-1}\sum_{j=1}^{N_{u}}W_{uj}^{2}/(\lambda_{j}+s^{2})$

Algorithm 1 Calculation of the local response function $G_{uu}(s)$ on a
quantum computer. The first part estimates classically the register and sample
size required by the quantum circuit following in the loop over $t$. The
classical post-processing in the end is required to extract the response
function from the measurements.

In Appendix D, we discuss a more general method which uses a modified Hadamard
test to estimate the non-local response function $G_{uv}(s)$.

##  IV Solving the random glued-trees problem

In this section, we show that a straightforward adaptation of the algorithm
presented in Sec. III can be used to solve the random glued-trees problem
first studied in Ref. [26] in the context of quantum walks. This problem is a
well-known example for which quantum algorithms give a provable, oracular
exponential speedup compared to classical methods. Recently, Ref. [24] also
formulated the problem in terms of coupled harmonic oscillators showing that
it can be solved in time polylogarithmic in the problem size, by using a
particular encoding and Hamiltonian simulation. We here complement this result
by showing that the problem can also be solved in polylogarithmic time using
QPE. In what follows, we describe the problem and give our main results, while
we provide the details in Appendix E.

Figure 4: Example of random glued-trees with $n_{c}=3$.

The basic setup of the random glued-trees problem is shown in Fig. 4. It
consists of a graph $G=(V,E)$ with two balanced trees each with $n_{c}$
columns, where the $n_{c}$-th column of the right tree is glued randomly to
the $n_{c}$-th column of the left tree. In this way, each vertex in the
$n_{c}$-th column of the right tree is connected to exactly two vertices of
the $n_{c}$-th column of the left tree, and vice versa. Let us denote by $A$
the adjacency matrix of a random glued-trees graph. As we see in Fig. 4, we
denoted the root of the left tree as the ENTRANCE vertex, and the root of the
right tree as the EXIT vertex. Access to the adjacency matrix $A$ allows us to
identify the ENTRANCE and EXIT vertices, since they are the only two vertices
with two neighbors. This allows us to define the following random glued-trees
problem.

######  Problem 1.

(Random glued-trees problem) Consider a random glued-trees graph with
$N=2(2^{n_{c}}-1)$ vertices and $n_{c}$ the number of columns in each tree.
Let $A$ denote its adjacency matrix. Let $\mathrm{deg}(i)$ be the function
that given the row $i$ returns the degree of the associated vertex, i.e., the
number of vertices connected to it. We assume that the vertices are randomly
ordered, but we are given the index of the ENTRANCE vertex, i.e., the row
associated with it. Assuming oracle access to the adjacency matrix $A$ and to
the function $\mathrm{deg}(i)$ the goal is to find the index of the EXIT
vertex.

Note that the function $\mathrm{deg}(i)$ is merely used to identify the
ENTRANCE and EXIT vertices, since they are the only vertices of degree $2$. 1
has an immediate mapping to a system of coupled oscillators as discussed in
Ref. [24]. In fact, we can associate with each vertex $u\in V$ a unit mass so
that the mass matrix is simply the identity matrix $I_{N}$. Moreover, we add a
unit spring constant to each edge in the graph and an additional unit spring
constant to the ENTRANCE and EXIT masses attached to a wall. In this way,
according to Eq. 4, the Hamiltonian $H$ coincides with the stiffness matrix
and is given by $H=3I_{N}-A$. Thus, the eigenvalues and eigenvectors of $H$
are in one-to-one correspondence with those of the adjacency matrix $A$. In
what follows, we assume that we have a block-encoding of $A$, which is
$3$-sparse and has $\lVert A\rVert_{\mathrm{max}}=1$.

The mapping of the problem to a quantum computer is the same as the one
described in Sec. II.2 for any coupled oscillator problem. In order to solve
it using QPE, we start by preparing the computational basis string associated
with the ENTRANCE vertex (known to us) that we denote as
$\ket{\mathrm{ENTRANCE}}$. The idea is now to perform QPE as in Fig. 3, using
a block-encoding of the adjacency matrix $A$. The glued-trees problem enjoys
an effective, exponential dimensionality reduction and under the action of the
adjacency matrix the evolution of the system is confined to the _column
subspace_ $\mathcal{H}_{\mathrm{col}}$ (see Appendix E for more details).
Following Ref. [26], one can obtain an estimate for the minimum eigenvalue gap
$\Delta_{\lambda}$ of the adjacency matrix projected onto the column subspace
$A_{\mathrm{col}}$ defined in Eq. E.8 666Ref. [26] gives the result in a
slightly different form, namely stating that the eigenvalue spacing is larger
than
$\frac{2\pi^{2}}{(1+\sqrt{2})n_{c}^{3}}+\mathcal{O}\left(\frac{1}{n_{c}^{4}}\right)$.
However, adapting the derivation one can show that Eq. 48 also holds.:

$\Delta_{\lambda}=\frac{8\pi^{2}}{n_{c}^{3}}+\mathcal{O}\left(\frac{1}{n_{c}^{4}}\right).$
(48)

Eq. 48 holds for any random glued-trees graph. Thus, we see that, despite the
fact that the size of the adjacency matrix increases exponentially with the
number of columns $n_{c}$, i.e., with the number of qubits $n$ needed to
encode it, the relevant, minimum gap scales only inverse polynomially. This is
exactly one of the features needed for a possible exponential quantum
advantage. The idea now is to take the number of qubits $m$ in the phase
register large enough to resolve the minimum eigenvalue gap
$\Delta_{\lambda}$. From Eq. C.8, we require

$m=\gamma\left\lceil\log_{2}\left(\frac{3\pi}{\Delta_{\lambda}}\right)\right\rceil\approx\gamma\left\lceil\log_{2}\left(\frac{3n_{c}^{3}}{8\pi}\right)\right\rceil,$
(49)

with $\gamma$ a constant, integer factor of our choice such that $\gamma\gg
1$. Accordingly, the number of oracle queries for a single run of QPE will be
of order

$N_{\mathrm{queries}}=\mathcal{O}(n_{c}^{3}).$ (50)

If Eq. 49 is satisfied, the QPE subroutine described in Sec. III, followed by
measurement of the ancilla register with post-selection on the
$\ket{0}^{\otimes a}$ outcome, behaves de facto as a measurement in the
eigenbasis of $A$ projected onto the column subspace. Thus, after one run of
QPE followed by post-selection on $\ket{0}^{\otimes a}$, we have prepared an
eigenstate of the adjacency matrix $A$ projected onto the column subspace in
the $n$-qubit register. A final measurement in the computational basis on this
register gives a vertex $u\in V$. We are interested in the probability
$\mathrm{Prob}(\mathrm{EXIT})$ that this procedure outputs the EXIT vertex. In
Appendix E, we obtain an analytical lower bound for this probability (see Eq.
E.19), from which we conclude that

$\mathrm{Prob}(\mathrm{EXIT})=\mathcal{O}\left(\frac{1}{n_{c}}\right),$ (51)

which finally proves the overall polynomial scaling of the algorithm in the
number of columns $n_{c}$, i.e., polylogarithmic in the size of $A$,
complementing the results obtained in Ref. [24] and Ref. [26] with time-
evolution-based algorithms.

##  V Conclusions

We have presented a quantum algorithm for obtaining response functions on
systems of coupled harmonic oscillators. These kinds of problems also
effectively emerge in the study of mechanical, vibrating structures, after
discretization, as well as other areas of interest. The particular analytical
structure of these functions makes them suitable for a quantum computer. The
algorithm is based on the QPE subroutine, that we use as an eigenvalue solver.
Importantly, the state preparation at the input of the algorithm just amounts
to the preparation of a product state. This allows us to analyze the scaling
of the algorithm in terms of space and queries of the matrix oracles. While
the algorithm achieves an exponential saving in space, the scaling of the
total number of oracles queries is problem-dependent and can be polynomial in
the number of qubits only when very specific conditions are met. We identify
such a list of sufficient conditions leading to a rigorous worst case
guarantee, but also leave open the possibility for improved practical scaling
when a priori information about the structure of the problem is available. In
addition, the scaling is mostly limited by statistical error coming from the
estimation of the weights of the response functions. Possible improvements of
the algorithm can be obtained by combining the algorithm with quantum
eigenvalue filtering techniques, as well as quantum amplitude estimation.
Finally, we have shown that the algorithm can be adapted to efficiently solve
the random-glued trees problems, similarly to other approaches based on
Hamiltonian simulation.

While we performed a basic scaling analysis, we believe that a deeper
understanding of our algorithm and its potential can be obtained by studying
its performances in practice. This could be simply understood via basic
classical simulations that aim at obtaining the minimum, relevant eigenvalue
gap for some class of oscillator problems. Such analysis would also shed light
on the kinds of problems for which large advantages are possible.
Additionally, a practical study of the algorithm performances using methods
developed for the estimation of multiple eigenvalues [20, 64, 70] tailored to
our problem, would be beneficial. Finally, an open question is whether
estimates of the response functions can also be obtained using time-evolution
based algorithms for coupled harmonic oscillator as in Ref. [24].

## Acknowledgements

SD is funded by the Bundesministerium für Wirtschaft und Klimaschutz (BMWK,
Federal Ministry for Economic Affairs and Climate Action) in the quantum
computing enhanced service ecosystem for simulation in manufacturing project
(QUASIM, Grant number. 01MQ22001A). AC acknowledges funding from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – Cluster of Excellence Matter and Light for Quantum
Computing (ML4Q) EXC 2004/1 – 390534769. AC acknowledge funding from the
German Federal Ministry of Education and Research (BMBF) in the funding
program “Quantum technologies – from basic research to market” (Project QAI2,
contract number 13N15585). MB is supported by the EPSRC Grant number
EP/W032643/1 and the Excellence Cluster - Matter and Light for Quantum
Computing (ML4Q).

## Data availabiltiy

The algorithm presented in this manuscript were tested for the example of
linearly coupled oscillators. For this case, it wasimplemented and emulated
from end-to-end in Qiskit [72], including compilation of the matrix oracles.
The corresponding repository containing all the required packages and a
demonstration notebook can be found at https://go.fzj.de/resp-QPE.

## References

  * Abrams and Lloyd [1999] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162 (1999).
  * Bauer  _et al._ [2020] B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Chemical Reviews 120, 12685 (2020).
  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * Dalzell  _et al._ [2023] A. M. Dalzell, S. McArdle, M. Berta, P. Bienias, C.-F. Chen, A. Gilyén, C. T. Hann, M. J. Kastoryano, E. T. Khabiboulline, A. Kubica, G. Salton, S. Wang, and F. G. S. L. Brandão, arXiv e-prints , arXiv:2310.03011 (2023).
  * 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, Nature Communications 5, 4213 (2014).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, Physics Reports 986, 1 (2022).
  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_, 10th ed. (Cambridge University Press, Cambridge, 2011).
  * Kitaev  _et al._ [2002] A. Y. Kitaev, A. H. Shen, and M. N. Vyalyi,  _Classical and Quantum Computation_ (American Mathematical Society, Boston, MA, USA, 2002).
  * Cleve  _et al._ [1998] R. Cleve, A. Ekert, C. Macchiavello, and M. Mosca, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 339 (1998).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Quantum 3, 163 (2019).
  * Low and Chuang [2017] G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118, 010501 (2017).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, 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–204.
  * Martyn  _et al._ [2021] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, PRX Quantum 2, 040203 (2021).
  * Tang and Tian [2023] E. Tang and K. Tian, arXiv e-prints , arXiv:2302.14324 (2023).
  * Poulin  _et al._ [2018] D. Poulin, A. Kitaev, D. S. Steiger, M. B. Hastings, and M. Troyer, Phys. Rev. Lett. 121, 010501 (2018).
  * Berry  _et al._ [2018] D. W. Berry, M. Kieferová, A. Scherer, Y. R. Sanders, G. H. Low, N. Wiebe, C. Gidney, and R. Babbush, npj Quantum Information 4, 22 (2018).
  * Grover and Rudolph [2002] L. Grover and T. Rudolph, arXiv e-prints , quant-ph/0208112 (2002).
  * Ozols  _et al._ [2013] M. Ozols, M. Roetteler, and J. Roland, ACM Trans. Comput. Theory 5, 10.1145/2493252.2493256 (2013).
  * McArdle  _et al._ [2022] S. McArdle, A. Gilyén, and M. Berta, arXiv e-prints , arXiv:2210.14892 (2022).
  * O’Brien  _et al._ [2019] T. E. O’Brien, B. Tarasinski, and B. M. Terhal, New Journal of Physics 21, 023022 (2019).
  * Dutkiewicz  _et al._ [2022] A. Dutkiewicz, B. M. Terhal, and T. E. O’Brien, Quantum 6, 830 (2022).
  * Ding and Lin [2023] Z. Ding and L. Lin, Quantum 7, 1136 (2023).
  * Somma [2019] R. D. Somma, New Journal of Physics 21, 123025 (2019).
  * Babbush  _et al._ [2023] R. Babbush, D. W. Berry, R. Kothari, R. D. Somma, and N. Wiebe, Phys. Rev. X 13, 041041 (2023).
  * Lee and Kanno [2023] Y. Lee and K. Kanno, arXiv e-prints , arXiv:2307.07478 (2023).
  * Childs  _et al._ [2003] A. M. Childs, R. Cleve, E. Deotto, E. Farhi, S. Gutmann, and D. A. Spielman, in  _Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing_, STOC ’03 (Association for Computing Machinery, New York, NY, USA, 2003) p. 59–68.
  * Costa  _et al._ [2019] P. C. S. Costa, S. Jordan, and A. Ostrander, Phys. Rev. A 99, 012323 (2019).
  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Phys. Rev. Lett. 110, 250504 (2013).
  * Montanaro and Pallister [2016] A. Montanaro and S. Pallister, Phys. Rev. A 93, 032324 (2016).
  * Jung and Langer [2013] M. Jung and U. Langer,  _Method of finite elements for engineers. An introduction to the numerical foundations and computer simulation._ , Vol. 2nd revised and extended ed (Wiesbaden : Springer Fachmedien Wiesbaden, 2013).
  * Trefethen and Bau [1997] L. N. Trefethen and D. Bau,  _Numerical Linear Algebra_ (SIAM, 1997).
  * Press  _et al._ [1992] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes in C: The Art of Scientific Computing_ , 2nd ed. (Cambridge University Press, USA, 1992).
  * Arbenz [2016] P. Arbenz, Lecture notes on solving large scale eigenvalue problems (2016).
  * Klocke [2011] F. Klocke,  _Manufacturing Processes 1: Cutting_ (Springer Berlin Heidelberg, 2011).
  * Maslo  _et al._ [2020] S. Maslo, B. Menezes, P. Kienast, P. Ganser, and T. Bergs, Procedia CIRP 93, 850–855 (2020).
  * Munoa  _et al._ [2016] J. Munoa, X. Beudaert, Z. Dombovari, Y. Altintas, E. Budak, C. Brecher, and G. Stepan, CIRP Annals 65 (2) , 785–808 (2016).
  * Quintana and Ciurana [2011] G. Quintana and J. Ciurana, International Journal of Machine Tools and Manufacture 51 (5) , 363–376 (2011).
  * Budak  _et al._ [2012] E. Budak, L. Tunç, S. Alan, and H. Nevzat Özgüven, CIRP Annals 61 (1) , 339–342 (2012).
  * Altintas  _et al._ [2008] Y. Altintas, G. Stepan, D. Merdol, and Z. Dombovari, CIRP Journal of Manufacturing Science and Technology 1 (1) , 35–44 (2008).
  * Schröder  _et al._ [2024] S. Schröder, S. Danz, P. Kienast, V. König, P. Ganser, and T. Bergs, Procedia CIRP 121C, 13–18 (2024).
  * Rudel  _et al._ [2022] V. Rudel, P. Kienast, G. Vinogradov, P. Ganser, and T. Bergs, Frontiers in Manufacturing Technology 2,1021029 (2022).
  * Note [1] The response function is also called the Green’s function for the system of coupled oscillators.
  * Triverio  _et al._ [2007] P. Triverio, S. Grivet-Talocia, M. S. Nakhla, F. G. Canavero, and R. Achar, IEEE Transactions on Advanced Packaging 30, 795 (2007).
  * Note [2] The number of qubits in the first register does not need to be equal to the one in the second register.
  * Ruiz-Perez and Garcia-Escartin [2017] L. Ruiz-Perez and J. C. Garcia-Escartin, Quantum Information Processing 16, 152 (2017).
  * Cuccaro  _et al._ [2004] S. A. Cuccaro, T. G. Draper, S. A. Kutin, and D. Petrie Moulton, arXiv e-prints , quant-ph/0410184 (2004).
  * Draper [2000] T. G. Draper, arXiv e-prints , quant-ph/0008033 (2000).
  * Vedral  _et al._ [1996] V. Vedral, A. Barenco, and A. Ekert, Phys. Rev. A 54, 147 (1996).
  * Note [3] The smaller the error on the eigenvalues, the smaller the error on the weights, and vice versa. However, for our analysis it is easier to think about them as independent tolerances.
  * Christandl  _et al._ [2004] M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, Phys. Rev. Lett. 92, 187902 (2004).
  * Gilles Brassard and Tapp [2002] M. M. Gilles Brassard, Peter Hø yer and A. Tapp, in  _Quantum Computation and Information_ , Vol. 305, edited by J. Samuel J. Lomonaco and H. E. Brandt (Contemporary Mathematics, 2002) pp. 53–74.
  * Lin and Tong [2020] L. Lin and Y. Tong, Quantum 4, 361 (2020).
  * Parrish and McMahon [2019] R. M. Parrish and P. L. McMahon, arXiv e-prints , arXiv:1909.08925 (2019).
  * Wall and Neuhauser [1995] M. R. Wall and D. Neuhauser, The Journal of Chemical Physics 102, 8011 (1995).
  * Chen and Guo [1996] R. Chen and H. Guo, The Journal of Chemical Physics 105, 1311 (1996).
  * Mandelshtam and Taylor [1997] V. A. Mandelshtam and H. S. Taylor, The Journal of Chemical Physics 106, 5085 (1997).
  * Ikegami  _et al._ [2010] T. Ikegami, T. Sakurai, and U. Nagashima, Journal of Computational and Applied Mathematics 233, 1927 (2010).
  * Sanders  _et al._ [2019] Y. R. Sanders, G. H. Low, A. Scherer, and D. W. Berry, Phys. Rev. Lett. 122, 020502 (2019).
  * Note [4] In fact, the oracles are ultimately used to prepare the state in Eq. B.5.
  * Giovannetti  _et al._ [2008] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. Lett. 100, 160501 (2008).
  * Aaronson [2015] S. Aaronson, Nature Physics 11, 291 (2015).
  * Gustavsen and Semlyen [1999] B. Gustavsen and A. Semlyen, IEEE Transactions on Power Delivery 14, 1052 (1999).
  * Kitaev [1995] A. Y. Kitaev, arXiv e-prints , quant-ph/9511026 (1995).
  * Svore  _et al._ [2014] K. M. Svore, M. B. Hastings, and M. Freedman, Quantum Info. Comput. 14, 306–328 (2014).
  * Rall [2021] P. Rall, Quantum 5, 566 (2021).
  * Childs [2010] A. M. Childs, Commun. Math. Phys 294, 581 (2010).
  * Berry and Childs [2012] D. W. Berry and A. M. Childs, Quantum Information and Computation 12, 29–62 (2012).
  * Lin [2022] L. Lin, arXiv e-prints , arXiv:2201.08309 (2022).
  * Note [5] If $U_{H}$ is not Hermitian, one can add one ancilla qubit and define the Hermitian unitary $U_{H}^{(\mathrm{new})}=\outerproduct{+}{-}\otimes U_{H}+\outerproduct{-}{+}\otimes U_{H}^{\dagger}$, which is a $(\alpha,a+1,0)$-block-encoding of $H$.
  * Wiebe and Granade [2016] N. Wiebe and C. Granade, Phys. Rev. Lett. 117, 010503 (2016).
  * Note [6] Ref. [26] gives the result in a slightly different form, namely stating that the eigenvalue spacing is larger than $\frac{2\pi^{2}}{(1+\sqrt{2})n_{c}^{3}}+\mathcal{O}\left(\frac{1}{n_{c}^{4}}\right)$. However, adapting the derivation one can show that Eq. 48 also holds.
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).
  * Newcomb [1966] R. Newcomb,  _Linear Multiport Synthesis_, McGraw-Hill electronic sciences series (McGraw-Hill, 1966).
  * Russer [2006] P. Russer,  _Electromagnetics, Microwave Circuit and Antenna Design for Communications Engineering_ (Artech House, Boston, 2006).
  * Foster [1924] R. M. Foster, Bell System Technical Journal 3, 259–267 (1924).
  * Nigg  _et al._ [2012] S. E. Nigg, H. Paik, B. Vlastakis, G. Kirchmair, S. Shankar, L. Frunzio, M. H. Devoret, R. J. Schoelkopf, and S. M. Girvin, Phys. Rev. Lett. 108, 240502 (2012).
  * Ciani  _et al._ [2024] A. Ciani, D. P. DiVincenzo, and B. M. Terhal,  _Lecture Notes on Quantum Electrical Circuits_ (TU Delft OPEN Publishing, Delft, 2024).
  * Note [7] We set $\Theta(0)=0$ and $\mathrm{sgn}(0)=1$.
  * Note [8] We use the standard definition for the $y$-rotation in here $R_{y}(\vartheta)=e^{-i\vartheta\sigma_{y}/2}$, where $\sigma_{y}$ is the second Pauli matrix.
  * Note [9] The first bit in the phase register just represents the sign of $\pm y\mod M$, so $y$ is represented with $m-1$ bits.
  * Note [10] An analogous derivation is possible if we assume $\lfloor\tilde{y}_{j}\rfloor=\lfloor y_{j}\rfloor$.
  * Qi  _et al._ [2010] F. Qi, S. Guo, and B.-N. Guo, Journal of Computational and Applied Mathematics 233, 2149 (2010).
  * Note [11] In our version of QPE we have two peaks associated with each eigenvalue, but we can simply sum the probability of each peak to get $W_{uj}^{2}$.
  * Hoeffding [1963] W. Hoeffding, Journal of the American Statistical Association 58, 13 (1963).
  * Dembo and Zeitouni [1998] A. Dembo and O. Zeitouni,  _Large Deviations Techniques and Applications_ , 2nd ed., Applications of Mathematics (Springer, New York, 1998).

##  Appendix A Analytical form of the response functions

In this appendix, we show that for a system of coupled harmonic oscillators
the matrix elements of the response function $G_{uv}(s)$ in Laplace domain can
be written as

$G_{uv}(s)=\frac{1}{\sqrt{m_{u}m_{v}}}\sum_{j=1}^{N}\frac{W_{uj}W_{vj}}{\lambda_{j}+s^{2}}.$
(A.1)

We remark that this is a well-known fact in the study of harmonic systems, in
particular in the theory of electrical circuits, where it appears in various
forms [73, 74, 75]. Also, it is the base for the black-box quantization method
for quantum electrical circuits [76]. Here we provide a derivation for our
problem adapted from Ref. [77].

Let us consider the case in which we apply to oscillator $v$ an impulse
$f_{v}(t)=\tilde{f}_{v}\delta(t)$, while no force is applied to the other
oscillators, i.e., $f_{v^{\prime}}(t)=0$ for $v^{\prime}\neq v$. In terms of
the normal mode variables in Eq. 6, the dynamical equation in time domain Eq.
9 reads

$\frac{d^{2}y_{j}}{dt^{2}}=-\lambda_{j}y_{j}+\frac{W_{vj}}{\sqrt{m_{v}}}\tilde{f}_{v}\delta(t),$
(A.2)

which in Laplace domain gives

$Y_{j}(s)=\frac{1}{\sqrt{m_{v}}}\frac{W_{vj}}{\lambda_{j}+s^{2}}\tilde{f}_{v},$
(A.3)

for $j\in\\{1,\dots,N\\}$. From Eq. 8, we get that the Laplace transform of
$x_{u}(t)$ can be written as

$X_{u}(s)=\frac{1}{\sqrt{m_{u}}}\sum_{j=1}^{N}W_{uj}Y_{j}(s)\\\
=\frac{1}{\sqrt{m_{u}m_{v}}}\sum_{j=1}^{N}\frac{W_{uj}W_{vj}}{\lambda_{j}+s^{2}}\tilde{f}_{v},$
(A.4)

from which Eq. A.1 follows.

##  Appendix B Hermitian block-encoding from the matrix oracles

Here we detail how a Hermitian $(\alpha,a,0)$-block-encoding $U_{H}$ of a
Hamiltonian $H$ can be constructed using the matrix oracles $O_{P}$ and
$O_{\vartheta}$ defined in Eq. 15 and Eq. 17, respectively. In particular, we
limit ourselves to the case when $H$ is real, i.e., symmetric, and all its
diagonal elements are positive, which is the case for the matrix in Eq. 4 we
are interested in. We start by considering a $2(n+1)$-qubit register, where
the first $a=n+2$ qubits will be the ancilla qubits needed for the block-
encoding. We take the block-encoding unitary $U_{H}$ to be of the following
form

$U_{H}=U_{T}^{\dagger}\mathrm{SWAP}_{n+1}U_{T},$ (B.1)

where $\mathrm{SWAP}_{n+1}$ swaps two $(n+1)$-qubit registers, i.e.,
$\mathrm{SWAP}_{n+1}\ket{\psi}\ket{\varphi}=\ket{\varphi}\ket{\psi}$ for any
$(n+1)$-qubit states $\ket{\psi},\ket{\varphi}$. Thus, Eq (30) and Eq. (B.1)
imply that

$\alpha H=\bra{0}^{n+2}U_{T}^{\dagger}\mathrm{SWAP}_{n+1}U_{T}\ket{0}^{\otimes
n+2}.$ (B.2)

We further assume that the map $U_{T}\ket{0}^{\otimes n+2}$, which maps an
$n$-qubit quantum state to a $2(n+1)$-qubit quantum state can be written as
($N=2^{n}$)

$U_{T}\ket{0}^{\otimes n+2}=\sum_{u=1}^{N}\ket{\psi_{u},0,u}\bra{u},$ (B.3)

where $\ket{\psi_{u}}$ denotes an $(n+1)$-qubit state to be determined,
$\ket{0}$ the single-qubit zero state and $\ket{u}$ an $n$-qubit computational
basis state. Plugging Eq. B.3 into Eq. B.2, we obtain that the states
$\ket{\psi_{u}}$ need to satisfy

$\alpha H_{uv}=\innerproduct{\psi_{u},0,u}{0,v,\psi_{v}}.$ (B.4)

A possible solution is to set $\alpha$ as in Eq. 40 and take

$\ket{\psi_{u}}=\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}(i\mathrm{sgn}(u-v))^{\Theta(-H_{uv})}\\\
\times\left(\sqrt{\frac{|H_{uv}|}{\|H\rVert_{\mathrm{max}}}}\ket{0}+\sqrt{1-\frac{|H_{uv}|}{\|H\rVert_{\mathrm{max}}}}\ket{1}\right)\ket{v},$
(B.5)

where $\mathrm{sgn}(\cdot)$ is the sign function 777We set $\Theta(0)=0$ and
$\mathrm{sgn}(0)=1$.. Note that Eq. B.5 requires the diagonal elements to be
non-negative in order to give a valid block-encoding. If this is not the case,
one can always enforce this property by adding a suited matrix proportional to
the identity to $H$.

Given the form of the block-encoding in Eq. B.1, all we need to show is how to
construct a unitary $U_{T}$ that satisfies Eq. B.3 with states
$\ket{\psi_{u}}$ given in Eq. B.5. The implementation of $U_{T}$ consists of
three parts. First preparing a superposition of all states $\ket{v}$ with non-
zero matrix entries $H_{uv}$. Second, encoding the matrix entries in the
amplitudes, and last, preparing the corresponding signs. All these steps are
shown in Fig. 5 as a quantum circuit that involves the matrix oracles $O_{P}$
and $O_{\vartheta}$.

[ wire types=q,b,q,b,b,q, classical gap =0.07cm] $\ket{0}$&
[5,steps=5,style=dashed,rounded corners,fill=teal!20, inner xsep=2pt,
background, label style=label position=above,
anchor=south,yshift=-0.2cm]matrix encodingR_y
$\ket{0}^{\otimes r}$ [4]O_ϑ $\rightarrow$ -1 [4]O_ϑ^† $\rightarrow$
$\ket{0}$ $\rightarrow$ S[4,steps=3,style=dashed,rounded corners,fill=teal!40,
inner xsep=2pt, background, label style=label position=below,
anchor=north,yshift=-0.2cm]complex sign 3$\rightarrow$
$\ket{0}^{\otimes n}$ [2]O_P [2,steps=1,style=dashed,rounded
corners,fill=teal!35, inner xsep=2pt, background, label style=label
position=below, anchor=north,yshift=-0.15cm]superposition$\rightarrow$
[3,label style=yshift=4mm]SUB [3,label style=yshift=4mm]ADD
$\ket{u}$ $\rightarrow$ $\rightarrow$
$\ket{0}$ $\rightarrow$ Z $\rightarrow$

Figure 5:  Circuit for the state transformation $U_{T}$. It consists of three
parts. The preparation of the superposition with $O_{P}$, the matrix encoding
in the amplitude and the multiplication with the complex sign. The order of
the registers in here is chosen for a compact circuit layout. The first and
fourth register combine into $\ket{\psi_{u}}$ of Eq. B.3. The second last
register is $\ket{u}$ in Eq. B.3. The second, third and last registers (in
blue) are of auxiliary nature and will be reinitialized after the full block-
encoding of $H$. In the second register we store the intermediate angle (Eq.
16). The third register is used for the binary sign $\Theta(-H_{uv})$ and the
last for $\Theta(\mathrm{sgn}(u-v))$.

We now describe in detail each transformation. As we see from Fig. 5 we have a
total of six different registers. The three register highlighted in blue are
additional ancilla registers needed for the implementation of the oracles,
namely:

  * •

an $r$-qubit ancilla register to store the value of the angles
$\vartheta_{uv}$ given in Eq. 16;

  * •

an ancilla qubit to store the value of the sign of the matrix elements
$H_{uv}$;

  * •

an ancilla qubit (bottom one) to store the value of the sign of $u-v$.

The remaining registers form a $(2n+1)$-qubit register that together with an
additional qubit (not shown in the figure) will be used to implement the
$2(n+1)$-qubit block-encoding we described at the beginning of this appendix.

The system starts in the state $\ket{0}^{\otimes r+2+n}\ket{u}\ket{0}$. After,
the application of the position oracle $O_{P}$ defined in Eq. 15 the state
becomes

$\xrightarrow{O_{P}}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}\ket{0}^{\otimes
r+2}\ket{v}\ket{u}\ket{0},$ (B.6)

and after the oracle $O_{\vartheta}$ given in Eq. 17

$\xrightarrow{O_{\vartheta}}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}\ket{0}\ket{\vartheta_{uv},\Theta(-H_{uv})}\ket{v,u}\ket{0}.$
(B.7)

We continue by reading the angle $\vartheta_{uv}$ in the ancilla register with
a series of controlled $y$-rotations
$\mathrm{c}R_{y}=\prod_{k=1}^{r}\mathrm{c}^{l}R_{y}(2^{-k}\pi)$ with target
the ancilla qubit on the top of Fig. 5 and control the $k$-th qubit in the
binary representation of $\ket{\vartheta_{uv}}$ in Eq. 19 888 We use the
standard definition for the $y$-rotation in here
$R_{y}(\vartheta)=e^{-i\vartheta\sigma_{y}/2}$, where $\sigma_{y}$ is the
second Pauli matrix. . This gives

$\xrightarrow{\mathrm{c}R_{y}}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}\left(\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{0}+\sqrt{1-\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{1}\right)\\\
\otimes\ket{\vartheta_{uv},\Theta(-H_{uv})}\ket{v,u}\ket{0}.$ (B.8)

It remains to implement the coefficient $i\text{\,sgn}(u-v)$ for negative
$H_{uv}$. The phase $i^{\Theta(-H_{uv})}$ can be realized with an S-gate to
the sign-qubit $\ket{\Theta(-H_{uv})}$:

$\xrightarrow{S}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}i^{\Theta(-H_{uv})}\\\
\times\left(\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{0}+\sqrt{1-\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{1}\right)\\\
\otimes\ket{\vartheta_{uv},\Theta(-H_{uv})}\ket{v,u}\ket{0}.$ (B.9)

For the relative sign, we subtract the $\ket{v}$ register from the $\ket{u}$
for which we need to increase the target register by one qubit to store the
sign

$\ket{v}\ket{u,0}\xrightarrow{\text{SUB}}\ket{v}\ket{u-v\mod 2^{n+1}}.$ (B.10)

Notice that the most significant bit of $u-v\mod 2^{n+1}$, which is a
$(n+1)$-bit number, encodes the sign of the regular subtraction $u-v$, since
both $u$ and $v$ are $n$-bit numbers. Thus, we can write

$\ket{u-v\mod 2^{n+1}}=\ket{u-v\mod 2^{n}}\ket{\Theta(v-u)},$ (B.11)

where the register associated with $\ket{u-v\mod 2^{n}}$ is an $n$-qubit
register. The subtraction gives the state

$\xrightarrow{\mathrm{SUB}}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}i^{\Theta(-H_{uv})}\\\
\times\left(\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{0}+\sqrt{1-\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{1}\right)\\\
\otimes\ket{\vartheta_{uv},\Theta(-H_{uv})}\ket{v,u-v\mod
2^{n}}\ket{\Theta(v-u)},$ (B.12)

A controlled-$Z$ gate with control the qubit in the state
$\ket{\Theta(-H_{uv})}$ and target the one in the state $\ket{\Theta(u-u)}$
(represented by the qubit at the bottom in Fig. 5) yields the desired
coefficient

$\xrightarrow{\mathrm{c}Z}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}(i\mathrm{sgn}(u-v))^{\Theta(-H_{uv})}\\\
\times\left(\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{0}+\sqrt{1-\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{1}\right)\\\
\otimes\ket{\vartheta_{uv},\Theta(-H_{uv})}\ket{v,u-v\mod
2^{n}}\ket{\Theta(v-u)},$ (B.13)

At last we invert the subtraction SUB using the adder ADD

$\xrightarrow{\mathrm{ADD}}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}(i\mathrm{sgn}(u-v))^{\Theta(-H_{uv})}\\\
\times\left(\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{0}+\sqrt{1-\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{1}\right)\\\
\otimes\ket{\vartheta_{uv},\Theta(-H_{uv})}\ket{v,u}\ket{0},$ (B.14)

and the oracle $O_{\vartheta}$ using $O_{\vartheta}^{\dagger}$

$\xrightarrow{O_{\vartheta}^{\dagger}}\frac{1}{\sqrt{s}}\sum_{v\in\tilde{N}_{G}[u]}(i\mathrm{sgn}(u-v))^{\Theta(-H_{uv})}\\\
\times\left(\sqrt{\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{0}+\sqrt{1-\frac{\absolutevalue{H_{uv}}}{\norm{H}_{\text{max}}}}\ket{1}\right)\\\
\otimes\ket{0}^{\otimes r+1}\ket{v,u}\ket{0},$ (B.15)

and thus, we have effectively prepared the state $\ket{\psi_{u}}$ in Eq. B.5.
Notice that in this procedure the $r+1$ ancilla qubits needed for the
implementation of $O_{\vartheta}$ and the ancilla qubit that stores the sign
of $u-v$ return to the zero state, and this is also the case for the full
block-encoding in Eq. 30. This is why we should not count them among the $a$
ancilla qubits of the block-encoding, which are thus $a=n+2$.

The circuit implementation of $U_{H}$ in the block-encoding Eq. B.1 is shown
in Fig. 6. Note that the total number of qubits was increased by one. This is
the last one of the $a=n+2$ ancilla qubits necessary for the Hilbert-space
extension. The operator in front is reflecting around $\ket{0}^{\otimes n+2}$
and transforms $U_{H}$ into the quantum walk operator defined in Eq. 35. It
consists of a multi-controlled NOT-gate and a regular NOT-gate applied to an
ancilla qubit initialized in $\ket{-}$. The QPE demands a controlled version
which is why we have another qubit which controls the CNOT gate in the
reflection and the SWAP gates in $U_{H}$.

[wire types=q,q,q,q,b,b,b,q, classical gap = 0.07cm] Control:& 1 1 5 2
$\ket{-}$ H
$\ket{0}$-1[ 6,steps=4, style=dashed,rounded corners,fill=teal!20, inner
xsep=2pt, background, label style=label position=below, anchor=north,
yshift=-0.2cm]block-encoding $U_{H}$1
$\ket{0}$ -1[5]U_T [5]U_T^†
$\ket{0}^{\otimes r+1}$
$\ket{0}^{\otimes n}$ -2 1
$\ket{u}$
$\ket{0}$

Figure 6:  Controlled-$V$ operator with $V$ defined in Eq. 35. The circuit
starts with a controlled reflection around $\ket{0}^{\otimes n+2}$ that is
implemented by a multi-controlled and a normal CNOT gate to an ancilla qubit
initialized in $\ket{-}$. After this we see the controlled version of the
block-encoding $U_{H}=U_{T}^{\dagger}\mathrm{SWAP}_{n+1}U_{T}$.

###  B.1 Qubit requirements

The full QPE estimation is applied to three registers. The $m$-qubit phase
register, the $n$-qubit state register and the $(a+l)$-ancilla register used
by the controlled walk operator $\mathrm{c}V$. We will focus on the $n+a+l+1$
qubits required for the implementation of $\mathrm{c}V$. For this we start
with $U_{T}$ described in Fig. 5 . The total number of qubits $n_{U_{T}}$
required for the state transformation is given by

$n_{U_{T}}=2n+1+\max(a_{P},r+1+a_{\vartheta},r+2+a_{\mathrm{ADD}}).$ (B.16)

where $(r+1)$ is the number of qubits used by $O_{\vartheta}$ to store the
angle $\vartheta_{uv}$ and the sign of $H_{uv}$. Similar to $m$ it effects the
accuracy of the results. Moreover, the number of ancilla qubits necessary for
the position and angle oracles, and the quantum adder are given by $a_{P}$,
$a_{\vartheta}$ and $a_{\mathrm{ADD}}$, respectively. In Eq. B.16, we used
that the ancilla qubits used in $O_{P}$ can be reused in $O_{\vartheta}$,
while the ancilla qubits used there can be used again in $O_{\mathrm{ADD}}$.

The state transformation $U_{T}$ appears twice in the circuit for
$\mathrm{c}V$ (cf. Fig. 6). The circuit is further extended by a control qubit
and the last qubit necessary for the block-encoding. Hence, we have a total
number of $n_{\mathrm{c}V}$ qubits necessary for the controlled quantum walk,
which is

$n_{\mathrm{c}V}=2n+3+\max(a_{P},r+1+a_{\vartheta},r+2+a_{\mathrm{ADD}},n+2).$
(B.17)

Here we see additional $n+2$ ancilla qubits which can be reused. The first
$n+1$ of them are required for the Toffoli-gates in the mutlicontrolled NOT
gate in Fig. 6 and the last one is the target of the same. In total we have
$m$ repetitions of $cV$ with varying control qubit which leads to the scaling
described in Eq. 23. There we replaced $m$ with the scaling Eq. C.34 described
in Appendix C.

##  Appendix C Detailed error analysis

In this appendix, we discuss the effect of discretization and statistical
errors on the estimated eigenvalues $\tilde{\lambda}_{j}$ and weights
$\tilde{W}_{uj}^{2}$ and how one can control the maximum error by increasing
the size of the phase register. This analysis proves the statement of 1.

In particular, we analyze the effect of the discretization of the phase when
it is stored as binary number in the $m$-qubit register with focus on the
eigenvalue in Sec. C.1.1 and the weights in Sec. C.1.2. The size of the phase
register determines the number of oracle queries necessary to implement the
$\mathrm{controlled}$-$V^{2^{k}}$ operations in Algorithm 1 shown in Fig. 3.
We discuss this connection in Sec. C.1.3. We point out that the number of
oracle queries constitutes the main limitation of the algorithm in most
_practical_ cases (see discussion in Sec. II.4). At last, in Sec. C.2, we
describe the connection between the sample size and the statistical accuracy
of the estimated weights. In the following, for notational simplicity, we
focus on the non-degenerate systems, but the results carry over to the
degenerate case as discussed in Sec. II.2.

###  C.1 Phase register size for given tolerances

####  C.1.1 Eigenvalue approximation

The following analysis is more easily understood by considering the case in
which we perform QPE using the operator $V$ as in Fig. 3, but taking as
initial state in the $n$-qubit register a single eigenstate of $H$
$\ket{\lambda}$. When we measure the $m$-qubit phase register at the end of
QPE we get a bitstring that represents an integer number between $0$ and $M-1$
with $M=2^{m}$. This gives an estimate of $\pm y\mod M\in[0,M)$ with $y$ (cf.
Eq. 45) given by

$y=\frac{M}{2\pi}\arccos\left(\frac{\lambda}{s\lVert
H\rVert_{\mathrm{max}}}\right).$ (C.1)

Accordingly, there exists a $\tilde{y}\in\\{0,1,\dots,2^{m-1}-1\\}$ such that
999The first bit in the phase register just represents the sign of $\pm y\mod
M$, so $y$ is represented with $m-1$ bits.

$|y-\tilde{y}|\leq\frac{1}{2},$ (C.2)

to which we can associate an estimated eigenvalue $\tilde{\lambda}$ by
inverting Eq. C.1. The phase stored in the ancilla register corresponds to the
estimated $\tilde{\lambda}$ as follows

$\tilde{y}=\frac{M}{2\pi}\arccos\left(\frac{\tilde{\lambda}}{s\lVert
H\rVert_{\mathrm{max}}}\right).$ (C.3)

Thus, we get

$\frac{M}{2\pi}\absolutevalue{\arccos\left(\frac{\lambda}{s\lVert
H\rVert_{\mathrm{max}}}\right)-\arccos\left(\frac{\tilde{\lambda}}{s\lVert
H\rVert_{\mathrm{max}}}\right)}\leq\frac{1}{2}$ (C.4)

The mean value theorem for the arc-cosine function states that for
$x>\tilde{x}$, $\exists x_{0}\in(\tilde{x},x)$ such that

$\frac{|\arccos
x-\arccos\tilde{x}|}{|x-\tilde{x}|}=\frac{1}{\sqrt{1-\left(x_{0}\right)^{2}}}.$
(C.5)

Applying the mean value theorem to our case we get that
$\exists\lambda_{0}\in(\tilde{\lambda},\lambda)$ if $\lambda>\tilde{\lambda}$
or $\exists\lambda_{0}\in(\lambda,\tilde{\lambda})$ if
$\lambda<\tilde{\lambda}$, such that

$s\lVert H\rVert_{\mathrm{max}}\frac{\left|\arccos\left(\frac{\lambda}{s\lVert
H\rVert_{\mathrm{max}}}\right)-\arccos\left(\frac{\tilde{\lambda}}{s\lVert
H\rVert_{\mathrm{max}}}\right)\right|}{|\lambda-\tilde{\lambda}|}\\\
=\frac{1}{\sqrt{1-\left(\frac{\lambda_{0}}{s\lVert
H\rVert_{\mathrm{max}}}\right)^{2}}}\geq 1,$ (C.6)

where the last inequality holds since $0\leq\lambda_{0}\leq s\lVert
H\rVert_{\mathrm{max}}$. Combining Eq. C.4 and Eq. C.6, we get

$|\lambda-\tilde{\lambda}|\leq\frac{\pi s\lVert H\rVert_{\mathrm{max}}}{M}.$
(C.7)

which tells us that the error of the estimated eigenvalue is below
$\varepsilon$ if the number of qubits in the phase register satisfies

$m\geq m_{\text{min}}^{(1)}=\left\lceil\log_{2}{\left(\frac{\pi
s\norm{H}_{\text{max}}}{\varepsilon}\right)}\right\rceil.$ (C.8)

####  C.1.2 Weight approximation

We analyze the error on the estimated weights $W_{uj}^{2}$ of the local
response function defined in Eq. 47 due to the finite size of the phase
register. In our analysis we assume that our estimated peak $\tilde{y}_{j}$
satisfies $\lceil\tilde{y}_{j}\rceil=\lceil y_{j}\rceil$ with $y_{j}$
associated with the true eigenvalue $\lambda_{j}$ as in Eq. 45 101010An
analogous derivation is possible if we assume
$\lfloor\tilde{y}_{j}\rfloor=\lfloor y_{j}\rfloor$. Accordingly, our estimator
for the weight $W_{uj}^{2}$ is

$\tilde{W}_{uj}^{2}=2\sum_{q=-Q}^{Q-1}P(\lceil y_{j}\rceil+q),$ (C.9)

with the range $Q$ to be determined.

The parameter $Q$ should be taken large enough to guarantee

$\absolutevalue{\tilde{W}_{uj}^{2}-W_{uj}^{2}}\leq\delta$ (C.10)

An increase in $Q$ prevents two things. First, it limits the probability we
miss by adding only over the $2Q$ probabilities for values closest to the
center of a peak. This is connected to the lower bound of
$\tilde{W}_{uj}^{2}-W_{uj}^{2}$. The second is that we need to increase the
distance between two peaks to fit $2Q$ between them which prevents probability
“leakage” between them. This is connected to the upper bound of
$\tilde{W}_{uj}^{2}-W_{uj}^{2}$.

For the lower limit we need to study the probability to measure $x$ in the
phase register defined in Eq. 46 and focus only on one contribution

$P(x)\geq\frac{W_{uj}^{2}}{2}\absolutevalue{a_{x+}^{(j)}}^{2}.$ (C.11)

We will now further develop $|a_{x+}^{(j)}|^{2}$ by using its definition Eq.
44 and the geometric series

$\sum_{z=0}^{M-1}x^{z}=\frac{1-x^{M}}{1-x}.$ (C.12)

Thus, we have

$\absolutevalue{a_{x+}^{(j)}}^{2}=\frac{1}{M^{2}}\absolutevalue{\frac{1-e^{2\pi
i(y_{j}-x)}}{1-e^{2\pi i(y_{j}-x)/M}}}^{2}\\\
=\frac{1}{M^{2}}\left[\frac{\sin(\pi(y_{j}-x))}{\sin(\pi(y_{j}-x)/M)}\right]^{2}\geq\left[\frac{\sin(\pi(y_{j}-x))}{\pi(y_{j}-x)}\right]^{2}.$
(C.13)

where we used $\absolutevalue{\sin x}\leq\absolutevalue{x}$. Using Eq. C.9,
Eq. C.11 and Eq. C.13, we get

$\tilde{W}_{uj}^{2}\geq
W_{uj}^{2}\sum_{q=-Q}^{Q-1}\left[\frac{\sin(\pi(y_{j}-\lceil
y_{j}\rceil-q))}{\pi(y_{j}-\lceil y_{j}\rceil-q)}\right]^{2},$ (C.14)

which is minimal for $\lceil y_{j}\rceil-y_{j}=1/2$. Hence,

$\tilde{W}_{uj}^{2}\geq
W_{uj}^{2}\sum_{q=-Q}^{Q-1}\left[\frac{1}{\pi(1/2+q)}\right]^{2}\\\
=2W_{uj}^{2}\sum_{q=1}^{Q}\left[\frac{1}{\pi(q-1/2)}\right]^{2}$ (C.15)

Now, we introduce the trigamma function

$\Psi^{(1)}(z)=\sum_{q=0}^{+\infty}\frac{1}{(z+q)^{2}},$ (C.16)

which for half integer values takes the form

$\Psi^{(1)}\left(Q+\frac{1}{2}\right)=\frac{\pi^{2}}{2}-\sum_{q=1}^{Q}\frac{1}{\left(q-1/2\right)^{2}},$
(C.17)

with $Q\in\mathbb{N}$. $\Psi^{(1)}(z)$ satisfies

$\Psi^{(1)}\left(z\right)\leq\frac{z+1}{z^{2}}\leq\frac{2}{z},\qquad\forall
z>0.$ (C.18)

More details about the trigamma function and a derivation of Eq. C.18 can be
found in [82]. Combining Eqs. C.15, C.17 and C.18 leads to

$\tilde{W}_{uj^{2}}\geq
W_{uj}^{2}\left(1-\frac{2}{\pi^{2}}\psi^{(1)}\left(Q+\frac{1}{2}\right)\right)\\\
\geq W_{uj}^{2}\left(1-\frac{4}{\pi^{2}\left(Q+\frac{1}{2}\right)}\right)$
(C.19)

Since, $W_{uj}^{2}\leq 1$, this tells us that our estimation
$\tilde{W}_{uj}^{2}$ is at most $\delta$ less than $W_{uj}^{2}$ if we choose a
$Q$ such that

$\frac{4}{\pi^{2}\left(Q+1/2\right)}\leq\delta.$ (C.20)

or

$Q\geq
Q_{\mathrm{min}}^{(1)}=\left\lceil\frac{4}{\pi^{2}\delta}-\frac{1}{2}\right\rceil.$
(C.21)

For the upper limit, let us consider

$\tilde{W}_{uj}^{2}-W_{uj}^{2}\\\
=\sum_{i=0}^{N-1}W_{ui}^{2}\sum_{q=-Q}^{Q-1}\left(\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(i)}}^{2}+\absolutevalue{a_{\lceil
y_{j}\rceil+q,-}^{(i)}}^{2}\right)-W_{uj}^{2}.$ (C.22)

Since $\sum_{q=-Q}^{Q-1}\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(j)}}^{2}\leq 1$, we obtain

$\tilde{W}_{uj}^{2}-W_{uj}^{2}\\\ \leq\sum_{i\neq
j}W_{ui}^{2}\sum_{q=-Q}^{Q-1}\left(\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(i)}}^{2}+\absolutevalue{a_{\lceil
y_{j}\rceil+q,-}^{(i)}}^{2}\right)\\\
+W_{uj}^{2}\sum_{q=-Q}^{Q-1}\absolutevalue{a_{\lceil
y_{j}\rceil+q,-}^{(j)}}^{2}.$ (C.23)

Let us consider the coefficients $\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(i)}}^{2}$. We can rewrite them as

$\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(i)}}^{2}=\frac{1}{M^{2}}\left[\frac{\sin(\pi(\gamma_{i}+\lceil
y_{i}\rceil-\lceil y_{j}\rceil-q))}{\sin(\pi(\gamma_{i}+\lceil
y_{i}\rceil-\lceil y_{j}\rceil-q)/M)}\right]^{2}$ (C.24)

with $\gamma_{i}=y_{i}-\lceil y_{i}\rceil\in(-1,0]$. We assume that $|\lceil
y_{i}\rceil-\lceil y_{j}\rceil>2Q$ $\forall i\neq j$, i.e., that the peaks are
separated by at lest $2Q$. Using the fact that $\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(i)}}^{2}$ decreases by increasing $|\lceil
y_{i}\rceil-\lceil y_{j}\rceil-q|$ for integer values for any fixed
$\gamma_{i}\in(-1,0]$, we get

$\absolutevalue{a_{\lceil
y_{j}\rceil+q,+}^{(i)}}^{2}\leq\frac{1}{M^{2}}\left[\frac{\sin(\pi(\gamma_{i}+Q+1))}{\sin(\pi(\gamma_{i}+Q+1)/M)}\right]^{2}=\\\
|a_{\gamma_{i}+Q+1}|^{2}.$ (C.25)

The function $|a_{\gamma+Q+1}|^{2}$ has exactly one local maximum for
$\gamma\in(-1,0]$, which coincides with the global maximum in the interval. We
simply denote by $\bar{\gamma}\in(-1,0]$ the point where
$|a_{\gamma+Q+1}|^{2}$ is maximal. The same analysis can be repeated for
$\absolutevalue{a_{\lceil y_{j}\rceil+q,-}^{(i)}}^{2}$. Thus, we get

$\tilde{W}_{uj}^{2}-W_{uj}^{2}\leq\sum_{i\neq
j}2W_{ui}^{2}\sum_{q=-Q}^{Q-1}\absolutevalue{a_{\bar{\gamma}+Q+1}}^{2}\\\
+W_{uj}^{2}\sum_{q=-Q}^{Q-1}\absolutevalue{a_{\lceil
y_{j}\rceil+q,-}^{(j)}}^{2}.$ (C.26)

We will further assume that $Q<\absolutevalue{\pm y_{j}}<M/2-Q$, which
guarantees that the two peaks $y_{j}$ and $-y_{j}$ are also $2Q$ away from
each other and Eq. C.25 is satisfied for $i=j$. The previous condition fails
only for $\lambda_{j}\approx s\norm{H}$, where $\pm y_{j}\approx 0$. This can
be prevented by further reducing $H$ with the renormalization factor
$\cos(2\pi QM^{-1})$.

Now, we replace the amplitude in the last term of Eq. C.26 as well and have

$\tilde{W}_{uj}^{2}-W_{uj}^{2}\leq\left(\sum_{i\neq
j}2W_{ui}^{2}+W_{uj}^{2}\right)\sum_{q=-Q}^{Q-1}\absolutevalue{a_{\bar{\gamma}+Q+1}}^{2}.$
(C.27)

At last, we replace the amplitudes using Eq. C.13 and use the normalization of
the eigenvalues $\sum_{i}W_{ui}^{2}=1$ to find

$\tilde{W}_{uj}^{2}-W_{uj}^{2}\leq\frac{4Q}{M^{2}}\absolutevalue{\frac{\sin(\pi(\bar{\gamma}+Q+1))}{\sin(\pi(\bar{\gamma}+Q+1)/M)}}^{2}.$
(C.28)

Now we use $\absolutevalue{\sin x}\leq 1$ and $\absolutevalue{\sin
x}\geq\absolutevalue{2x/\pi},\forall\absolutevalue{x}\leq\pi/2$ and obtain
that our estimation is less than $\delta$ larger than the exact weight
$W_{uj}^{2}$ if

$\tilde{W}_{uj}^{2}-W_{uj}^{2}\leq\frac{Q}{(\bar{\gamma}+Q+1)^{2}}\leq\frac{1}{Q}\leq\delta,$
(C.29)

if $2Q\leq M-1$. This yields another lower limit for $Q$, which is

$Q^{(2)}_{\mathrm{min}}=\left\lceil\frac{1}{\delta}\right\rceil.$ (C.30)

A comparison between Eqs. C.21 and C.30 shows that the latter is larger and
therefore sufficient to fulfill both conditions. Thus, we set

$Q\geq\left\lceil\frac{1}{\delta}\right\rceil$ (C.31)

The previous analysis depends on the fact that we have a sufficient number of
qubits $m$ in the phase register, so that we can distinguish between two
neighboring phases. Mathematically, this translates into the condition

$\min_{i}\absolutevalue{y_{i}-y_{j}}\geq\min_{i}\frac{M\absolutevalue{\lambda_{i}-\lambda_{j}}}{2\pi
s\norm{H}_{\text{max}}}\geq 2Q,$ (C.32)

where we used the results obtained in Sec. C.1.1. We can guarantee this if we
take the number of qubits in the phase register $m$ such that

$m\geq m_{\text{min}}^{(2)}=\left\lceil\log_{2}\frac{4\pi
s\norm{H}_{\text{max}}}{\delta\Delta_{\lambda}^{(u)}}\right\rceil,$ (C.33)

with $\Delta_{\lambda}^{(u)}$ is defined in Eq. 22. Combining the results of
Eq. C.8 and Eq. C.33, we obtain

$m\geq\mathrm{max}\left(m_{\mathrm{min}}^{(1)},m_{\mathrm{min}}^{(2)}\right)\\\
=\left(\left\lceil\log_{2}\frac{\pi
s\norm{H}_{\text{max}}}{\varepsilon}\right\rceil,\left\lceil\log_{2}\frac{4\pi
s\norm{H}_{\text{max}}}{\delta\Delta_{\lambda}^{(u)}}\right\rceil\right),$
(C.34)

which justifies Eq. 23 in the main text.

####  C.1.3 Oracle queries for a single run of QPE

The number of qubits $m$ in the phase register determines the number of
oracles queries in the algorithm in Fig. 3. As we can see from Fig. 6, the
implementation of a controlled-$V$ operator requires to use once $U_{T}$ and
once $U_{T}^{\dagger}$. The circuit for $U_{T}$ is shown in Fig. 5 and uses
once the position oracle $O_{P}$, once the angle oracle $O_{\vartheta}$ and
once its inverse $O_{\vartheta}^{\dagger}$. Thus, for each $U_{T}$ we have $3$
oracle calls and a controlled-$V$ requires $6$ oracle calls. In order, to
perform a single run of QPE with $m$ qubit on the phase register we would need

$\sum_{k=0}^{m-1}2^{k}=2^{m}-1,$ (C.35)

controlled-$V$ operations, which yields the number of oracle queries

$N_{\mathrm{queries}}=6(2^{m}-1).$ (C.36)

Using the result of Eq. C.34 yields the scaling result in Eq. 25 in the main
text.

###  C.2 Sample size for given tolerances

We now turn our attention to the scaling of the number of samples $N_{S}$. Our
aim here is to perform a simplified analysis that should capture the
asymptotic behaviour of the scaling, while the details will eventually depend
on the procedure used to infer eigenvalues and weights from measured data in
the phase register (see Ref. [20] for instance). In particular, we work with
the assumption that the number of qubits $m$ in the phase register is taken
large enough to resolve all the peaks. From Eq. C.8 this means

$m\gg\left\lceil\log_{2}\frac{\pi s\lVert
H\rVert_{\mathrm{max}}}{\Delta_{\lambda}}\right\rceil.$ (C.37)

If this holds we can unambiguously identify the eigenvalue peaks and neglect
possible errors in the assignment of the samples to the eigenvalues. In this
regime, QPE de facto behaves as a measurement in the eigenbasis of $H$, where
each eigenvalue $\lambda_{j}$ appears with probability $p(j)=W_{uj}^{2}$
111111In our version of QPE we have two peaks associated with each eigenvalue,
but we can simply sum the probability of each peak to get $W_{uj}^{2}$. By
performing QPE many times, we obtain samples from which we can reconstruct the
probability distribution $W_{uj}^{2}$.

Let $\Lambda_{u}$ be the set of eigenvalues $\lambda_{j}$ that have “support”
on the oscillator $u$, i.e.,

$\Lambda_{u}=\\{\lambda_{j}\,|\,W_{uj}^{2}>0\\},$ (C.38)

and $N_{u}=|\Lambda_{u}|$ the number of elements of $\Lambda_{u}$ ($N_{u}\leq
N$). $W_{uj}^{2}$ is a probability distribution over the finite alphabet
$\Lambda_{u}$.

Let $p(\lambda)$ be a generic probability distribution over the finite
alphabet $\Lambda_{u}$. With each $\lambda_{j}\in\Lambda_{u}$ we can associate
a Bernoulli random variable $X_{\lambda_{j}}(\lambda)$ defined
$\forall\lambda\in\Lambda_{u}$ as

$X_{\lambda_{j}}(\lambda)=\begin{cases}1,\quad\lambda=\lambda_{j},\\\
0,\quad\lambda\neq\lambda_{j}.\end{cases}$ (C.39)

Notice that $\mathrm{E}[X_{\lambda_{j}}]=p(\lambda_{j})$. By drawing $N_{S}$
samples from $p(\lambda)$, we get samples $x_{\lambda_{j}}^{(i)}$,
$i\in\\{1,\dots,N_{S}\\}$ for $X_{\lambda_{j}}$. By taking the arithmetic
average of these samples we can construct the empirical estimator for
$p(\lambda_{j})$:

$\tilde{p}(\lambda_{j})=\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}x_{\lambda_{j}}^{(i)}.$
(C.40)

Our aim is to understand how many samples $N_{S}$ we have to draw from
$p(\lambda)$ to have, with probability $1-\zeta$, that each
$\tilde{p}(\lambda_{j})$, $\lambda_{j}\in\Lambda_{u}$ estimated
$p(\lambda_{j})$ with accuracy $\delta$. Mathematically, this means

$\mathrm{Prob}\left(\bigcup_{\lambda_{j}\in\Lambda_{u}}\\{|p(\lambda_{j})-\tilde{p}(\lambda_{j})|\geq\delta\\}\right)\leq\zeta.$
(C.41)

Thus, our aim is to obtain bounds for the left-hand side of Eq. C.41. First,
using the union bound we get

$\mathrm{Prob}\left(\bigcup_{\lambda_{j}\in\Lambda_{u}}\\{|p(\lambda_{j})-\tilde{p}(\lambda_{j})|\geq\delta\\}\right)\\\
\leq\sum_{\lambda_{j}\in\Lambda_{u}}\mathrm{Prob}\left(|p(\lambda_{j})-\tilde{p}(\lambda_{j})|\geq\delta\right).$
(C.42)

Now, an application of Hoeffding’s inequality [84, 85] to each random variable
$X_{\lambda_{j}}$ yields

$\mathrm{Prob}\left(|p(\lambda_{j})-\tilde{p}(\lambda_{j})|\geq\delta\right)\\\
=\mathrm{Prob}\left(\left\lvert\mathrm{E}[X_{\lambda_{j}}]-\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}x_{\lambda_{j}}^{(i)}\right\rvert\geq\delta\right)\leq
2e^{-2N_{S}\delta^{2}}.$ (C.43)

Plugging Eq. C.43 into Eq. C.42 and enforcing Eq. C.41 we get

$\mathrm{Prob}\left(\bigcup_{\lambda_{j}\in\Lambda_{u}}\\{|p(\lambda_{j})-\tilde{p}(\lambda_{j})|\geq\delta\\}\right)\\\
\leq 2N_{u}e^{-2N_{S}\delta^{2}}\leq\zeta,$ (C.44)

which holds if

$N_{S}\geq\frac{1}{2\delta^{2}}\ln\left(\frac{2N_{u}}{\zeta}\right).$ (C.45)

Eq. C.45 justifies Eq. 26 in the main text.

##  Appendix D Modification for non-local response

In this appendix, we describe how to modify our algorithm to compute non-local
entries $G_{uv}(s)$ of the response function whose analytical form is given in
Eq. A.1.

[wire types= q,q,q,q,b,b,classical gap=0.07cm] [wires=4, label style=anchor=
east, xshift=-0cm]$\text{Phase:}\ket{0}^{\otimes m}$ & H 5 … [4,disable auto
height]QFT^†
H 4 …

H … 2
[label style=anchor= east, xshift=-0.2cm]$\text{Ancilla:}\ket{0}^{\otimes l}$
[2, disable auto height]V[2, disable auto height]V^2 … [2,disable auto
height]V^2^m-1
[label style=anchor= east, xshift = -0.2cm ]$\text{State:}\ket{0}^{\otimes n}$
∏_i=0^n-1X_i^u_i ∏_i=0^n-1X_i^v_i⊕u_i …
[label style=anchor= east, xshift=-0.2cm]$\text{Hadamard:}\ket{0}$ H-1H

Figure 7:  Modified circuit for the computation of the non-local response
function. One can see two additions compared to Fig. 3. The first is an
additional qubit and the second is a modified Hadamard test which we use to
prepare the state $\ket{u}$ next to $\ket{v}$.

For the local contributions the preparation of the computational basis state
$\ket{u}$ is crucial to extract the weights $W_{uj}^{2}$. The disadvantage of
this method is that we can only extract the absolute square of those weights,
which does not contain the relative sign and makes this method not suited to
evaluate the terms $W_{uj}W_{vj}$ that appear in Eq. A.1.

A way around this is to add another qubit and still prepare $\ket{u}$ in the
state register as shown in Fig. 7. We will now use the modified Hadamard test
to “add” $\ket{v}$ to this state. At first we apply a Hadamard gate to the
additional ancilla qubit starting in $\ket{0}$ (bottom qubit in Fig. 7). Thus,
we have

$\xrightarrow{H}\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})\ket{u}.$ (D.1)

We now apply a series of CNOT gates with the ancilla qubit as control and the
qubits in the $n$-qubit register as target. The CNOT will be only applied to
target qubit $i$ for which $v_{i}\oplus u_{i}=1$. This brings the state
$\ket{u}$ to $\ket{v}$ when the control qubit is in $\ket{1}$. Accordingly,
after this series of CNOTs we get

$\xrightarrow{\mathrm{CNOTs}}\frac{1}{\sqrt{2}}(\ket{0}\ket{u}+\ket{1}\ket{v}).$
(D.2)

At last we apply another Hadamard gate to the ancilla qubit and obtain

$\xrightarrow{H}\frac{1}{2}\left(\ket{0}\left(\ket{u}+\ket{v}\right)+\ket{1}\left(\ket{u}-\ket{v}\right)\right).$
(D.3)

Measuring the additional qubit results in either the sum or difference of
those two states

$\frac{1}{\sqrt{2}}\left(\ket{u}\pm\ket{v}\right)=\frac{1}{\sqrt{2}}\sum_{j=1}^{N}\left(W_{uj}\pm
W_{vj}\right)\ket{\lambda_{j}}.$ (D.4)

If we proceed now as before with the local version of the algorithm, depending
on the ancilla being in $0$ or $1$, we would extract the square of those
amplitudes

$\frac{1}{2}(W_{uj}\pm
W_{vj})^{2}=\frac{1}{2}W_{uj}^{2}+\frac{1}{2}W_{vj}^{2}\pm W_{uj}W_{vj},$
(D.5)

where we used that $W_{uj}$ and $W_{vj}$ are real. Let $P(0,x)$, $P(1,x)$ the
probability of measuring the bitstring $x$ in the phase register given that
the ancilla is measured in $0$ and $1$, respectively. By adapting Eq. 47 with
the new amplitudes we find

$\frac{1}{2}(W_{vj}+W_{uj})^{2}\approx 2\sum_{q=-Q}^{Q-1}P(0,\lceil
y_{j}^{\prime}\rceil+q),$ (D.6a) $\frac{1}{2}(W_{vj}-W_{uj})^{2}\approx
2\sum_{q=-Q}^{Q-1}P(1,\lceil y_{j}^{\prime}\rceil+q),$ (D.6b)

or, equivalently,

$W_{vj}W_{uj}\approx\sum_{q=-Q}^{Q-1}\left(P(0,\lceil
y_{j}^{\prime}\rceil+q)-P(1,\lceil y_{j}^{\prime}\rceil+q)\right).$ (D.7)

##  Appendix E Details of the glued-trees problem

In this appendix, we provide the necessary details of the glued-trees problem.
We consider a random glued-tree graph $G=(V,E)$ shown in Fig. 4 and described
in Sec. IV. We work directly with the quantum formalism so that, with each
vertex $u\in V$, we associate an element in the computational basis $\ket{u}$.

First, we note that, by definition, the action of an adjacency matrix $A$ on a
vertex $\ket{u}$ reads

$A\ket{u}=\sum_{v\in N_{G}[u]\setminus\\{u\\}}\ket{v}.$ (E.1)

As we see from Fig. 4 the vertices in the glued-trees graph can be organized
in $2n_{c}$ columns. Let us denote by $\mathrm{col}(j)$ the set of vertices in
the $j$-th columns starting from the left of the random glued tree graph with
$j\in\\{1,\dots,2n_{c}\\}$. We define the “column” states

$\ket{\mathrm{col}(j)}=\frac{1}{\sqrt{2^{j-1}}}\sum_{v\in\mathrm{col}(j)}\ket{v},\quad
j\in\\{1,\dots,n_{c}\\},$ (E.2a)
$\ket{\mathrm{col}(j)}=\frac{1}{\sqrt{2^{2n_{c}-j}}}\sum_{v\in\mathrm{col}(j)}\ket{v},\quad
j\in\\{n_{c}+1,\dots,2n_{c}\\}.$ (E.2b)

Notice that ENTRANCE and EXIT vertices correspond to $\ket{\mathrm{col}(1)}$
and $\ket{\mathrm{col}(2n_{c})}$, respectively. We denote by
$\mathcal{H}_{\mathrm{col}}$ the column subspace:

$\mathcal{H}_{\mathrm{col}}=\mathrm{span}\\{\ket{\mathrm{col}(j)}:j\in\\{1,\dots,2n_{c}\\}\\}.$
(E.3)

The projector onto $\mathcal{H}_{\mathrm{col}}$ is

$\Pi_{\mathrm{col}}=\sum_{j=1}^{2n_{c}}\outerproduct{\mathrm{col}(j)}{\mathrm{col}(j)}.$
(E.4)

The crucial fact is that if the system starts in the column subspace, an
application of the adjacency matrix $A$ will give a state in the column
subspace. In fact, the action of the adjacency matrix $A$ of the random glued-
trees graph on the column states reads

$A\ket{\mathrm{col}(j)}=\sqrt{2}\ket{\mathrm{col}(j+1)}+\sqrt{2}\ket{\mathrm{col}(j-1)},$
(E.5)

for $j\in\\{2,\dots,n_{c}-1\\}\cup\\{n_{c}+1,\dots 2n_{c}-1\\}$, and

$A\ket{\mathrm{col}(1)}=\sqrt{2}\ket{\mathrm{col}(2)},$ (E.6a)
$A\ket{\mathrm{col}(2n_{c})}=\sqrt{2}\ket{\mathrm{col}(2n_{c}-1)},$ (E.6b)
$A\ket{\mathrm{col}(n_{c})}=\sqrt{2}\ket{\mathrm{col}(n_{c}-1)}+2\ket{\mathrm{col}(n_{c}+1)},$
(E.6c)
$A\ket{\mathrm{col}(n_{c}+1)}=2\ket{\mathrm{col}(n_{c})}+\sqrt{2}\ket{\mathrm{col}(n_{c}+2)}.$
(E.6d)

For notational simplicity we will denote the column states
$\ket{\mathrm{col}(j)}$ as $\ket{j}$ from now on:

$\ket{\mathrm{col}(j)}\mapsto\ket{j}.$ (E.7)

Since in 1 the system starts at the ENTRANCE, i.e., in
$\ket{\mathrm{ENTRANCE}}=\ket{j=1}$, we can effectively restrict to the column
subspace and consider the projected adjacency matrix

$\frac{A_{\mathrm{col}}}{\sqrt{2}}=\frac{\Pi_{\mathrm{col}}A\Pi_{\mathrm{col}}}{\sqrt{2}}=\sum_{j=1}^{n_{c}-1}(\outerproduct{j}{j+1}+\mathrm{h.c.})\\\
+\sqrt{2}(\outerproduct{n_{c}}{n_{c}+1}+\mathrm{h.c.})\\\
+\sum_{j=n_{c}+1}^{2n_{c}-1}(\outerproduct{j}{j+1}+\mathrm{h.c.}),$ (E.8)

as the effective Hamiltonian of the system. Eq. (E.8) is the Hamiltonian of a
quantum walk on one line with $2n_{c}$ sites with a “defect” in the middle.
Note that from a matrix point of view $A_{\mathrm{col}}$ is a tridiagonal
symmetric matrix.

As shown in Ref. [26], $A_{\mathrm{col}}/\sqrt{2}$ has eigenvalues of the form

$\lambda_{z}=2\cos z,$ (E.9)

where $z=p+iq\in\mathbb{C}$, $p,q\in\mathbb{R}$, satisfies the quantization
condition

$\frac{\sin(z(n_{c}+1))}{\sin(zn_{c})}=\pm\sqrt{2}.$ (E.10)

The corresponding eigenvector is of the form

$\ket{\lambda_{z}}=\frac{1}{\sqrt{\mathcal{N}(z,n_{c})}}\biggl{(}\sum_{j=1}^{n_{c}}\sin(zj)\ket{j}\\\
\pm\sum_{j=n_{c}+1}^{2n_{c}}\sin(z(2n_{c}+1-j))\ket{j}\biggr{)},$ (E.11)

with $\mathcal{N}(z,n_{c})$ a normalization factor.

Figure 8: Quantization condition Eq. E.10 for $n_{c}=8$.

We denote by $\mathcal{Z}$ the set of solutions of Eq. E.10 that are
associated with a distinct eigenpair. Thus, $|\mathcal{Z}|=2n_{c}$. The
quantization condition Eq. E.10 has exactly $2n_{c}-2$ real solutions $z=p$
(see Fig. 8 for an example). We denote by $\mathcal{P}$ the set of these
solutions. The two remaining eigenpairs can be obtained by setting $z=iq$ and
$z=\pi+iq$, but we do not need them to show our result.

In terms of the eigenvectors $\ket{\lambda_{z}}$, the
$\ket{\mathrm{ENTRANCE}}$ state reads

$\ket{\mathrm{ENTRANCE}}=\sum_{z\in\mathcal{Z}}\frac{\sin
z}{\sqrt{\mathcal{N}(z,n_{c})}}\ket{\lambda_{z}}.$ (E.12)

We now assume that the number of qubits in the phase register is taken to
satisfy Eq. 49 in the main text. With this assumption, the probability of
obtaining the eigenstate $\ket{\lambda_{z}}$ after applying QPE to
$\ket{\mathrm{ENTRANCE}}$, as described in Sec. III, and measurement on the
ancilla register with post-selection on the $\ket{0}^{\otimes a}$ state, is

$\mathrm{Prob}(\lambda_{z}|\mathrm{ENTRANCE})=\frac{|\sin
z|^{2}}{2\mathcal{N}(z,n_{c})}.$ (E.13)

Once we obtained one of the $\ket{\lambda_{z}}$, the probability of finding
the EXIT after measurement in the computational basis is

$\mathrm{Prob}(\mathrm{EXIT}|\lambda_{z})=\frac{|\sin
z|^{2}}{\mathcal{N}(z,n_{c})},$ (E.14)

Thus, the total probability of finding the EXIT is

$\mathrm{Prob}(\mathrm{EXIT})=\sum_{z\in\mathcal{Z}}\frac{|\sin
z|^{4}}{2\mathcal{N}(z,n_{c})^{2}}.$ (E.15)

We now want to study analytically the scaling of this probability
$\mathrm{Prob}(\mathrm{EXIT})$. First of all, we restrict the sum over the
real solutions $z=p$ of Eq. E.10. In this case, the normalization factor
$\mathcal{N}(p,n_{c})$ reads

$\mathcal{N}(p,n_{c})=\sum_{j=1}^{n_{c}}\sin^{2}(pj)+\sum_{j=n_{c}+1}^{2n_{c}}\sin^{2}(p(2n_{c}+1-j))\\\
=\frac{1}{2}\biggl{[}1+2n_{c}-\frac{\sin(p(1+2n_{c}))}{\sin
p}\biggr{]}<2n_{c}.$ (E.16)

Thus, we get

$\mathrm{Prob}(\mathrm{EXIT})>\sum_{p\in\mathcal{P}}\frac{\sin^{4}p}{2\mathcal{N}(p,n_{c})^{2}}>\frac{1}{8n_{c}^{2}}\sum_{p\in\mathcal{P}}\sin^{4}p.$
(E.17)

As discussed in Ref. [26] the real solutions concentrate around
$\pi\ell/n_{c}$ with $\ell\in\\{1,\dots,n_{c}-1\\}$, with corrections of order
$\mathcal{O}(1/n_{c}^{2})$. This is also visible in the example shown in Fig.
8 for $n_{c}=8$. In particular, we see that there are two solutions around
$\pi\ell/n_{c}$. Thus, for large $n_{c}$, we can write

$\sum_{p\in\mathcal{P}}\sin^{4}p=2\sum_{\ell=1}^{n_{c}-1}\sin^{4}\left(\frac{\pi\ell}{n_{c}}\right)+\mathcal{O}(1)\\\
=\frac{3n_{c}}{4}+\mathcal{O}(1).$ (E.18)

We finally obtain

$\mathrm{Prob}(\mathrm{EXIT})>\frac{3}{32n_{c}}+\mathcal{O}\left(\frac{1}{n_{c}^{2}}\right),$
(E.19)

which justifies Eq. 51 in the main text.