Real time evolution for ultracompact Hamiltonian eigenstates on quantum hardware

In this work we present a detailed analysis of variational quantum phase estimation (VQPE), a method based on real-time evolution for ground and excited state estimation on near-term hardware. We derive the theoretical ground on which the approach stands, and demonstrate that it provides one of the most compact variational expansions to date for solving strongly correlated Hamiltonians. At the center of VQPE lies a set of equations, with a simple geometrical interpretation, which provides conditions for the time evolution grid in order to decouple eigenstates out of the set of time evolved expansion states, and connects the method to the classical filter diagonalization algorithm. Further, we introduce what we call the unitary formulation of VQPE, in which the number of matrix elements that need to be measured scales linearly with the number of expansion states, and we provide an analysis of the effects of noise which substantially improves previous considerations. The unitary formulation allows for a direct comparison to iterative phase estimation. Our results mark VQPE as both a natural and highly efficient quantum algorithm for ground and excited state calculations of general many-body systems. We demonstrate a hardware implementation of VQPE for the transverse field Ising model. Further, we illustrate its power on a paradigmatic example of strong correlation (Cr2 in the SVP basis set), and show that it is possible to reach chemical accuracy with as few as ~50 timesteps.


I. INTRODUCTION
In fulfilling the promise of quantum computation [1][2][3] and enabling the exact solution of many-body quantum systems, numerous algorithms of different resource requirements have been proposed [4][5][6][7], which require quantum and classical resources of different complexity.Many of these algorithms are focused on efficient eigenvalue extraction, important for solving problems in chemistry [8], physics [9], and materials science [10] and limited classically by the exponential scaling of Hilbert space with system size.Though immense progress has been made in the development of quantum algorithms for eigenvalue estimation, the resource requirements remain prohibitively high with regards to noisy intermediate-scale quantum (NISQ) hardware [10,11].
For example, quantum phase estimation (QPE) [12][13][14] is considered an algorithm that will need considerable quantum resources to run, but will ultimately be a highly accurate approach to quantum simulation.Adiabatic state preparation [15,16] allows the ground state of a particular Hamiltonian to be reached by preparing an initial ground state of a simpler system and slowing changing the Hamiltonian to the desired system, requiring long coherence times and low gate errors.Particularly in the current era of NISQ quantum computers, the variational quantum eigensolver (VQE) framework [17], and its non-orthogonal variant NOVQE [18], are promising approaches for the exact solution of many-body quantum systems.However, the common formulation of this family of methods relies on the solution of a highly complex variational optimization problem on classical computers, which remains an open challenge [19,20].
The methods described above are generally used to solve the time-independent Schrödinger equation to determine the Hamiltonian eigenvalues and eigenstates.However, time evolution is a more natural operation on a quantum computer and thus simulation of the time-dependent Schrödinger equation is a more ideal framework to implement.Given the intrinsically quantum-mechanical relation between the time and energy domains [21], a different family of quantum algorithms focuses on using a time-dependent perspective to solve timeindependent problems [22,23].These methods propose a linear wave function Ansatz expanded in time evolved states and solve the thereby defined generalized eigenvalue equation classically, while the Hamiltonian and overlap matrix elements are measured on quantum hardware.Exploiting this, quantum computers hold unique potential to outperform their classical counterparts with algorithms based on real-time evolution.In this work we focus on that advantage and use realtime evolution to generate a basis of states to extract Hamilto-arXiv:2103.08563v3[quant-ph] 7 Apr 2021 nian eigenvalues.
Using real time evolution to generate a basis of states to solve a Hamiltonian is not a new idea [24][25][26], but it is not widely used in classical simulation due to the computational limitations of simulating real time evolution.Approximate imaginary time evolution [27,28] and Krylov diagonalization methods [29] are far more widely used in classical simulation, and the intuition behind such approaches is simple to understand: Each new state generated in these approaches has a larger overlap with the true ground state.This of course is not how real time evolution works, as the expectation value of the energy remains constant, and one never gets closer to the ground state during the evolution.However, the states that can be generated through real time evolution in fact do provide a basis from which one can extract ground and excited states.This is not only highly efficient, but in some cases may be faster than other quantum methods that use time evolution, such as QPE.Thus the main goal in this work is to develop an approach for computing ground and excited states, using states generated by real time evolution, that is as fast and efficient as possible.With this in mind, we analyse the theoretical underpinnings of a class of algorithms we term variational quantum phase estimation (VQPE) which is based on real time expansion methods [22,23].We use the term VQPE because of its relationship to both QPE and VQE, as we detail below (see sections II A and II G).
Our VQPE algorithm and analysis goes beyond previous proposals in the following ways: We reduce the number of quantum measurements needed to be linear instead of quadratic in the number of expansion states.We introduce the phase cancellation conditions, providing the underlining theory for why this approach works and deriving a direct link between time steps and band gaps in the spectrum.We also analyze the effects of noise on our convergence properties, providing significantly improved intuition for ill-conditioning of VQPE methods.Further, we demonstrate the method classically for several weakly and moderately correlated molecules as well as a strongly correlated transition metal dimer, Cr 2 .We show that for all the systems, regardless of the level of electronic correlation, less than 50 real-time evolved expansion states are needed to reach agreement within chemical accuracy for ground state energies obtained with state-of-theart classical methods requiring O(10 6 ) variational parameters [30].Additionally, we describe and implement the algorithm on quantum hardware for the transverse field Ising model.Since real-time evolution is natural to implement on quantum hardware, this approach holds immense promise for NISQ implementation.
The paper is structured as follows: in Section II, we analyze the theoretical structure of VQPE.This starts with a brief overview and intuition of existing methods corresponding to the VQPE family (II A, II B).Subsequently, we present a novel examination of the theoretical underpinnings of the method based on the phase cancellation picture in II C and we propose a procedure for choosing optimal time step sizes.With the practical implementation of VQPE on NISQ devices in mind, we present the details of our unitary formulation of VQPE in Subsections II D and II E. We conclude the theory section with a careful analysis of the effects of noise II F, a comparison between VQPE and QPE II G, and a discussion of the inclusion of other time evolved states into VQPE II H. Section III provides details of the systems studied as well as the classical and quantum simulation methods.Section IV summarizes and discusses the results of the simulations.Concluding remarks are found in Section V.

II. THEORY
A. Landscape of Existing Variational Algorithms VQE approaches are highly relevant for the NISQ era of quantum computation [4].They comprise fairly simple quantum circuit implementations at the price of relying on the solution of a high dimensional, classical optimization problem in the presence of noise.Despite the optimization challenge, many of the currently existing examples of actual quantum simulations for many-body physics correspond to implementations of this algorithm.The basic premise of VQE relies on the variational approximation: A parametrized wave function Ansatz |Φ( α j ) , where α j are the variational parameters, is chosen such that it can be efficiently implemented on a quantum computer.The energy expectation value of this Ansatz is evaluated on quantum hardware, and then the optimization problem ∇E( α) = 0 is solved on a classical computer.As in any variational approach, the efficacy of the approximation depends on the flexibility of the Ansatz.A way to increase this flexibility is to choose a more general ground state estimate, namely as a linear combination of several parametrized expansion states |Φ j ( α j ) .The total wave function Ansatz thus becomes where c j are the expansion coefficients and α j the optimization parameters of the j-th expansion states.Applying the variational principle to the coefficients c alone leads to the secular equations [31] j where ε I is an estimate for the I-th Hamiltonian eigenvalue E I .Equation 2 is a standard generalized eigenvalue equation, which can be solved classically.The Hamiltonian, which we assume to be time independent, and overlap matrix elements are measured on quantum hardware in the "basis" of expansion states following Thus, the expansion coefficients c can be determined by classically solving the noisy, generalized eigenvalue problem in Eq. ( 2).The expansion states |Φ j ( α j ) themselves can be then optimized with a classical minimization method, further improving the energy estimates.One example in which this has been used recently is with chemically motivated Ansätze, such as unitary coupled cluster expansions [32], which had nonetheless some difficulties related to the optimization of parameters.
The optimization of the α j parameters is an open field due to (i) the high dimensional nature of the optimization problem and (ii) the presence of noise, which is necessarily part of any approach on quantum hardware.To alleviate these problems, an alternative framework has been recently explored [5,22,23,[33][34][35][36][37], which completely bypasses the need for optimization routines.In this family of methods, one does not employ a set of parameterized expansion states |Φ j ( α j ) , but instead generates a set of expansion states |Φ j systematically from one or several reference states |Ψ I .The only variational parameters left are thus the expansion coefficient c, and consequently the only task to be performed by a classical computer is solving the (noisy) generalized eigenvalue problem in Eq. ( 2).
The way of generating these expansion states should balance ease of implementation on quantum hardware with creating a flexible expansion set, in a variational sense, to obtain accurate energies.While a priori, by changing from parametrized |Φ j ( α j ) to systematically generated |Φ j , we are reducing the variational flexibility of our Ansatz, the expansion state generation can still be performed in such a way that it is natural to both the description of ground and excited states, as well as to the implementation on a quantum computer.One such approach is the quantum subspace expansion (QSE) method [33,35,36], which, after optimizing a ground state estimate with regular VQE, generates expansion states by applying single excitation operators on top of this VQE reference.A recent variation introduced a more general multireference Ansatz, targeting ground and exited states simultaneously at the optimization step [34].
Alternatively, the expansion states can be formed by applying the time evolution operator U (t) = e −iHt to the reference states.This is the approach followed in the QLanczos method [5,37], where the time evolution is performed along imaginary time (i.e.t → iτ ), and the quantum filter diagonalization and quantum Krylov approaches [22,23], where the time evolution is performed along real time.In the case of the QLanczos method, it is clear that evolving to large enough imaginary times will provide an expansion set which is well suited to describe the ground state, provided that the reference state is not orthogonal to it.In the next subsections, we discuss the formalism behind using a set of the real-time states in the expansion, which we refer to as VQPE.We also analyze VQPE's robustness to noise, which is critical to asses its applicability in NISQ devices.

B. VQPE -Basic Intuition
In the VQPE approach [22,23], the expansion states |Φ j,I are generated from the reference states |Ψ I through time evo-lution as If there are N R reference states, and N T time steps are considered, this produces a "basis" of N R (N T + 1) expansion states, where |Φ 0,I = |Ψ I .For simplicity, we will consider a single reference state |Ψ 0 for now, and will discuss the use of multiple reference states further below.Throughout the paper, we set the reduced Planck constant = 1.On a first glance, it seems counter-intuitive that the set of expansion states in Eq. ( 4) would improve the ground state estimate given by the reference state |Ψ 0 .After all, the |Φ j,0 states all have the same energy expectation value.Stair et.al. [23] suggest an interpretation based on short time evolution, in which Taylor expanding Eq. ( 4) shows that the expansion states span the same space as the Krylov vectors H j |Ψ 0 , making it equivalent to power methods such as the Lanczos algorithm [29].This justifies why a set of expansion states concentrated in a time grid over a short time scale should produce a good variational Ansatz for the ground state.Including multiple reference states then should provide for good and stable approximates for the first few excited states as well, in the spirit of the band Lanczos method [38].This interpretation suggests that the VQPE method should work best for short time evolution.We want to complement this interpretation with a more general one, not limited to short time steps, although this still remains the most interesting regimes from an implementation perspective.In particular, the VQPE approach is reminiscent of the computation of the autocorrelation function g(t) = Ψ 0 |e −iHt |Ψ 0 , which contains the full spectral information of H [21].In g(t), the ground state information is encoded in the long time limit, rather than the short time one, since the ground state evolves with the slowest frequency ω 0 = E 0 in units where = 1.Thus, there should be nothing particular about the short time evolution limit.Instead we articulate the precise requirements for the implementation that will lead to accurate energy estimates.This is related to the linear independence of the expansion states which poses a lower bound to the optimal time step size.We also note that VQPE is closely related to the classical filter diagonalization method [24][25][26]39], as pointed out by Parrish et.al.[22].
We thus want to depart from the Krylov intuition, and provide a different, hopefully more general heuristic, which we will then rigorously formalize, to understand why VQPE approaches should provide good ground and excited state approximations.We begin by writing out the decomposition of the reference state |Ψ 0 into Hamiltonian eigenstates |N , such that H |N = E N |N .This can be written out explicitly as where ψ 0 N = N |Ψ 0 are the coefficients of the reference state |Ψ 0 in the eigenbasis of H.We will refer to those Hamiltonian eigenstates |N for which ψ 0 N is beyond some nonnegligible threshold as the support space of state |Ψ 0 with respect to (w.r.t.) H.This decomposition gives for the expansion states The equation above just states the obvious: each component of |Ψ 0 in the support space w.r.t.H evolves with its own frequency.This, however, makes transparent why the VQPE method can work: Choosing the time grid {t j } accordingly, it is possible to make linear combinations of the expansion states |Φ j,0 such that the different phases e −iE N tj cancel out targeted components of |Ψ 0 along the support space state |N .In this way it is possible to "extract" eigenstates in the support space of |Ψ 0 by including enough expansion states |Φ j,0 .Note that this is not exclusive to the ground state, nor is this limited to short time scales t j .The only requirement is given by the number of eigenstates of H in the support space of |Ψ 0 , defining how many time steps are needed for perfect state extraction, and by the energy gaps (i.e.relative frequencies) of those states, which govern the phase cancellation conditions.
In essence, VQPE allows one to extract "the most out of the reference state", in the sense that if there are Q states in its support space, it should be possible to produce Q time evolved states from which to reconstruct the corresponding Q Hamiltonian eigenstates, by solving the secular equation Eq. (2).Of course, this presumes that it is possible to produce Q linearly independent time evolved states, and that our time evolution is noiseless and performed at arbitrary numerical precision.For general reference states, the size of the support space will be too large in general to recover all eigenstates, but a modest amount of these should be enough to approximate the lowest lying energy eigenstates in it.

C. Phase Cancellation Conditions and Relation to
Filter-Diagonalization We now formalize the phase cancellation heuristic on a solid mathematical footing.To this end, we derive a set of equations, the phase cancellation conditions, which set sufficient conditions to exactly extract the Hamiltonian eigenstates from the support space.These conditions embody the intuition in terms of auto-correlation functions described before, and are effectively discrete versions of the main relations at the heart of the classical filter-diagonalization approach [24][25][26].
We consider the overlap matrix S j,k in Eq. (3) for the expansion states in Eq. (6).It is convenient to write the overlap matrix in operator form in the span of the expansion states |Φ j,0 , and it is easy to verify that This means that the operator corresponding to the overlap matrix projects onto the span of the expansion states [40].
Substituting Eq. ( 6) into the overlap operator gives, after some minor reordering of terms In the equation above, Q is the number of Hamiltonian eigenstates in the support of |Ψ 0 , and we can ignore Hamiltonian eigenstates outside the support space due to their small coefficients ψ 0 N .Now, we can define the phase cancellation conditions (PCCs) as These are the Q(Q − 1)/2 conditions for the N T + 1 time steps in the time grid.Given that the support space is spanned by just Q vectors, it seems that the PCCs impose stricter conditions on the time grid than absolutely necessary to recover the full support space.Still, they embody mathematically the phase cancellation heuristic which we have discussed above.Indeed, the condition in Eq. ( 9) enforces the cancellation of the time evolved phase between all Hamiltonian eigenstates in the support of |Ψ 0 , and can be represented graphically as a sum of phases in the unit circle.When the phase cancellation conditions are fulfilled, the overlap operator simplifies into a weighted projector into the support of |Ψ 0 w.r.t.H, namely weighted by the coefficients of the reference state |Ψ 0 on the support space, c.f. Eq. (5).In this case, the expansion states span exactly the same space as the Hamiltonian eigenstates in the support space, and solving the secular equation (2) returns the exact eigenstates and eigenvalues of H.The PCCs in Eq. ( 9) can be understood as the discrete limit of the eigenstate extraction through Fourier transform of a time evolved state exploited in the classical filter diagonalization literature [24].Further, in the limit where long-time evolutions are used the phase cancellation conditions will also be approximately satisfied with high probability as N T tends to infinity.To see this, let t j = (j + 1)/ω where ω −1 is a uniform random variable on [0, 1/ min(E N − E M )] := [0, ∆E −1 min ] (where M = N ).First we have that When solving the secular equation, we choose the same threshold for the SVD decomposition of the overlap matrix.The vertical dashed line marks the 15-th time step, after which we have as many expansion states as vectors in the support space.The large subplot corresponds to the smallest time step which fulfills the phase cancellation condition Eq. ( 9), which reduces to a single condition for this Hamiltonian.The insets in each subfigure correspond to a geometric representation of the phase cancellation condition, with each phase e −it j ∆E a point in the unit circle on the complex plane.See text for details.
Intuitively it is reasonable to expect that if the mean is small then with high probability the PCCs should hold approximately.In order to demonstrate such a concentration for the oscillating functions that we use here we, however, need to also bound the variance.
Thus from Chebyshev's inequality we have that with high probability j e −ij(E N −E M )ω −1 will be within of the expectation value.Thus the phase cancellation condition's error for the N, M component is in Thus the value of N T needed to ensure that the PCC holds within error at most (with high probability) obeys Here O(•) denotes an asymptotic upper bound with multiplicative polylogarithmic factors neglected.Thus an approximate solution to the phase cancellation conditions will generically hold for a gapped system.We exemplify the previous theory on the example of a Hamiltonian of linear spectrum E N = N ∆E, akin to a harmonic oscillator, in Fig. 1.In this case, the PCCs in Eq. ( 9) can be fulfilled exactly by a linear time grid t j = j∆t P with the perfect time step size ∆t P defined as Indeed, it is easy to check that in the case of a linear spectrum, a linear time grid with time step size given by Eq. ( 16) fulfills the PCCs exactly after N T = Q − 1 time steps.This can be accomplished with a single time step size since in the case of a linear spectrum the PCCs effectively reduce to a single condition.This can be seen in the rightmost panel of Fig. 1, where exactly after 15 time steps the first four eigenvalues of the secular equation match the exact eigenvalues to the maximal precision.This precision is determined by the singular value threshold, s SV , introduced into the general eigenvalue problem, which truncates the singular values of the overlap matrix.In Fig. 1 this precision corresponds to s SV = 10 −12 , i.e. midway between double and single machine precision.This threshold also determines the support space size Q.The reference state in all examples in Fig. 1 is defined as |Ψ 0 ∝ N e −E N |N , excited states being exponentially suppressed.The support space is then defined as the Hamiltonian eigenstates with squared coefficients in |Ψ 0 above s SV .The inset in this panel represents the PCCs graphically, as the phases of all eigenstates in the support space perfectly span the unit circle, thus cancelling each other.The smaller panels on the left of Fig. 1 show the outcome of choosing a time step size differing from ∆t p .Small time step sizes are shown in the upper two subfigures, and would be the natural choice from the Krylov interpretation of VQPE [23].These clearly show a significantly slower convergence than the perfect time step derived from the PCCs, which is easy to explain from the phase distribution on the unit circle in the insets.Only once we cover the unit circle close to homogeneously, thus approximately fulfilling the PCCs, can we extract all eigenstates essentially exactly after the minimal number of time steps (see lower left panel in Fig. 1).Due to the periodic nature of the complex phases in Eq. ( 9), large time steps can result in as poor approximations as short ones, as shown in the lower right panel in Fig. 1.The worst such longer time step sizes correspond to particular integer multiples of ∆t P , namely zQ∆t P , where z is an integer.For these time step sizes, the PCCs in Eq. ( 9) cannot be fulfilled, even approximately.These are very particular time steps, and thus, for the linear spectrum, a randomly chosen time step ∆t ≥ ∆t p is still likely to provide good results for the Hamiltonian eigenstates.
For general spectra, a single time step size ∆t in a linear grid is unlikely to fulfill all Q(Q − 1)/2 PCCs exactly.From the above analysis, a valid strategy would be to choose a time step size ∆t and number of timesteps N T such that we sample a full period of the slowest oscillation in the support space.This is given, e.g., by the minimal energy gap ∆E min if we are interested in all excited states contained in the support space, and the ground state gap ∆E 1,2 = E 2 − E 1 if we only need an estimate of the ground state.However, in practical implementations it is advantageous, and sometimes necessary, to limit the total simulation time in order to minimize the error.At the same time, we have to choose a time step size large enough such that each new state is linearly independent from the previous ones.Otherwise no new information is added and the variational Ansatz is not improved (see Fig. 1, upper left hand panel, where the energy decreases in a step-like fashion).Reconciling these two notions, we propose the following systematic approach: 1. Choose a small enough time step size such that the energy convergence is step-like.
Step-like convergence refers to the situation in which adding a new expansion state, i.e. propagating for an additional time step, does not improve the variational Ansatz, resulting in the same (or slightly worse) energy estimates as before including the new step (c.f.upper left panel in Fig. 1).This happens when the inclusion of the new expansion state produces an overlap matrix which has no additional singular value over the threshold s SV .
2. Perform the VQPE algorithm using the previously identified small time step size, until the first expansion state resulting in an improvement of the energy estimates is produced.
3. Plotting the lowest eigenvalue of the previous VQPE simulation as a function of the propagation time will result in a nearly horizontal line.This plateau ends after the addition of the final expansion state, which does improve the energies.The length of this plateau defines a new, larger time step size, which can be used in a new VQPE simulation.
4. Repeat simulation with the new time step size.If there is still step-like convergence, go back to 3. Otherwise use this as your simulation time step.
This procedure underlies the fact that for devising practical implementations of VQPE, the guiding principle should be to generate linearly independent expansion states rather than to exactly fulfill the PCCs, which nonetheless are a useful perspective for the theoretical analysis of the algorithm.This is the strategy we adopt in the results section below.

D. Towards an Optimal Implementation -Toeplitz Structure of S
Besides being natural to implement on quantum hardware, and presenting the interesting phase cancellation structure described above, VQPE approaches show a further theoretical advantage: when using a linear time grid t j = j∆t, the Hamiltonian and overlap matrices in Eq. ( 3) have a restrictive structure, which formally reduces the number of measurements that should be needed to solve the generalized eigenvalue problem.Indeed, as pointed out by Parrish et.al. in Ref. [22], using the real time expansion set these matrices become Toeplitz, meaning that e.g. S j,k = S j+1,k+1 .In particular, the concrete expressions read where we have only used the fact that a time-independent Hamiltonian commutes with itself at all times, relying thus exclusively on time translational symmetry.From Eq. ( 17), it follows that we can reconstruct the 2 (N T + 1) × (N T + 1) matrices by measuring only 2(N T + 1) overlaps in total.
Unfortunately, as pointed out in Ref. [22], the Toeplitz property of the Hamiltonian matrix is lost in actual quantum hardware implementations, if the time evolution operator is Trotterized.In those cases, the commutativity of U (t j ) with H is lost, and thus we either need to evaluate all (N T + 1) 2 Hamiltonian matrix elements separately or transition to a higher-order Trotter formula that better approximates the commutation relations.Nevertheless, the same is not the case for the overlap matrix.As long as the expansion states |Φ j,0 are constructed using a linear grid with a unitary time evolution approximation the Toeplitz condition will prevail.Here, U a (∆t) is the Trotterized time evolution operator and ∆t is the common time step size of the linear time grid, such that the approximate expansion states obey which in the limit of an exact time evolution operator recovers Eq. ( 6).Now, if U a (∆t) is unitary, the overlap matrix of the approximated expansion states |Φ a (∆t) will clearly be Toeplitz, since Given that gate operations in quantum hardware are naturally unitary, this means that it is always possible to guarantee the Toeplitz condition of the overlap matrix, simply by choosing a linear time grid.This is true of course, for the often invoked first order Trotterization [41] approximation to the time step evolution U (∆t), which is indeed unitary.

E. Towards an Optimal Implementation -Unitary Formulation
It is possible to rewrite this generalized eigenvalue problem in a simpler form, exploiting the particular relationship between the Hamiltonian and overlap matrices in Eq. ( 17), essentially formulating it equivalently to the classical filter diagonalization problem found in signal processing [26].This proves to be the ideal formulation of VQPE for quantum computation.
The main insight relies on substituting the Hamiltonian in the secular equation Eq. ( 2) by the time evolution operator U (∆t) = e −iH∆t .This operator is effectively isospectral with the Hamiltonian, indeed the eigenstates |N of H fulfill It is important to note that unlike the Hamiltonian, the time evolution operator is not Hermitian, but unitary, having thus complex eigenvalues of unit modulus.We can therefore write a secular equation for the time evolution operator U (∆t) where the overlap matrix, eigenvalues ε I , and expansion coefficients c I j are the same as in Eq. ( 2), and the time evolution matrix elements follow, in the single reference implementation, (23) To transform from eigenvalues of Eq. ( 21) to Eq. ( 2), ∆t must also be small enough that we can distinguish a physical E N value from its unphysical periodic images E N ± 2π/∆t.
From Eq. ( 2) to Eq. ( 22), we have simply reformulated the VQPE problem into an equivalent generalized eigenvalue problem with a unitary matrix.The key simplification for the implementation on quantum hardware relies on the realization that the time evolution matrix elements in Eq. ( 23) have the same structure as the overlap matrix elements in Eq. (17).Thus, choosing again the time grid {t j } to be linear, i.e. t j = j∆t, the time evolution matrix elements coincide with the overlap matrix elements as The last equality is a manifestation of the Toeplitz structure.Thus, according to Eq. ( 24), for linear time grids there is no need to measure the time evolution matrix explicitly, since it can be recovered from the measurements for the overlap matrix plus an additional measurement involving an extra expansion state |Φ N T +1,0 .In this way, exploiting Eq. ( 24) and the Toeplitz structure of the overlap matrix, the number of measurements reduces from 2(N T +1) 2 to just N T +2.Further, as shown in the previous subsection, the Toeplitz structure prevails when implementing the time evolution operators with a unitary approximation, such as first order Trotterization, making the reduction in number of measurements applicable for real implementation on quantum hardware.Intuitively, the unitary formulation of VQPE is the quantum algorithm equivalent to measuring the autocorrelation function g(t) = Ψ 0 |e −iHt |Ψ 0 and analyzing its Fourier spectrum.The overlap matrix elements are essentially sampling g(t) at different points, and one can approximate the underlying spectrum once enough samples are obtained.Thus, we are fundamentally expressing the VQPE algorithm in its most natural language, that of autocorrelation functions.In this work, we implement both the traditional and unitary formulations of the VQPE secular equation.

F. Diminishing the Effect of Noise through Singular Value Decomposition
In the previous subsections, we have briefly reviewed the theoretical formalism of VQPE approaches, in the new light of the phase cancellation interpretation, but without considering the effects of noise.We turn our attention now to how the presence of noise, comprising both finite numerical precision on the classical computer and measurement uncertainty from the quantum hardware, limits the final accuracy of the VQPE results.We consider systematic noise, due to a priori uncontrollable or unavoidable sources, plus any remaining statistical uncertainty after repeated measurements.To this end, the phase cancellation formalism will simplify the analysis.We will restrict ourselves for simplicity to the single reference implementation, but generalizing our conclusions to the multi-reference case is straightforward.
As shown in Eq. ( 4), VQPE generates a series of N T states from a reference |Ψ 0 defined by a time grid {t j }.In a noiseless simulation, any given time grid is more likely than not to produce a set of N T linearly independent vectors.For them to be linearly dependent requires the following determinant to vanish exactly This equation is one constraint on N T unknowns {t j }, which is generically satisfied by an (N T − 1)-dimensional manifold of {t j } embedded in R N T .For a linear grid t j = j∆t, the choice of time step size ∆t will generically cause linear dependencies on a subset of R with measure zero.For example, ∆t = 2πn E2−E1 causes a linear dependency in the case of N T = 2. Thus, with the exception of Hamiltonians with a restricted spectrum such as E j = j∆E, it seems safe to assume that in almost any time grid chosen, a noiseless simulation will generate N T linearly independent vectors.Since all the expansion states share the same support space [42], a noiseless simulation with Q steps, Q being the size of the support space, should recover all eigenstates exactly.We will assume a linear time grid with time step size ∆t henceforth.This ideal notion stops holding the moment we consider noise, both from numerical and measurement origins.Noise can for example make states close to linearly dependent, and thus introduce errors in the eigenvalues ε I of the secular equation.We will quantify noise introducing the parameter .Two measured or computed values α, β are only distinguishable if |α − β| > .Noise becomes important, for example, in the small time step size limit.When ∆t∆E min is small, where ∆E min is the minimal spectral gap in the support space, the first expansion steps will produce states that are only marginally different to the reference |Ψ 0 .These will not improve the variational Ansatz if where we have recovered Planck's constant to make the units clear.If Eq. ( 26) is fulfilled, the magnitude of the difference between the expansion state and the reference will fall below the noise threshold, making the new expansion state U (∆T ) |Ψ 0 useless from a variational perspective.This is the reason behind the step-like decreasing behavior in the small time step panels of Fig. 1.A finite thus determines a minimal time step ∆T > ∆Emin .
There will be cases where it is hard to generate precise expansion states with the minimal time step size ∆T required to offset a given noise level.It becomes thus important to prune the Hamiltonian and overlap matrices of numerical and measurement noise.This can be done by means of a singular value decomposition (SVD) of the overlap matrix: Neglecting all singular values bellow some threshold, which should be larger than the magnitude of the noise.In the case of measurement error, this noise scales as 1/ √ M where M is the number of samples.In the majority of this work, we have thus conservatively chosen a threshold of 10 −1 , corresponding to > 100 samples.We used such a truncation already in the results shown in Fig. 1.This singular value truncation produces effectively a new but smaller expansion basis of N SV D elements.As a consequence, the number N T is not the significant measure of how much information is collected in the expansion set, and instead N SV D ≤ N T becomes the measure to follow.Only when N SV D = Q will the secular equations recover the exact support space spectrum.
We exemplify this on the harmonic spectrum in Fig. 2, where in the upper panels we show the relative noise error for the first four eigenstates as a function of number of expansion states, introducing Gaussian noise N (0, ) of standard deviation on the Hamiltonian and overlap matrix elements.We choose the singular value truncation threshold s SV to be between 10 −1 and 1, and in the lower panels we plot the singular values of S j,k as a function of the number of expansion states, marking s SV as a dashed line.In each simulation, we choose as time step size ∆t the optimal time step in Eq. ( 16), considering a possible support space of 16 elements, regardless of s SV [43].In a noiseless simulation, this choice of time step would result in an optimally compact number of effective expansion vectors N SV D , which equals the number of actual expansion vectors until the maximal number N max SV D = Q is reached, after which all eigenstates would be resolved accurately.The presence of statistical noise has two consequences: on the one hand, the asymptotic accuracy decreases with increasing noise variance.As mentioned above, this type of statistical noise can be reduced by sampling.
As a second effect of statistical noise, not all Q eigenstates in the support space are resolved after exactly Q steps, since the corresponding singular values of the overlap matrix fall bellow the truncation threshold s SV (see lower panels of Fig. 2).As long as a given singular value falls bellow s SV , the corresponding eigenvalue cannot be extracted from the generalized eigenvalue equation, which is represented in the upper panels of Fig. 2 by horizontal straight lines.Thus the noise limits what states can be extracted from the reference state |Ψ 0 , by setting the minimal singular value truncation threshold s SV .Those Hamiltonian eigenstates with smaller absolute coefficient squared than s SV cannot be resolved.How- ever, examining Eq. ( 10), we observe that the singular values of the overlap matrix are enhanced linearly with increasing number of expansion states N T , i.e. with increased number of time steps in the VQPE approach, once the PCCs are reasonably fulfilled.Thus, it should be possible to extract eigenstates with reference state components below the error threshold by increasing the number of time steps.We exemplify this in Fig. 3, again on the harmonic spectrum example with exponentially suppressed initial state.Because of this choice of starting state, it takes an exponentially large number of extra time steps to resolve every new eigenstate, but it is in principle possible.The ideal strategy is of course to propose a reference state with large overlap with the eigenstate of interest, but this discussion shows that it is possible to extract states beyond the dominant one accurately.Once enough time steps have been produced, any singular value of the overlap matrix can be made to increase above s SV , the horizontal dashed line in the lower panel of the figure.Of course, this is not unexpected, since the initial starting signal on each excited state is exponentially small by construction.When a uniform overlap is used as a starting state, then all the SVD values are the same.
We want to address another notion that has been brought up with regards to the effect of noise in VQPE simulations: the overlap matrix condition number n(S) [23].As shown in [23], typical eigenvalue problems arising in VQPE have extremely large condition numbers, which may suggest high sensitivity to noise.This is expected for Hartree-Fock starting states, which ideally have exponentially small overlap with most of the Hilbert space, leading to very large condition numbers.Still, this can be dealt with by performing an SVD of the overlap matrix and truncating the singular values below s SV .For our application, this reduces the number of linearly independent states we can resolve within our error threshold, thus reducing the best possible accuracy of the results.This can be remedied be improved by including additional time evolved states, as discussed above.
In the Supplementary Information, we investigate the role of n(S) for the eigenvalue accuracy in the VQPE approach.Our results show that, upon pruning the expansion space from the singular values of the overlap matrix below the threshold s SV , we consistently obtain accurate eigenvalues in the presence of noise even with matrices of large condition number.This is in itself not surprising, since the singular value truncation is effectively proposing an auxiliary generalized eigenvalue problem with smaller condition number.
At this point, we can return our attention to the notion of multi-reference implementations of the VQPE, such as the multi-reference quantum Krylov method of Stair et.al. [23].In this work the use of several reference states is proposed in order to reduce the condition number of the Hamiltonian and overlap matrices in the expansion space, at the cost of requiring a larger number of expansion states for the same ground state accuracy.Using the phase cancellation picture, we argue that the worsened ground state energy convergence is due to two distinct, cooperating factors: the increased size of the total support space, and the smaller number of expansion vectors in each individual support space.In the multi-reference formulation of VQPE, each state |Ψ I has its own support space w.r.t.H, which we refer to as individual support spaces, the union of these forming the total support space of the implementation.The individual support spaces will be in general distinct from each other.Clearly, the larger total support space allows for a more flexible variational ansatz, from which it is possible to extract more Hamiltonian eigenstates than in the single reference case.However, this comes at the price of requiring more expansion states to perform the phase cancellation procedure to purify individual eigenstates.From our results, performed in classical simulations with noise, and on actual noisy quantum hardware, the larger condition numbers do not result in large errors in the eigenvalue estimates ε I , and thus we conclude the condition number alone should not be a reason to employ multi-reference VQPE implementations.However, in some excited states simulations a multi-reference approach might accelerate covergence, in the same way that band Lanczos improves normal Lanczos in this regard [29,38].
The previous considerations hold for statistical errors but a more careful analysis needs to be performed for systematic errors in the implementation of the Hamiltonian dynamics.For example, in the case of a Trotterized Hamiltonian, the reference states are created under the evolution of a Hamiltonian different than the one of interest, which will result in errors in the eigenvalue estimation which cannot be reduced through sampling.

G. Comparison between unitary VQPE and QPE
Here we compare the unitary formulation of VQPE in Eq. ( 21) to conventional QPE in the general case of a multidimensional support space.For the special case of a 1dimensional support space, there are adaptive variants of QPE that use one ancilla qubit and achieve Heisenberg-limit measurement [44].Recently developed methods [45,46] devise QPE variants that use one ancilla qubit and are suitable for larger support spaces (which is also the case for VQPE).In an adaptive approach where |Ψ 0 is the ground state, a different t would be chosen for each measurement to maximize the extraction of information about E 0 rather than performing multiple measurements to estimate the expectation value of g(t) = Ψ 0 |e −iHt |Ψ 0 = e −iE0t for a single choice of t.
Quantum phase estimation (QPE) is a natural algorithm to compare VQPE against; however, there are a wealth of different phase estimation algorithms known in the literature and further some applications of phase estimation can even be used in concert with VQPE.Our aim in this section is to compare and constrast different flavors of phase estimation to VQPE and also show how QPE can be used to accelerate learning the expectation values of the VQPE circuit through amplitude estimation.
There are broadly two categories of phase estimation algorithms, iterative phase estimation and Fourier-based phase estimation.Fourier-based phase estimation is perhaps the best understood approach to performing phase estimation.An advantage of this approach is that it is known to precisely achieve optimal scaling of the uncertainty with the number of applications of the underlying unitary (i.e. it saturates the Heisenberg limit [47,48]).The optimal approach to Fourierbased phase estimation deviates slightly from traditional approaches by using an optimized initial state which deviates from the Fourier state typically used in older approaches.Specifically, the input state is taken to be an m-qubit state of the form This state is chosen to minimize an estimate of the circular variance, known as the Holevo variance, of the eigenphases of the unitary e −iHt that results from the phase estimation protocol.Next let us define notation for a controlled directional evolution below † This notation is useful here because in quantum phase estimation for the simulation of electronic structure, the cost of performing a controlled-directional evolution as seen above is approximately the same as performing an ordinary controlled evolution.Applying phase estimation with a circuit of this  form provides twice as much phase per application of the unitary as a controlled-U application would experience.The uncertainty from this phase estimation circuit, as quantified by the Holevo variance, is approximately , where m is the number of qubits used in the QPE register.This is precisely the Heisenberg limit.
Of course, error in the simulation dynamics needs to be considered in any such simulation.For simplicity we will consider using the lowest-order Trotter-Suzuki formula within the phase estimation protocol.In the limit where the uncertainty in the eigenphase is small relative to π, the Holevo variance approximately corresponds to the variance.Then since the Trotter-Suzuki error is uncorrelated with the phase estimation error, the total error in the phase estimate that arises from using a timestep of duration t is ≈ π This suggests that if we choose both error sources to be t/(2 √ 2π) then the overall error will be .This corresponds to The error in the Trotter-Suzuki formula is difficult to estimate a priori and for the purposes of simplifying the comparison we will choose the symmetric Trotter formula.Using the results of Proposition 16 from [49] we have that if H = j α j H j then the error in the symmetric Trotter formula is upper bounded by This bound can be further simplified using the fact that each of the Hamiltonian terms H j in quantum chemistry is of norm 1. Therefore let us denote S to be the set of all tuples (p, q, j) such that [H q , [H p , H j ]] = 0. Then we can further bound T S (t) ≤ 2t 3 3 j p>j q>j |α p α q α j |δ (q,p,j)∈S Under the assumption that the error from QPE is equal to the Trotter error we then find that a sufficient choice for the Trotter-step is Therefore the value of m needed in the phase estimation step is Thus if M is the number of terms in the Hamiltonian then the maximum number of operator exponentials that need to be simulated in the circuit to achieve Holevo variance is This process needs to be repeated a number of times to find an eigenstate j which occurs with a probability of |ψ 0 j | 2 .This implies that the total number of repetitions needed to reconstruct the entire spectrum with high probability is in O(1/ min j |ψ j | 2 ) and thus Note that the performance of this algorithm can be improved using optimized simulation methods and the use of amplitude amplification, but we forgo these optimizations here since they are less appropriate for experiments in the NISQ era.Now let us consider the analogous problem for estimating the corresponding eigenvalue using VQPE.The first step involves learning the matrices H j,k and S j,k .In order to do so, we need to apply e −iHtj for all j considered.There are of course a host of quantum simulation algorithms that can be employed to perform this simulation.For simplicity, we will consider the second-order Trotter-Suzuki formula.Consider (32).From this expression we have that we can perform e −iHtj using O(M (Γt max ) 3/2 / √ T S ) operator exponentials, where |t j | ≤ t max for all t j .Thus we can use this circuit and the Hadamard test to compute the real and imaginary components of φ j |φ k within error at most with high probability (using the Chernoff bound to justify the high-probability statement) using Similarly, using the arguments of [50] we have that the number of operator exponentials needed to approximate the mean energy within the same error requirements is in Next, using the fact that the induced 2-norm is at most N T times the max-norm, we have that if S and H are the approximated matrices then we need to take → /N T in order to ensure that the total error in the reconstructed matrix (as measured by the operator norm) is at most .This implies that the number of exponentials needed to reconstruct both matrices within the required error budget is in Then from the Hellman-Feynman theorem we have that the eigenvalues of both matrices lie within an -neighborhood of the original eigenvalues.Next using standard matrix inequalities for the error in the matrix inverse [51] (under the assumption that is less than the minimum eigenvalue gap) yields that the error in any eigenvalue reconstructed from (2) has error in O( H S −1 2 ).Thus if we wish the error in the reconstructed eigenvalues to be at most then we finally need to take → / H S −1 2 .This, combined with the bound that H ≤ j |α j | leads to a final complexity scaling of the number of exponentials needed to reconstruct the eigenvalues within error at most with high probability is in .
(41) Asymptotically, this analysis suggests that VQPE may be favorable in cases where due to the potential benefits that could arise from inferring having to sample low probability events using phase estimation.However, such asymptotic analysis does not conclusively show which method will be superior as both are upper bounds that are potentially loose.For this reason, numerical studies such as the one undertaken in this paper, are essential for gauging the actual performance of these algorithms.
For Toeplitz matrices such as in Eqs. ( 17) and ( 24) to be diagonal in a Fourier basis, they must have additional circulant matrix structure satisfying which corresponds to an aliasing condition whereby all energy eigenvalues E N are an integer multiple of a base energy, An accurate approximation of this condition requires a very small value of ∆t and thus a large number of ancilla qubits in the QPE circuit and large N T .There is a formal equivalence between QPE and unitary VQPE when this condition is exactly satisfied, although this occurs beyond the intended small-N T operating regime of any VQPE variant.Practically, VQPE and QPE remain distinct in that VQPE uses measurements of one ancilla qubit to statistically estimate overlap matrix elements, while QPE uses measurements of many ancilla qubits to sample from approximate eigenpairs.Besides requiring fewer ancilla qubits, the main benefit of VQPE over QPE in the small-N T regime is its exact diagonalization of a generalized eigenvalue problem rather than relying on approximate diagonalization within a Fourier basis.The overlap matrix of this eigenvalue problem has a numerical rank corresponding to the size of the support space, and this reduced rank is not reliably exposed by the diagonal elements of the overlap matrix in the Fourier basis.As a result, VQPE achieves a higher effective energy resolution than QPE for small N T values and an even smaller support space by exploiting an SVD-based projection of the eigenvalue problem into a numerically relevant subspace.
In Fig. 5, we plot two different resource (for LiH discussed further in Section III) estimates that exemplify the near-term benefits of VQPE and the long-term benefits of QPE, here shown as the Heisenberg limit.VQPE is able to achieve substantially higher accuracy than idealized QPE for the same number of time steps per circuit, thus it is better suited for utilizing near-term hardware with shorter circuit depths and fewer qubits available for use as ancillae.However, QPE utilizes its deeper circuits to achieve Heisenberg-limited energy resolution, which is more efficient in overall run time for achieving high accuracy.But, as with modern parallel algorithms, it is clear that each matrix element can be calculated independently and thus can be run as a parallel algorithm over multiple quantum computers.The expected QPE run time is also amplified by | Ψ 0 |0 | −2 ≈ 1.02 to account for the probability of collapsing into a state other than the ground state, which is a negligible overhead in this example.Other than finite-sampling costs and errors, this example does not consider other sources of error such as Trotterization or quantum circuit noise, which can further degrade the accuracy of both QPE and VQPE.For example, Fig. 6 depicts a ∆t-dependent error floor induced by Trotterization errors (for the transverse field Ising model defined in Eq. ( 45) and discussed further in Section III) that remains overwhelmed by sampling errors in VQPE.With shorter circuit depths and fewer ancilla qubits, VQPE is less susceptible than QPE to these extrinsic sources of error.

H. Inclusion of Other Time-Evolved States
In this work we focus on real-time evolution.However, in practice any time-evolved state shown in Fig 7, real or imaginary, could be included in the expansion set for solving the generalized eigenvalue equation.Moreover, any unitaries f (H; j) that commute with the Hamiltonian would be viable candidates to produce an expansion set with the following important property: all the states in the expansion set share an

Accuracy
Figure 5. Accuracy comparison between VQPE and the Heisenberg limit for the ground state energy of LiH.Varying numbers of samples per expectation value are shown for VQPE alongside the idealization of VQPE without finite-sampling errors.We use SVD cutoffs of 2, 0.9, and 0.3 respectively corresponding to 10 3 , 10 4 , and 10 5 samples and 10 −6 for the ideal case.All methods use ∆t = 0.5.The top and bottom plots show the the maximal evolution time (i.e.how many time steps need to be included for energy convergence), related to the circuit depth, and the total evolution time, i.e. the sum of all time segments needed.identical support space w.r.t the Hamiltonian.Thus any expansion |Φ j,I = f (H; j) |Ψ I formed by any such unitary will eventually cover the whole support space without introducing more eigenstates.Recent algorithms such as QLanczos [5,37], quantum filter diagonalization [22], quantum Krylov approaches [23], and quantum power methods [53,54] have all proposed using real or imaginary time evolution in various ways to generate states for creating a wave function Ansatz.The real/imaginary time evolution plane is illustrated in Figure 7, along with the states used in these various algorithms.As long as the support space is retained in the time evolution and the states are linearly independent, efficient eigenvalue extraction will be possible.There are sev- eral benefits derived from using only the real-time axis which include only a linear number of measurements for setting up the generalized eigenvalue problem, almost no optimization of parameters for determining the time evolution, and a well conditioned problem for extracting information about the ground and low lying excited states.

III. METHODOLOGICAL DETAILS
We turn now to presenting concrete numerical data for Hamiltonians from condensed matter physics and quantum chemistry.For practical applications, the most important quality of the VQPE approach is doubtlessly its compactness.By this, we refer to the number of variational parameters required to obtain accurate eigenvalues.The number of variational parameters is (N T + 1)N R , where N R is the number of reference states (for all examples shown below, we use a single reference N R = 1), and we choose chemical accuracy (∼ 1.6 mHa) as our target accuracy.Thus, we apply VQPE to several many-body systems of different complexity, taken from the fields of quantum chemistry and condensed matter physics: several small molecules, including weakly correlated < l a t e x i t s h a 1 _ b a s e 6 4 = " l 6 a l c m g h k 8 n U 6 j 5

Ground State
< l a t e x i t s h a 1 _ b a s e 6 4 = " 5 E P i A d f P 2 f l 8 8 M 1 S 2 x c t y w q t

Filter Diagonalization
< l a t e x i t s h a 1 _ b a s e 6 4 = " A P R R p d J 8 c t 6 6 M T A + w 0 Q S A 8 j q e n 4 = "

Eigenstate Projection
< l a t e x i t s h a 1 _ b a s e 6 4 = " y d L i 6 A p m 6 3 W r z Z L 8 9 D 7 a T k n Z X c 2 9 N i + T K r I 0 c O y C E 5 J h 4 5 J 2 V y T S q k S j h J y T N 5 J W / O k / P i v D s f 0 9 E F J 9 v Z I 3 / k f P 4 A O X i V p Q = = < / l a t e x i t > Initial State Figure 7.A schematic picture of recent eigenvalue algorithms relying on time evolution proposed for quantum computers.The plane corresponds to the complex time plane, horizontal axis being real time, the vertical axis imaginary time.The colored dots denote expansion states |φj,I .The operator f (H; j) denotes any unitary that commutes with the Hamiltonian, and thus the expansion states formed by applying this unitary will remain in the support space and can be used to extract eigenstates.The real time axis is where the phase cancellation picture applies.On the imaginary time axis, the coefficients of the excited states are suppressed as time increases, leading to the ground state.The dotted lines indicate that short time imaginary and real time states will be very similar, as argued in the short time Krylov method [52].We emphasize that combinations of states generated on the entire real/imaginary time plane have essentially the same support space and can be used as a basis to solve the generalized eigenvalue problem.ones, LiH (3-21g basis set [55], 11 orbitals and 4 electrons), H 2 O (cc-pVDZ basis set [56], 24 orbitals and 10 electrons), N 2 (cc-pVDZ basis set [56], 28 orbitals and 14 electrons), two moderately correlated, linear H 6 (STO-6g basis [57], bond length 1.5 Å, 6 orbitals and 6 electrons) and C 2 (cc-pVDZ basis set [56], 28 orbitals and 12 electrons), and one strongly correlated Cr 2 (def2-SVP basis set [58], 30 orbitals and 24 electrons), and the transverse field Ising model (2 site on hardware, 10 site on simulator).We simulate the VQPE calculation for the molecular systems classically, while running the Ising model calculations on IBM's quantum hardware and simulator.Further details, such as molecular geometries, can be found in the Supplementary Material.
For the classical simulations, we perform time evolution using exact diagonalization (ED) dynamics based on the Lanczos algorithm [59][60][61].In the case of all molecular systems except for Cr 2 , the small basis sets allow us to perform the time evolution including all electrons and all orbitals explicitly, while for Cr 2 we restrict the simulation to the widely studied 30 orbital active space.Still, only for LiH and H 6 is it computationally feasible to perform exact dynamics.For all other systems, we must truncate the Hilbert space, considering only a subset of all Slater determinants in the corresponding active spaces.This presents an approximation, and we perform a finite-size effect study by comparing the dynamics of progressively larger truncations, from one thousand to one million determinants, which we show for Cr 2 , the most complicated   of the systems considered.The determinants of each truncation are chosen using the adaptive sampling configuration interaction (ASCI) algorithm [30,62,63].This is an iterative selected configuration interaction approach, which explores the Hilbert space and identifies the most important determinants for a ground state approximation, providing highly accurate yet moderately sized truncations.Thus, by optimizing the truncated spaces with respect to the ground state description, we can reliably infer the full Hilbert space limit from our finite size calculations.In practice, we determined a one million determinant Hilbert space truncation with ASCI for Cr 2 (10 For all molecular systems, we determined the optimal time step size ∆t using the strategy described in the theory section (see Fig. 9), and use the Hartree-Fock state as reference state |Ψ 0 .For these classical simulations, we did not exploit the unitary formulation of VQPE, instead solving the generalized eigenvalue equation 2.
For the simulations on IBM's quantum hardware and simulator, we consider the transverse field Ising Model (TFIM) with open boundary conditions, defined by the following Hamiltonian where X, Y , and Z above and henceforth denote the Pauli operators, J = 1 corresponds to the spin coupling and h = 2 the external field.In these simulations, we chose a time step size ∆t = 0.05 and approximated the time evolution operator via first order Trotterization [41].Each run was performed with 8192 shots.For all hardware data shown, |Ψ 0 readout error mitigation was performed.
As discussed in the theory section, the Toeplitz property of the overlap matrix means that the number of measurements is linear with the number of timesteps N T + 2. Further, we take advantage of the unitary formulation of VQPE to measure only the real and imaginary components of the overlap matrix explicitly, not the Hamiltonian matrix elements.To measure the overlap matrix on quantum hardware, we follow the approach described in [64], where one additional ancilla qubit is required, as shown in Fig 8 .The price of this simplicity is the requirement for an efficient construction of a controlled time evolution operator.We do this with the use of QSearch [65] for circuit synthesis, where a unitary, U , is read in and a circuit with minimal CNOT depth is constructed while keeping the similarity with the synthesized unitary (U s ) below some threshold σ (here σ = 10 −10 is the Hilbert-Schmidt inner product between the conjugate transpose of U and U s ).For the 2 qubit transverse field Ising model with one ancilla, a total of 6 CNOTs are needed, regardless of time step [66].The overlap matrix elements are obtained by measuring S 0j = X j + i Y j on the ancilla qubit for each timestep, where j denotes the number of timesteps after which we make the measurement.As noted in the theory section, because this has the structure of a time correlation function, the other rows of the matrix can be filled out accordingly, with one extra measurement needed for the unitary formulation of the generalized eigenvalue equation.

A. Finding the Optimal Time Step Size
Here we exemplify the strategy for finding optimal time step sizes on LiH. Figure 9 shows the convergence of the ground state (solid lines) and first excited state included in the support space (dashed lines) for several time step sizes ∆t, as a function of the number of expansion states (left) and the total simulation time (middle).The left panel in Fig 9 shows how the step-like behavior appearing for the smallest time step size (∆t = 0.05) decreases with increasing ∆t.The convergence as a function of simulation time (middle) clearly shows the implicit balance to be made between simulations with short total time and small number of expansion states.The ∆t = 2.0 simulation (green curves) achieves chemical accuracy with the smallest number of expansion states, i.e. with the most compact variational Ansatz.However, this comes at the cost of a significantly larger total simulation time, at least twice as large as the simulations with smaller time step sizes.For practical implementations on NISQ devices, the marginal decrease in number of expansion states is not worth the drastic increase in total simulation time.This underlies the importance of finding a time step size that is long enough, but not too long.
We further note that the excited state convergence is worse than the ground state in general, requiring for LiH about twice as many expansion states and total simulation time.This is likely related to the magnitude of the expansion coefficients of the Hartree-Fock reference vector in the basis of Hamiltonian eigenstates.
The rightmost plot of Fig. 9 shows the total simulation time for ground state convergence to chemical accuracy using ∆t = 0.2, for different singular value cutoffs s SV .Clearly, the larger the value of s SV , the longer the simulation time becomes, which relates to the notions of linear independence, which we elaborated upon in the theory section (see Eq. ( 26)).

B. Classical Simulation of Molecular Systems
Figure 10 shows the power of the VQPE method on molecular systems with varying degrees of electronic correlation, using a singular value cutoff of 10 −1 .As discussed in III, for all systems except for H 6 and LiH, we use a truncated Hamiltonian obtained with ASCI.In panel (b) we show that this approximation has no effect on our conclusions, since vastly different truncation sizes present the same convergence even though the actual ground state energies are very different.Panel (a) in this figure shows the number of variational parameters to reach chemical accuracy using ASCI (a classical selected configuration interaction approach), and the VQPE algorithm with different s SV values.These parameters correspond to determinants in ASCI, and to expansion states in VQPE.Given the differences in complexity with respect to the classical algorithm, it is remarkable how many of these systems converge to chemical accuracy at essentially identical rates, as shown in panel (c).Note that the differences in convergence times are likely connected to the different levels of correlation in these systems: the existence of low lying excited states (i.e. more strongly correlated) means longer convergence time because of a smaller energy gap.Even with these differences, for all systems studied we reach chemical accuracy with as little as 50 variational parameters and about 30 a.u. of total simulation time.This exceptionally low number of variational parameters is particularly striking when considering that classically, e.g. with ASCI, it is necessary to include 10 5 − 10 6 parameters (determinants) to reach the same energies.
The inset in panel (c) shows the convergence of the first excited state in the support space.Here it was necessary to use a significantly smaller SVD cutoff (s SV = 10 −6 ) to converge in reasonable simulation times.Even with this drastic restriction in the allowable noise, the excited state energies converge significantly slower than the ground state energies.Still, they should be possible to access efficiently with reference states tailored to describe excited states [67][68][69][70][71].

C. Quantum Simulation of Transverse Field Ising Model
Finally, Fig. 11 shows our results using IBM's QPUs (2 qubit TFIM) and qasm simulator (10 qubit TFIM), compared to the exact ground state energy.Our initial state corresponds to a ferromagnetically ordered chain where all spins point down in the Z direction, which corresponds to the |000... state in the computational basis.For processing the hardware data we use s SV = 2 to account for larger noise.Both the simulator and hardware runs used the same Trotter and VQPE step size ∆t = 0.05.Error bars on the hardware data correspond to the variance of the average ground state energy at each time point over 10 individual hardware runs on IBM's Paris machine, using qubits 0, 1, and 2. Each hardware run used 8192 shots per matrix element.Figure 11 shows convergence to the exact ground state energy for both simulator and hardware data.The saw-tooth behavior is indicative of linear dependence in the expansion states, as discussed and shown in Sections II C and IV A.  V. CONCLUSIONS Here we analyze the theoretical underpinnings of VQPE algorithms as well as the effects of statistical noise on its performance.We have presented evidence for VQPE as a particularly compact and natural algorithm for the NISQ era, showing that with the unitary formulation of the generalized eigenvalue equation and the prevailing Toeplitz structure of the overlap matrix, only a linear number of measurements are needed.We provided a heuristic for choosing an optimal time step size, which balances two opposing factors: On the one hand, generating linearly independent expansion states with each new time evolution, such that the variational Ansatz be-comes compact, which sets a lower bound to the time step size, and on the other hand minimizing the total simulation time, as required by NISQ hardware.We have also demonstrated the effect of statistical noise on the optimal time step and final accuracy of the energies, showing that simple regularization techniques suffice to mitigate the effects of noise.
We have exemplified the power of the VQPE approach on a wide range of molecules of different complexities, simulating the algorithm classically, as well as the transverse field Ising model on IBM's quantum simulator and hardware.On a traditional example of strong correlation, Cr 2 , our results suggest that VQPE achieves comparable accuracy to state of the art classical simulations with orders of magnitude fewer varia- tional parameters (∼ 50 vs. 10 6 ).This compactness, together with their NISQ compatibility (with respect to total simulation time, number of measurements, and noise resilience), marks VQPE approaches as some of the most promising platforms to realize the long sought quantum computation goal of performing many-body simulations beyond the reach of classical computations.

Figure 1 .
Figure 1.Relative error of the first four eigenvalues from the real time NOVQE secular equation for a Hamiltonian of linear spectrum EN = N ∆E (∆E = 0.75 here) and different time steps ∆t as a funciton of the number of time steps.The reference state follows |Ψ0 ∝ N e −E N |N , such that only the 16 Hamiltonian eigenstates of lowest energy are part of the support space, choosing a coefficient threshold of 10 −12 .When solving the secular equation, we choose the same threshold for the SVD decomposition of the overlap matrix.The vertical dashed line marks the 15-th time step, after which we have as many expansion states as vectors in the support space.The large subplot corresponds to the smallest time step which fulfills the phase cancellation condition Eq. (9), which reduces to a single condition for this Hamiltonian.The insets in each subfigure correspond to a geometric representation of the phase cancellation condition, with each phase e −it j ∆E a point in the unit circle on the complex plane.See text for details.

Figure 2 .
Figure 2. Relative error and corresponding singular values of the overlap matrix for Hamiltonian of linear spectrum with different noise values.Upper panels: Relative error of the first four eigenvalues from the VQPE secular equation for a Hamiltonian of linear spectrum EN = N ∆E (∆E = 0.75 here) and perfect time step ∆tP , including Gaussian noise N (0, ) to the Hamiltonian and overlap matrix elements.We choose the singular value truncation threshold sSV to be at least two orders of magnitude larger than the noise standard deviation .The reference state follows |Ψ0 ∝ N e −E N |N , the effective support space size determined using the singular value truncation threshold N max SV D = − 1 2∆E ln(sSV ) − 1. Lower panels: Corresponding singular values of the overlap matrix as a function of time.The horizontal dashed line represents the threshold sSV .Note that, until a given singular value is larger than sSV , we cannot extract the corresponding eigenvalue from the generalized eigenvalue problem.This is represented by the horizontal lines in the upper panels.See text for details.

Figure 3 .
Figure 3. Relative error and corresponding singular values of the overlap matrix for Hamiltonian of linear spectrum for a large enough number of time steps to extract eigenstates below the error threshold.Upper Panel: Relative error of the first five eigenstates from the VQPE secular equation for a Hamiltonian of linear spectrum EN = N ∆E (∆E = 0.75 here) and perfect time step ∆tP , including Gaussian noise N (0, 10 −2 ) to the Hamiltonian and overlap matrix elements.We choose the singular value truncation threshold sSV = 9•10 −1 .The reference state follows |Ψ0 ∝ N e −E N |N .Lower Panel: Corresponding singular values of the overlap matrix as a function of time.The horizontal dashed line represents the singular value threshold sSV .The time step at which each singular value becomes larger than sSV is marked with a vertical dashed line, connecting upper and lower panels.See text for details.

Figure 4 .
Figure 4. Four qubit version of asymptotically optimal phase estimation circuit.
v g q s y I o s u i C C 4 r 2 A e 0 p W T S 2 z Y 2 k x m S O 2 I Z u n P j r 7 h x o Y h b f 8 G d f 2 P 6 W G j r g c D h n H u 5 O S e I p T D o e d 9 O Z m F x a X k l u 5 p b W 9 / Y 3 H K 3 d y o m S j S H M o 9 k p G s B M y C F g j I K l F C L N b A w k F A N + p c j v 3 o P 2 o h I 3 e I g h m b I u k p 0 B G d o p Z a 7 3 0 B 4 w z Z y P O k c l z w T w v e z U m + e D G t I 0 v 2 y A E 5 I j 4 5 I 0 V y T U q k T D h 5 J M / k l b w 5 T 8 6 L 8 + 5 8 T E Y z z n R n l / y B 8 / k D K f a Z e g = = < / l a t e x i t > g y X o x 3 4 2 P a m j O y m V 3 0 B 8 b n D x b 1 l E s = < / l a t e x i t > f (H; t j ) = e iHt j < l a t e x i t s h a 1 _ b a s e 6 4 = " H j 2 L d 8 4 + N k U r a m 9 L O c a + / e I b P zE = " > A A A B / n i c b Z B L S w M x F I U z P m t 9 V c W V m 2 A R X J U Z U X R Z d K O 7 i v Y B 7 V A y 6 W 0 b m s k M y R 2 x D A X / i h s X ir j 1 d 7 j z 3 5 i 2 s 9 D W A 4 G P c + 4 l y Q l i K Q y 6 7 r e z s L i 0 v L K a W 8 u v b 2 x u b R d 2 d m s m S j S H K o 9 k p B s B M y C F g i o K l N C I N b A w k F A P B l f j v P 4 A 2 o h I 3 e M w B j 9 k P S W 6 g j O 0 V r u w 3 0 J 4 t e x i t s h a 1 _ b a s e 6 4 = "+ f Z x f 0 v O G 6 p U d U 6 R Q e M Q q H + W e q 4 = " > A A A B 6 n i c b V B N S 8 N A E J 3 4 W e t X1 a O X x S J 4 q o k o e i x 6 8 V j R f k g b y m a 7 a Z d u N m F 3 I p T Q n + D F g y J e / U X e / D d u 2 x y 0 9 c H A 4 7 0 1 b L 2 w U N 7 e 2 d 3 Z L e / s N E 6 e a 8 T q L Z a x b A T V c C s X r K F D y V q I 5 j Q L J m 8 H w Z u I 3 n 7 g 2 I l Y P O E q 4 H 9 G + E q F g F K 1 0 3 z p 9 7 J b K b s W d g i w S L y d l y F H r l r 46 v Z i l E V f I J D W m 7 b k J + h n V K J j k 4 2 I n N T y h b E j 7 v G 2 p o h E 3 f j Y 9 d U y O r d I j Y a x t K S R T 9 f d E R i N j R l F g O y O K A z P v T c T / v H a K 4 Z W f C Z W k y B W b L Q p T ST A m k 7 9 J T 2 j O U I 4 s o U w L e y t h A 6 o p Q 5 t O 0 Y b g z b + 8 S B p n F e + i 4 t 6 d l 6 v X e R w F O I Q j O A E P L q E K t 1 C D O j D o w z O 8 w p s j n R f n 3 f m Y t S 4 5 + c w B / I H z + Q P P 4 4 1 8 < / l a t e x i t >X/Y< l a t e x i t s h a 1 _ b a s e 6 4 = " 4 6 k S 4 2 x Q / p t N q B I s r p X L L j F Q b h U = " > A A A B 8 X 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 0 W P R i 8 c K 9 g P b U D b b S b t 0 s w m 7 G 6 H E / g s v H h T x 6 r / x 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S A T X x n W / n c L K 6 t r 6 R n G z t L W 9 s 7 t X 3 j 9 o 6 j h V D B s s F r F q B 1 S j 4 B I b h h u

Figure 8 .
Figure 8. Hadamard test circuit, which computes the real and imaginary parts of the overlap matrix, where |Ω denotes the system and the top qubit is the ancilla.Here U = e −itH .

Figure 9 .Figure 10 .
Figure 9. Ground (solid) and excited state (dashed) energy of LiH as a function of timestep (left) and total simulation time (middle).The initial reference vector is the Hartree-Fock state.We use the 3-21g basis set which has 4 electrons and 11 orbitals.The SVD cutoff is 10 −1 for both the left and the middle plot and the legend is the same for both.The right figure shows the total time to reach chemical accuracy as a function of SVD.

Figure 11 .
Figure 11.The top figure shows the ground state energy as a function of timestep for the 10 qubit transverse field Ising model run on IBM's qasm simulator with sSV = 10 −1 .The bottom figure shows results for hardware runs for a 2 qubit transverse field Ising model.For both the top and bottom plots, VQPE and Trotter step size are the same, ∆t = 0.05.The QPU results were run on IBM's Paris machine.All classical processing in the bottom plot (for both QPU and simulator data) used SVD cutoff sSV = 2.
5determinants for H 2 O, C 2 and N 2 ), and built all smaller truncations of size N from this one, by picking the N determinants with the largest coefficients in the one million (10 5 for H 2 O, C 2 and N 2 ) determinant ground state wave function.