Exponential quantum speedup in simulating coupled classical oscillators

We present a quantum algorithm for simulating the classical dynamics of $2^n$ coupled oscillators (e.g., $2^n$ masses coupled by springs). Our approach leverages a mapping between the Schr\"odinger equation and Newton's equation for harmonic potentials such that the amplitudes of the evolved quantum state encode the momenta and displacements of the classical oscillators. When individual masses and spring constants can be efficiently queried, and when the initial state can be efficiently prepared, the complexity of our quantum algorithm is polynomial in $n$, almost linear in the evolution time, and sublinear in the sparsity. As an example application, we apply our quantum algorithm to efficiently estimate the kinetic energy of an oscillator at any time. We show that any classical algorithm solving this same problem is inefficient and must make $2^{\Omega(n)}$ queries to the oracle and, when the oracles are instantiated by efficient quantum circuits, the problem is BQP-complete. Thus, our approach solves a potentially practical application with an exponential speedup over classical computers. Finally, we show that under similar conditions our approach can efficiently simulate more general classical harmonic systems with $2^n$ modes.


I. INTRODUCTION
The efficient simulation of quantum dynamics is among the most promising areas of quantum computing [1].This is due to having practical applications as well as universality providing strong complexity theoretic evidence for exponential quantum advantage [2,3].Several other applications with similarly favorable qualities are known (e.g., factoring [4]) but the discovery of more compelling use cases remains critical for understanding and motivating the value proposition of quantum computers.
One approach for leveraging the power of Hamiltonian simulation would be to consider the space of problems that reduces to it.Hamiltonian dynamics is an example of a homogeneous, first order partial differential equation that gives rise to unitary evolutions.Thus, it is natural to explore whether other differential equations can be mapped to Hamiltonian simulation.The goal would be a more natural (albeit perhaps more narrow) alternative to solving differential equations compared with methods [5][6][7] leveraging quantum linear systems algorithms [8][9][10].Past work has sought to develop Hamiltonian simulation approaches for a limited set of other differential equations but thus far, has only succeeded in obtaining polynomial speedups [11][12][13].
In this paper, we discuss a rich set of problems in classical dynamics that can be mapped to Hamiltonian simulation and solved with exponential quantum advantage.As a prominent example, our approach can simulate the dynamics of exponentially many coupled classical oscillators in polynomial time.Such systems describe a variety of physical phenomena including: networks of masses and springs [14], circuits with capacitors and inductors [15], models of neuron activity [16], and vibrations in molecules [17], materials, and mechanical structures.More generally, the harmonic approximation (defined by a quadratic potential) arises as the first order correction to equilibrium in bound systems.
That the dynamics of coupled classical oscillators can be studied via Hamiltonian simulation is perhaps not too surprising.Solutions to the Schrödinger equation contain oscillatory terms and interference, and we are effectively using these properties to simulate the same phenomena in the classical system.This correspondence was also presaged by the finding that Grover search can be implemented (in the absence of errors) using mechanical wave interference [18].Similarly, the problem of simulating the (discrete) classical wave equation is a special case of the oscillator dynamics problem considered here, and a quantum algorithm to solve the classical wave equation via Hamiltonian simulation has been studied in Ref. [11].However, despite the seemingly natural connection, our mapping involves a number of subtle technical features and modern quantum algorithms are required to efficiently simulate the resultant Hamiltonians.
Besides providing a classical-to-quantum reduction, another key to obtain exponential quantum advantage for this problem is the way we encode physical quantities in quantum states.In our case, this encoding is motivated by energy conservation (being different to the one in Ref. [11]).Specifically, we encode quantities related to the displacements and momenta of the classical system in the amplitudes of a quantum state.Thus, extracting the full configuration of the classical system would scale polynomially in the Hilbert space dimension, precluding a large quantum speedup.Nevertheless, interesting global properties, like estimating the kinetic or potential energies at any time, can still be computed efficiently.Likewise, our methods only provide an exponential speedup for the dynamics of sparsely coupled oscillator networks when masses, spring constants, and initial states can be efficiently computed.Fortunately, many interesting systems meet those conditions.
Finally, we provide strong evidence that such simulations cannot also be performed efficiently on a classical computer.In particular, we are able to show that estimating the kinetic energy of simple instances of these systems at any time that is poly(n) requires 2 Ω(n) queries to the oracles in the worst case and, when the oracles are instantiated by quantum circuits of poly(n) size, the problem is BQP-complete.Thus, all problems efficiently solved by a quantum computer can be reduced to an instance of simulating exponentially many coupled classical oscillators, which can also be solved efficiently by our approach.
The rest of the paper is organized as follows.In Sec.II we describe classical systems of oscillators, define the associated simulation problems of interest, and state our main results.In Sec.III we show how these problems reduce to instances of Hamiltonian simulation, giving rise to an efficient quantum algorithm.We discuss some applications of this algorithm in Sec.V, making emphasis on the problem of estimating the time-dependent kinetic (or potential) energies of the oscillators.In Sec.VI A we show the exponential lower bound for classical algorithms for this problem in the oracle setting, and in Sec.VI B we show that this problem is BQP-complete when the oracles can be accessed via efficient quantum circuits.Finally, in Sec.VII we generalize our approach to efficiently simulating classical systems under the harmonic approximation.Detailed proofs of our claims are provided in the Appendix.We also provide a comparison of our results with those of related work on quantum algorithms for differential equations in Appendix G. V S E 1 7 5 E y 6 T 1 K B k i 0 V h K o i J y e x v 0 u c K m R F j S y h T 3 N 5 K 2 J A q y o x N p 2 B D 8 J Z f X i X N S t m 7 K F f v q q X a d R Z H H k 7 g F M 7 B g 0 u o w S 3 U o Q E M B v A M r / D m C O f F e X c + F q 0 5 J 5 s 5 h j 9 w P n 8 A E C K N q g = = < / l a t e x i t > w H P 8 A p v j n J e n H f n Y 9 a 6 4 s x n j u A P n M 8 f 1 j 6 O t w = = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " l N V + p 7 x Z e X S e u i 4 l 1 W q g / V c u 0 2 j 6 M A x 3 A C Z + D B F d T g H h r Q B A I S n u E V 3 h z j v D j v z s e 8 d c X J Z 4 7 g D 5 z P H 3 I I k L g = < / l a t e x i t > m N 1 = m N < l a t e x i t s h a 1 _ b a s e 6 4 = " E 7 9 S Q e b 6 7 1 8 U f p O s y V a V 4 f o j y a M = " > A A A B 8 n i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I v g q e y K q M e i F 4 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " q 2 1 + s s b

II. SIMULATING COUPLED OSCILLATORS: MAIN RESULTS
We consider a classical system of coupled harmonic oscillators, i.e., N = 2 n point (positive) masses m 1 , . . ., m N that are coupled with each other by springs.At any time t ≥ 0, the displacements (with respect to their rest position) and velocities of the masses are given by ⃗ x(t) = (x 1 (t), . . ., x N (t)) T ∈ R N and ⃗ x(t) = ( ẋ1 (t), . . ., ẋN (t)) T ∈ R N , respectively, where ȧ(t) = d dt a(t) and ä(t) = d 2 dt 2 a(t).We let κ jk = κ kj ≥ 0 be the spring constants that couple the j th and k th oscillator, and κ jj ≥ 0 is the spring constant that connects the j th oscillator to a "wall".We have described it for dimension D = 1 for simplicity (i.e., we assigned one spatial coordinate to each oscillator), but the same formulation holds for coupled harmonic oscillators in D spatial dimensions, for arbitrary D, by using D coordinates to represent the position of a single oscillator (see Fig. 1 for an example).
In the harmonic approximation, the dynamics of the oscillators can be determined from the initial values ⃗ x(0) and ⃗ x(0) and Newton's equation (for all j ∈ [N ] := {1, . . ., N }): , where M is a N × N diagonal matrix with entries m j > 0 and F is the N × N matrix whose diagonal and offdiagonal entries are f jj = k κ jk and f jk = −κ jk , respectively.(Observe that F is the discrete Laplacian of a weighted graph.)Solutions to Eq. ( 1) are well understood and can be expressed in terms of normal modes [14], which is essentially a way of describing the system as N uncoupled harmonic oscillators in a different basis.Classical algorithms that compute ⃗ x(t) and ⃗ x(t) have complexity poly(N ) or, equivalently, exp(n).(In this paper we use exp(m), poly(m) and polylog(m) to mean O(k m ), O(m k ), and O(log k m) respectively, for some constant k. (We use "log" to mean logarithm to base 2.) Our goal is to provide a quantum algorithm for simulating the dynamics of the classical oscillators efficiently, in time poly(n).Doing so requires a particular notion of "simulation", since any algorithm that outputs the full vectors ⃗ x(t) or ⃗ x(t) would necessarily have complexity at least linear in N .Specifically, we consider a problem formulation where the output is a quantum state with some amplitudes proportional to √ m j ẋj (t) and others to or √ κ jj x j (t), and the input masses and spring constants are provided through oracles in the usual way; see App. A. (Here and throughout this paper, √ x is the principal square root of x ≥ 0 and √ X is the principal square root of a positive semidefinite matrix X ≽ 0.) The formal problem is as follows: Problem 1.Let K be the N × N symmetric matrix of spring constants κ jk ≥ 0 and assume it is d-sparse (i.e., there are at most d non-zero entries in each row).Let M be the N × N diagonal matrix of masses m j > 0 and define the normalized state where E > 0 is a constant, and ⃗ µ(t) ∈ R M (M := N (N + 1)/2) is a vector with N entries √ κ jj x j (t) and N (N − 1)/2 entries √ κ jk (x j (t) − x k (t)), with k > j.Assume we are given oracle access to K and M, and oracle access to a unitary W that prepares the initial state, i.e., W |0⟩ → |ψ(0)⟩.Given t ≥ 0 and ϵ > 0, the goal is to output a state that is ϵ-close to |ψ(t)⟩ in Euclidean norm.
Our main result is a quantum algorithm that prepares |ψ(t)⟩ efficiently: Theorem 1. Problem 1 can be solved with a quantum algorithm that makes Q = O(τ + log(1/ϵ)) queries to the oracles for K and M, uses two-qubit gates, and uses W once, where τ := t √ ℵd ≥ 1, ℵ := κ max /m min , m max ≥ m j ≥ m min > 0 and κ max ≥ κ jk for all j, k ∈ [N ] are known quantities.
In the O notation above, the asymptotically large parameters are N , τ , 1/ϵ, and m max /m min .This complexity has explicit dependence on n = log(N ) that is polynomial.If τ , log(1/ϵ), and log(m max /m min ) are poly(n), and all oracles (including W) can be performed with cost poly(n), then the entire algorithm has complexity poly(n).An interesting feature of this complexity is that it scales as the square root of the sparsity d, whereas Hamiltonian simulation complexities are usually linear in d.The complexity is also linear in τ , which bounds the maximum number of oscillations of the normal modes after time t.
Our algorithm can be used to determine properties of the system at time t.For example, the state |ψ(t)⟩ encodes the velocities and displacements of the oscillators in a way that makes it easy for the estimation of the kinetic or potential energies, as we discuss below.We discuss in App.F that other encodings are also possible, and the choice of encoding may determine which initial states and properties (observables) are efficient to prepare and measure.Nevertheless, our main goal is to establish results for this particular encoding.
The constant in Eq. ( 2) is E = K(t) + U (t), where and U (t) = 1 2 ⃗ µ(t) T ⃗ µ(t) is the potential energy, i.e., ( Hence, E is the total energy of the system, which is time-independent as the system is closed. where K V (t) := 1 uses of the quantum algorithm from Thm. 1.We also show a related result for the estimation of the potential energy stored on a subset V ⊆ [N ] × [N ] of springs at time t.Problems 1 and 2 are formulated using a specific input and output format, and one might wonder if a classical algorithm could also solve these problems efficiently.The answer is no: we show that any classical algorithm solving Problem 2 must make poly(N ) or exp(n) queries to the oracles in general.
To provide more evidence of the impossibility of efficient classical simulations and prove a stronger hardness result, we define a non-oracular version of Problem 2 and show that it is BQP-complete.Thus, no classical algorithm can solve this non-oracular problem efficiently (in time poly(n)) unless it can also solve every problem solved by quantum computers in polynomial time.The same observation essentially applies to Problem 2 as a result.Formally, to be BQPcomplete a problem must be a decision problem (i.e., a problem where there are only two possible answers), so we define a decision version of Problem 2 with |V | = 1 and a sparse initial state: Problem 3. The setup is as in Problem 2, but we are given efficient quantum circuits (i.e., circuits with poly(n) gates) to implement the oracles for K and M, an explicit description of |ψ(0)⟩, which is required to have a constant number of nonzero entries, and a label v ∈ [N ] of a single oscillator.We additionally require that τ , 1/ϵ, and m max /m min are bounded by poly(n).The problem is to decide if K v (t)/E, the kinetic energy of oscillator v as a fraction of total energy, is at least promised that one of these holds.
Finally, we show that our results can be extended to address more general classical systems that arise naturally under the harmonic approximation [14].In these cases, for example, the matrix M resulting from Eq. 1 is not necessarily diagonal and the off-diagonal entries of F are not necessarily non-positive.A simple change of variables takes Eq. 1 to the standard form ⃗ y(t) = −A⃗ y(t), where ⃗ y(t) = √ M⃗ q(t), and ⃗ q(t) ∈ R N are the generalized displacements.We then consider the following generalization of Problem 1: where E > 0 is a constant and ⃗ µ(t) := √ A⃗ y(t) ∈ C N .Assume we are given oracle access to A and oracle access to a unitary W that prepares the initial state |ψ(0)⟩.Given t and ϵ, the goal is to output a state that is ϵ-close to |ψ(t)⟩ in Euclidean norm.
The constant E in Eq. ( 7) is also the energy of the system.The encoding is different from the one used for Problem 1; for example √ A is of dimension N ×N and might not be sparse even if A is.However, the support of |ψ(t)⟩ on the subspace spanned by the first N basis states is still K(t)/E, where K(t) is the kinetic energy.We then show the following: Theorem 4. Problem 4 can be solved with a quantum algorithm that makes ) queries to the oracles, uses O(Q×polylog(N/ϵ)) twoqubit gates, and uses W once.
The Õ notation hides logarithmic factors in several parameters that specify the input.

III. REDUCTION TO QUANTUM EVOLUTION
We show how to reduce Problem 1 to time evolution of a quantum system.A change of variables where ⃗ y(t) := √ M⃗ x(t) allows us to write Eq. ( 1) as where ≽ 0 is PSD and realsymmetric, and M ≻ 0 is the diagonal matrix with entries m j .Any solution to Eq. ( 9) satisfies This is simply Schrödinger's equation induced by the Hamiltonian − √ A. Hence its solution is Unfortunately, we do not have direct access to − √ A and casting the problem as quantum evolution with simpler Hamiltonians requires additional steps.It is possible, however, to do Hamiltonian simulation with − √ A given oracle access to A using quantum phase estimation, as we describe in Sec.VII, but that is less efficient than what we describe here.The phase estimation approach also gives rise to a different encoding, and some properties of the system are not as easy to access.For example, in the encoding used in Problem 1, the support of |ψ(t)⟩ in basis states are either proportional to the kinetic energies or potential energies of the springs, but this might not be the case in other encodings.
Let B be any N ×M matrix that satisfies BB † = A and H be the block Hamiltonian: where 0 are matrices of all zeroes.(The dimension of each 0 is clear from context.)Note that H acts on the space C N +M and functions as the "square root" of A, since the first block of H 2 is A. Schrödinger's equation induced by H is where |ψ(t)⟩ ∈ C N +M is the state of a quantum system at time t.
It can be verified by direct substitution that is a valid solution to Eq. ( 13).Hence we have which lets us compute ⃗ y(t) and B † ⃗ y(t) using Hamiltonian simulation.This generalizes Eq. ( 11) since B is not necessarily − √ A. (This formulation also stores the positions and velocities in different components of the vector, but this is a minor difference.) To match the state representation of Problem 1, we need to choose B such that ⃗ µ }, which is of size M .Since F is a graph Laplacian, we can obtain B from the incidence matrix of F given by This choice satisfies √ MB( √ MB) † = F, and hence BB † = A, because we can write F = √ M( j≤k B|j, k⟩⟨j, k|B † ) √ M and each term in this sum describes the interaction between the j th and k th oscillators: and These reproduce F once we perform the sum.To obtain the action of B we apply √ M −1 to Eq. ( 16).
Note that a similar approach based on the incidence matrix was taken in Ref. [11] to provide a quantum algorithm that simulates the wave equation.The setting considered in that paper corresponds to the special case of our problem where the graph is spatially local, and the masses and spring constants are uniform.

IV. QUANTUM ALGORITHM
Our quantum algorithm solves Problem 1 by preparing |ψ(0)⟩ and simulating H for time t.
Access model: The Hamiltonian H is built from the matrix of spring constants K (or F) and M. As we wish to avoid complexities that are polynomial in N , we require succinct representations of these matrices.For maximum generality, we assume that access to M and the d-sparse K is provided by a black box, similar to that used in prior quantum algorithms [19,20].The black box is a unitary S that computes the masses m j on input j and the nonzero entries of K, i.e. κ jk , on input (j, k), as well as their locations.This black box readily gives query access to H in Eq. (12), which is also d-sparse.Its entries are ± κ jk /m j and these factors can be applied using the values of m j and κ jk ; see App. A.
Simulation of H: Time evolution with the sparse H can be simulated using one of the methods in Refs.[21][22][23][24], which apply an approximation of the exponential e −itH .Such methods are tailored to our access model and achieve almost optimal scaling in the parameters ϵ, d, t, and ∥H∥ max , the largest entry of H in absolute value.Here ∥H∥ max is at most √ ℵ, which is defined in Thm. 1.
Complexity: The complexity of our approach is determined by the simulation of H for time t from the initial state.We consider the query complexity Q of the number of calls to S, and we allow inverses and controlled forms of these oracles.This query complexity can also be taken to include the unitary for preparing the initial state, but only one call to that preparation is needed.We also consider the gate complexity G, which is the number of additional elementary gates, which could, for example, be singlequbit gates and CNOTs.We describe these as "twoqubit gates" in Thm. 1 for simplicity.
Using Ref. [23], the Hamiltonian evolution may be simulated with calls to a block encoding of H, where the block encoding gives a factor of 1/Λ.In App.A we show how to block encode H with Λ = √ 2ℵd using O(1) calls to S, and so Eq. ( 19) gives the total query complexity Q given in Thm. 1 since τ := t √ ℵd.The reason for the square root dependence on the sparsity is then because the algorithm involves simulation of H, which is effectively the square root of the d-sparse operator A. Note that ℵd is an upper bound on ∥A∥, implying that √ ℵd is an upper bound on the largest frequency of the normal modes.See App.A for a more technical explanation.
To achieve final error ϵ, each block encoding should be given within error which will determine the gate complexity.In Lemmas 8 and 9 of App.A we also show that it suffices to represent the relevant inputs, i.e., masses and spring constants, with a number of bits that is , respectively, to achieve error ϵ ′ in the block encoding.
The gate complexity for the arithmetic is essentially O(r m r κ ), and there is also another contribution O(n) to the gate complexity for other operations in the block encoding (e.g, for inequality tests and other state preparations).Hence, we can bound the overall gate complexity as Substituting the value of ϵ ′ from Eq. ( 20), this gives which is the expression provided in Thm. 1.To simplify this expression we used the fact that log((τ + log(1/ϵ))/ϵ) ≤ log((τ + (1/ϵ))/ϵ).Because we consider complexity for large τ and 1/ϵ, we can upper bound that by log τ /ϵ 2 = O(log(τ /ϵ)).
In Thm. 1 we have not quantified the state preparation complexity (i.e., the complexity of W), since Problem 1 assumes W is given as an input.However, in many situations this preparation can often be performed efficiently.One example is when we have oracles to separately prepare states with amplitudes proportional to ⃗ x(0) and ⃗ x(0).In App.E, we show the query complexity to prepare |ψ(0)⟩ is , where is the energy of N uncoupled oscillators of mass m max and spring constants κ max with similar initial conditions.We also establish the gate complexity

V. EXAMPLE APPLICATION
Our quantum algorithm can simulate classical oscillators in time polynomial in n under some assumptions.As we are concerned with dynamics and seek to preserve the speedup, we focus on computing certain time-dependent properties of the system.In particular, we explain how to use the quantum algorithm to solve Problem 2, i.e., to estimate the kinetic energy of all or some oscillators V ⊆ [N ] at any time.
Theorem 5. Let V be an oracle for V , i.e., a unitary that performs the map V |j⟩ = − |j⟩ if j ∈ V and V |j⟩ = |j⟩ otherwise, and δ > 0. Then Problem 2 can be solved with probability at least 1 − δ by a quantum algorithm that makes O(log(1/δ)/ϵ) uses of V and the quantum circuit that prepares |ψ(t)⟩ (along with its inverse and controlled version).
, and a solution to Problem 2 provides an estimate of U (t)/E as well.As usual, we include controlled applications of the oracle V, which means that the global phase is distinguishable (this way we have that Although here we focus on estimating kinetic energies, a similar idea can be used to estimate the potential energies stored on a subset of springs at any time, as we discuss below. To prove Thm. 5, we note that K V (t)/E = ⟨ψ(t)|P V |ψ(t)⟩, where P V = (1l − V)/2 is the projector onto V .This expectation, which is the support of |ψ(t)⟩ on the subspace spanned by the basis states whose labels correspond to V (i.e., a subspace of that spanned by the first N basis states), can be obtained by making measurements on many copies of |ψ(t)⟩, but that approach is not optimal (i.e., the scaling is quadratic in 1/ϵ).The problem reduces to estimating ⟨ψ(t)| V |ψ(t)⟩ with additive error at most 2ϵ and error probability δ.A method based on high-confidence amplitude estimation [25] then provides the result in Thm. 5.
Hence, when |ψ(t)⟩ can be efficiently prepared and V can be efficiently implemented, we can obtain an estimate of K V (t)/E with error ϵ and high probability in time polynomial in n.If in addition E is known or is efficiently computed, this translates to an efficient estimation of K V (t) with error ϵE.
A related result shows that our quantum algorithm can be applied to estimate the potential energy of a subset of oscillators: Problem 5. Given the same inputs as Problem 1 and an oracle for where is the potential energy of the springs in V at time t.
By noting that U V (t)/E is the support of |ψ(t)⟩ on the subspace spanned by the basis states |j, k⟩ such that (j, k) ∈ V , the problem reduces to estimating the expectation of another unitary V on |ψ(t)⟩.Like in the prior example for the kinetic energy, this can be done efficiently using our quantum algorithm: Theorem 6.Let V be an oracle for V , i.e., a unitary that performs the map V |j, k⟩ = − |j, k⟩ if (j, k) ∈ V and V |j, k⟩ = |j, k⟩ otherwise, and δ > 0.Then, Problem 5 can be solved with probability at least 1 − δ by a quantum algorithm that makes O(log(1/δ)/ϵ) uses of W and the quantum circuit that prepares |ψ(t)⟩ (along with its inverse and controlled version).

VI. IMPOSSIBILITY OF EFFICIENT CLASSICAL SIMULATIONS
An important question is whether our quantum algorithm results in an exponential quantum speedup, or whether it can be efficiently simulated by classical algorithms.We address this question in two ways: (i) by showing that our approach can solve an oracular problem -the "glued trees" problem of Ref. [26] -in poly(n) time, while Thm. 2 is implied by Ref. [26], and (ii) by showing that our approach can simulate any quantum circuit of size L acting on q qubits in poly(L, q) time.More precisely we show that Problem 3 is BQP-complete, thereby proving Thm. 3.Although our results use physical systems that may seem artificial (e.g., resulting interactions between oscillators are not spatially local), they are strong evidence that no polynomial-time (in n) classical algorithm for simulating systems of coupled classical oscillators exists.Our results also add a new problem to the list of "natural" BQP-complete problems [27].

A. Oracle lower bound
At a high level, in the glued trees problem of Ref. [26], we are given oracle access to the adjacency matrix of a sparse graph with an ENTRANCE vertex, and the goal is to find the EXIT vertex.We map this problem to a system of masses and springs by putting a unit mass at each vertex and a spring of unit spring constant at every edge.Then we show that if the system starts at rest except with the ENTRANCE vertex having unit velocity, after polynomial time, the EXIT vertex will have inverse polynomial velocity, and thus it is a special case of Problem 2. Formally, we study the following problem: Problem 6.Consider a network of N = 2 n+1 − 2 masses coupled by springs, where M = 1l and K coincides with the adjacency matrix of a graph consisting of two balanced binary trees of depth n "glued" randomly as in Fig. VI A (i.e., the spring constants are κ jk = 1 if (j, k) is an edge of the graph or κ jk = 0 otherwise).Each mass (vertex) of the network is labeled randomly with a bit string of size 2n.The network contains two special and verifiable masses, ENTRANCE and EXIT, which correspond to the roots of both trees.Given oracle access to K and the label of the ENTRANCE mass, the problem is to find the label of the EXIT mass.
This problem is equivalent to that studied in Ref. [26] where the authors presented an efficient f F e X c + 5 q 0 F J 5 8 5 h D 9 w P n 8 A y J 2 N e w = = < / l a t e x i t > l = 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " E f z i v s Y K j m F a 5 l y + quantum walk based algorithm while showing that no classical algorithm can solve the problem efficiently in this oracular setting.In that work the root vertices do not have the extra edge that connects them to a "wall", so they can be easily verified as they are the only vertices of degree two.We add these extra edges to simplify the analysis below; this small change allow us to use some results on the spectral properties of the adjacency matrix already studied in Ref. [26].With this change, the oracle to access K allows us to verify whether a certain vertex is a root or not, since they are the only vertices (or masses) with nonzero diagonal spring constants; i.e., κ jj = 1 only if j corresponds to ENTRANCE or EXIT and κ jj = 0 otherwise.Since our problem is a minor modification of the glued trees problem, and any oracle query to our problem can be answered using a single oracle query to the glued trees problem, the original classical lower bound of Ref. [26] implies that our problem requires 2 Ω(n) queries to solve classically.
We then show that our quantum algorithm for simulating coupled classical oscillators can also be used to solve Problem 6 efficiently.This is a different approach from that of Ref. [26] as we are not simulating the quantum walk but rather the classical dynamics obtained from Newton's equation with a quantum algorithm.This shows that our quantum algorithm cannot be simulated classically efficiently in the oracular setting in general.Our main claim is the following lemma, which we prove in App.B. at which the kinetic energy of the EXIT mass is K EXIT (t) := 1 2 ( ẋEXIT (t)) 2 = Ω(1/poly(n)).Given this lemma, it is easy to see that Problem is a special case of the problem of estimating the kinetic energy, Problem 2. The quantum algorithm can prepare |ψ(t)⟩ efficiently under the conditions of Lemma 7. Because Lemma 7 proves that the kinetic energy of the EXIT mass will be inverse polynomially large after time t = O(poly(n)), a projective measurement on |ψ(t)⟩ will return the label of the EXIT with probability K EXIT (t)/E = Ω(1/poly(n)), where the energy is E = 1 2 ( ẋENTRANCE (0)) 2 = 1 2 in this case.Since the EXIT mass can be verified with the oracle by assumption, the probability of finding the label of EXIT can be made a constant close to 1 after O(poly(n)) repetitions of the previous procedure, giving an overall query and gate complexity O(poly(n))).Since Problem 6 requires 2 Ω(n) queries to solve classically, this yields the exponential lower bound in Thm. 2.

B. BQP-completeness
To prove this, we start from the standard BQPcomplete problem of simulating universal quantum circuits, i.e., the problem of determining some property of the output state |ϕ⟩ = U L . . .U 1 |0⟩ ⊗q , where L = poly(q).The unitary gates U l belong to a universal set such as {H, X, Toff}, where H and X are single-qubit Hadamard and Pauli X (bit flip) gates, and Toff is the three-qubit Toffoli gate.(Without loss of generality, all Hadamard gates can act on the last qubit, a property that is not necessary but we use it to simplify the presentation of the proof.)The specific problem is that of deciding if |ϕ⟩ is essentially the all zero state or has almost no overlap on it; see Problem 7. In App.C we show that this problem reduces to one in which the kinetic energy of a specific oscillator is either exponentially close to zero or at least 1/poly(q) large.This is essentially Problem for n = O(q).
We use the standard Feynman-Kitaev construction [28,29] to encode the U l 's into a Hamiltonian with an extra "clock" register.One property of this construction is that evolution under the Hamiltonian is known to apply the sequence of gates U L . . .U on a given initial state (with 1/poly(q) amplitude).We want this Hamiltonian to be the matrix A of a corresponding system of coupled oscillators.This way we could use the connection between classical and quantum systems discussed in Sec.III, which suggests that the evolved state of the oscillators encodes the amplitudes of |ϕ⟩.However, this has two problems: 1.The off-diagonal entries of the Hamiltonian would not be real negative numbers, i.e., they cannot be related to spring constants; 2. The evolution of the oscillators is induced by the √ A rather than A, so the evolution property of the Hamiltonian does not apply.
The first problem is due to the fact that Hadamards have negative matrix entries.To fix this, we observe that the same effect can be obtained using an operator with non-negative entries and an ancilla in the [30,31].This operator is no longer unitary, but this poses no additional problems as it is only being encoded in a Hamiltonian that, at this stage, corresponds to a system of oscillators coupled by springs.The second problem is addressed by showing that the spectral properties of √ A are such that there is an appropriate evolution time t = poly(q) after which the evolved state of the oscillators indeed encodes the output of the quantum circuit with 1/poly(q) amplitude.
Formally, we consider problem instances, i.e., systems of coupled oscillators, where N = (L + 1)2 q+1 , M = 1l, and The matrix 1l N is the N × N identity matrix.The 2 q+1 × 2 q+1 matrices W l are real-symmetric and act on one more qubit than do the U l matrices.We define them as W l = U l ⊗ 1l 2 if U l is the X or Toff gate acting on the space of q qubits.When U l is a Hadamard gate, which acts on the last qubit (the q th qubit), we replace it by the following 4 × 4 matrix acting on qubits q and q + 1: This matrix acts like Hadamard on qubit q when qubit q + 1 is set to |−⟩.Thus the entries of W l are non-negative in this case as well, the off-diagonal entries of A are {0, −1/ √ 2, −1}, and all diagonal entries are 4. Hence, the off-diagonal entries of the matrix K are {0, 1/ √ 2, 1} and one can show that κ jj ≥ 0 for all j ∈ [N ].These spring constants can be efficiently accessed using Eq. ( 25).The gates U l are 2-sparse, implying that the W l 's are also 2-sparse, and A and K are 4-sparse.These properties are outlined in Thm. 3.
For our proof, we consider the initial condition for the N oscillators where ⃗ x(0) = (0, 0, . . ., 0) T , ⃗ x(0) = (+1, −1, 0, . . ., 0) T , so that the energy is E = 1.Labeling each oscillator j ∈ [N ] by (l, r), where l ∈ [L + 1] and r ∈ [2 q+1 ], the oscillator under consideration is the one with l = L + 1 and r = 1.In App.C we show that computing the kinetic energy of this oscillator after time t = poly(q) = polylog(N ) solves the above BQP-complete problem.In addition, this classical system satisfies all of our conditions for the quantum algorithm to be efficient (i.e., the problem is in BQP).

VII. GENERALIZED COORDINATES
In general, when modeling a classical system using the harmonic approximation, Newton's equation can be simply written as ⃗ y(t) = −A⃗ y(t), where ⃗ y(t) ∈ R N encodes the generalized coordinates and A ≽ 0 is of dimension N × N .To prove Thm. 4 we use a standard approach based on quantum phase estimation [32][33][34] to first estimate the eigenvalues of H (2) := −X ⊗ A, which are γ η,j := (−1) η λ j , within certain precision that gives the eigenvalues of H within precision O(ϵ/t).Here, λ j ≥ 0 are the eigenvalues of A, and η ∈ {0, 1} determine the eigenvalues of −X as (−1) η , with eigenstates |0 X ⟩ = |−⟩ and |1 X ⟩ = |+⟩.We specifically use a standard QFTbased approach to phase estimation that eschews the need for measurement and in turn maps (within small error) each eigenvector of H (2) where |ϕ η,j ⟩ is some unspecified error state of the ancillas of unit norm, x |b x | 2 = 1, δ η,j ≤ δ PE , and ϵ PE and δ PE need to be set according to the error requirements.The states |x⟩ encode the eigenvalue estimates of H (2) .From each estimate x, we obtain the eigenvalue estimate of H by taking its sign and computing a square root, i.e., sign(x) |x|.Then we apply a phase factor yielding Last, we invert quantum phase estimation and other operations to (approximately) uncompute the estimated eigenvalues.The result is, for a value of the error tolerances where δ PE = O(ϵ 2 ) and a unitary approximating e −itH within error ϵ.The result is implied by the complexity of quantum phase estimation that requires O(∥A∥ max d log(1/δ PE )/ϵ PE ) queries to the oracle.See App.D for details.

VIII. CONCLUSION
We introduced an approach to simulating systems of coupled classical oscillators using a quantum computer.We showed how to map Newton's equations for these dynamics (coupled second order ordinary differential equations) to the Schrödinger equation (a first order partial differential equation), with polylogarithmic overhead in the number of oscillators under certain assumptions.Using this mapping, we demonstrated that any classical algorithm that simulates these classical dynamics requires at least 2 Ω(n) queries to the oracle, and that they can provide solutions to BQP-complete problems, and would thus require exponential time to solve classically under reasonable complexity theoretic assumptions.We also generalized our approach to the simulation of other classical harmonic systems.
While providing a large quantum advantage in certain contexts, these techniques also have significant limitations.For example, the approach is only efficient for computing particularly large or global properties and when masses and spring constants can be computed in time polylogarithmic in system size.Another feature of our algorithm is that its complexity is (almost) linear in the evolution time t, a scaling that might not be avoided in general [19,35,36], being efficient only if t is also polylogarithmic in system size.This would discourage applications where, for example, N = poly(t) (and when fast forwarding is not possible).This feature is expected to arise when simulating physical systems with geometrically-local interactions, and for initial states that are locally supported.In these examples, the relevant system size is determined by the "light cone", whose size in D spatial dimensions would scale as N ∼ t D .Nevertheless, even for these examples our quantum algorithm would still result in significant quantum speedups (e.g., super-quadratic), suggesting a new application area for quantum computers [37].
As many systems from molecular vibrations, to structural mechanics, to electrical grids, to neuronal activation can be modelled within the harmonic approximation, and since interesting dynamical features can appear at relatively short times (e.g., t independent of N ), we expect that there exist specific applications that meet all the requirements for our quantum algorithm to be efficient.Such applications will then benefit from this speedup with appreciable real-world impact; identifying them is one important next step in this line of research.
Last, we note that our classical-to-quantum reduction provides yet another way to think about quantum algorithms.For example, once mapped to a one dimensional system, the glued-trees problem of masses and springs becomes a simple wave propagation problem.Since waves propagate ballistically, it is now clear why our quantum algorithm, whose complexity is mainly dominated by evolution time but not by the number of masses, works in this case.Another example results from Grover's classical system of coupled pendulums [18] that, after using our reduction, provides another way to solve the unstructured search problem with O(t) queries, where t = O( √ N ).
acting on 2n + r + 4 qubits with r = O(log(1/ϵ ′ )) that provides a block encoding of H to within additive error ϵ ′ as follows: The quantum circuit that implements U H uses a controlled version of U B and its inverse once, and additional O(n) two-qubit gates (e.g., CNOTs and other single qubit gates).

State preparation using inequality testing
Before presenting the proofs of the lemmas, we revisit how inequality testing can be used for state preparation [38], since our constructions use this approach, which is more efficient than controlled rotation of an ancilla qubit as in Ref. [8].Let {β 1 , . . ., β N } be such that β j ≥ 0 can be expressed using r bits (see below) and β ≥ β j for all j ∈ [N ].Assume we have access to an oracle that, on input |j, z⟩ outputs j, z ⊕ βj , where z and βj ∈ [2 r ] are given as bit strings of size r (we will fix z = 0 . . .0).That is, in this notation, | βj ⟩ is a basis state, and the bits are obtained from the binary fraction β j = β[.b 1 . . .b r ].The method in Ref. [38] can be used to perform: where |ω j ⟩ is a state orthogonal to |0⟩ ⊗r |0⟩.The basic steps of the method are as follows: 1. Apply the oracle to compute βj , i.e., perform the map |j⟩ → j, βj .
3. Apply inequality testing to prepare j, βj 5. Apply the inverse of the oracle to reverse the computation of βj , i.e., perform the map j, βj → |j, 0⟩.
After these operations there will be amplitude βj /2 r = β j /β on |0⟩ ⊗r |0⟩.This is because the amplitude can be found by taking the inner product of the state in step 3 with x=1 |x⟩ |0⟩.The method requires one use of the oracle and its inverse, in addition to O(r) Hadamard gates for steps 1 and 4, and O(r) gates for the inequality test in step 3.It is also possible to compute more complicated functions of the oracle by rearranging the inequality in step 3, which is the method we will use here.

Proof of Lemma 8
We consider a block encoding of B † for simplicity, and later use it to provide a block encoding of B by taking the conjugate transpose.Using the padding described above, the dimension of B † is N 2 × N and its entries are of the form κ jk /m j or zero; specifically, the definition of B in Sec.III (obtained from Eq. ( 16)) implies We will construct a unitary U † B that is a block encoding of B † following simple steps and then explain how these can be simulated with quantum circuits.The steps require using a work register of qubits for some computations that we discard at the end.
1. Apply a unitary that performs the map 2. Apply the oracle S for the positions of non-zero entries of K to map |ℓ⟩ to |a(j, ℓ)⟩ according to Eq. (A1).
4. Prepare the equal superposition state of r ancilla qubits, i.e., 1 x=1 |x⟩, where r is given below.5. Using coherent arithmetic, compute the square of x and multiplications needed for the inequality test where r κ and r m are the numbers of bits of κjk and mj , respectively, also determined below.This gives the factor in the amplitude approximately κ jk /(m j ℵ).
6. Reverse the computation of the arithmetic for Eq.(A7) and the oracles for mj and κjk .This transforms the working registers back to an all-zero state.The previous sequence of unitaries define U † B .To show that the method is correct, consider steps 1-6.If we discard the working registers (used to store mj , κjk and perform the arithmetic), these steps combined implement approximately for some state |ω j ⟩ that is orthogonal to |0⟩ ⊗r |0⟩ on the ancilla qubits.We want the the error in this approximation to be O(ϵ ′ ), and we show how to achieve this below.Step 7 is then needed to rearrange the sum depending on whether k < j or k ≥ j, according to Eq. (A5).This step transforms Eq. (A8) to where ω ′ j is still orthogonal to |0⟩ ⊗r |0⟩ on the ancilla qubits.Applying step 8 to Eq. (A9), and considering the part of the state where the ancillas are in |0⟩ |0⟩ ⊗r |0⟩ only, we obtain ) is an edge and κ jk = 0 otherwise).The ENTRANCE and EXIT masses are additionally connected to a wall each with a spring of constant also 1.We write ⃗ x(t) ∈ R N and ⃗ x(t) ∈ R N for the positions and velocities the masses, where we assume their motion is constrained to one spatial dimension, such as the "horizontal" direction.Although we do not know the actual names of the masses a priori, our labels are such that j = 1 refers to the ENTRANCE, j = 2, 3 refer to the masses in the second column, and so on, until j = N = 2 n+1 − 2 represents the EXIT mass.The energy of the system is E = T (t) + U (t), where ( ẋj (t)) 2 , (B1) are the kinetic and potential energies, respectively.Newton's equation gives then a set of N coupled second-order differential equations: where A is the N × N -dimensional symmetric matrix Note that A = 3I N − A, where A is the adjacency matrix of the graph constructed from the two binary trees randomly glued, if we disregard the edges that connect the roots to their respective walls, and I N is the N × N identity matrix.Also, A is positive semi-definite because it is symmetric and diagonally dominant.This property also follows from the potential being U (t) = 1 2 ⃗ x(t) T A⃗ x(t) ≥ 0 for all ⃗ x(t) ∈ R N , which implies A ≽ 0. In fact, A is positive definite: The matrix A ′ = A − |1⟩⟨1| − |N ⟩⟨N | represents the above network of oscillators where the masses at the roots are not connected to any wall.Hence, A ′ contains a non-degenerate eigenvector of eigenvalue zero corresponding to the translations of the system, i.e., the eigenvector |u 1 ⟩ := 1 √ N j |j⟩.The only way for A to contain an eigenvalue zero is if this eigenvector |u 1 ⟩ was also an eigenvector of |1⟩⟨1| + |N ⟩⟨N |, which is not the case.Nevertheless, while A ≻ 0, its smallest eigenvalue is exponentially small in n since the expectation ⟨u 1 | A |u 1 ⟩ = 2/N is an upper bound on the smallest eigenvalue.
We would like to show the following property.If ⃗ x(0) = (0, 0, . . ., 0) T and ⃗ x(0) = (1, 0, . . ., 0) T , then ⃗ x(t) is such that the magnitude of its N th entry, corresponding to EXIT, is at least polynomially small in n for a time t that is at most polynomial in n (i.e., the kinetic energy of the N th oscillator is Ω(1/poly(n))).For these initial conditions, the solution to Eq. (B3) implies

B5)
The matrix A, and hence √ A, possesses a symmetry that allows one to simulate the dynamics of the N oscillators by considering that of 2n oscillators in one spatial dimension, instead.(A similar idea was used in coincides with the average probability of projecting |ψ(t)⟩ into a basis state that correspond to the EXIT mass.Thus, the goal reduces to showing that there exists T = O(poly(n)) such that P EXIT (T ) = Ω(1/poly(n)).If we prove this, then Eq. (B15) automatically implies the existence of t = O(poly(n)) such that K N (t) = Ω(1/poly(n)), which is our main goal.Finding such t can be done classically by simulating Eq. (B12) in time polynomial in n.
We can write Let {|λ 1 ⟩ , . . ., |λ 2n ⟩} be the 2n eigenvectors of the adjacency matrix Ã of eigenvalues λ 1 < λ 2 < . . .< λ 2n .These are also eigenvectors of Ã of eigenvalues γ l := √ 3 − λ l > 0. This implies To obtain this we used lim T →∞ T 0 dt e it(γ l +γ l ′ ) = 0 since the eigenvalues are positive, |⟨λ l | 1⟩| = |⟨λ l | 2n⟩| for all l ∈ [2n] due to a reflection symmetry of the network (i.e., the adjacency matrix of the graph in the line is invariant under the transformation l ↔ 2n − l + 1), and also the Cauchy Schwarz inequality for the last line as there are 2n different eigenvectors.This would already prove the desired result, but we are interested in finite times T < ∞ and, more precisely, showing a similar bound for T = O(poly(n)).
Using again the reflection symmetry of the network, for T < ∞ the average success probability can be written as where we used the identity cos 2 (α) = 1 2 (1 + cos(2α)).Then, the correction to where ∆ is the minimum spectral gap between any pair of eigenvalues γ l , that is, ∆ := min l |γ l+1 − γ l |. (Also ) Also, if we order the eigenvalues so that γ l+1 > γ l for all l = [2n − 1], for l ≥ 2 we have γ l ≥ ∆ and For l = 1, the eigenvalue γ l is exponentially small in n and the corresponding average can be O(1).Combining these equations and using the reflection symmetry of the network we obtain where we also used l |⟨1| λ l ⟩| 2 = 1.We need to show that both these terms are small.
In Ref. [26] it has been shown that the eigenvalues λ l 's are separated by spectral gaps bounded by ∆ ′ = Ω(1/n 3 ).The same follows for the γ l 's.More precisely, suppose that the two closest eigenvalues of Ã are 3 − λ ≤ 6 and 3 − (λ + ∆ ′ ).Then, the two closest eigenvalues of Ã are separated as Then, since ∆ ′ = Ω(1/n 3 ), we have ∆ = Ω(1/n 3 ).This implies that we can choose T = O(n 4 ) so that, for example, It is also possible to show that |⟨1| λ 1 ⟩| 4 is exponentially small in n.Consider the vector The length of this vector is ) and that |v 1 ⟩ /∥ |v 1 ⟩ ∥ is exponentially close to the normalized eigenvector of lowest eigenvalue as the spectral gaps are at least ∆ = Ω(1/n 3 ); that is  O(1/ exp(n)).Since the first entry of Combining these bounds, and for the choice of T above, we obtain Together with Eq. (B17), this gives which is the desired result: we can choose T = O(n 4 ) so that the average probability P EXIT (T ) is Ω(1/n), implying the existence of t ∈ [0, T ] with the desired property | ẋN (t)| = Ω(1/poly(n)).We note that numerical solutions show a stronger result, where t ∝ n suffices.See Fig. 4 for an example.
measuring the first qubit of the output state, we have to decide if it is 1 with probability at least 2/3 or 1 with probability at most 1/3, promised that one of these is the case.By known results, we can assume our gate set is real and contains only Hadamard and Toffoli, since this is a universal gate set for real quantum computation, which is computationally as powerful as complex quantum computation [39,40].We also add the single-qubit Pauli X gate so that we can create the state |1⟩ from |0⟩, as we want the computation to start with all qubits in state |0⟩.Also, the X and Toffoli gates generate swap gates, and we can then assume that all Hadamards act on the same qubit.This property will be useful to simplify the presentation of the proof, but it is not necessary, and one could in principle allow Hadamards to act on any qubit.Given such a circuit, we can amplify the constants 2/3 and 1/3 to 1 − 1/ exp(q) and 1/ exp(q) by running O(q) copies of the circuit in parallel and taking the majority vote of the first output qubits of all circuits.Now we have a new circuit on O(q 2 ) qubits whose first qubit is almost certainly (with probability 1 − 1/ exp(q)) |1⟩ when the answer is yes and almost certainly |0⟩ when the answer is no.We then use an additional qubit and copy this answer to that qubit.Since the qubit we are copying is exponentially close to being |0⟩ or |1⟩, let us assume it is one of these, and this will only introduce an error of 1/ exp(q).We then run the circuit's inverse on the remaining qubits to restore them to the all zero state.The resulting output state is now exponentially close to all zeros if the additional qubit was also zero, which happens when the input was a no instance.If it was a yes instance, the additional qubit we added will be exponentially close to |1⟩, and hence the overall state will have almost no amplitude on the all zeros state.
Thus deciding if the output state of a O(q 2 ) qubit circuit of size poly(q) is 1 − 1/ exp(q) close to the all zeros state or has at most 1/ exp(q) overlap on the all zeros state is BQP-complete.This yields the stated result by renaming q to q 2 .(Our proof below does not require the closeness to be exponentially small, and it would suffice to have inverse polynomial closeness, as long as this polynomial was smaller than all the other polynomials appearing in the problem.)

From a circuit to a network of oscillators
We now show that our quantum algorithm for simulating coupled classical oscillators allows us to solve the BQP-complete problem above.We begin with the given quantum circuit on q qubits with L = poly(q) gates.If the gates are U 1 to U L , the output of this circuit is U L . . .U 1 |0⟩ ⊗q .Our goal is to decide if this state is essentially |0⟩ ⊗q or has essentially no overlap with |0⟩ ⊗q .The unitaries U 1 to U L are either the single-qubit gates Hadamard H or Pauli X (tensored with identity on all other qubits) or the three-qubit Toffoli gate Toff (tensored with identity on all other qubits).These gates are 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 in the corresponding one-qubit and three-qubit subspaces.It will be very useful in our context that these are real as the resulting Hamiltonian will have real entries, which we require because the κ jk are real.
From this circuit we will create a network of N oscillators, where N = (L + 1)2 q+1 .It is convenient to label each oscillator j ∈ [N ] using two indexes j → (l, r), where l ∈ [L + 1], r ∈ [2 q+1 ], and j = (l − 1)2 q+1 + r.We will call the oscillator with l = L + 1 and r = 1 the output oscillator.In our construction, the N × N matrix of spring constants K that describes the couplings between oscillators, which depends on the U l 's, is 4-sparse and each entry is bounded as κ jk ≤ 4. We will have all masses be m j = 1.Note that and ⃗ y(t) = √ M⃗ x(t) = ⃗ x(t) in this case.This simplifies our analysis.In the following, we show that for our oscillator network, there exists a time t = polylog(N ) such that, for the initial conditions where ⃗ x(0) = (0, . . ., 0) T , ẋ1 (0) = 1, ẋ2 (0) = −1, and ẋj (0) = 0 for all j > 2, the kinetic energy of the output oscillator is either exponentially close to 0 (when the original circuit had almost no overlap with |0⟩ ⊗q ), or 1/polylog(N ) (when the original circuit's output is very close to |0⟩ ⊗q ).Since our algorithm allows us to estimate the kinetic energy of a given oscillator to additive 1/polylog(N ) precision, we can distinguish these two cases with complexity O(polylog(N )).The initial quantum state |ψ(0)⟩ is simply |−⟩ tensored with the all zeros state, and hence easy to prepare.This proves that Problem 2 is BQP-complete even with these constraints on the initial conditions, t, ϵ, M, and K, which is the desired result.
We now provide the details of this construction.We use the bra-ket notation for specifying vectors and matrices, which is convenient for relating properties of the classical system with properties of quantum states.We will specify our oscillator network using the A matrix, instead of the K matrix.This real matrix is PSD and has non-negative entries on the diagonal and non-positive entries on the off-diagonal.The standard Feynman-Kitaev [28,29] circuit-to-Hamiltonian construction gives us a Hamiltonian of the form where W l is the l th gate in the circuit, i.e., U l .We cannot let A equal this Hamiltonian, since the matrix is not PSD and there are off-diagonal entries of both signs.The first issue is easily fixed by adding a multiple of the identity, but it will require some work to fix the second issue.We will adjoin an additional qubit to work in a slightly larger space.Formally, we consider the following Hilbert space: Here, H clock is the (L + 1)-dimensional Hilbert space that is used to describe the state of the clock (i.e., corresponding to the first register in Eq. (C2)), which is used to track the progress of a simulated quantum computation on a state in H comp ⊗ H 2 where H comp is 2 q -dimensional and H 2 is two-dimensional.The purpose of H comp is to store the state of the given quantum circuit at a time given by the clock.The purpose of H 2 is to address the issue on the off-diagonal entries raised above.The off-diagonal entries of A need to have the same sign, but our gate set includes a gate, H, with entries of both signs.We address this by providing a resource state |−⟩ = 1 √ 2 (|0⟩ − |1⟩) in the last register which we use to effectively create negative signs in a subspace although the operators will only have positive entries.Similar constructions that create Hamiltonians with non-negative entries have been used within the context of Hamiltonian complexity; see Refs.[30,31] for example.
We will map our circuit to a system of oscillators that satisfy the assumptions made in Problem 1 by choosing the couplings to implement an "encoded" sequence of gates {W l } l as where each W l corresponds to an encoded version of the l th gate in the original circuit, which is either a Hadamard, X, or a Toffoli, and 1l N is the identity matrix on the system of dimension N .Notice that W l lives in H comp ⊗ H 2 , whereas the original unitaries U l lived in H comp .So we encode each gate into a larger space with 1 additional qubit.We want all the entries in the encoded versions of Hadamard, X, and Toffoli to be non-negative, which will make the off-diagonal entries non-positive in A due to Eq. (C4).Let's start with the Hadamard matrix.Without loss of generality, because we can swap qubits, we assume the Hadamard always acts on qubit q, the last qubit of H comp .We describe our encoded Hadamard gate via the following positive-valued real-symmetric block matrix acting on the last qubit of H comp and H 2 : Hence, when the first register in the tensor product is l = L + 1, the second register is U L . . .U So all we have to show is that there is exists a t = polylog(N ) such that |α L+1 (t)| 2 = Ω(1/polylog(N )) = Ω(1/poly(q)), which will allow us to distinguish the two cases by using our algorithm O(polylog(N )) = poly(q) times.

Establishing inverse polynomial overlap
It is not too hard to establish inverse polynomial overlap |α L+1 (t)| 2 = Ω(1/polylog(N )) for t = polylog(N ).(In fact, it is possible to get perfect overlap as we discuss in the next section.)To prove this result, we use some known results on the spectral properties of X ′ .Its eigenvectors are and the eigenvalues are where l ∈ [L + 1].Note that 6 > γ l > 2. The |ϕ l ⟩'s are also the eigenvectors of √ X ′ , whose eigenvalues are γ ′ l := √ γ l and √ 6 > γ ′ l > √ 2. Let ∆ l = γ l+1 − γ l be the spectral gaps of X ′ and note that ∆ l > 0 and ∆ l /γ l > 0. The spectral gaps of √ X ′ satisfy (for l ∈ [L]) This derivative is zero at x = 0 and π, which is outside the region of values of x.The derivative will therefore take its smallest (nonzero) values at the nearest allowed values of x: , and x = π(L + 1) We have (L ≥ 1) Thus we have (C23) Hence, the smallest spectral gap of √ X ′ , ∆ ′ = min l (γ ′ l+1 − γ ′ l ), satisfies ∆ ′ = Ω(1/L 2 ).The amplitude α L+1 (t) that we are interested in can be written as Recall that the probability of measuring |L + 1⟩ after evolving with √ X ′ for time t, when the initial state is |1⟩, is where we used that cos(θ) = 1 2 (e iθ + e −iθ ).
The result gives us the bound on the error in the signed square root as follows: ∆ η,j := max Hence, from the estimate x of γ η,j within error at most ϵ PE , we can compute an estimate of (−1) η λ j , the eigenvalue of H, within error at most ∆ η,j as above.
To implement e itX⊗ √ A , we apply a phase factor proportional to the estimated eigenvalue from QPE.Then the state after applying the phase factor is transformed as The states |ϕ η,j ⟩ and |ϕ ′ η,j ⟩ are states of the ancillas of unit norm, supported in the space orthogonal to {|x⟩ : x ∈ S η,j }.After the phase was applied, we need to invert QPE.The inverse phase estimation and projection onto the zero state of the ancillas that were used to estimate the eigenvalues corresponds to applying to give The error can be bounded from the Euclidean distance between the above state and the correct state e −it(−1) η √ λj |η X , λ j ⟩.Using δ η,j ≤ δ PE , this distance is upper bounded by where the expectation value is on the probability distribution given by |b to satisfy max η,j ∆ η,j t ≤ ϵ/2, we can take Using these expressions in Eq. (D2) gives us an overall query complexity for the phase estimation approach.This is the stated result.
There are further elementary two-qubit gates arising from three main areas.
1. Computing sign(x) |x|, the estimate of the eigenvalue of H, multiplying by t, and then applying the phase factor.The dominant complexity is that from computing the square root, which scales using Newton iteration and textbook multiplication as O(log 2 (1/ϵ)).The phase rotation then requires a linear number of controlled rotation gates based on the target, which in turn requires at most O(log(1/ϵ)) two-qubit gates to implement.Thus the former cost dominates.
2. The gates for the implementation of the block encoding will be logarithmic in N and the allowable error of the block encoding.Because the allowable error in the block encoding needs to be ϵ divided by the number of block encodings in the phase estimation, there will be logarithmic factors in many of the parameters here.
3. The gate complexity of preparing the control states states for phase estimation with optimal confidence intervals is at most linear in the dimension of this control register [47], which is the same as the number of oracle calls.This is the least significant contribution to the complexity and gives no logarithmic factors.
In quoting the complexity we give the logarithmic factor in N coming from implementing the block encoding, but for simplicity use O and do not explicitly give the logarithmic factors in other parameters.
separately create superpositions over the initial positions and initial velocities.Our construction is related to that in Lemma 8 of App. A. As in Sec.IV, we consider a setting where the N × N matrices M and K can be queried through the use of a unitary S that gives the positions and non-zero entries of these matrices.Lemma 10.Let K be the N × N symmetric and d-sparse matrix of spring constants, κ max ≥ κ jk ≥ 0, M ≻ 0 be the N × N diagonal matrix of masses, and m max ≥ m j > 0, where κ max and m max are known.Assume we are given access to K and M through an oracle S, and access to a unitary U that performs the map x j (0) |j⟩ , are normalized states that encode the initial states of the oscillators in their amplitudes, and α and β are known norms of the vectors ⃗ x(0) and ⃗ x(0), respectively.Then, there exists a quantum algorithm that prepares a state that is ϵ-close in Euclidian norm to where E > 0 is a known constant (the energy of the classical oscillators) and ⃗ µ(0) ∈ R M for M = N (N + 1)/2 is a vector whose entries are √ κ jj x j (0) or √ κ jk (x j (0) − x k (0)), k > j.The quantum circuit makes ) uses of U, S, and its inverses, in addition to is the energy of a system of N uncoupled oscillators where all masses are m max and all individual spring constants are κ max , and for the same initial conditions.
Proof.For ease of implementation, we can encode the state on 2n + 1 qubits, where n = log(N ).In this space, the initial state is represented as First we rotate a qubit and apply the state preparation oracle as The way we have described the oracle for ⃗ x(0) and ⃗ x(0) allows for the case where α or β may be zero, in which case U can be arbitrary on that subspace because it has no effect on the state above.Then our goal is to apply √ M to the | ⃗ x(0)⟩ portion, and B † √ M to the |⃗ x(0)⟩ portion.The most challenging part is for applying B † √ M. To do this, we essentially apply the same construction in Lemma 8 in App.A, but for block encoding B † √ M rather than B † .First consider the operation where we prepare a superposition over d values in an ancilla to give which is the correct state (Eq.(E3)), but subnormalised indicating that it is not produced deterministically and we need to use amplitude amplification.
The number of amplitude amplification rounds scales as the inverse of the amplitude, so the complexity in terms of α, β, E is We assume we know the constants for this Lemma to simplify the amplitude amplification.It is also possible to perform amplitude amplification when the amplitude is unknown but bounded [48].Note that which is the kinetic energy of the system of oscillators at t = 0 if all masses were m j = m max .Also, is the potential energy of a system of N uncoupled oscillators at t = 0 if all spring constants were κ jj = κ max and is the total energy of such a system, it satisfies The number of amplitude amplification rounds is then O( E max d/E).Each amplitude amplification round makes two uses of S, S † , and also one use of U and U † .Hence, the overall query complexity of our algorithm that prepares |ψ(0)⟩ is To find the accuracy of the block encoding we need to account for the error due to applying the factors of κ jk /κ max and m j /m max by inequality testing.If we give the maximum errors in these factors as δ (corresponding to log(1/δ) bits in the inequality testing), then the error in Eq. (E12) is upper bounded as 1 The inequality on the third line is obtained by noting that we can move from the state on the third line to that on the second line by the procedure described above, with an inequality test between j and k, a controlled swap and phase, then projection onto |+⟩ on the ancilla.Thus, δ corresponds to the error before amplitude amplification, and it is amplified by a factor of O( E max d/E).
If we are aiming for error ϵ in the state preparation, we should therefore take Because the state preparation by inequality testing uses squares, we have a gate complexity that is the square of log(1/δ) for each step of amplitude amplification, and therefore the gate complexity from this source is Moreover, there are O(log(N )) gates needed to implement an inequality test and controlled swap for each step of amplitude amplification.Using d ≤ N , we can give the overall gate complexity as which are the stated results.
The factor E max /E in the complexity can become large when the masses m j or spring constants κ jk lack uniformity.In this case, the state |ψ(0)⟩ is more complicated, and one has to "pay extra" for its preparation.In addition, the query complexity in Lem. 10 matches a lower bound for some instances.For example, consider the case where x j (0) = 0 for all j ∈ [N ], which implies β 2 = 0, and ẋj (0) = 1/ √ N for all j ∈ [N ], which implies α 2 = 1.Let the masses be m j = 1, for some unknown j, and m j ′ = 0 otherwise.Then, |ψ(0)⟩ is simply a state |j⟩ in a computational basis, and our state preparation algorithm would output the unknown j with high probability, as in the unstructured search problem [49].A lower bound on the query complexity for this case is Ω( √ N ) [50].Since E max = 1 2 m max α 2 = 1 2 , E = 1 2 ( ẋj (0)) 2 = 1 2N , and d = 1 for this instance, the query complexity of our approach is O( √ N ), matching the lower bound.In Refs.[5,51,52] this is done by approximating the solutions to the differential equations as solutions to systems of linear equations.The complexity of this approach depends on several parameters, in particular the condition number of the matrix to be inverted.In Ref. [13] a solution is found by giving the solution to the differential equation in exponential form, and then using the LCU approach to approximate the exponential; a related approach that approximates the exponential operator is given in Ref. [55].That approach is not applicable here, because it requires the matrix to be normal.Here the matrices -the 2×2 block matrices including A -are not normal, so the approach cannot be used.Nevertheless, even disregarding such difficulties, we show that this way of encoding the coordinates of the oscillators as in Eqs.(G1) or (G3) is readily problematic.
To observe this, it suffices to consider the evolution of a single normal mode of frequency ω k , i.e., the eigenvector of A of eigenvalue (ω k ) 2 .We obtain where ϕ k ∈ R is the initial phase and ⃗ a k ∈ R N is the eigenvector of A. Computing the time average we obtain and (The latter is also the time average of the kinetic energy.)In a case where A gives rise to normal modes of low frequency, we have |ω k | ≪ 1 (in the corresponding units), and then The implication is that the support of |⃗ v(t)⟩ is, on average, mostly concentrated on the subspace spanned by basis vectors that encode ⃗ y(t).These do not encode the terms appearing in the kinetic energy or ⃗ y(t).Hence, the complexity of solving Problem 2 using this encoding, where |⃗ v(t)⟩ is prepared rather than |ψ(t)⟩, is at least linear in 1/|ω k |, which is the factor needed to increase the desired amplitudes to a constant.This can become unbounded if |ω k | → 0.
In the second case we note that the support of | ⃗ w(t)⟩ is, on average, mostly concentrated on the subspace spanned by basis vectors that encode ⃗ y(t) instead.That is, the case ⟨∥A⃗ y(t)∥ 2 ⟩ ≪ ⟨∥ ⃗ y(t)∥ 2 ⟩ can occur when there are normal modes of low frequency.To solve Problem 1 or Problem 4, and especially Problem 5, the second component of the state must be proportional to B † ⃗ y(t) or √ A⃗ y(t), which requires the application of the pseudo-inverse of B or √ A to | ⃗ w(t)⟩.The condition number of these pseudo-inverses is also polynomial in 1/|ω k |; for example ∥(⃗ a k ) T B∥ = |ω k | ≪ 1 for low frequencies.This implies that the complexity of solving Problem 1, Problem 4, or Problem 5 using this encoding, where | ⃗ w(t)⟩ is prepared rather than |ψ(t)⟩, is at least linear in 1/|ω k |, which is the complexity of the most efficient methods to implement the pseudo-inverse [53].As in the previous case, this can also become unbounded if |ω k | → 0.
Similar observations follow directly from the results for the complexity in Refs.[5,51].There the complexity is polynomial in the condition number of the matrix that diagonalises the matrix appearing in the differential equations, which is here that with the blocks −A and 1l N in Eq. (G1), or the negative of these for ⃗ w(t) in Eq. (G3).The matrix that diagonalises this matrix has condition number max(|ω k |, |ω k | −1 ) .
(G11) Thus it will become unbounded when |ω k | → 0, as described above.Technically, we cannot use the result as given in Ref. [5] directly, because the stability condition it uses requires the matrices appearing in the differential equations to not have eigenvalues exactly on the imaginary axis, as is the case here.
An alternative approach to avoid dependence on the condition number of the diagonalising matrix is that using the log-norm [52].There we would need to subtract the identity times the log-norm from the original matrix.Here, the log-norm is This means it is always positive, and so subtracting a multiple of the identity would result in a norm that exponentially decreases.Hence the approach of Ref. [52] would yield a complexity that is exponential in time.At a high level, our encoding is motivated by the conservation of energy in time, and places equal emphasis on those terms that encode the kinetic energy ( ⃗ y(t)) and the potential energy (B † ⃗ y(t)) of the oscillators.The other two encodings do not have this feature: |⃗ v(t)⟩ underrepresents the kinetic energy terms and | ⃗ w(t)⟩ underrepresents the potential energy terms, bringing the issues discussed above.
Last, we provide a comment on the relation between our results and a closely related result in Ref. [11] for simulating the wave equation.In Ref. [11] it is shown that the wave equation, a second-order and homogeneous differential equation, can also be mapped to a Schrödinger equation.Their construction is related to ours in that it uses a factorization of the (discrete) Laplacian as L = BB † , and indeed it is well-known that the wave equation is one example of a classical system of coupled oscillators where all masses and springs are uniform, and where the couplings between oscillators are geometrically local (e.g., on the square grid).However, the encoding used in Ref. [11] is different to ours; for example, N amplitudes are reserved to encode the intensity of the wave, which corresponds to ⃗ y(t) in our case.(The component of ⃗ y(t) on the eigenvector of eigenvalue 0 of L is projected out.)The other amplitudes are proportional to B + ⃗ y(t), where B + is the pseudo-inverse of B. That is, the state prepared by the algorithm of Ref. [11] is of the form where the constant of proportionality is also time-independent.This is essentially the same encoding as the one discussed in Appendix F, since P⃗ y(t) = ⃗ y(t) by assumption.The length of this vector is preserved in time since the evolution is unitary.Nevertheless, one implication of using this encoding to solve Problem 1 or Problem 2 for this example is the need to invert B + initially, to obtain amplitudes proportional to B + ⃗ y(0).The condition number can be large, i.e., it is O(poly(N )) for the example of the wave equation, and hence the initial state preparation step (t = 0) can be inefficient.Indeed, Ref. [11] manages to give evidence of a polynomial quantum speedup only.In addition, when considering the wave equation, the system is geometrically local and implies that N is polynomial in the evolution time t (i.e, t is exponential in n).In our problem, however, we can treat a significantly larger set of instances: we do not impose uniform masses, uniform spring constants, or even geometrically-local interactions, and we can allow for times t = O(poly(n)).This generality together with our improved encoding are key to obtaining our claimed exponential quantum speedups for the problems defined.
< l a t e x i t s h a 1 _ b a s e 6 4 = " A R 8 b w 1 i L V + h / T J 1 W t H I G v C Z C 5 c 0 = " > A A A B 6 n i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y J U Y 9 E L x 4 x y i O B D Z k d e m H C 7 O x m Z t Z I C J / g x Y P G e P W L v P k 3 D r A H B S v p p F L V n e 6 u I B F c G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q e N U M W y w W M S q H V C N g k t s G G 4 E t h O F N A o E t o L R z c x v P a L S P J Y P Z p y g H 9 G B 5 C F n 1 F j p / q l X 6 R V L b t m d g 6 w S L y M l y F D v F b + 6 / Z i l E U r D B N W 6 4 7 m J 8 S d U G c 4 E T g v d V G N C 2 Y g O s G O p p B F q f z I / d U r O r N I n Y a x s S U P m 6 u + J C Y 2 0 H k e B 7 Y y o G e p l b y b + 5 3

x 2 < 1 < 4 <
l a t e x i t s h a 1 _ b a s e 6 4 = " h 5 e K r k j e G s j m p a + b W d e U V y U R G 5 I = " > A A A B 7 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i y e p Y N p C G 8 p m u 2 2 X b j Z h d y K W 0 N / g x Y M i X v 1 B 3 v w 3 b t s c t P X B w O O 9 G W b m h Y k U B l 3 3 2 y m s r K 6 t b x Q 3 S 1 v b O 7 t 7 5 f 2 D h o l T z b j P Y h n r V k g N l 0 J x H w V K 3 k o 0 p 1 E o e T M c 3 U z 9 5 i P X R s T q A c c J D y I 6 U K I v G E U r + U / d 7 G 7 S L V f c q j s D W S Z e T i q Q o 9 4 t f 3 V 6 M U s j r p B J a k z b c x M M M q p R M M k n p U 5 q e E L Z i A 5 4 2 1 J F I 2 6 C b H b s h J x Y p U f 6 s b a l k M z U 3 x M Z j Y w Z R 6 H t j C g O z a I 3 F f / z 2 i n 2 r 4 J M q C R F r t h 8 U T + V B G M y / Z z 0 h O Y M 5 d g S y r S w t x I 2 p J o y t P m U b A j e 4 s v L p H F W 9 S 6 q 5 / f n l d p 1 H k c R j u A Y T s G D S 6 j B L d T B B w Y C n u E V 3 h z l v D j v z s e 8 t e D k M 4 f w B 8 7 n D / 9 F j t I = < / l a t e x i t > xN < l a t e x i t s h a 1 _ b a s e 6 4 = " Q R l F K T d Z F T H 2 T Z m C L u q 5 i s p B s x g = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i 8 e K 9 g P a U D b b T b t 0 s w m 7 E 7 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j m 6 n f e u T a i F g 9 4 D j h f k Q H S o S C U b T S / V P P 6 5 U r b t W d g S w T L y c V y F H v l b + 6 /Z i l E V f I J D W m 4 7 k J + h n V K J j k k 1 I 3 N T y h b E Q H v G O p o h E 3 f j Y 7 d U J O r N I n Y a x t K S Q z 9 f d E R i N j x l F g O y O K Q 7 P o T c X / v E 6 K 4 Z W f C Z W ky B W b L w p T S T A m 0 7 9 J X 2 j O U I 4 t o U w L e y t h Q 6 o p Q 5 t O y Y b g L b 6 8 T J p n V e + i e n 5 3 X q l d 5 3 E U 4 Q i O 4 R Q 8 u I Q a 3 E I d G s B g A M / w C m + O d F 6 c d + d j 3 l p w 8 p l D + A P n 8 w c O n o 2 p < / l a t e x i t > x l a t e x i t s h a 1 _ b a s e 6 4 = " h t b A d 2 s 2 V 0 A C G 7 D p 8 H 7 4 9 J 3 d p 8 4 = " > A A A B 7 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 V T F t o Q 9 l s t + 3 S z S b s T s Q S + h u 8 e F D E q z / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 M J H C o O t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U N H G q G f d Z L G P d D q n h U i j u o 0 D J 2 4 n m N A o l b 4 X j 2 5 n f e u T a i F g 9 4 C T h Q U S H S g w E o 2 g l / 6 m X 1 a a 9 c s W t u n O Q V e L l p A I 5 G r 3 y V 7 c f s z T i C p m k x n Q 8 N 8 E g o x o F k 3 x a 6 q a G J 5 S N 6 Z B 3 L F U 0 4 i b I 5 s d O yZ l V + m Q Q a 1 s K y V z 9 P Z H R y J h J F N r O i O L I L H s z 8 T + v k + L g O s i E S l L k i i 0 W D V J J M C a z z 0 l f a M 5 Q T i y h T A t 7 K 2 E j q i l D m 0 / J h u A t v 7 x K m h d V7 7 J a u 6 9 V 6 j d 5 H E U 4 g V M 4 B w + u o A 5 3 0 A A f G A h 4 h l d 4 c 5 T z 4 r w 7 H 4 v W g p P P H M M f O J 8 / 1 8 O O u A = = < / l a t e x i t > x l a t e x i t s h a 1 _ b a s e 6 4 = " z j N o D 0 H a A d t k X Z G a d e 0 n o N d t 8 n T b i H F r S B w A S e 4 R X e n M R 5 c d 6 d j 2 V r w c l n T u E P n M 8 f d P a P A Q = = < / l a t e x i t > m 1 = m 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " s + n 8 o X P 8 m Q e n a / C o z x 4 7 V M o I 4 z a P o w g n c A r n 4 M E V N O A e m t A C A m N 4 h l d 4 c x L n x X l 3 P h a t B S e f O Y Y / c D 5 / A H s O j w U = < / l a t e x i t > m 3 = m < l a t e x i t s h a 1 _ b a s e 6 4 = " r J b d / B t E x c C j T 6 U M S W i c H X r o a w U = " > A A A B 8 n i c b V B N S w M x E J 3 1 s 9 a v q k c v w S J 4 s e x K U S 9 C 0 Y u n U s F + w H Y p 2 T T b h i a b J c k K Z e n P 8 O J B E a / + G m / + G 9 N 2 D 9 r 6 Y O D x 3 g w z 8 8 K E M 2 1 c 9 9 t Z W V 1 b 3 9 g s b B W 3 d 3 b 3 9 k s H h y 0 r / 4 a b / 4 b 0 3 Y P 2 v p g 4 P H e D D P z Y i 2 4 B d / / 9 g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 5 q m S Z 0 T I Y s d F S S h N k o m 5 8 8 x W d O 6 e O B M q 4 k 4 L n 6 e y I j i b W T J H a d C Y G R X f Z m 4 n 9 e m M L g J s q 4 1 C k w S R e L B q n A o P D s f 9 z n h l E Q E 0 c I N d z d i u m I G E L B p V R y I Q T L L 6 + S 1 k U 1 u K p e P l x W a r d 5 H E V 0 g k 7 R O Q r Q N a q h e 9 R A T U S R Q s / o F b 1 5 4 L 1 4 7 9 7 H o r X g 5 T P H 6 A + 8 z x 9 W Q J F N < / l a t e x i t >  NN < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 3 g I Q a 2 h D t Y 8 Z C 1 O j 9 Q b T V H x i O Y = " > A A A B 7 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g x b A r Q T 0 G v X i S C O Y B S Q i z k 9 5 k y O z s M j M r h i U f 4 c W D I l 7 9 H m / + j Z N k D 5 p Y 0 F B U d d P d 5 c e C a + O 6 3 0 5 u Z X V t f S O / W d j a 3 t n d K + 4 f N H S U K I Z 1 F o l I t X y q U X C J d c O N w F a s k I a + w K Y / u p n 6 z U d U m k f y w Y x j 7 I Z 0 I H n A G T V W a j 7 1 0 r s z b 9 I r l t y y O w N Z J l 5 G S p C h 1

1 FIG. 1 .
FIG. 1.An example system of N/2 oscillators in D = 2 dimensions.We use x1(t) for the first coordinate of the first mass and x2(t) for the second coordinate of the first mass, and so on.Since the first and second mass now correspond to the same original mass, m1 = m2.
FIG.2.A network of N = 2 n+1 − 2 coupled oscillators obtained by randomly gluing two binary trees.Each mass is mj = 1 for all j ∈ [N ] and is labeled by a random bit string of size 2n.Edges denote springs of constant κ jk = 1.The goal is to find the label of the EXIT mass given the label of the ENTRANCE mass and given oracle access to the network.

5 <
l a t e x i t s h a 1 _ b a s e 6 4 = " 4 m S R i A O C 1 H P b U s b y d 7

)
This can be seen to act as a logical Hadamard when acting on a state of the form |ψ⟩ |−⟩ because of the In general, we can write cos √ X ′ t |l = 1⟩ = L+1 l=1 α l (t) |l⟩ , (C16) with α l (t) ∈ R, because the Taylor expansion of the cosine function has only even powers of √ X ′ .Then ⃗ y