Calculating transition amplitudes by variational quantum deflation

Variational quantum eigensolver (VQE) is an appealing candidate for the application of near-term quantum computers. A technique introduced in [Higgot et al., Quantum 3, 156 (2019)], which is named variational quantum deflation (VQD), has extended the ability of the VQE framework for finding excited states of a Hamiltonian. However, no method to evaluate transition amplitudes between the eigenstates found by the VQD without using any costly Hadamard-test-like circuit has been proposed despite its importance for computing properties of the system such as oscillator strengths of molecules. Here we propose a method to evaluate transition amplitudes between the eigenstates obtained by the VQD avoiding any Hadamard-test-like circuit. Our method relies only on the ability to estimate overlap between two states, so it does not restrict to the VQD eigenstates and applies for general situations. To support the significance of our method, we provide a comprehensive comparison of three previously proposed methods to find excited states with numerical simulation of three molecules (lithium hydride, diazene, and azobenzene) in a noiseless situation and find that the VQD method exhibits the best performance among the three methods. Finally, we demonstrate the validity of our method by calculating the oscillator strength of lithium hydride, comparing results from numerical simulations and real-hardware experiments on the cloud enabled quantum computer IBMQ Rome. Our results illustrate the superiority of the VQD to find excited states and widen its applicability to various quantum systems.


I. INTRODUCTION
We are now in an era where quantum computing practitioners can regularly use noisy quantum computers with tens of qubits [1]. Although the existing quantum computers have no immediate practical applications, we believe that there is a possibility that such a device outperforms existing classical algorithms in specific tasks. Various quantum algorithms for near-term devices have been suggested recently [2][3][4][5][6][7][8][9][10][11], and among such, the variational quantum eigensolver (VQE) [12] is considered to be an appealing candidate for applications of near-term quantum computers.
The original VQE [12] is a method for constructing an approximate ground state of a Hamiltonian on a programmable quantum device based on the variational principle of quantum mechanics. The VQE constructs the approximate ground state by iteratively tuning a quantum circuit to minimize an energy expectation value of the generated state. Because the near-term quantum devices are believed to be capable of generating a wavefunction that is not classically achievable, the VQE has the potential to explore a variational space that has not been investigated before. To expand the potential application of the VQE other than for the ground state, a lot of works have extended the method to evaluate properties * ibe@qunasys.com of excited states of a target Hamiltonian [13][14][15][16][17][18]. Those methods generally inherit the iterative and variational feature of the VQE, i.e., they also iteratively optimize a quantum circuit relative to some cost function.
The major and perhaps popular algorithms among such extensions are the subspace-search VQE (SSVQE) [13], the multistate contracted VQE (MCVQE) [14], and the variational quantum deflation (VQD) [15]. There are pseudo-eigenvalue extensions of VQE such as quantum subspace expansion (QSE) [19] or quantum equation of motion (qEOM) [20], but these algorithms are not fully variational methods and are therefore not included in our comparisons. While the SSVQE and the MCVQE can readily evaluate the transition amplitude | ψ 1 |A|ψ 2 | 2 of an observable A with respect to two approximate eigenstates |ψ 1 and |ψ 2 in a hardware-friendly manner (i.e., using less-costly quantum gates and circuits), the VQD, to the best of our knowledge, lacks such a method for evaluating the transition amplitude [21]. Since the transition amplitude is related to properties of the system such as the absorption and emission spectrum of photons [22,23], this severely limits the application range of the VQD method.
In this work, we fill this gap by providing a technique to evaluate the transition amplitude without using costly quantum circuits such as the Hadamard test [24]. Our technique is not only for the VQD, but can also be applied in a general setting where we have two orthogonal states |ψ 1 and |ψ 2 and a means to evaluate the overlap | ψ 1 |ψ 2 | 2 , and wish to evaluate | ψ 1 |A|ψ 2 | 2 . To support the significance of the proposed technique, we present a comprehensive comparison of the SSVQE, the MCVQE, and the VQD by conducting noiseless numerical simulations, where we use exact energy expectation values in the optimization routine of the parametrized circuit. In this test, we use molecular Hamiltonians of LiH and two azo compounds: diazene and azobenzene (AB). We find that, under this setting, the VQD generally exhibits better performance than the other two, which validates the importance of our proposed technique. Finally, as a demonstration of the technique, we conduct a proof-of-principle calculation both on a noisy simulator, i.e., expectation values of observables are simulated with the realistic hardware noise, including shot noise and environmental noise such as T1/T2 and readout errors, and on IBM's real quantum hardware.

II. EVALUATION OF TRANSITION AMPLITUDES
First, we propose a method to evaluate the transition amplitude of Hermitian operators between two quantum states. More concretely, we present how to evaluate | ψ 1 |A|ψ 2 | 2 for a Hermitian operator A and quantum states |ψ 1 and |ψ 2 such that ψ 1 |ψ 2 = 0 and where P i ∈ {I, X, Y, Z} ⊗n . I, X, Y, and Z are Pauli operators and a i ∈ R. We assume that, for any given two states |ψ 1 and |ψ 2 , we can evaluate the overlap | ψ 1 |ψ 2 | 2 . This evaluation can be performed by, e.g., the so-called swap test [25]. Note that the requirements on the hardware of quantum computers to evaluate the overlap can be relaxed when we know quantum circuits U 1 and U 2 that generate |ψ 1 and |ψ 2 , that is, we can evaluate the overlap by | ψ 1 |ψ 2 | 2 = | 0|U † 1 U 2 |0 | 2 . This is the case for calculating the transition amplitude for approximate eigenstates obtained by the VQD.
Let us consider unitary gates which can be realized as a product of Pauli rotation gates U ij,± = e ±i π 4 Pi e ±i π 4 Pj . We can show the following equality holds under the assumption of ψ 1 |ψ 2 = 0, which can be employed to evaluate | ψ 1 |A|ψ 2 | 2 . More concretely, we can measure each term in the right-hand side of Eq. (3) on a quantum device by regarding it as a overlap between two states because P i , P j , U ij,± are unitary and P i |ψ , P j |ψ , U ij,± |ψ can be realized on the device. We then combine the results of measurement according to the equation. Note that the assumption ψ 1 |ψ 2 = 0 should always be satisfied.
Equation (3) is one of the main results of this work. It reduces the evaluation of the transition amplitudes to a sequence of measurements of overlaps of two states. On the other hand, if we allow more complicated circuits to be executed on a device, we can construct an ancilla based technique to evaluate ψ 1 |A|ψ 2 as shown in Appendix B.

III. COMPARISON OF ALGORITHMS FOR EXCITED STATES
In this section, we compare the accuracy and capability of three algorithms for obtaining excited states of a given Hamiltonian on near-term quantum computers, namely the SSVQE, the MCVQE, and the VQD, by noiseless numerical simulations. By "noiseless", we mean all of the expectation values required in the algorithms are exactly computed. This situation is ideal compared with the real near-term devices but still appropriate for discerning the capability of the algorithms. This section aims to support the significance of our method proposed in the previous section by showing that the VQD gives the best performance among the three, and given a sufficiently low error rate and large enough number of sampled shots, one can reasonbly expect that these results will remain consistent when run on an actual quantum processor.
For the comparison, we use electronic Hamiltonians of LiH, diazene, and AB molecule. LiH is considered and employed as a simple "benchmark" molecule by a variety of studies on quantum computational chemistry [26][27][28][29]. Diazene and AB, on the other hand, are more relevant to applications of quantum chemistry to industry. In particular, AB is one of the most representative organic molecules which show cis-trans photoisomerization in photochemistry. AB has been attracting significant in-terests from the viewpoints of its photophysics and photochemistry associated with its various applications of photo-functional materials, and its derivatives are widely used as important photo-functional dyes in the industry [30]. For photo-functional molecules such as AB and its derivatives, it is crucial to theoretically predict their photophysical properties such as absorption and emission spectra or emission quantum yields and to elucidate their photochemical reaction mechanisms. Although the elucidation of its photoisomerization mechanism has been made theoretically so far, it remains controversial whether it proceeds with rotation or inversion or others. The simulations presented here do not only support the significance of the proposed method, but can also be viewed as a first step toward the real-world application of near-term quantum devices.

A. Settings of numerical simulation
Here we describe setups of numerical simulations for comparing the SSVQE, the MCVQE, and the VQD (the details of three algorithms are explained in Appendix A). As a variational quantum circuit for trial wavefunctions, we adopt an ansatz shown in Fig. 1 for all three algorithms. We call this ansatz real-valued symmetrypreserving (RSP) ansatz, U RSP (θ), where θ are classical parameters to be optimized. It is a slightly modified version of a heuristic ansatz introduced in Ref. [31], which preserves the number of particles of a reference state, so that the generated wavefunction is always a real-valued vector in the computational basis. As reference states for trial wavefunctions, we use the spin-restricted closedshell Hartree-Fock state and singly excited states. The BFGS method [32] is employed for classical optimization of the ansatz quantum circuit. The convergence criterion is set so that the optimization terminates when the relative difference of energy expectation value between iterations becomes lower than 10 −8 . The electronic Hamiltonian is computed by PySCF [33,34], an open-source quantum chemistry library, and mapped to the qubit Hamiltonian by the Jordan-Wigner transformation [35]. Note that the RSP ansatz only works for the Jordan-Wigner transformation. For all simulations in this section, the weight vector for the SSVQE (see Appendix A 2) is set as w = (k, k − 1, . . . , 1), where k is the number of quantum states to be calculated.
In the simulations, we compute several singlet and triplet eigenstates in the low-energy spectrum of the electronic Hamiltonians of LiH, diazene, and AB. To circumvent the need to find all three degenerate eigenstates in each triplet subspace, we slightly modify the original electronic Hamiltonian H to where S z is an operator representing the z-component of the total electron spin, and α > 0. When α is sufficiently large, all eigenstates of H which have non-zero S z are projected out from the low-energy subspace of H . In the following subsections, we use this H with α = 4 in atomic units as the target Hamiltonian in the optimization process presented above. This approach of adding "penalty terms" to a Hamiltonian can be found in Refs. [36][37][38]. All simulations in this section are performed using Qulacs [39].

B. Simple benchmark molecule: LiH
As for LiH molecule, we take STO-3G minimal basis set and all molecular orbitals into consideration, resulting in the number of simulated qubits being 12. The RSP ansatz mentioned above with D = 10 is used as the variational quantum circuit, where the total number of parameters is 110. We calculate three energy levels S 0 , T 1 , and S 1 [40] for 36 points of the interatomic length of LiH and compared with the full configuration interaction (full-CI) calculations. As the initial values of ansatz parameters θ, we use uniform random numbers drawn from [0, 2π] for the first point of the potential energy curve, and after that, we employ the optimized parameters of an adjacent point of the potential energy curve.
Calculated energies and their errors of LiH molecule by the SSVQE, the MCVQE, and the VQD are shown in Fig. 2 along with the result of the full-CI calculations. As seen in Fig. 2 (bottom), the VQD (the green lines) gives the most accurate energies than the other two, keeping the "chemical accuracy" throughout the plot.
C. More complex systems: diazene and azobenzene

Diazene
For diazene, we present VQE simulations using molecular structures along its minimum energy path (MEP) between trans and cis isomers. First, we obtain the structures by MEP calculations using the state-averaged complete active space self-consistent field (SA-CASSCF) method implemented in Molpro 2015 [41] with the 6-31G* basis set. The MEP calculations are done with the Gonzalez-Schlegel method [42]. We use an active space consisting of 4 orbitals (2×(lone pair on N) + π (HOMO) + π * (LUMO)) with 6 electrons. The S 0 -S 3 and T 1 -T 3 states are averaged in the SA-CASSCF calculations. In the simulations, we use several points on S 2 MEP, from S 2 Franck-Condon (FC) state of trans/cis isomers to the S 2 minimum, which proceeds with a rotation about the N-N bond associating the disruption of the double bond. Then, we perform the SSVQE, the MCVQE, and the VQD simulations of the Hamiltonians constructed in the same active space along the MEP and compare them with complete active space configuration interaction (CASCI) calculations, i.e., energies obtained by the exact diagonalization within the active space. The number of qubits to be simulated is eight. Here we utilize the RSP ansatz mentioned above with D = 20, where the total number of parameters is 140. Again, the initial values of ansatz parameters are taken as uniform random numbers drawn from [0, 2π] for the first point on the MEP. For other points on the MEP, optimized parameters at an adjacent point are used as initial values of parameters.
In addition to the energy of each eigenstate, we also calculate oscillator strength f ij between each pair of spinsinglet states |S i and |S j . It is defined as where E(S) is the energy of |S , R α = N l=1 r l,α is the electric dipole moment operator, and r l,α is the αcoordinate of the l-th electron. The oscillator strength gives the normalized strength of the absorption and emission spectrum of molecules [23], so it is fundamental for studying photochemical dynamics and reactions of molecules in quantum chemistry. Note that the oscillator strength involves the transition amplitude of R α operator, so it is impossible to evaluate it by the VQD on a quantum device in a hardware-friendly manner without using our proposed method in Sec. II. In this subsection, we evaluate it by exact values since we know all components of wavefunctions |S i , |S j and the matrix elements of R α by virtue of numerical simulations.
Results for diazene along the MEP are shown in Fig.  3. Calculated energies and deviations from the exact ones ( Fig. 3 (a-f)) show that the VQD method gives more accurate results than the other two in our simulation settings. Moreover, from Fig. 3 (g, h, i), the VQD gives the most accurate oscillator strengths, which indicates that the method enables precise calculations of the wavefunctions as well as the energy spectrum.

Azobenzene
For AB, we perform VQE simulations using two structures: trans/cis isomers. First, we obtain optimized trans/cis isomers using SA-CASSCF calculations, using an active space consisting of 3 orbitals ((lone pair on N) + π (HOMO) + π * (LUMO)) with 4 electrons. The S 0 -S 4 and T 1 -T 3 states are averaged in the SA-CASSCF calculations. The structures used in the VQE simulations are the minimum energy structures of S 0 state. Then, we perform the VQE simulations for the Hamiltonian in the same active space to obtain the eigenenergies and their deviations from the exact energies obtained by the exact diagonalization of the Hamiltonian. In this case, the number of qubits is six and we utilize RSP ansatz with D = 10, where the total number of parameters is 50. The oscillator strength is calculated as well by using the converged states. As the initial values of ansatz parameters, we use uniform random numbers within [0, 2π].
Results are shown in Fig. 4. Inheriting the trend from the previous results, the VQD gives more accurate results than the other two. Moreover, as evident in Fig. 4 (g, h,  i), the VQD gives the most accurate oscillator strengths, which indicates that the method enables us to generate precise wavefunctions.

D. Discussion
We have observed that the VQD has a better performance compared to the SSVQE and the MCVQE in this section. This result is probably because the requirement for the ansatz quantum circuit U (θ) is looser for the VQD than the SSVQE and the MCVQE; the ansatz with optimal parameters must make all reference states reside in the low-energy subspace simultaneously in the SSVQE and the MCVQE, while the VQD can do it separately for each reference state by the ansatz with different optimized parameters.
Moreover, the performance of the SSVQE seems worse than the MCVQE in our simulations. The part of the reason is because we employ the "weighted-sum" version of the SSVQE in the simulation [13]. If one used the "equal-weight" version of the SSQVE, which is equivalent to the MCVQE, the difference will vanish as long as the optimization of the ansatz quantum circuit goes well.

IV. IMPLEMENTATION OF TRANSITION AMPLITUDE EVALUATION FOR VQD
The previous section demonstrates the significance of the technique presented in Sec. II. That is, we provide the hardware-friendly way to evaluate the transition amplitude between two eigenstates with the most accurate method for generating excited states among the previously proposed three methods. In this section, to verify the correctness of the technique, we run it both on a more realistic simulator and the real hardware, which contain noise in expectation values of observables.
We employ the molecular Hamiltonian of LiH and take the active space of (2e, 2o) with STO-3G minimal basis set. We run the VQD simulation to calculate the potential energy curves of two singlet states, S 0 and S 1 , and calculate the oscillator strength between them by the method in Sec. II. We take 36 points of the bond distances from 0.5 to 4.0Å. The parity-mapping [43] is used to map a fermionic (electron) Hamiltonian to a qubit Hamiltonian and to reduce the number of qubits from 4 to 2 utilizing the symmetry of particle number and S z [44]. As the ansatz for the VQD, we use the "Two Local" with R y rotation blocks, CZ entanglement blocks, and a depth D = 4, which is available in Qiskit Circuit library in Qiskit Terra v0.17.0 [45] as qiskit.circuit.library.TwoLocal (see Fig. 5). For classical minimization of the cost function, the SLSQP method [32] is employed for statevector or QASM simulations and the SPSA [46] optimizer is used for realhardware experiments. Similar to the previous section, we set the convergence criterion so that the optimization terminates when the relative difference of energy expectation value between iterations becomes lower than 10 −8 . Energy derivatives concerning the circuit parameters are calculated with the so-called parameter-shift rule [47,48] to mitigate the shot noise explained in the next paragraph. To obtain only spin-singlet states, we add the expectation value of the total spin S 2 as a penalty term to the cost function of the VQD. The simulation in this section is performed using Qiskit [45].

A. Sampling and Hardware simulations
To elucidate where accuracy limitations arise in estimation of the transition amplitude, we evaluate simulating energy expectation values and overlaps in the optimization routine of the VQD in multiple approaches. First, we run the VQD optimization routine to obtain parameters using either an exact noiseless (statevector or SV) simulator or a backend noise model (QASM) simulator and then use these parameters in the statevector simulation (SV+SV or QASM+SV, respectively, in Fig. 6) to obtain physical quantities. Following this, we use the parameters obtained from simulation with a backend noise model and do the final evaluations on a quantum hard-ware ibmq rome, or with the evaluations completely run on the ibmq rome (QASM+HW or HW, respectively, in Fig. 7).
The statevector simulator is the same as explained in the previous section. For the noisy simulations, the expectation value of the energy is obtained by sampling 122,880 shots for each Pauli term in the Hamiltonian throughout the simulation. This sampling introduces the fluctuation in the energy expectation values, i.e., shot noise, as well as a realistic model of the quantum backend that includes T1,T2, gate and readout errors. Similarly, 122,880 shots measurements are sampled for evaluating | 0|U † 1 U 2 |0 | 2 overlaps with our method in Sec. II. For hardware evaluations, we use the shot count of 8192.
In the VQD parameter optimization using either QASM simulator or the real quantum hardware, we utilized the so-called state purification technique. Due to the environmental noise, the quantum state generated by the ansatz circuit can be a mixed state instead of a pure state. This state can be purified by diagonalizing the density matrix of the state obtained from quantum state tomography technique and selection of the pure state with the maximum eigenvalue resulting from this process. The energy of the pure state is then refined by computing the expectation value of the Hamiltonian with respect to the obtained pure state. The values of parameters in ansatz corresponding to the purified state is determined by fitting the vector of the ansatz to the purified maximum eigenstate. For a thorough discussion on the purification approach, we refer the reader to the literature of Ref. [49].
In Fig. 6, we show the result of limitations in the parameter optimizations due to finite sampling and backend noise. We run experiments with the VQD iterative loop with the exact statevector simulations, as well as with a noise model to obtain the ansatz parameters of Fig. 5, and then evaluate the energies and oscillator strength using those parameters in a noise free statevector simulation as shown in Fig. 6. When using the statevector simulator in the whole process (blue dots in Fig. 6), calculated energies achieve the exact solution. However, oscillator strengths deviate from the exact values when using the parameters obtained from the noisy simulations especially at bond lengths in which the energy spacing is small. This is, in part, because of the imperfection of the optimizations, as well as the imperfection of state purification, in which the ansatz generates a mixture of ground and excited state parameters.
In Fig. 7, we show the result of the parameters obtained from noisy simulations with the final evaluation on cloud enabled quantum backend ibmq rome, corresponding to the black dots for 36 different bond lengths. For two bond lengths, we carry out the simulation of energies and oscillator strength entirely on the quantum hardware, correspond to the green dots. We reduce the impact of state preparation and measurement (SPAM) errors by preparing each of possible 2 2 states and measuring the outcome to create an inversion matrix [50,51]  (top) Potential energy curve of S0 and S1 states and its deviation from exact CASCI calculations. (bottom) Oscillator strength between S0 and S1 states. For the VQD calculations, "QASM+SV (SV+SV)" indicates that the numerical simulation optimizes the circuit parameters with (without) the shot+environmental noise. In both cases, we calculate oscillator strengths without the shot noise after the optimization. For the VQD optimization routine in "QASM+SV", we utilize the state purification technique described in the text. Error bars are calculated in a way explained in Appendix C.
after energy convergence is achieved (see Appendix D). The initial parameters are selected at random when minimizing the ground state energy. These parameters are then used for the ground state parameters when minimizing the excited energy. From Fig. 7 we can see that the finite sampling on the real quantum hardware when determining the ansatz parameters (labeled HW in Fig. 7), leads to a more significant deviation from the exact answer when compared to the results that are obtained with a realistic backend noise model simulation, but with a significantly large number of shots (labeled QASM+HW in Fig. 7). This deviation of parameters on the real hard- (top) Potential energy curve of S0 and S1 states and its deviation from exact CASCI calculations. (bottom) Oscillator strength between S0 and S1 states. Here, QASM+HW indicates that we first run the VQD optimization on a numerical simulator with shot+environmental noise, and then calculate the energy and the oscillator strength on the real device. On the other hand, simply "HW" means that we run the whole process on the hardware. In both cases, we utilize the state purification technique described in the text for the VQD optimization routine. Error bars are calculated in a way explained in Appendix C.
ware can be relieved with the use of error mitigation techniques such as Richardson extrapolation [52,53].

V. CONCLUSION
In this work, we propose a general technique to evaluate transition amplitudes between two orthogonal states in a hardware-friendly manner on a quantum device. Its immediate application is the evaluation of transition amplitudes between the approximate eigenstates of the Hamiltonian obtained by the VQD method. The significance of the proposed method is supported by the comprehensive comparison of the three method, namely the SSVQE, the MCVQE and the VQD, in noiseless simulations which show the advantage of using the VQD for generating approximate excited states. Finally, we also verify the proposed method by running it in on a real near-term quantum devices. This work enlarges the possibility of the VQD and greatly advances the field of excited states calculations on a quantum device. In this section, we provide a review of the algorithms used in the main text.

VQE
The VQE [12] is a variational algorithm for finding the ground state of a system of n-qubits whose Hamiltonian is in the form of H = P ∈{I,X,Y,Z} ⊗n h P P. (A1) where I, X, Y, Z are single-qubit Pauli operators and h P ∈ R is a coefficient. If the number of terms with h P = 0 in the summation is not too large, we can evaluate the expectation value of the Hamiltonian, or energy of the system, by evaluating expectation values of each P and then summing them on a classical computer. To approximate the ground state of the Hamiltonian, the VQE uses parameterized quantum circuit U (θ), and iteratively optimize the parameter θ so that the energy expectation value E(θ) := 0|U † (θ)HU (θ)|0 , where |0 is an initialized state of the quantum computer, is minimized. The VQE algorithm proceeds as follows: 1. Define a quantum circuit U (θ) with parameters θ.
2. Repeat the followings until the convergence of E(θ).
When the convergence is reached, we expect that |ψ(θ) and E(θ) is an approximate ground state and its energy from the variational principle of the quantum mechanics.

Subspace-search VQE
The SSVQE [13] uses multiple initial states to search low-energy subspace of a Hamiltonian. The SSVQE algorithm can be summarized as follows: 1. Define an ansatz quantum circuit U (θ) and mutually orthogonal initial states (reference states) The reference states must be chosen so that one can readily make superpositions of them on a quantum device such as the computational basis.
2. Repeat the following steps until the convergence.
where the weight vector w is chosen such that w 1 > w 2 > · · · > w k > 0.
(c) Update parameter θ to decrease L.
The weight vector w has the effect of choosing which |ϕ i converges to which excited state. The cost function L w (θ) reaches its global minimum when the ansatz circuit U (θ) maps |ϕ i to the i-th excited state |E i of the Hamiltonian.
We note here that the assumed ability to create the superposition of {|ϕ i } i enables us to evaluate transition amplitudes of an operator between two eigenstates. It can be performed by creating two different superpositions of two eigenstates, measuring the operator of the interest, and postprocessing on a classical computer [13].

Multistate contracted VQE
The protocol of MCVQE [14] is similar to that of the SSVQE. It works as follows: 1. Perform Step 1 and 2 of the SSVQE, using a cost function where the weight vector is omitted: 2. Using the converged θ * , evaluateH ij := ψ i (θ)|H|ψ j (θ) for all i and j.
Energy spectrum ofH approximates that of the original H, and approximate eigenstates are obtained by superposing {|ψ i (θ) } i with coefficients determined by the eigenvectors ofH. The evaluation of transition amplitudes between the approximate eigenstates can be performed in the same manner as the SSVQE.

Variational quantum deflation
The VQD algorithm [15] is probably the most straightforward way to construct approximate eigenstates of a Hamiltonian H. The algorithm for finding the k-th excited state is as follows.
2. Set j = 1 and repeat the following until j = k.
(a) Define a Hamiltonian where {β i } is a set of sufficiently large realvalued coefficient.
(b) Perform the VQE to find an approximate ground state of H j .
The above algorithm works because H j has the j-th excited state of the original H. To evaluate the expectation value of H j , we need to evaluate the overlap between two states |ψ(θ) and |ψ(θ ) . It is suggested in [15] that we can either employ so-called destructive swap test [54] or measure them by 0|U † (θ)U (θ )|0 2 exploiting the knowledge of the circuit. We note that, while the previous two methods, namely the SSVQE and the MCVQE, can measure the transition amplitudes by creating the superposition of the initial states, there has been no efficient method for the VQD. Here, we describe a method to evaluate ψ 1 |A|ψ 2 using an ancilla qubit. We assume that we have descriptions of circuits U 1 and U 2 which generates |ψ 1 and |ψ 2 respectively from the initialized state |0 . Let, Λ(U 1 ) = |0 0| ⊗ U 1 + |1 1| ⊗ I, (B1) Λ(U 2 ) = |0 0| ⊗ I + |1 1| ⊗ U 2 .
(B2) Then, we have the following equality, We can recover ψ 1 |A|ψ 2 by combining them according to ψ 1 |A|ψ 2 = i a i ψ 1 |P i |ψ 2 . However, the above method uses expensive controlled-U gates, which might make it unfeasible on a near-term device.
Appendix C: Statistical errors in shot noise simulator The error bars of the energy of the S 0 and S 1 states in the top panel of Fig. 6 and 7 are drawn by the outputs of evaluate_with_result method of qiskit.aqua.operators.WeightedPauliOperator of Qiskit Aqua v0.7.3 [45]. The output is calculated as follows. For a Hamiltonian H = N i=1 c i P i with P i being Pauli operator, its statistical sampling error is where (∆P i ) 2 is the variance of measurement outcome of P i (= ±1).
The error bars for the oscillator strength (bottom panel of Fig. 6 and 7), on the other hand, are calculated by propagation of the error for energies of S 0 and S 1 states described in the above and error in transition amplitude | S 0 |R α |S 1 | 2 , according to the definition of the oscillator strength (Eq. (5)). The latter error is estimated by five realized values of the transition amplitude obtained by using Eq. (3), where each term in the right-hand side is computed as | 0|U |0 | 2 (return probability of |0 state after some circuit U ).

Appendix D: Measurement error mitigation
State preparation and measurement errors can significantly impact the accuracy of hardware based results and, importantly, these error will not necessarily remain stable over the course of an experiment. To mitigate the impact of these errors we use the readout error mitigation technique described in [50], in which 2 n state preparation circuits are run for n qubits. Specifically, we calibrate the A Full : y|A Full |x = m(y, x) N cal where x, y ∈ {0, 1} n are the input and measured bitstring respectively, m(x, y) is the number of shots for each input state x and N cal is the total number of shots used. During the experiments on the real hardware, this matrix is calibrated every thirty minutes with 8192 shots by preparing 2 2 basis states C = {{0, 0}, {0, 1}, {1, 0}, {1, 1}}. The matrix A Full is then inverted using the Qiskit function qiskit.ignis.mitigation.CompleteMeasFitter and applied to the measured experimental outcomes to mitigate the impact of readout errors.
Appendix E: Backend Properties