Practical quantum computation of chemical and nuclear energy levels using quantum imaginary time evolution and Lanczos algorithms

Various methods have been developed for the quantum computation of the ground and excited states of physical and chemical systems, but many of them require either large numbers of ancilla qubits or high-dimensional optimization in the presence of noise. The quantum imaginary-time evolution (QITE) and quantum Lanczos (QLanczos) methods proposed in Motta et al. (2020) eschew the aforementioned issues. In this study, we demonstrate the practical application of these algorithms to challenging quantum computations of relevance for chemistry and nuclear physics, using the deuteron-binding energy and molecular hydrogen binding and excited state energies as examples. With the correct choice of initial and final states, we show that the number of timesteps in QITE and QLanczos can be reduced significantly, which commensurately simplifies the required quantum circuit and improves compatibility with NISQ devices. We have performed these calculations on cloud-accessible IBM Q quantum computers. With the application of readout-error mitigation and Richardson error extrapolation, we have obtained ground and excited state energies that agree well with exact results obtained from diagonalization.


I. INTRODUCTION
Noisy intermediate scale quantum (NISQ) computers have recently become workhorse platforms for the study of codesign and the design of near-term quantum algorithms.Thus far, the variational quantum eigensolver has proved to be one of the most useful applications for these devices.Variational methods have been used to solve problems in chemistry, nuclear physics, quantum field theory, high energy physics, and others [2][3][4][5][6][7][8].While these small-scale applications show promise for using NISQ devices to sample from distributions and calculate expectation values, short coherence times make calculations involving time evolution exceedingly difficult on NISQ devices.Time evolution calculations hold promise for calculating scattering amplitudes [9] and, excited [10,11], and non-equilibrium states [12].One approach to the problem of short coherence times is quantum imaginary time evolution (QITE) [13], in which nonunitary evolution can be calculated variationally.Combining QITE with the Lanczos optimization method (re-ferred to as QLanczos in the context of quantum computing), one can obtain time-evolved phenomena of various many body systems [1].
Here, we demonstrate the practical application of QITE and QLanczos on current cloud-based NISQ hardware in order to calculate ground and excited states in different fields of study.We use the method to obtain the ground state of the deuteron nucleus in one instance, and we calculate both the ground and excited states of the H 2 molecule in another.The quantum computations were done on several cloud-accessible IBM Q Experience devices, i.e. 20-qubit Johannesburg, 20-qubit Poughkeepsie, 53-qubit Rochester, and 5-qubit Yorktown hardware.The results obtained from the quantum computations were compared with the classical calculations obtained from exact diagonalization.Despite the fact that we used a simplified version of the deuteron Hamiltonian we were able to obtain the ground state energy of deuteron without the need for any non-linear optimization or ancillae.We also obtained the energy spectrum of H 2 molecule very close and even within chemical accuracy (1.6 × 10 −3 Hartree).These demonstrations show great promise for scaling up time evolution as a solution method on near-term quantum hardware, and they illustrate that the approaches have practical, near-term applicability to an array of fields from high energy physics to chemistry.
Quantum imaginary time evolution addresses the problem of exponentially increasing resource requirements for computation as a function of the number of interacting particles.It replaces the real time in the time-dependent Schrödinger equation with imaginary time (t → −iβ).The solution to this equation involves an imaginary-time evolution operator, U = e −βH .This operator leads to the decay of all states except for the ground state.Therefore, the normalized imaginary-time evolution of a state can be expressed as where β is the imaginary time [14] and Quantum computation of the ground state energy of many-body systems using the imaginary-time evolution can be thought of as a natural alternative as quantum computers provide exponential speed ups.The basic idea behind QITE [13] is to approximate the non-unitary imaginary-time evolution in small steps with unitary updates on a set of qubits including data qubits and ancilla qubits.By tracing over the ancilla, the data qubits effectively see non-unitary evolution, which allows us to approximate imaginary time evolution and calculate the decay to the ground state via (1).However, the algorithm of Motta et al. [1] eliminates ancillae as a requirement, considerably simplifying the algorithm.On a quantum computer, the unitary evolution utilizes Trotterization.Current quantum computers are incapable of simulating long time evolution, or a large number of Trotter steps, due to short coherence times and excessive gate noise that further reduces coherence time.However, since QITE seeks to approximate non-unitary evolution with a unitary operator, we can reduce the number of Trotter steps by calculating a specific unitary that corresponds to the largest possible steps in imaginary time that yield a given desired accuracy.This amounts to solving a linear system of equations that provide coefficients of expansion, in terms of Pauli operators, for the unitary evolution operators.In the case of the deuteron, we found that solving this system of equations for the largest timesteps provided a unitary evolution operator that corresponded to the familiar unitary coupled cluster (UCC) ansatz [15].
While this was a large simplification of the QITE algorithm, a key advantage over variational methods is the ability to use the method in a QLanczos algorithm to calculate excited states.The basic idea behind the QLanczos algorithm is to fill the Krylov space with vectors in powers of e −2∆τ H , which is done using QITE, and then these vectors are used to calculate Hamiltonian matrix elements, which leads to a generalized eigenvalue equation, yielding a computation of ground and excited states.Using the single-step method in QITE, we reduced the depth of the quantum circuit, which makes these algorithms more compatible with NISQ [16] devices.
The rest of the paper is organized as follows.In section II we present how the quantum computations are conducted, including the algorithms (sections II A 1 and II A 2).Following this, in section II C 1 we present the results of our quantum computations conducted using several IBM Q quantum processors.Finally, in section III, we provide a summary and future work.In the appendix we provide information on the model Hamiltonians that we used for deuteron and molecular Hydrogen (section A), error mitigation strategies used (section B) and some details of the calculations (section C).

A. Algorithms
Here, we present a brief review of the QITE and QLanczos algorithms that were proposed in [1].

Quantum Imaginary Time Evolution (QITE)
To be able to simulate the dynamics of many-body systems we need to break down the Hamiltonian of these systems into local components such that H = M m h m where h m are non-commuting local terms of the system [17].For many-body systems, the number of terms in the Hamiltonian scales polynomially with the number of particles in the system.For example, the N = 2 deuteron Hamiltonian in (A4) can be decomposed into Because of the non-commuting terms in the Hamiltonian the decomposition of the evolution into small time steps and decomposing these steps into local gates can be done using the first order Lie-Trotter-Suzuki decomposition formula [18] which gives where n = β ∆τ is the number of steps in the evolution.For two non-commuting operators the matrix exponential can be written as following the Baker-Campbell-Hausdorff lemma.This formula is given for two operators only, but it can be generalized to n operators.In our calculations assuming that ∆τ is small we can approximate the imaginarytime evolution up to an order of O(∆τ ) as follows. where is the normalization constant.The s-th step of the imaginary-time evolution can be written as where s = 1, 2, . . ., n.The purpose of the QITE algorithm is to approximate (7) with unitary updates such that where A s can be written in terms of Pauli operators (defined in (A5)) up to D + 1 qubits and can be expressed as For our two (three)-qubit systems we used D = 1 (D = 2).To be able to approximate the imaginary-time evolution with these unitary updates we need to calculate the coefficients a[s].For small ∆τ , up to an order of O(∆τ ), the coefficients are found by solving a linear system of equations Sa[s] = b at every step of the imaginary-time evolution, where with The solution to this equation minimizes the operator norm ||c More detailed discussion on the calculation of the coefficients a[m] can be found in the supplementary information of ref. [1].
The calculation of the unitary updates for our deuteron and molecular Hydrogen examples gave us interesting results.For N = 2 case the unitary updates have the form of which are in the same form as UCC (unitary coupled cluster) ansätze that were proposed for molecular Hydrogen in [2] and deuteron in [5].This means that the unitary updates recover the UCC ansatz.
Using QITE it is possible to obtain the excited state energies since the system does not necessarily converge to the ground state, but rather depends on the initial state, |Ψ 0 , choice.In general, the system converges to the eigenvalue of the Hamiltonian whose eigenvector is non-orthogonal to the initial state, |Ψ 0 .

Quantum Lanczos (QLanczos) Algorithm
The QLanczos algorithm is based on the QITE algorithm, but provides the advantages of faster convergence in special cases, and it can be used to calculate excited state energies.The basic idea behind the QLanczos algorithm is to fill in the Krylov subspace with vectors in powers of e −2∆τ H at each Lanczos iteration such that K : {|Φ , e −2∆H |Φ , e −4∆H |Φ , . . .}.The vectors in the Krylov subspace are obtained using the QITE algorithm as for 0 ≤ l < L max assuming l is an even number.Here, |Ψ t = c t t s=1 e −i∆τ As |Ψ 0 = |Φ 0 is the initial QLanczos state which is obtained from QITE subroutine.After building the Krylov subspace we need to calculate the overlap matrix elements (T l,l ) and Hamiltonian matrix elements (H l,l ) in terms of the expectation values since they are the only experimentally accessible values.The calculations give overlap and Hamiltonian matrix elements as where r = l+l 2 .The normalization constants can be recursively calculated in terms of expectation values using The next step of the QLanczos algorithm is to utilize the calculated overlap and Hamiltonian matrix elements and solve the generalized eigenvalue equation The ground and excited states can then be found from the eigenvectors of the generalized eigenvalue equation.For example, the normalized ground (g) (excited (e)) state approximation is where the coefficients x l g(e) are obtained from the eigenvector that corresponds to the ground (excited) state energy such that x 0g(e) x 2g(e) . . .x Lmax g(e) T .Then the energy expectation values are calculated from which then leads to calculation of the ground and excited state energies using QLanczos algorithm.In the exact calculations the energy values obtained from the eigenvalues of the generalized eigenvalue equation ( 16) match with the values obtained from (18).Our quantum computation shows that using ( 18) is numerically more stable and gives much better results than using the eigenvalues of ( 16) as seen in Table I.The QLanczos method converges much faster than the QITE algorithm but one needs to do measurements at each imaginary-time projection of the Krylov subspace vectors to obtain the corresponding overlap and Hamiltonian matrix elements from the expectation values.The more vectors in the Krylov subspace the more QITE measurements with an increasing quantum circuit depth are required.At this point, the single-step method we proposed that is explained in section II B plays an important role in terms of reducing the circuit depth and possible noise that will arise due to the gates in the circuit.

B. Quantum Program
As mentioned in section II A 1, the imaginary-time evolution in QITE algorithm is provided by unitary updates of the form U s = e −i∆τ a[s](X0Y1−X1Y0) for our two-qubit examples.One way to obtain the ground state energy using QITE is to start with an initial product state, say |Ψ 0 = |10 and apply the unitary updates while calculating the coefficients a[s] that give the state |Ψ s at every step of the imaginary-time evolution.At the end of the n−th step of the imaginary-time evolution one expects to reach the ground state energy.This version of QITE would require a quantum circuit as seen in Fig. 1, which only shows the first two steps of the imaginary-time evolution; the depth of the quantum circuit increases as the number of steps increases.At every step of the imaginary-time evolution, the quantum circuit in the shaded area is repeated such that θ s = 2∆τ a[s].Naturally, large depth circuits are very noisy, and not necessarily amenable to error mitigation techniques.To make QITE more compatible with NISQ devices we reduce the number of time steps.In the single step version, instead of building the quantum circuit that combines each unitary update which gives |Ψ s ≈ e −i∆τ A[s] |Ψ s−1 we build the quantum circuit based on the calculated coefficient A s that gives |Ψ s ≈ e −i∆τ sA s |Ψ 0 .In this case, the quantum circuit is given in Fig. 2   In addition to our single-step QITE approach we also applied the error mitigation strategies to improve results (explained in Sec.B in detail).In what follows, we applied these error mitigation strategies to obtain the energy expectation values.

C. Results and Discussion
Here, we present the experimental results from IBM Q hardware for QITE and QLanczos algorithms.Information on the experiments and the hardware used can be found in Table IV of the supplementary information.

Deuteron
Using the QITE algorithm we were able to calculate the ground state energy of deuteron for both N = 2 and N = 3 cases.Fig. 3  Readout error mitigation suffice for β = 0 data points for both N = 2 and N = 3 case since they don't involve any CNOT gates.To obtain the energy measurements in Fig. 3(a) only readout error mitigation was conducted, except for β = 0.30 where both readout and extrapolation were used.Each experimental point in Fig. 3(b) is the result of post processing with readout error mitigation and Richardson extrapolation.

Molecular Hydrogen
Although we would expect QITE algorithm to converge to the ground state energy only, we found that depending on the choice of the initial state, |Ψ 0 , the excited state energy of the system can also be calculated.
In Fig. 5(a) we plotted the ground and first excited state energies as a function of bond length, R, that are obtained using QITE on hardware and compared with the values obtained from exact diagonalization.Because of the availability of devices we used two separate processors for calculation of the ground (on IBM Q 5 Yorktown) and first excited (on IBM Q Poughkeepsie) state energies.The ground (first excited) state energy values are calculated with an initial state of choice |Ψ 0 = |00 (|Ψ 0 = |10 ).In the case of chemical systems we would like to calculate energy values within chemical accuracy which, is 1.6 × 10 −3 Hartree.Therefore, in the inset of Fig. 5(a) we show the relative error in energy (∆E(R)) as a function of bond length compared with chemical accuracy.QITE was able to obtain chemical accuracy for one or two steps depending on the trial state.

QLanczos Results
As mentioned earlier, QLanczos can also be used for quantum computation of both the ground and excited state energies.The choice of the initial state, |Ψ 0 , is the one that determines which energies are being calculated.Here, we present our quantum computation of the ground (for deuteron and molecular Hydrogen) and excited state energies (for molecular Hydrogen only -note that the deuteron does not have a bound excited state) using QLanczos.
Quantum computation of the ground and excited state energies using QLanczos might require stabilization of the algorithm as the generalized eigenvalue equation in ( 16) might be numerically ill-conditioned.In our particular deuteron problem, due to the linear dependence of the vectors, |Φ l , in Krylov subspace, we had to perform the stabilization process explained in the supplementary
We ran QLanczos on two different devices: IBM Q 20-qubit Poughkeepsie (for N = 2 deuteron and first and second excited state energies of molecular Hydrogen) and IBM Q 53-qubit Rochester (for N = 3 deuteron and ground and third excited state energies molecular Hydrogen).The statistical error is calculated for N runs = 5 for deuteron and N runs = 3 for molecular Hydrogen, each run having 8192 shots.Results of our quantum computation of the ground state energies for N = 2 and N = 3 deuteron Hamiltonian are summarized in Table I.
In Table I and Fig. 5(b) we present the results for QLanczos with and without readout error mitigation (indicated as ROEM) for the deuteron and molecular Hy-drogen, repsectively.The results obtained using (18) are in good agreement with the values obtained from exact diagonalization, while energies obtained from the stabilized generalized eigenvalue equation do not agree well with the exact values due to stability issues in the case of molecular Hydrogen.Choosing a smaller regularization parameter would make these values closer to the exact values with a cost of adding more vectors to the Krylov subspace.In our example, a Krylov subspace with two vectors out of {|Φ 0 , |Φ 2 , |Φ 4 } subspace were sufficient to obtain the ground and excited state energies for the deuteron and molecular Hydrogen examples.For molecular Hydrogen, we used two different initial states (|Ψ 0 = |00 and |Ψ 0 = |10 ) which helped us to calculate the energy spectrum as a function of the bond length, R.
We found that using (18) gives very close values to exact diagonalization with or without readout error mitigation, meaning that QLanczos is potentially noise resilient.Combined with fast convergence the algorithm has a few advantages that make it useful for quantum computation of the ground and excited state energies of many-body systems.Since our QLanczos results are in good agreement with the exact values from diagonalization, we did not perform Richardson extrapolation.This would require 3 more measurements at every QITE step to build the Krylov space.
Although the computational limits of the quantum computers require us to truncate the harmonic oscillator basis, different schemes were proposed for extrapolating the bound state energies to infinite basis.We will follow the scheme that is based on the Lüscher's formula [19] N that was used in [5].The extrapolation of the bound state energy values to the infinite basis is listed in Table II.More information on the extrapolation of the ground state energy to the infinite harmonic oscillator basis can be found in the supplementary material in appendix C.

III. CONCLUSION
In this study, we presented a practical alternative for calculation of the ground and excited state energies of the many-body systems by using single-step version of the QITE and QLanczos algorithms presented in [1] using deuteron and molecular Hydrogen as specific examples.This approach may be a good low-depth circuit alternative to other contemporary methods.We also noted that QITE can be used to calculate the excited state energy whose eigenvector is non-orthogonal to the initial state |Ψ 0 .We also presented examples of the applications of readout error mitigation and Richardson extrapolation with these algorithms.On the other hand, QLanczos gave results that are good agreement with the exact diagonalization calculations therefore, it did not require additional error mitigation procedures.
We obtained the bound state energy of the deuteron at the next-to-leading order with a 0.5% (0.9%) error for N = 2 (N = 3) using QITE and with a 2.2% (1.6%) error for N = 2 (N = 3) case using QLanczos, compared to its experimental value of −2.22 MeV.We also showed the ground and excited state energies of the two-qubit molecular Hydrogen can be calculated within chemical accuracy using the QLanczos algorithm for a few bond lengths.

ACKNOWLEDGMENTS
We acknowledge useful discussions with C. W. Johnson, T. Morris, and E. Dumitrescu.The quantum circuits were drawn using Q-circuit package [20].This work was supported by the Quantum Information Science Enabled Discovery (QuantISED) for High Energy Physics program at ORNL under FWP number ERKAP61 and used resources of Oak Ridge Leadership Computing Facility located at ORNL, which is supported by the Office of Science of the Department of Energy under contract No. DE-AC05-00OR22725.The authors acknowledge use of the IBM Q for this work.The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Q team.
where σ ∈ {X, Y, Z} is in the j-th position with j = 0, . . ., N − 1, ⊗ indicates tensor product and I is the identity matrix.

Hydrogen Molecule
We will use the two-qubit molecular Hydrogen Hamiltonian [2] where coefficients h i (R) for i ∈ {0, 1, . . ., 5} are realvalued and functions of bond length, R, of the molecule.For calculation of the binding and excited state energies of the Hydrogen molecule we will use the coefficients calculated in STO-3G basis given in Table I of supplementary information of [7].

Appendix B: Error Mitigation
The noise due to the nature of the quantum simulators requires the application of the error mitigation strategies.Although there are various error mitigation strategies proposed in the literature, for our purposes, we used readout error mitigation and Richardson extrapolation techniques to reduce the noise involved in our calculations.
Out of the different sources of errors in a quantum circuit the readout errors are the errors associated with the final measurements in the quantum circuit.Therefore, we start by mitigating these errors in our quantum computation.To this end, we use the readout error mitigation scheme proposed in [25].In that scheme, the expectation values of the operators in the Hamiltonian are calculated using the following formula.
where p(x) is the probability of each qubit outcome and it takes 2 N values.For example, for N = 2, x ∈ {00, 01, 10, 11}.The symmetric and anti-symmetric combinations of the probability of i-th qubit flipping from 0 to 1 (p i (0|1)) or from 1 to 0 (p i (1|0)) is defined as Although p i (0|1) and p i (1|0) values are provided by IBM's Qiskit library, to get the most up-to-date values we obtained the readout error probabilities by preparing each qubit in computational basis 0 and 1 and then performing a measurement on each qubit in each case which gives us p(1|0) = # of states prepared in |1 measured in |0 # of shots (B3) or vice versa for p(0|1).We ran the simulations using 8192 number of shots.To propagate the error due to the statistical error in the readout errors for N = 2 deuteron case we did the readout error measurements 10 times and propagated the statistical error in measurements and statistical error in readout measurements in our results.As a result of our experimental measurements the statistical error in measurements is not different than the statistical error in readout error measurements therefore, we calculated the statistical error only for our N = 3 deuteron and molecular Hydrogen calculations.
Although we were able to reduce the depth of the quantum circuit using the single-step method, the decoherence effects became apparent in the expectation value measurements.Therefore, in addition to the readout error mitigation we also used the Richardson extrapolation ( [15], [26]) technique for the short-depth quantum circuits [27] to mitigate the errors associated with the noise produced by the gates used in the quantum circuit.The basic idea in this technique is to increase the error rate deliberately by a constant factor of r which is followed by an extrapolation to obtain the noise free expectation value.In this particular study, we increase the error rate by adding pairs of CNOT gates.The process of adding CNOT pairs is not expected to change the result of measurements since it corresponds to an identity matrix but it will contribute to the noise produced by CNOT gates.Our results showed that for two-qubit systems the expectation values of the observables scale linearly as and for N = 3 deuteron system they scale quadratically as where the coefficients A, B, and the extrapolated noiseless expectation value O(0) are found from the linear and quadratic fit to the data points of the expectation values of the operators for each case.We did not apply Richardson extrapolation technique to the QLanczos measurements since the results obtained using the QLanczos algorithm were in good agreement with the exact diagonalization results.
Appendix C: Extrapolation to the infinite harmonic oscillator basis The finite-size corrections to the infinite size harmonic oscillator basis based on the Lüscher's method can be stated as  where The values and definitions of the variables in (C1) are given in Table III.The terms in right-hand side of (C1) refer to leading order (LO), next-to-leading order (NLO) and N2LO, respectively.Curve fitting the LO and NLO terms gives the binding momentum, k ∞ and the asymptotic normalization coefficient, γ, for each order by using E 1 and E 2 .Fitting to N2LO term adding E 3 data helps calculating an effective range parameter, w 2 .

FIG. 1 :
FIG. 1: Two-qubit QITE quantum circuit with initial state |Ψ0 = |10 .The quantum circuit in the box is repeated at each QITE step after the second step of the algorithm for convergence.
(a) which only includes one CNOT gate for a specific initial state of |Ψ 0 = |10 .The rotation angle is now defined as θ s = 2s∆τ a [s] such that β = s∆τ is the imaginary-time corresponding to a specific expectation value, and at β = n∆τ the energy converges to the ground (or excited) state energy.We run the same quantum circuit with different calculated a [s] coefficients until the energy expectation value converges to the ground (or excited) state energy.Applying the same strategy to our three-qubit deuteron example with an initial state of |Ψ 0 = |100 gives the unitary updates of the form U s ≈ e −i∆τ sa 1 [s](X0Y1−X1Y0) e −i∆τ sa 2 [s](X0Z1Y2−X2Z1Y0) .(19) with θ s i = 2s∆τ a i [s] for i = 1, 2 which can be approximated with the quantum circuit in Fig. 2(b).

FIG. 4 :
FIG. 4: (Color online) Richardson extrapolation of the expectation values of the Pauli operators (on the right axis) and Hamiltonian operator (on the left axis) from their plots as a function of the CNOT gates corresponding to each CNOT gate in quantum circuit in Fig. 2(b) for N=3 qubit Hamiltonian at β = 0.30.This simulation was run on IBM Q 20-qubit Johannesburg hardware using the qubit layout [q0, q1, q2] = [8, 7, 9].

FIG. 5 :
FIG. 5: (Color online) Energy expectation values of two-qubit molecular Hydrogen as a function of bond length, R. We compared the values from exact diagonalization with the values obtained from hardware.The inset shows the relative errors of the quantum computed energy values compared to chemical accuracy.(a) The ground state energy (GSE) calculations (with |Ψ0 = |00 ) were done on IBM Q 5 Yorktown and the first excited state energy (1st ESE) calculations (with |Ψ0 = |10 ) were done on IBM Q Poughkeepsie hardware using the QITE algorithm.ROEM and Richardson extrapolation were applied.(b) The GSE and third excited state energy (3rd ESE) calculations (with |Ψ0 = |00 ) were done on IBM Q Rochester and the first and second excited state energy (2nd ESE) calculations (with |Ψ0 = |01 ) were done on IBM Q Poughkeepsie hardware using the QLanczos algorithm.The values with and without ROEM are presented.

TABLE III :
The values and definitions of the variables in (C1).

TABLE IV :
Information about the experimental runs on hardware.