Witnessing eigenstates for quantum simulation of Hamiltonian spectra

We introduce the concept of an eigenstate witness and use it to find energies of quantum systems with quantum computers.


INTRODUCTION
The simulation of quantum mechanical systems, using conventional classical methods, requires resources that make the problems rapidly intractable when the size of the system grows (1). Since Feynman's seminal proposal, several algorithms for quantum simulation have followed (2)(3)(4), and many demonstrations have been reported in different physical systems (5)(6)(7)(8)(9)(10)(11)(12)(13)(14). Calculating the spectrum of a given Hamiltonian is a fundamental problem of widespread applicability, necessary, for example, to understand reaction rates or optical spectra in quantum chemistry. In particular, characterization of excited states is required to study energy and charge transfer processes, such as those in bulk heterojunction solar cells or photosynthetic light-harvesting complexes (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16). These kinds of problems are hard for classical computers and, in the most general case, for quantum computers as well. However, quantum devices are expected to provide scalable solutions to some instances of interest (2,3). Furthermore, also in those cases where classical methods can successfully describe the ground state (for example, for weakly interacting Hamiltonians), excited states are often hard to access (17,18), increasing the interest toward quantum methods able to address the problem of finding an efficient description of excited states. Here, we show promising progress in this direction by introducing the concept of an "eigenstate witness," a quantity that has no known efficient analog in classical algorithms. This witness is based on the entropy acquired by a quantum register, whose time evolution is controlled by an ancillary qubit.
Given an approximate eigenstate, the quantum phase estimation algorithm (QPEA) can efficiently estimate the corresponding eigenvalue. A more practical version, the iterative phase estimation algorithm (IPEA) (19), has been demonstrated using different quantum hardware, such as nuclear magnetic resonance, photonic, and superconducting systems (9)(10)(11). To prepare the input eigenstates, adiabatic state preparation has been proposed as a potentially scalable solution (3), at the cost of expensive state preparation and deep circuits, making it unsuitable for near-term implementations on quantum computers.
Variational quantum eigensolvers (VQEs), using a hybrid quantumclassical approach, were designed to address these shortcomings (10,12,(20)(21)(22)(23). These methods prepare states described via a chosen parameterization, known as ansatz, leveraging pre-existing knowledge about the system. Different types of ansatz have been proposed for the variational search, such as the unitary coupled cluster (UCC), which is among the most promising ones to tackle quantum chemistry problems (12,20,24). In addition, VQE methods are believed to have unique robustness to certain errors in estimating the ground state and its eigenvalue (10,25). They are, however, quadratically less precise than QPEA, because they rely on sampling for the energy estimation in the original formulation. Crucially, variational methods could only target ground states to date.
A linear response methodology and a spectrum folding method have been proposed as possible solutions (12,25,26). However, although the linear response methodology maintains the low coherence time advantages of the original VQE, it requires additional sampling measurements and cannot refine approximate excited states. Instead, the folded spectrum (FS) method requires a quadratic increase in the number of terms of the effective Hamiltonian and, consequently, in the computational cost of the procedure. Thus, experimental demonstrations have been limited to ground states, despite the practical importance of excited states.
Here, by introducing the concept of an eigenstate witness, we develop a new method that also targets excited states. A crucial limitation for the solution of the eigenvalue problem is that no method for eigenstate preparation is expected to be scalable in general (3). It remains unanswered, whether variational methods can solve particular classes of this problem. However, it is widely conjectured that eigenstates of physically relevant Hamiltonians often can be efficiently represented within an ansatz (3,4,12,20). In these cases, we estimate that the number of applications of controlled operations required to perform our algorithm increases polynomially with the size of the system.
We demonstrate this method in a proof-of-principle experiment and test its performance via numerical simulations on higher-dimensional Hamiltonians. For the implementation of the algorithm, we developed a two-qubit quantum photonic processor on a silicon photonic platform (27). This device embeds the key functionalities of on-chip entangled state generation (28)(29)(30), tomography (28,31), and arbitrary controlled unitary operations (CÛ). To perform the latter, we exploited an entanglement-based scheme (32,33).

The WAVES protocol
The approach proposed here, witness-assisted variational eigenspectra solver (WAVES), is divided into three main steps (Fig. 1A): (i) an ansatzbased variational search for the ground state; (ii) a witness-assisted variational search for excited states, starting with an initial guess obtained from the ground-state reference as outlined below; and (iii) IPEA for the accurate energy estimate of the eigenstates found.
The quantum logic circuits for WAVES are shown in Fig. 1 (B and C). The search (Fig. 1B) proceeds by preparing trial states |Y〉 T in the target register, according to the ansatz, and setting the control qubit to |+〉 C . The combined state |+〉 C ⊗ |Y〉 T is then evolved through a controlled unitary (CÛ) operation that embeds the unitary Û = e −iĤt for the evolution of |Y〉 T according to the Hamiltonian Ĥ, for a time t. The emerging control qubit state r C = Tr T (r) is then analyzed by singlequbit state tomography, from which it is possible to calculate the von Neumann entropy Sðr T Þ ¼ Sðr C Þ. The entropy acts as an eigenstate witness: It is zero if the target state is an eigenstate of the Hamiltonian. In particular, for small t, the von Neumann entropy and the linear entropy are upper-and lower-bounded by monotonic functions of the support of |Y〉 in the eigenbasis of Ĥ; that is, they are sensible measures of such support (see section S1).
This measurement of the entropy enables us to variationally target excited states as well as the ground state. The control qubit also provides an energy estimator E = −Arg[ T 〈Y|e − iĤt |Y〉 T ]/t, evaluated using the offdiagonal elements of r C . The variationally optimal ground state simultaneously minimizes the entropy Sðr C Þ and the energy estimate E. Computationally, the task of finding the ground state can be interpreted as an optimization problem using a physically motivated objective function. Here, we adopt F obj ¼ E þ TSðr C Þ ¼ E À T Tr½r C ln ðr C Þ, analogous in form to a free energy, where T is a parameter that trades off between energy optimization and entropy optimization. We call P ¼ Tr½r 2 C the purity of the reduced density matrix of the control qubit, which is used to measure the linear entropy 1 À P, here chosen as an approximation of the von Neumann entropy. We can therefore introduce the more practical objective function up to negligible constants. The optimization of F obj also permits one to identify excited states, because they occupy local minima in the high-T limit for almost all evolution times t > 0 (section S1). Defining an initial reference state |F〉 (usually obtained by mean-field approximations) and the complex vector q → as the list of parameters describing the ansatz-based state preparationÂðq → Þ, that is, jY〉 T ¼Âðq → ÞjF〉, our algorithm proceeds as follows: (1) Variationally search for the state parameters q → g that minimize the objective function F obj , thus obtaining the unitary for the ground stateÂ g ¼Âðq → g Þ.
(2) Construct a unitary for an approximate ith target excited state viaÊ piÂ ðq → g Þ, withÊ pi being a system-dependent perturbation. Variationally search for the q → e i that minimizes F obj in the high-T limit (entropy), obtaining the unitary for the target excited stateÂ ei ¼Ê piÂ ðq → ei Þ.
(3) Using the Â g for the ground state orfÂ e i gfor the excited ones in the state preparation, perform the IPEA, which further projects each state onto the closest eigenstate (34) and refines the energy estimate.
Here, we adopted a swarm optimization method for the experimental variational searches, where, for each iteration,F obj is measured for a swarm of trial states (particles), randomly sampled from a prior distribution, and used to infer a posterior with lower F obj (for more details on the optimization method, see section S2). We call each of the iterations an "epoch." The computational complexity of using our variational method to learn eigenvalues of the Hamiltonian can be quantified by the number of controlled unitary operations performed in the simulation, which depend on the optimization method used. For the particle swarm gradient-free optimization, the computational cost of sampling from the eigenspectrum of Ĥ is described by Theorem 1. A version for gradient-based methods is instead reported in the "Computational cost of WAVES for gradient-based methods" section (both demonstrations can be found in section S3).
Theorem 1. LetĤ ∈ ℂ 2 n Â2 n be Hermitian and assume that after where Ĥ|l i 〉 = l i |l i 〉 for l i ≥ 0, and that the sequence of sets of particles fXðkÞ :¼ fq Þ‖ max ≤ x max and dim(X(k)) = N ∀k. Then, the number of applications of controlled e −iĤt , for [0, p/(2‖Ĥ‖) ∋ t ∈ Q(‖Ĥ‖ −1 ), required for our particle swarm optimization algorithm to learn an eigenvalue within error D, with a probability of at least where d is the maximum error in the evaluation of F obj allowed and fD 2 m ðkÞg (fD 4 S ðkÞg) is the corresponding tolerance in the (variance of the) trace of the covariance matrix of the sample mean. Finally, we define G :¼ max k ðx max ðkÞ=D m ðkÞ; x 2 max ðkÞ=D 2 S ðkÞÞ. The above theorem implies that, in this regime, the relevant scaling parameter for iteration cost is the dimension of the parameter space. The problem of finding an appropriate ansatz is beyond the scope of this work: It is expected though to be polynomial in the number of spin orbitals for many physically relevant systems (3,4,20,21,24,35). Similarly, the number of swarm particles required (N) and the number of variational steps (N iter ) depend on both the dimension of the relevant parameter space and prior knowledge about the solution. Because the particles are moved toward the true model as the algorithm learns, N is expected to scale polynomially (36) for problems, such as chemistry, where a good ansatz and a high degree of prior knowledge are possible. These considerations lead to the implicit scaling of the number of controlled unitary applications, which is expected to increase with the number of spin orbitals (n) in the system. The number of variational parameters, together with the number of swarm particles required for these specific ansätze to achieve chemical accuracy, will likely require empirical studies to be precisely estimated. Further breakdown in the cost estimates can be considered by decomposing the controlled unitary into fundamental gates using Trotter-Suzuki or linear combination-based methods (37), but here, we ignored these issues for simplicity. If Trotter-Suzuki methods are also taken into account for the simulation, then there is a factor of roughly n 5.5 multiplied by the above costs (38).
Another fundamental point is how to choose the excitation operators used in the excited-state variational search. Consistent choices for the system-and state-specific perturbing unitariesÊ pi , required to construct the excited states, can be inferred from readily computable properties of the simulated system (25). General many-body Hamiltonians for interacting particles decompose intoĤ ¼Ĥ 0 þV , whereĤ 0 ¼ iâ i is a one-particle term andV is an interaction term. Because Ĥ 0 dominates Ĥ, a transition from the ground state to an ex-cited state can be approximated by the action of a sequence of singleexcitation operatorsâ i †â

If excitation operators based on the
Hartree-Fock approximation do not provide sufficient accuracy, then alternative approximations can be used. Advanced methods, such as multiconfiguration self-consistent field approximations (39), may ultimately be needed for hard instances with substantial electron correlations. In the "Excitation operators for chemical Hamiltonians" section, we further discuss how post-Hartree-Fock methods can be used to tackle these problems through the use of natural orbitals.

Silicon quantum photonic chip and experimental setup
The experimental demonstration of WAVES was performed on a twoqubit silicon quantum photonic processor schematically described in Fig. 2. The device is fabricated via deep-ultraviolet lithography on a silicon-on-insulator wafer with 450 nm × 220 nm single-mode waveguides. A continuous-wave (CW) pump laser in the telecom C band with an off-chip power of approximately 10 mW is coupled into the chip using polymer spot-size converters and lensed fibers, with a facet loss of approximately 8 dB. Pairs of single photons are generated in two 1.2-cm-long waveguide spiral sources through spontaneous four-wave mixing (SFWM) (27). Output photons are filtered using arrayed waveguide gratings (AWGs) with a 0.9-nm bandwidth, spectrally selecting signal photons (blue) for the control qubit and idler ones (red) for the target. The photons are detected by superconducting nanowire single-photon detectors (SNSPDs) with approximately 85% efficiency, obtaining a maximum photon coincidence rate of ≈150 Hz. Optical interferometers consisting of thermo-optic phase shifters and multimode interferometer (MMI) beam-splitters are used for photonic qubit manipulation and analysis, driven by an electronic controller with 12-bit digital-toanalog converters. The automation for the WAVES algorithm, including the control of quantum gates, the data collection, and real-time analysis, is realized by a classical computer interfaced with the quantum photonic chip. More experimental details are reported in section S5. Because of the low-power CW pump used in our experiment, multiphoton terms can be safely neglected. The use of the two spiral sources generates the Fock state ð 20 〉 þ j j02 〉 Þ= ffiffi ffi 2 p . High-visibility two-photon quantum interference (V Qua ¼ 1:00 ± 0:02) was observed in this experimental setup, as shown in the inset of Fig. 2. The photons are probabilistically split by two MMIs and then swapped by a waveguide crossing, yielding the maximally path-entangled state in the Fock basis (29,31). The state of the target photon is then expanded by adding two optical spatial modes. These additional modes represent the two components of the target qubit, which is prepared in |Y〉 T for both paths and undergoes an arbitrary Û. That is, the operation performed on the target qubiteither Î or Û-depends on which path the photon is traveling on, indicated by |0〉 P or |1〉 P . Path, in turn, is controlled by the state of the control qubit (the two qubits are entangled), |0〉 C or |1〉 C , which yields a superposition of Î and Û gates in the circuit By erasing the path information with the use of an additional waveguide crossing and two final MMIs and by detecting the signal photon and idler photon, we obtain the controlled unitary CÛ operation (32, 33, 40, 41) Note that this approach implements the CÛ gate without decomposing it into multiple two-qubit gates (11). The state preparation is realized byÂ ¼ e if a e if bŝ z =2 e if cŝ y =2 operations, whereas the Û used to map the Hamiltonian is obtained byÛ ¼ e if cŝ z =2 e if dŝ y =2 e if eŝ z =2 , where f i are phases applied by on-chip thermal phase shifters. The control qubit undergoes another single-qubit operation that can be used to perform the required operations both for tomography and for the IPEA.

Experimental results
We used the quantum photonic chip to perform WAVES, calculating the eigenspectrum of a simplified exciton transfer Hamiltonian of two chlorophyll units in the 18-mer ring of the LHII complex. We remark that this simplified model is not intended to provide an accurate description of the LHII system, yet it serves as a useful demonstration and test for our algorithm. The Hamiltonian is parameterized by a ≃ 1.46 eV, the exciton energy of a single chlorophyll unit, and b ≃ 0.037 eV, the coupling strength between the two units (42). This 2 × 2 Hamiltonian is written asĤ ¼ ða À ℓÞÎ þ bŝ x , on the basis of Pauli operators (20), where ℓ is a reference energy that can be chosen arbitrarily (see "The single-exciton Hamiltonian: Hamiltonian parameters, mapping, and eigenvalues" section). For this Hamiltonian, the perturbing unitary for the excited state corresponds tô In Fig. 3, we show the experimental results of the WAVES approach for the ground-and excited-state search. The minimization of the objective function was performed in both cases, adopting the particle swarm method outlined above. In the experiment, the energy E and purity P used to evaluate F obj were obtained by performing singlequbit tomography of the control photon. In Fig. 3 (A and B), we show the color-coded evolution of the swarm, achieving rapid convergence of the particle distribution toward the expected eigenstates of the Hamiltonian: the ground state |−〉 and the first excited state |+〉. For the ground-state search, pessimistically assuming no pre-existing knowledge of the system, the prior is initialized to uniformly span the subsection of the Hilbert space identified by the ansatz. For the excited state, instead, the search is initialized with the guessed state obtained by applying Ê p to the ground state.
As shown in Fig. 3 (C and D), within 10 to 13 search steps, F obj converges well to its minimal value, which corresponds to the ground state and excited state, respectively. Figure 3 (E and F) shows that the mean of the particle distribution achieves a high overlap with the eigenstate targeted: fidelities of 99.48 ± 0.28% with the ground state and 99.95 ± 0.05% with the excited state. All uncertainties are given by the variance of the prior distributions: A well-motivated error bar is among the amenable features derived from the adoption of a swarm optimization method.
The successful convergence of |Y〉 T is achieved by optimizing the F obj function. In particular, for the ground-state search, we used a small value of T (T = 1.25) in F obj , whereas for the excited state case, we used F obj ≡ À P (that is, high-T limit). Imperfect measurements of   Silicon quantum photonic processor. The quantum device enables one to produce maximally path-entangled photon states, perform arbitrary single-qubit state preparation and projective measurements, and, more importantly, perform any CÛ operation in the two-dimensional space. Photons are guided in the silicon waveguides and controlled by thermo-optical phase shifters. Photon pairs are directly generated inside the silicon spiral sources through SFWM, off-chip-filtered and postselected by AWG filters (not shown), and measured by SNSPDs. The generated signal (blue) and idler (red) photons are different in wavelength and form the control and target qubits, respectively. The quantum chip is interfaced with a classical computer. Inset: High-visibility quantum (blue) and classical (green) interference fringes obtained in the device using the photon sources part and configuring the top final interferometer. The high visibility is essential to verify the high-performance and correct characterization of the device. Step Step

S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
F obj , more evident in the regime close to convergence, are mainly due to experimental noise in the phase settings, given by residual thermal cross-talk among the phase shifters. The fast convergence of the algorithm, however, indicates a good robustness of the protocol to this kind of experimental noise. After the eigenstate search, the WAVES algorithm embeds the IPEA to improve the accuracy of eigenvalue estimates (3,11). In our implementation, we took advantage of circuit reconfigurability, mapping each Û 2 k directly into the chip parameters (Fig. 1C). However, in universal quantum computers, Û k can be efficiently achieved without classical precompilation by cascading k copies of Û (10). The IPEA estimated the binary fraction expansion of the eigenphase ϕ(mod 2p) for both the ground-and excited-state energies up to 32 bits (that is, a precision of 2.9 × 10 −9 eV). The normalized photon counts are shown in Fig. 3 (G and H) for all the 32 bits. This precision is higher than what is typically achievable by spectroscopic methods.
Numerical results for higher-dimensional systems We complement these proof-of-principle experimental results with a set of numerical simulations, providing insight into the performance of our approach when applied to more complex Hamiltonians. For our numerical tests, we chose a set of molecular hydrogen systems (H 2 ,H þ 3 , H 3 , and H 4 ). The corresponding Hamiltonians (up to 8 qubits) are represented in a Slater-type orbital basis (43) in the Jordan-Wigner representation (see the "Hydrogen molecules: Mapping and ansatz" section) (20) and exhibit several degeneracies in the spectrum. For each set of degenerate excited states, we will refer generically to the excited subspace they span. Figure 4 (A and B) shows the simulation results of the groundstate search and some exemplary excited-state variational searches, respectively, addressing the latter ones with a set of excitation operators of the formÊ p i . Note that this is only the first (variational) part of WAVES and that the second part (IPEA) will further project the state found into the eigenstate with a higher overlap. This feature is absent in previous VQE implementations. For the different cases, we increased the number of particles to 8, 16, 30, and 50 for H 2 , H þ 3 , H 3 , and H 4 , respectively, which follows approximately linearly the number of parameters involved in the "parameterized Hamiltonian" (PH) ansatz provided in section S7 and adopted for these simulations.
In all the cases investigated here, WAVES is able to consistently find both the ground and excited states with high fidelities (≃99% in average). The insets of Fig. 4 show that the final fidelities achieved by each variational search do not decrease (within error tolerance) when increasing the size of the Hilbert space. Although these study cases do not imply scalability of the approach, they provide an encouraging result, suggesting that to keep constant the algorithm performances with the dimensionality of the problem, a subexponential increase in number of particles and iterations is enough, provided that a polynomial parameterization applies.

DISCUSSION
We have introduced the concept of eigenstate witness and used it to develop WAVES, a new quantum method for targeting both ground and excited states of a physical Hamiltonian. We showed its proof-ofprinciple implementation on a silicon quantum photonic chip for a simplified exciton transfer Hamiltonian, obtaining its eigenstates with high fidelities and estimating the eigenvalues up to spectroscopic accuracy. Additional analysis of WAVES performances is provided by numerical simulations, where the protocol yields eigenstate estimates with high fidelity for Hamiltonians of up to 8 qubits.
All states found using the variational search, both in the experiments and in the numerical simulations, exhibited high fidelities with the target eigenstates. This preliminary refinement provides IPEA with an improved approximation of the target eigenstate, leading to an exponentially higher success probability in estimating the corresponding eigenvalue and reducing the overall complexity. Using IPEA, in addition to the variational search, allows the projection onto the eigenvectors, which is not guaranteed by the solely variational methods using a polynomial-sized ansatz. As the size of the system simulated increases, the shrinking of some energy gaps may lead to  On the x axis, we refer to the cumulative number of trial states probed (that is, the number of particles in the swarm times the variational steps). For ease of comparison, the x-axis origin has been shifted in (A) for the various cases to have equivalent fidelity for the average initial guess. Dashed lines denote average fidelities, with the shaded areas indicating a 67.5% confidence interval. The average fidelities achieved by the particle swarm optimization for both ground and excited states are calculated for 100 independent runs of WAVES. In all simulations, a binomial noise model has been taken into account when performing projective measurements. Insets: Bar charts summarizing final fidelities obtained by each search. All the simulations converged to the same high fidelity within errors, as indicated by the dashed black line in the inset. eigenstates close in energy being sensed as effectively degenerate by the variational witness, within the precision achievable by the experimental platform of choice. In these cases, the VQE-refined guess will exhibit consistent overlap with more than one eigenstate. Nevertheless, careful modifications to the phase estimation procedure may allow one to learn exponentially quickly either one of the eigenvalues belonging to almost-degenerate eigenstates (section S4). In summary, the variational search acts as a state preparation stage for the phase estimation, whereas the IPEA step addresses the shortcomings present in the variational ansatz.
From our study, for a particular ansatz (for example, UCC) and using low-order Trotter formulas, the time required by WAVES is expected to explicitly scale with the problem size as O(Mn 5.5 ), where n is the number of spin orbitals and M is the number of variational parameters, in accordance with simulation results reported by Reiher et al. (37) for systems difficult to simulate classically, such as nitrogenase. This does not automatically imply that WAVES or any other eigenstate preparation method is efficient for any Hamiltonian, because that would imply that QMA = BQP, where QMA stands for "quantum Merlin Arthur" and BQP for "bounded-error quantum polynomial time," respectively (44), which is widely believed to be false. However, an optimization based on an eigenstate witness allows the variational algorithm to address the problem of efficiently estimating an eigenvalue in the vicinity of a generic targeted state, in those cases where a polynomialsized ansatz can be provided. This problem is generally expected to be hard on classical machines, and it is challenging to solve using traditional variational methods.
WAVES offers key improvements over previous protocols. For those instances where a good ansatz is found, it can be used to locate excited states with a quantum method in a purely variational manner, in contrast to quantum-classical linear response methods (25). These methods avoid the need for additional nonlinear optimization, but this may limit their accuracy, and they do not yet use quantum phase estimation to improve the final accuracy and readout precision as in WAVES. Furthermore, one can speculate how the eigenstate witness provides an independent test of the protocol's success, detecting failure cases of convergence to local optima that do not represent a single eigenstate or excited subspace (see also section S7). These advantages come at the cost of controlling the evolution of the target register with an ancillary qubit, which is avoidable in previous VQE proposals. In addition, in WAVES, the ability to find specific eigenstates relies on the quality of the excitation operators. Further optimization on the objective function, for example, including the use of an energy penalty, can, in principle, overcome some of these limitations.
In terms of resource costs, the use of IPEA gives a quadratic speedup compared to standard VQE in estimating the energy of an eigenstate within a chosen precision. These advantages are significant given the high accuracy required in quantum chemical applications (24). Moreover, WAVES does not require lengthy adiabatic preparation of targeted eigenstates (3,9,20) or an increase in the number of terms of the implemented Hamiltonian and a precise knowledge of the spectral gap, unlike the FS method. In particular, the FS method reduces to variational optimization of a shifted and squared Hamiltonian H′ = (H − l) 2 , thus squaring both its norm and number of terms. The choice of the shift parameter l can also result in accidental degeneracies in the spectrum and dramatic closing of small spectral gaps. In the case of a problem, such as quantum chemistry, this can lead to O(n 8 ) terms formally in the Hamiltonian, drastically increasing the cost and making it cumbersome for even small instances (12,20). Moreover, the ultimate accuracy of FS methods matches the accuracy of the variational ansatz used. However, WAVES corrects essentially all of these difficulties. It retains the original norm, spectrum, and number of terms in the Hamiltonian O(n 4 ); does not depend on a shift parameter; and exceeds the accuracy of the variational ansatz used through projective phase estimation. In sections S6 and S7, the interested reader can find numerical simulations comparing WAVES with previous VQE implementations and with the FS method. The results show, in particular, that the FS method finds states with poor overlap with any true eigenstate in cases exploiting its weaknesses, whereas WAVES provides estimates exceeding 99% fidelity with a correct eigenstate in all cases tested. This direct comparison indicates higher reliability for WAVES, adding to the improvements in terms of resource costs.
In conclusion, WAVES is a new approach to tackling the search for both ground and excited states of physical Hamiltonians. The analysis performed shows that the method is expected to be scalable, under the assumption that a good ansatz can be found. The experimental demonstration on a quantum photonic chip and numerical simulations show the method performance on small-scale scenarios, indicating good noise resilience properties and better performances if compared to previous approaches. Our algorithm is, in principle, amenable to short circuit depths and leverages methods known to exhibit error robustness, thus enabling near-term experiments on non-fault-tolerant machines. By introducing new objective functions for variational algorithms, this protocol opens the way to the investigation of new methods for computing Hamiltonian spectra and represents a promising tool for future developments of quantum simulation on quantum computers.

Computational cost of WAVES for gradient-based methods
In the following theorem, the computational cost for the case of gradient-based methods is reported. Proof is given in section S3.1.
Theorem 2. LetĤ ∈ ℂ 2 n Â2 n be Hermitian and assume that after k ∈ {1, …, N iter } epochs, the state |y T (k)〉 = ∑ i a i (k)|l i 〉, where Ĥ|l i 〉 = l i |l i 〉 for l i ≥ 0. Furthermore, assume that there exists a numerical differentiation formula that evaluates ∂ q q F obj using a constant number of function evaluations on a grid of spacing h > 0 within error at most k(hL) p for positive k and L for p ∈ Q(1). Then, the number of applications of controlled e −iĤt , for (0, p/(2‖Ĥ‖) ∋ t ∈ Q(‖Ĥ‖ −1 ), required in the algorithm is in where d is the maximum error in the two-norm of the gradient of F obj allowed and D is the maximum error allowed in phase estimation of the final system with a probability of 1/2. It is then clear from the analyses contained in Theorems 1 and 2 that the particle swarm method has the potential to outperform the gradient-based method in cases where many parameters are required to describe the ansatz state and G is modest. However, the rate at which the two learn can differ substantially, because the same number of iterations may provide more or less information than the other case. In practice, gradient-based methods may be more practical to find an optimal solution in the vicinity of local optima, whereas global methods, such as our particle swarm method, may provide a better method for approaching them. Because the scaling of the Bayesian optimization approach with the number of variational parameters is better than the bounds that we prove for gradient-based optimization, we assumed that these approaches would be better in high-dimensional problems. Moreover, being inspired by ideas from approximate Bayesian inference, the latter retains part of their noise robustness. For these reasons, we focused on the particle swarm method for experiment and simulations in this work.
Finally, although both methods scale quadratically with T, in practice, the scaling will not typically be so bad. If d is chosen to guarantee fixed relative error for the process, then the cost approaches T 2 /d 2 , which is constant. This means that the quadratic scaling of T is not necessarily problematic in cases where the WAVES algorithm is optimizing for purity.

Excitation operators for chemical Hamiltonians
Our method for locating excited states variationally used approximate excitation operators to enhance the rate at which excited states may be located. Quantum chemistry has a long history of using the theory of linear response to external perturbations to approximate excited states of the system (45). The accuracy of this approximation relies on the partitioning of the total Hamiltonian intoĤ ¼Ĥ 0 þV, where Ĥ 0 is a noninteracting Hamiltonian of the formĤ 0 ¼ ∑ ij h ij b † i b j andV is an interacting perturbation. In quantum chemistry, this partitioning is often taken to beĤ 0 ¼F, where F is the Fock operator that includes one-body and averaged two-body interactions, and V is the remainder. For many systems,V is small enough such that a perturbation treatment suffices (46).
As a noninteracting Hamiltonian, Ĥ 0 may be efficiently diagonalized by a unitary transformation such thatĤ 0 ¼ ∑ i D i a † i a i , where D i are the eigenvalues of the free-fermion Hamiltonian. In this model, excited states may be formed through excitation operators of the form a † i a j acting on the ground state, where j indexes sites currently occupied by electrons and i indexes unoccupied sites. IfV is comparatively small, these eigenstates will approximate eigenstates of the true Hamiltonian, and one may refine estimates within the single-particle approximation space by diagonalizing the Hamiltonian on the basis of vectors fa † i a j Y 〉 g j , where |Y〉 is some reference state. This method is called the configuration interaction singles method. The connection may also be seen in the context of the first-order time-dependent response to an external field. For quantum computers, variations of these states may be prepared by the unitary operatorsÊ ij ðqÞ ¼ exp ½qða † i a j À a † j a i Þ. In the case of weak interactions, classical methods, such as coupled cluster, have been successful in describing the ground state; however, even low-lying excited states in these systems may exhibit correlation structures and entanglement that prevent their efficient description. This is reflected in their difficulty of simulation by current classical methods (17,18) and represents a key motivation for quantum methods, such as WAVES, to study excited states. Moreover, we stress that the single excitations here represent initial guesses for WAVES to search through correlated states not accessible to classical simulation and that these single excitations may be derived from a reduced density matrix, using a procedure described below, which was not accessible classically due to quantum correlations in the ground state. The WAVES method refines these initial guesses through optimization and then projects to a dominant eigenstate by using phase estimation.
In quantum computing, one hopes to go beyond states that are well approximated by mean-field solutions through preparation of states with nontrivial entanglement. In the VQE approach, these states are defined by the parameterization of the ansatz; however, unlike classical approaches, we may not have efficient access to full knowledge of the wave function we are preparing. In these cases, a possible approach to generating excitation operators is to look for the "closest" one-body system. This problem defines the so-called natural orbitals in quantum chemistry (47), which are the orbitals that diagonalize the one-electron reduced density matrix (1-RDM) of the prepared state, given by that may be efficiently measured on any prepared quantum state, including those with entanglement. As a symmetric positive semidefinite matrix, it may be diagonalized to yield a set of excitation operators c † i c j to approximate the excited states of the interacting system. Note that in the case of an antisymmetric product state reference, such as that generated by Hartree-Fock, these orbitals are identical to those discussed above forĤ 0 ¼F, as the canonical Hartree-Fock orbitals diagonalize the 1-RDM of a single antisymmetric product state.
The single-exciton Hamiltonian: Hamiltonian parameters, mapping, and eigenvalues Previous demonstrations of digital quantum simulation have focused almost exclusively on systems of interacting fermions such as electronic structure in molecules or the Fermi-Hubbard spin lattice model. Here, we performed numerical simulations for several such cases in section S7, reporting performances of WAVES in correctly identifying the eigenstates for the molecules H 2 to H 4 .
However, physically interesting Hamiltonians are not restricted to interacting fermions and it is important to extend quantum simulation methodologies to general systems of interacting quantum particles and quasi-particles so that quantum simulation can have an impact on a broad range of problems relevant to physics, chemistry, biology, and materials science. The spectrum of a 2 × 2 bosonic Hamiltonian was adopted for the experimental demonstration of WAVES in the main paper. We therefore required a method to convert the bosonic Hamiltonian e −iĤt into a sequence of unitary operations that can be implemented on a quantum computer. This is significant because there is not a simple analog of the Jordan-Wigner transformation that maps bosonic occupation numbers to qubits. For example, if Ĥ had a concise Pauli decomposition, then Trotter-Suzuki formulas can be used to write e ÀiĤ t ≈ e ÀiP1 t e ÀiP2 t ⋯ for Pauli operatorsP 1 ;P 2 ; …. General-purpose simulation methods can be used to express Ĥ as a sum of (at most) O(N 6 ) one-sparse matrices, provided that Ĥ does not contain interactions higher than two-body (11,20). However, these methods are ill-suited for present-day experiments, because they require a coherent implementation of a graph coloring method, which requires additional qubits.
Notwithstanding this open challenge, we selected to demonstrate our WAVES approach on the exciton transfer between two chlorophyll units found in the light-harvesting complexes of purple bacteria.
On the basis of localized excitons on each chlorophyll unit, the exciton transfer Hamiltonian isĤ where a = 1.46 eV is the energy of the exciton on one of the chlorophyll units and b = 0.037 eV is the interaction between the excitons arising from the transition dipole between the two units. The qubit representation of this two-state Hamiltonian is obtained using compact mapping (20) and iŝ whereÎ andŝ x are the usual Pauli matrices in the computational basis. The WAVES approach sequentially performs a witness-assisted variational search to find the eigenstates and a QPEA to obtain an accurate energy. At this point, we recall a well-known property of eigenfunction equations: The HamiltonianĤ ′ ¼Ĥ À ℓÎ has the same eigenstates as Ĥ and has eigenvalues l′ = l − ℓ, where l are the eigenvalues of Ĥ and ℓ is a constant. The parameter ℓ simply redefines arbitrarily the energy zero and we were free to exploit this mathematical equivalence to improve the performance of our algorithm.
In many quantum simulation applications, the natural choice of energy zero results in Hamiltonians where the total energy is orders of magnitude larger than the energy differences relevant to the phenomena under investigation. This is particularly true, for example, for reaction energies in quantum chemistry and is also the case for our excitonic Hamiltonian where we are interested in the difference between the ground and excited state 2b ≪ a. Because QPEA requires a bitwise readout of the eigenvalue of each state of interest, any shift ℓ that reduces the magnitude of the corresponding energy increases the precision that can be obtained with a given length in the QPEA binary expansion. In practice, a reasonable choice for ℓ may be obtained, for example, from a mean-field calculation, which can be performed efficiently on a classical computer. That is, such an algorithm directly estimates the correlation energy rather than the ground-state energy. To mimic a realistic problem where mean-field theory provides a rather poor guess for the exact eigenvalues, we selected an arbitrary value of ℓ ≃ 1.24 eV in the experiment.
The energy estimation in WAVES adopts the form E = −Arg[ T 〈Y|e −iĤt |Y〉 T ]/t, and this imposes restrictions on the value of t to avoid issues due to the 2p periodicity of the Arg function. This is a limitation already known from QPEA, normally addressed by choosing t small enough to prevent the algorithm from providing any eigenvalues mod 2p (11).
However, in the WAVES protocol, additional boundaries for t emerge from considerations about the P estimator, as described in section S1.2. The span in purity within the accessible Hilbert space also dominated the choice of the evolution time t = 26. It is also easy to verify that 26(l g − l e ) ≠ 0 mod 2p; therefore, our choice satisfies all the conditions stated for t, concerning the value of P in the objective function.

Hydrogen molecules: Mapping and ansatz
In addition to the experimental verification described, we report numerical simulations of chemical Hamiltonians using classical com-puters in section S7, which are partially shown in Fig. 4. In particular, we had simulated ground and excited electronic states of H 2 , H þ 3 , and H 4 in a STO-3G basis (43) in the Jordan-Wigner representation (20) to assess the scalability of our proposal. These represent four-, six-, and eight-qubit Hamiltonians, respectively. In our investigations, we looked at two different ansatz. First, we used a PH ansatz, where we took and allowed variation of the terms t ij and t ijkl to define the ansatz. The a † i and a j represent creation and annihilation operators in the Hartree-Fock basis, respectively. Variation was performed after transformation to Pauli operators via the Jordan-Wigner transformation. In all cases, the reference state on which Uðt → Þ acts was taken to be the Hartree-Fock state with the correct number of particles. This is essentially a deformation of the original Hamiltonian, allowing one to preserve its symmetries and giving a natural connection to the original interaction structure of the problem. We also used an unrestricted UCC ansatz of the form The key difference between this ansatz and the previous one is that excitations between arbitrary orbitals are allowed, not just those found in the Hamiltonian. The consequence of this is that one may create or repair symmetry-broken states that have been produced by some other means, allowing additional flexibility in the description of the state at the cost of more parameters. In the following, we will refer synthetically to a parameterization of the ansatzÂðq → Þ, corresponding tô Similarly, the approximate excitation operators used were defined in this basis as where we take j to index the occupied orbitals of the Hartree-Fock reference and i to index the occupied orbitals of the reference.

SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/4/1/eaap9646/DC1 section S1. Notes on the objective function section S2. Swarm optimization algorithm section S3. Complexity analysis of the variational protocol section S4. Phase estimation without quantum collapse section S5. Experimental details section S6. Robustness against experimental noise of the variational search section S7. Numerical simulations of hydrogen molecules fig. S1. Median inference error in eigenvalue inference with a = 1/2 for a distribution with two randomly chosen eigenvalues and likelihood function (37) is used. fig. S2. Schematic representation of the experimental setup. fig. S3. Numerical simulations of the variational ground state search robustness for the bosonic one-qubit Ĥ against gate infidelities using F obj or E alone. fig. S4. Synopsis of numerical simulations of excited states searches for the molecules H 2 and H þ 3 . fig. S5. Numerical simulations of the WAVES variational search for the synthetically truncated PH ansatz are studied in the molecular hydrogen systems (H 2 , H þ 3 , H 3 , H 4 ). fig. S6. Convergence of the WAVES algorithm to a subspace of excited states for different hydrogen systems. fig. S7. Comparison between different ansaetze adopted in the search for excited states in the H þ 3 system. fig. S8. Behavior comparison of the first part of WAVES and an equivalent implementation of the FS method when applied to the initial guess provided by the Ê p 3 excitation operator for the H 2 system. table S1. Summary of ansätze used for simulations in the main paper and in the Supplementary Materials for the various systems investigated, along with the cardinality of their parameterization, dim ðq → Þ. table S2. Summary of possible situations occurring in numerical simulations when WAVES is performed with different ansätze and excitation operators, targeting a certain excited subspace E t from an initial guess |Y 0 〉. References (48)(49)(50)(51)(52)(53)(54)(55)(56)