Variational quantum algorithms for discovering Hamiltonian spectra

Calculating the energy spectrum of a quantum system is an important task, for example to analyse reaction rates in drug discovery and catalysis. There has been significant progress in developing algorithms to calculate the ground state energy of molecules on near-term quantum computers. However, calculating excited state energies has attracted comparatively less attention, and it is currently unclear what the optimal method is. We introduce a low depth, variational quantum algorithm to sequentially calculate the excited states of general Hamiltonians. Incorporating a recently proposed technique, we employ the low depth swap test to energetically penalise the ground state, and transform excited states into ground states of modified Hamiltonians. We use variational imaginary time evolution as a subroutine, which deterministically propagates towards the target eigenstate. We discuss how symmetry measurements can mitigate errors in the swap test step. We numerically test our algorithm on Hamiltonians which encode 3SAT optimisation problems of up to 18 qubits, and the electronic structure of the Lithium Hydride molecule. As our algorithm uses only low depth circuits and variational algorithms, it is suitable for use on near-term quantum hardware.


I. INTRODUCTION
Many physical properties of a quantum system are determined primarily by its energy spectrum. Diagonalisation of the Hamiltonian allows one to calculate various expectation values and correlation functions [2]. For example, the energy spectra of molecules inform their dynamics and therefore an understanding of such spectra is vital for molecular design [3]. But diagonalising the Hamiltonians of quantum systems on a classical machine is an often intractable task. The exponentially growing cost of storing and operating upon the quantum system makes diagonalising large systems prohibitively expensive. This precludes, for example, the study of complicated compounds [4].
It is widely believed that quantum computers will make these classically intractable molecular simulations possible [5]. This was formalised by Aspuru-Guzik et al., who suggested using the adiabatic state preparation and phase estimation algorithms to find the ground state energy of molecules [6]. Such a method necessitated deep quantum circuits and therefore long coherence times. The recently proposed variational quantum eigensolver (VQE) circumvents these requirements, exchanging them for an increased number of circuit repetitions [7,8]. To date, there have been several proof of principle experiments which have applied the VQE to find the ground state energy of small molecules [7,[9][10][11][12]. Other variational algorithms have been introduced which can simulate the real [13] or imaginary [14] time dynamics of quantum systems. In particular it was shown that imaginary time evolution can be used as an alternative to the VQE to find the ground state of molecular Hamiltonians.
While the ground state problem has received significant attention, the problem of finding the excited states of molecular systems has experienced comparatively less development. This is despite its particular importance in analysing chemical reactions, which is a vital ingredient in the quest to discover new drugs and industrial catalysts [15].
There have thus far been a handful of proposals for calculating excited states, all based on the VQE, such as the quantum subspace expansion method [16,17] and the von-Neumann entropy method [18]. These methods require either many measurements, or deep quantum circuits, for instance to perform quantum phase estimation.
In this work, we propose a variational algorithm which uses imaginary time evolution to sequentially calculate the energy levels of a Hamiltonian. The algorithm makes use of the shallow swap test [19,20] to evaluate the overlap of two input wavefunctions [1]. We first use imaginary time evolution to target the ground state of the Hamiltonian. Using the shallow swap test, we can energetically penalise the ground state wavefunction, then discover the other eigenstates through repeated evolution and penalisation. This method makes use of only shallow circuits, at the cost of additional measurements.
A recent work by Higgott et al. [1] introduces the use of the swap test with the VQE to discover the energy eigenstates of the diatomic Hydrogen molecule. Here, we contrast the performance of methods based on direct descent with our imaginary time approach, finding that the former is prone to becoming stuck in non-physical local minima of the parameter space [14,21]. This may render them unsuitable for probing the full spectra of bigger systems. We numerically demonstrate this for a 6-qubit molecular Hamiltonian.
Conversely, we find that when evolution is restricted to a submanifold of the full Hilbert space, variational imaginary time evolution tends to converge to energy eigenstates of the Hamiltonian, regardless of the initial state. This is a crucial mechanism exploited by our algorithm to reliably penalise and discover the physical energy spectrum. We test our method on Hamiltonians which encode the boolean satisfiability problem (3SAT), and to find the electronic spectrum of the Lithium Hydride (LiH) molecule. Our algorithm makes use of variational imaginary time evolution, to be performed by a hybrid quantum-classical machine. We briefly outline the procedure below. See Ref. [14] for a more detailed discussion.
For a time independent Hamiltonian, H, the normalised imaginary time evolution is given by which is the solution of the imaginary time Schrödinger equation where H(τ ) = ψ(τ )| H |ψ(τ ) . Forgoing normalisation, a general state |ψ = j c j |e j evolves in imaginary time like where the probability of energy eigenstates |e j decay exponentially with their energies E j . Provided that |ψ(τ ) has a non-zero overlap with the ground state |g , it can be verified that lim τ →∞ |ψ(τ ) = |g . While non-unitary imaginary time evolution cannot be directly implemented on a quantum computer, it can be simulated using a hybrid quantum-classical algorithm. The state |ψ(τ ) is approximated by a parametrized trial state |ϕ(θ 1 (τ ), ..., θ M (τ )) := |ϕ( θ(τ )) , and its evolution determined by the evolution of θ(τ ). The trial state is produced by an ansatz quantum circuit |ϕ( θ(τ )) = U (θ M )U (θ M −1 )...U (θ 1 ) |0 , where U (θ k (τ )) is in practice a single or two qubit gate. The evolution of the parameters θ(τ ) under imaginary time evolution is given by where These elements are obtained by the quantum processor using the shallow quantum circuit shown in Appendix A. The classical processor can then update the parameters using the Euler update rule If the ansatz is sufficiently powerful, repeatedly constructing and solving this linear equation will evolve the system to a state close to the ground, which we denote as |g . We monitor convergence by the change in the parameters, and halt when ∆ θ(τ ) = δτ M −1 V ≈ 0. The expected energy of a converged state is easily evaluated using a polynomial number of simple Pauli operators [7]. With a less powerful ansatz, imaginary time evolution may fail to reach the ground state, but tends to converge to a higher excited eigenstate of the Hamiltonian. We do not presently provide a proof of this, though this behaviour is seen consistently in our numerical simulations.

III. EVALUATION OF THE ENERGY SPECTRUM OF THE HAMILTONIAN USING IMAGINARY TIME EVOLUTION
We now describe how to evaluate the excited states of the Hamiltonian. Having found an approximate ground state |g , we can construct a modified Hamiltonian which, for sufficiently large α ∈ R, no longer has ground state |g . Instead, the first excited state |e 1 of H becomes the ground state of H , and |g is now an excited state of H with energy increased by α, which will decay exponentially faster in imaginary time. The rest of the spectrum, orthogonal to |g , is unaffected. A system evolving under H in imaginary time will then approach |e 1 instead. This state can in turn be excited, and the system evolved under Hamiltonian to reach the next excited state of the original Hamiltonian. We can repeat this process by preparing the effective Hamiltonian H +α |g g|+ N j=1 α |ẽ j ẽ j | to obtain the (N + 1) th excited state |ẽ N +1 . In principle, we can obtain the complete energy spectrum, including a count of the degeneracies, so long as α is kept greater than the gap between ground and the highest state sought. Note the order of the discovered and subsequently excited eigenstates is unimportant.
In practice we do not directly modify the Hamiltonian, as doing so would require full tomography of the state vector, which is exponentially costly. Instead, we modify the imaginary time evolution equations to describe the evolution under the modified Hamiltonian, H . We replace V by and so on to excite all discovered eigenstates by α. These additional terms to V i can be evaluated using the low depth swap test circuit [19,20], outlined in Appendix B. We use the swap test to evaluate terms | ϕ(θ i + δθ i )|g | 2 and | ϕ|g | 2 , and then approximate for some sufficiently small δθ i . There is no requirement that each discovered eigenstate is excited by the same amount in the modified Hamiltonian. We may vary α for each excited state. In that scenario, one can add to V i in Eq. (5) to emulate a Hamiltonian with energy eigenvalues

IV. ERROR MITIGATION
We raise the possibility of applying error mitigation to our algorithm, through a simple error detection routine. Instead of using the low depth swap test described above, we can also use the conventional swap test (Fig. 1 [22]). The depth of this circuit grows linearly with the number of qubits used. The conventional swap test calculates the overlap between the two states by measuring an ancilla. However, no measurements are performed on the register qubits, and so any information we gain from them is, in a sense, free. We consider the input states as |g , |e . After the conventional swap test circuit, the register is left in the state where the sign is determined by the measurement result of the ancilla qubit. The state |φ ± r will be invariant under a symmetry S, if both |g and |e are also invariant under S. If we make a measurement of this symmetry on the register, we will be able to detect errors which break this symmetry. We can then discard those results for which we detect an error.
In the case of molecular Hamiltonians, we are often interested in ground and excited states which conserve the number of electrons in the molecule. If we use an ansatz which conserves the number of electrons (such as the unitary coupled cluster ansatz [23]) then a measurement of the electron number operator,N e , on the output state of the swap test should give the total number of electrons. If it does not, then an error has occurred, and we can discard the measurement. This method can thus mitigate the effect of single qubit bit flip errors, and certain combinations of two qubit errors. This error mitigation method can also be applied to the method developed in Ref. [1]. Moreover, our algorithm is compatible with the other error mitigation techniques proposed in Refs. [24,25]. We do not test these strategies in the present work.

V. NUMERICAL SIMULATIONS
We numerically simulate our algorithm with several Hamiltonians and several ansätze. Each time, our initial parameter values are random, and parameters are re-randomised when we excite states in the Hamiltonian. We employ Tikhonov regularisation when updating the parameters to ensure smoothness. We elaborate on these details and further describe our numerical methods in Appendix D.
The choice of ansatz is very important in variational simulation, and in this work, we explore the use of two. We try the recently proposed low depth circuit ansatz (LDCA) [26] which is chemically motivated and was  found to outperform the unitary coupled cluster ansatz for the molecule cyclobutadiene. We also employ the ansatz recently used in Ref. [14] to find the ground-state with imaginary time, which we refer to as the compact ansatz. Both ansätze scale linearly with the number of qubits, and can be considered hardware efficient [14,26].
We task our algorithm with finding the spectra of simple Hamiltonians which encode the 3SAT optimisation problem, and more complicated Hamiltonians which encode the electronic structure of LiH. We describe their structure below -see Appendix C for a detailed descrition of their construction. The 3SAT Hamiltonians are diagonal in the classical basis; where n j is the number of 3SAT clauses violated by the j th classical state when treated as a boolean assignment. This yields equally-spaced, highly degenerate spectra. We consider 3SATs Hamiltonians of up to 18 qubits. The LiH Hamiltonian can be simplified by employing various physical approximations -we do this to reduce the full 12 qubit Hamiltonian to 10 and 6 qubit representations.

Molecular Hamiltonians can be written as a linear combination of products of local Pauli operators,
where σ j i represents one of I , σ x , σ y or σ z , i denotes which qubit the operator acts on, and j denotes which term in the Hamiltonian we apply. For example We compare the spectrum reported by our simulated method with the eigenvalues of these Hamiltonian as found by exact numerical diagonalisation.
In Fig. 2, we present a simulation of our method exploring the simple spectrum of some 3SAT problems. The vertical axis is the expected energy ϕ( θ(τ ))|H|ϕ( θ(τ )) of the ansatz state, which we note is not necessary to monitor experimentally. The expected energy monotonically decays under imaginary time evolution until the system converges into an eigenstate, which is subsequently excited. It is interesting to note that the ground state is not necessarily discovered first, as demonstrated by the 16 qubit (right) simulation in Fig. 2, where ground is the third state discovered. For the more complicated reduced 6-qubit LiH Hamiltonian, we show that the variational imaginary time evolution successfully discovers eigenstates in Fig. 3a. In contrast, Fig. 3b reveals gradient descent converging to non-eigenstates which when subsequently excited, modify the Hamiltonian in a non-trivial way. This leads to errors in the discovered spectrum, shown in Fig. 4.
We next task our method with finding several of the lowest lying states of the physical 10-qubit LiH Hamiltonian as a function of the bond length. The results for both ansätze are shown in Fig. 5. The compact ansatz with 70 parameters shows reasonable agreement with the true sepctrum, despite generating only a small submanifold of the full 2 10 Hilbert space. The smooth deviation of the lowest discovered energy with the true ground energy may result from the ansatz's inability to generate the ground state. The low depth ansatz with 145 parameters shows a marked improvement in accuracy, and a better discovery of the degenerate energy eigenvalues.
Both ansätze show decreased accuracy around bond length l ≈ 2.5Å. This was also seen in recent variational eigensolver experiments [9], and attributed to the insufficient power of the low depth ansatz to generate these particularly highly entangled eigenstates [7].

VI. DISCUSSION
In this article, we have proposed a variational algorithm for a hybrid quantum-classical computer to discover the spectra of Hamiltonians. Our algorithm can also offer a route to enhancing the performace of the ground state solver in Ref. [14] since it can eliminate low lying states once found, thus 'clearing the way' for a successful identification of the true ground state. We tested our method on SAT and LiH Hamiltonians, using two different ansätze, and successfully obtained estimates of the excitation spectra. In our simulations we rarely saw variational imaginary time evolution converging to non-eigenstates. In contrast, gradient descent was prone to becoming stuck in local minima which when excited, caused errors in the reported spectrum.
Our results suggest a number of directions for fruitful future research. For instance, how should the necessary ability to accurately generate the energy eigenstates inform the design of the ansatz? And, how faithfully must the variational evolution simulate the true imagi- There are also questions concerning the classical component of the hybrid algorithm. For example, we might seek a fuller understanding of the relationship between the numerical solving algorithm employed by the classical processor and the consequential convergence of variational imaginary time evolution. We elaborate upon this in Appendix D.
A final topic to mention is the challenge of predicting the number of iterations necessary to converge to an eigenstate; Fig. 6 demonstrates an anomalous simulation where, despite the energy stabilising to an eigenvalue, the parameters continued to change.
Our method is not limited to exciting eigenstates -we can excite states for which the generating parameters are a priori known. This could be applied to eliminate unwanted subspaces from the searched spectrum, such as those which break symmetries or indicate error. Furthermore, our algorithm can be adapted to modify Hamiltonians in real-time variational simulation [13]. Discovered eigenstates can be excited by different amounts to modulate their new energies, for instance to create or remove energy degeneracies, or create time-dependence in the spectrum. Updating the linear equations in variational simulation by the procedure outlined in this work will then effectively simulate the real-time dynamics under the modified Hamiltonian. We leave exploring these extensions for a future work.  Ref. [20]. This circuit evaluates the overlap of two input density matrices ρ and σ.

VII. ACKNOWLEDGEMENTS
the qubits of ρ. We then measure observable L n=1 C n where Here, |0 0| n ρ projects the nth qubit of ρ onto the classical 0 state, and similarly for σ and the 1 state.
Measuring L n=1 C n is accomplished with post processing, by first measuring each of C n . If C n is measured as |1 1| n ρ ⊗ |1 1| n σ , we assign c n = −1, assigning c n = 1 for all other outcomes. The result of observable L n=1 C n is then Π L n=1 c n . By repeating this process and averaging over the results, we can evaluate the overlap function Tr(ρσ). rameters under the variational imaginary time evolution described in Eqs. (5) and (6).
In general, Eq. (6) can be ill-posed, and direct inversion of M is numerically unstable. We instead, after populating M and V at every time-step, update the parameters under Tikhonov regularisation, which minimises where the Tikhonov parameter λ can be varied to tradeoff accuracy against keeping˙ θ small and the parameter evolution smooth. Our simulations estimate an ideal λ at each time-step by selecting the corner of a 3-point Lcurve [32,33], though force λ ∈ [10 −4 , 10 −2 ]. This is because too large a λ over-restricts the change in the parameters in an iteration and was seen to lead to eventual convergence to non-eigenstates. Meanwhile, no regularisation (λ = 0) saw residuals in M −1 disrupt the monotonic decrease in energy. Still, using Tikhonov regularisation affords us a larger time-step than other tested methods, which included LU decomposition, least squares minimisation, singular value decomposition (SVD) and truncated SVD. Our simulations typically employ a time-step of δτ = 10 −1 . We suspect the largest stable time-step possible relates to the greatest energy eigenstate with non-negligible probability in the initial ansatz state.
We continue simulating in imaginary time until detecting convergence by a change in the parameters smaller than some threshold for several iterations, typically ∆ θ < 10 −2 for 3. The parametrised state is then assumed an eigenstate and has its state-vector recorded, to be subsequently excited in the Hamiltonian through modifying V via Eq. (9) every time-step thereafter. At this point, we reset the parameters to their initial values, restoring the original excited state (or one now of greater energy), and then resume imaginary time evolution.

Populating M and V
To save time, our code simulates only the ansatz and Hamiltonian circuits, using several shorcuts to avoid direct simulation of the entire circuits involved in populating M and V. Firstly, we calculate each ∂ |ϕ( θ(τ )) /∂θ i term by a fourth-order central finite-difference approximation with a step-size of ∆θ i = 10 −5 , in lieu of simulating the circuits shown in Appendix A. Full simulation of these sub-circuits is performed in Ref. [14]. M is then populated by the inner product of these terms, and V via their inner product with the statevector produced by simulating the Hamiltonian circuit on the ansatz. Excitations in V are introduced merely by the inner product of these terms and the recorded statevectors of the discovered eigenstates, in lieu of simulating the swap test circuits described in Appendix B. Our simulations typically excited the eigenstates by α ∼ 10, comparable to the gap between the ground and the highest considered excited state of the system.