Calculating nonadiabatic couplings and Berry's phase by variational quantum eigensolvers

The variational quantum eigensolver (VQE) is an algorithm to find eigenenergies and eigenstates of systems in quantum chemistry and quantum many-body physics. The VQE is one of the most promising applications of near-term quantum devices to investigate such systems. Here we propose an extension of the VQE to calculate the nonadiabatic couplings of molecules in quantum chemical systems and Berry's phase in quantum many-body systems. Both quantities play an important role to understand the properties of a system beyond the naive adiabatic picture, e.g., nonadiabatic dynamics and topological phase of matter. We provide quantum circuits and classical post-processings to calculate the nonadiabatic couplings and Berry's phase. Specifically, we show that the evaluation of the nonadiabatic couplings reduces to that of expectation values of observables while that of Berry's phase also requires one additional Hadamard test. Furthermore, we simulate the photodissociation dynamics of a lithium fluoride molecule using the nonadiabatic couplings evaluated on a real quantum device. Our proposal widens the applicability of the VQE and the possibility of near-term quantum devices to study molecules and quantum many-body systems.


I. INTRODUCTION
Quantum computers currently available or likely to be available in the near future are attracting growing attention. They are referred to as noisy intermediate-scale quantum (NISQ) devices [1], comprising tens or hundreds of qubits without quantum error correction. While it remains unclear whether they have "quantum advantage" over classical computers, the fact that they work explicitly based on the principle of quantum mechanics motivates researches on finding applications and developing quantum algorithms for practical problems that are classically intractable [2][3][4][5][6][7][8][9][10][11][12][13][14]. In particular, investigating quantum many-body systems with the variational quantum eigensolver (VQE) [15] is believed to be one of the most promising applications for NISQ devices.
The VQE is an algorithm to obtain eigenenergies and eigenstates of a given quantum Hamiltonian. In the VQE, quantum and classical computations are separated appropriately, and interactive quantum-classical hybrid architecture eases the difficulty of implementing the algorithm in the NISQ devices [15][16][17][18][19][20]. The VQE, which was originally proposed for finding the eigenenergy of the ground state, has been extended to find the excited energies and states [18,[21][22][23][24][25][26], non-equilibrium steady states [27,28], derivatives of eigenenergies with respect to external parameters of the system [29][30][31], and the Green's function [32]. * tamiya@qi.t.u-toyko.ac.jp † koh@qunasys.com ‡ nakagawa@qunasys.com This study aims to add a new recipe to the catalog of the VQE-based algorithms for quantum systems. We propose a method to calculate the nonadiabatic couplings (NACs) [33,34] of molecules in quantum chemistry and Berry's phase [35][36][37] of quantum many-body systems by utilizing the results of the VQE. Both quantities are related to the variation of slow degrees of freedom of the system and play a crucial role in the study of quantum chemistry, condensed matter physics, optics, and nuclear physics [36][37][38][39][40][41][42].
The NACs in quantum chemistry are defined as couplings between two different electronic states under the Born-Oppenheimer approximation [43], which are induced nonadiabatically by motions of nuclei (vibrations). They are fundamental in the nonadiabatic molecular dynamics simulations to study various interesting dynamical phenomena such as photochemical reactions around the conical intersection and electron transfers [38][39][40][41]. On the other hand, Berry's phase is defined as a phase acquired by an eigenstate when external parameters of a system are varied adiabatically along a closed path in the parameter space. It reflects intrinsic information about a system such as topological properties of materials. For example, several symmetry-protected topological phases are characterized by Berry's phase [44][45][46][47]. Berry's phase has become influential increasingly in many fields of modern physics, including condensed matter physics and high-energy physics [36,37,42].
Mathematically, the NACs and Berry's phase are related to derivatives of eigenstates with respect to external parameters of a system. In this study, in order to evaluate the NACs and Berry's phase based on the VQE, we develop analytical formulas and explicit quantum circuits to calculate the inner products related to the derivatives of the eigenstates. A naive way of calculating the NACs based on the VQE requires the Hadamard test [48] with a lot of controlled operations In contrast, our proposed methods for the NACs are based on the measurements of expectation values of observables, which is tractable on NISQ devices, and do not require the Hadamard test. As for Berry's phase, there is a previous study [49] to calculate it by simulating adiabatic dynamics and performing the Hadamard test at each time step. That method cannot avoid the undesired time-and energy-dependent dynamical phase contribution in addition to Berry's phase. Our proposed method for Berry's phase can remove dynamical phase contribution by utilizing the definition of Berry's phase although it still requires the Hadamard test at most once. Finally, as a demonstration of our methods, we present the simulation of photodissociation dynamics of a lithium fluoride molecule with the value of the nonadiabatic couplings evaluated on the real quantum device, IBM Q Experience [50], by our proposed methods. Our results enlarge the possible scopes of the VQE algorithm and the NISQ devices for simulating various quantum systems.
The rest of the paper is organized as follows. We briefly review the definition of the NACs and Berry's phase in Sec. II. The VQE algorithm is also reviewed in Sec. III. Our main results are presented in Secs. IV and V, where we describe the ways to calculate the NACs and Berry's phase based on the VQE. The results of the experiment of estimating the nonadiabatic coupling using IBM Q hardware and the simulation of photodissociation dynamics with our methods are shown in Sec. VI. The discussion about the cost analysis for running our algorithms on quantum devices is provided in Sec. VII. We conclude our study in Sec. VIII. Appendices provide details of the experiments, mathematical proofs of the cost analysis, and further numerical demonstrations of our algorithms.

II. REVIEW OF THE NONADIABATIC COUPLINGS AND BERRY'S PHASE
In this section, we review definitions of the NACs [33,34] and Berry's phase [35]. Let us consider a quantum system which has external parameters R = (R 1 , . . . , R Nx ) ∈ R Nx . These parameters R characterize the system, e.g., coordinates of nuclei in the case of quantum chemistry, electromagnetic field applied to a system in the case of conducting metals. We call R as "system-parameters" and represent the Hamiltonian of the system which depends on R by H( R). The eigenvalues and eigenstates of H( R) are denoted by {E i ( R)} i and {|χ i ( R) } i . We assume that {E i ( R)} i and {|χ i ( R) } i depend on R smoothly and that there is no degeneracy in the eigenspectrum unless explicitly stated in the text. When there is a degeneracy in the spectrum, the NACs is not well-defined among degenerate eigenstates. Berry's phase is generalized to non-abelian one, i.e., SU(N) ma-trix for N -degenerate ground states [42], and the components of the matrix can be determined in a similar way for abelian Berry's phase for the non-degenerate ground state studied in this paper.

A. Nonadiabatic couplings
Here let us consider a molecular system and H( R) as the electronic Hamiltonian. Definitions of the first-order NAC (1-NAC) d I kl and the second-order NAC (2-NAC) D I kl are as follows, where k and l are different indices for eigenlevels and I = 1, . . . , N x denotes the index for the system-parameters. The Hellman-Feynman theorem [51,52] gives a simpler expression of the 1-NAC as which means that the 1-NAC becomes large when two eigenstates are close to degenerate (E k ∼ E l ). We take advantage of this expression when calculating the 1-NAC in Sec. IV. The 1-NAC lies in the heart of various nonadiabatic molecular dynamics algorithms such as the Tully's fewest switches method [38,39] and ab initio multiple spawning [53,54]. Equation (2) in the case of k = l is related to the diagonal Born-Oppenheimer correction (DBOC) defined as where k is the eigenlevel to be considered, M m is the mass of the nucleus m, and R mα is α-cordinate (α = x, y, z) of the nucleus m. It is argued that this correction sometimes brings out crucial differences in stability and dynamics of molecules [55][56][57][58].
In addition, we comment on the gauge invariance of the NACs. Overall phase factors of eigenstates are arbitrary in general, so there is a U (1) M degree of freedom in the definition of the NACs, where k = 0, . . . , M − 1, M is the number of eigenlevels to be considered, and Θ k ( R) ∈ R is an arbitrary smooth function of R. The 1-NAC (Eq. (1)) and the 2-NAC (Eq. (2)) are not invariant under the transformation (5). This dependence must be resolved in each algorithm utilizing the value of the NACs. For example, see Refs. [59][60][61]. We note that real-valued eigenfunctions are usually considered in quantum chemistry, but complex eigenfunctions may be obtained in the VQE in general.

B. Berry's phase
Berry's phase [35] is defined for a closed loop C in the parameter space R Nx as, where C . . . is the line integral along the closed loop C, |χ k ( R) is the k-th eigenstate of the Hamiltonian H( R).
If one prepares the k-th eigenstate of the system |χ k ( R 0 ) at some system-parameters R 0 and adiabatically varies them in time along C, the final state will obtain the phase e −iΠ C in addition to the dynamical phase. We note that Berry's phase is always real by definition because the normalization condition Finally, we point out the gauge invariance of Berry's phase. The eigenstates have U (1) gauge freedom stemming from arbitrariness of overall phases for them. Under U (1) gauge transformation (Eq. (5)), Berry's phase is invariant only up to an integer multiple of 2π. Since Berry's phase appears as e −iΠ C , this arbitrariness does not affect the physics, and we can consider Berry's phase as an observable property of the system [36,37,42].

III. REVIEW OF VARIATIONAL QUANTUM EIGENSOLVER
In this section, we review the VQE algorithm [15] to obtain a ground state and excited states of a given Hamiltonian. We also describe how to compute analytical derivatives of optimal circuit-parameters of the VQE with respect to system-parameters of the Hamiltonian. Methods described in this section are repeatedly used in Secs. IV and V to calculate the 1-and 2-NACs and Berry's phase.
Again, let us consider an n-qubit quantum system whose Hamiltonian is H( R). In the VQE, we introduce an ansatz quantum circuit U ( θ) and the ansatz state |ϕ 0 ( θ) in the form of where |ϕ 0 is a reference state and θ = (θ 1 , . . . , θ N θ ) ∈ R N θ is a vector of circuit-parameters contained in the ansatz circuit. We assume U ( θ) to be a product of unitary matrices each with one parameter, We also assume each unitary U a (θ a ) consists of nonparametric quantum gates and parametric gates in the form of e igaPaθa generated by a Pauli product P a ∈ {I, X, Y, Z} ⊗n with a coefficient g a ∈ R (a = 1, . . . , N θ ). Note that many ansätze proposed in previous studies fall into this category [15,17,23,[62][63][64][65][66]. We will represent U j (θ j ) · · · U i (θ i ) as U i:j for simplicity.
A. Variational quantum eigensolver for ground state and excited states The original VQE algorithm finds a ground state of a given Hamiltonian based on the variational principle of quantum mechanics. In the VQE, one optimizes the circuit-parameters θ variationally by classical computers so that the expectation value is minimized with respect to θ. When the ansatz circuit has sufficient capability of expressing the ground state of H( R) and the circuit-parameters θ converge to optimal ones θ * , we can expect the optimal state |ϕ 0 ( θ * ) will be a good approximation to the ground state. Since tasks of evaluation and optimization of quantum circuits are distributed to quantum and classical computers, it is easier to implement the algorithm on the near-quantum devices [15][16][17][18][19][20].
After the proposal of the original VQE algorithm, there are a variety of extensions of the VQE to find excited states of a given Hamiltonian [18,[21][22][23][24][25][26]. As we will see in Secs. IV and V, one has to compute (approximate) eigenenergies and transition amplitudes of several Pauli operators between obtained eigenstates to calculate the NACs. From this viewpoint, the most appropriate methods to calculate them are the subspace-search VQE (SSVQE) [22] algorithm and its cousin algorithm, the multistate contracted VQE (MCVQE) algorithm [23]. Here we briefly describe the SSVQE just for completeness, but formulas for the MCVQE are quite similar.
To obtain approximate eigenenergies and eigenstates up to i = 0, . . . , M − 1, the SSVQE algorithm uses M easy-to-prepare orthonormal states {|ψ i } M −1 i=0 (e.g. computational basis) as reference states. For our algorithms to work, the reference states also have to be chosen so that we can readily prepare the superpositions of them on quantum devices. The SSVQE proceeds so as to minimize the following cost function, where {w i } M −1 i=0 are positive and real weights which satisfy w 0 > w 1 > · · · > w M −1 > 0. When the cost function converges to the minimum at θ * ( R), it follows that are good approximations of the eigenstates and eigenenergies, respectively. One of the most distinctive features of the SSVQE and the MCVQE algorithms is that one can readily compute transition amplitudes ϕ k ( R)|A|ϕ l ( R) of any observable A between the (approximate) eigenstates obtained. Although evaluation of the transition amplitude between two quantum states requires the Hadamard test in general, which contains a lot of extra and costly controlled gates [67], the SSVQE and the MCVQE circumvent the difficulty by preparing superposition of two eigenstates. It is possible to evaluate the transition amplitude by lowcost quantum circuits without extra controlled gates as where |ϕ ± k,l ( R) = U ( θ * ( R))(|ψ k ± |ψ l )/ √ 2 and Since each term of the right hand sides of the equation is an expectation value of the observable, the evaluation of the transition amplitude is tractable on near-term quantum devices.

B. Derivatives of optimal parameters
To calculate the NACs with the result of the VQE on near-term quantum devices, we also need derivatives of the optimal circuit-parameters θ * ( R) with respect to the system parameters R. These derivatives are given by solving equations [29] where simultaneously for a = 1, . . . , N θ (with I, J = 1, . . . , N x fixed). Now we use notations as follows: These formulas (Eqs. (14) and (15)) can be derived by taking the derivative of ∂E0( θ * ( R), R) ∂θa = 0 with respect to R. For detailed derivation, see Appendix. A in Ref. [29].
The quantities appearing in Eq. (14) and Eq. (15), , can be evaluated quantum circuits on quantum devices using the method shown in Ref. [29]. Therefore one can solve Eq. (14) and Eq. (15) on classical computers and obtain the derivatives of the optimal circuit-parameters { ∂θ * a ( R)

IV. CALCULATING NONADIABATIC COUPLINGS WITH VARIATIONAL QUANTUM EIGENSOLVER
In this section, we explain how to calculate the 1-NAC and 2-NAC with the VQE.

A. First-order nonadiabatic coupling
Evaluation of the 1-NAC based on the VQE is simple by utilizing the formula (3). First, we perform the SSVQE or the MCVQE and obtain approximate eigenstates |ϕ i ( R) and eigenenergiesẼ i of H( R). Then we calculate the derivative of the Hamiltonian, ∂H ∂R I , on classical computers. Specifically, when we use the Hartree-Fock orbitals to construct the second-quantized Hamiltonian, the derivative ∂H ∂R I (more precisely, the derivatives of the one-and two-electron integrals in the molecular orbital basis) can be obtained by solving the coupled perturbed Hartree-Fock (CPHF) equation [29][30][31]. The solution of the CPHF equation can be obtained by the standard softwares for quantum chemistry.
Finally, evaluating the transition amplitude ϕ k ( R)| ∂H ∂R I |ϕ l ( R) on quantum devices by using the method of Eq. (13) and substituting it into Eq. (3) gives the value of the 1-NAC.

B. Second-order nonadiabatic coupling
Next, we introduce an analytical evaluation method of the 2-NAC on near-term quantum devices. After obtaining approximate eigenstates {|ϕ i ( R) } i by the SSVQE or the MCVQE, putting them into Eq. (2) yields where we denote ∂ ∂θa ∂ ∂θ b |ϕ j and ∂ ∂θc |ϕ j as |∂ a ∂ b ϕ j and |∂ c ϕ j , respectively. We note that plugging Eq. (18) when k = l into Eq. (4) gives the formula of the DBOC based on the VQE. The derivatives of the optimal circuit-parameters such as ∂θ * a ∂R I and can be calculated by the method reviewed in Sec. III. The terms ϕ k |∂ a ∂ b ϕ l and ϕ k |∂ c ϕ l can be evaluated with the Hadamard test [48] in a naive way, but its implementation is costly for near-term quantum devices. Therefore, in the following, we describe how to reduce the evaluation of ϕ k |∂ a ∂ b ϕ l and ϕ k |∂ c ϕ l to the measurements of the expectation value of observables, which is the standard process of the near-term quantum algorithms.
with |Φ being an arbitrary reference state.
When a = b, we assume 1 ≤ b < a ≤ N θ without loss of generality. By using the method in Ref. [67], the real and imaginary parts of Eq. (19) are evaluated separately in the following way.
The real part is calculated with quantum circuits containing projective measurements of the Pauli operator P b denoted by M P b , where is the conditional expectation value of P a when the projective measurement of P b yields ±1, and is the probability of getting the result ±1 for the projective measurement of P b . If P b is a single Pauli operator or even if P b is a multi-qubit Pauli operator, we expect that the projective measurement of it can be performed in near-term quantum devices [68]. The total circuit for evaluating Eq. (21) is shown in Fig. 1(a).
On the other hand, the imaginary part of Eq. (19) can be calculated as where is the expectation value of P a for the quantum state U b+1:a e ±iπP b /4 U 1:b |Φ . The circuit for calculation is shown in Fig 1(b). Then, to obtain ϕ k |∂ a ∂ b ϕ l we take advantage of the following equality where |ϕ ± k,l = U 1:N (|ψ k ± |ψ l )/ √ 2 and |ϕ i± k,l = U 1:N (|ψ k ± i |ψ l )/ √ 2. All terms in the right hand side of (25) can be evaluated by the method described above with taking |Φ appropriately, so the ϕ k |∂ a ∂ b ϕ l is also obtained.
2. Evaluation of ϕ k |∂cϕ l Next, we describe how to compute ϕ k |∂ c ϕ l . It follows that The term in the last line can be evaluated by the method of Eq. (13) by substituting |ϕ ± k,l ( R) by U 1:c 1 √ 2 (|ψ k ± |ψ l ) and |ϕ ±i k,l ( R) with U 1:c 1 √

Summary
In summary, calculation of the 2-NAC D I kl proceeds as follows: 1. Perform the SSVQE or the MCVQE and obtain approximate eigenstates |ϕ i ( R) and eigenenergies The main contribution of this paper is that we reduce the definition of the NACs (Eqs. (1) and (2)) to the formulas that we can evaluate on quantum devices by the existing techniques. Here we note that the procedure 2 follows the techniques in Ref. [29], the procedure 3 partially uses those in Ref. [67], and the procedure 4 basically follows those in Ref. [22].

V. CALCULATING BERRY'S PHASE WITH VARIATIONAL QUANTUM EIGENSOLVER
In this section, we describe a method for calculating Berry's phase with the VQE algorithm. From the results of the VQE, while we can access the density operators of the eigenstate ρ k ( θ * ) = |ϕ 0 ( θ * ) ϕ 0 ( θ * )| determined by the optimized circuit-parameters θ * , we cannot access the information about the phase of quantum state. Here we discuss how to calculate Berry's phase on quantum devices from the optimized circuit-parameters obtained by the VQE. In the following, without loss of generality, we only consider the ground state as the eigenstate. Let N 0 denote the set of normalized states in a complex Hilbert space H. We consider performing the VQE from one point R 0 := R(t 0 ) of the closed loop , R(t 0 ) = R(t 1 )} in the systemparameters space and continue doing it along C R , then we obtain a smooth curve C θ * := { θ * ( R(t)) | t ∈ [t 0 , t 1 ]} in the circuit-parameter space. For simplicity, let θ * s := θ * ( R(t 0 )) and θ * t := θ * ( R(t 1 )) denote the starting point and the end point of C θ * , respectively. We note that θ * s = θ * t may occur, i.e., the curve C θ * of the optimal parameters does not necessarily form the closed loop in the circuit-parameter space even when C R is the closed loop in the system-parameter space. This is because the VQE does not care about the overall phase of the ground state, and for most cases there is a redundancy in the ansatz |ϕ 0 ( θ) such that |ϕ 0 ( θ 1 ) = e iξ |ϕ 0 ( θ 2 ) , e iξ = 1 for some θ 1 = θ 2 . Next, we introduce the projective Hilbert space called Ray space. Ray space R is defined as the equivalent class R := N 0 / ∼ where the equivalence relation ∼ holds for two elements of N 0 which differ only by a global phase. We also define the projection map π : |ψ ∈ N 0 → ρ = |ψ ψ| ∈ R. For a given curve and this curve C ρ in R is determined uniquely according to the optimized circuit-parameters θ * .
Then let us describe and formulate the way to calculate Berry's phase based on the results of the VQE. Suppose that a curve C ρ = {ρ( θ * )} is given. Here, we consider a particular lift We assume that |ϕ 0 ( θ * ) is smooth, i.e., |ϕ 0 ( θ * ) is differentiable with respect to θ * . With this lift C N0 , Berry's phase can be defined as [69] We want to emphasize here that Berry's phase is a functional of the curve C ρ . Namely, for a given curve C ρ , though we can construct a new curve C N0 which differs only by U (1) phase degree of freedom from C N0 with a real smooth function Θ( θ * ), the value of Berry's phase Π Cρ calculated with the curve C N0 is identical to that with C N0 . As discussed above, by performing the VQE, we obtain the curves C θ * and C ρ . To calculate Berry's phase based on the Eq. (27), we have to choose some fixed lift C N0 from C ρ . Due to the arbitrariness of the lift, we can fix the phase freedom globally and choose C N0 so that the freedom does not depend on θ * . Therefore, given a ρ( θ * s ), we first choose |ϕ 0 ( θ * s ) up to phase, and then form a lift C N0 = {|ϕ 0 ( θ * ) } uniquely up to a phase degree of freedom in the starting point of the curve |ϕ 0 ( θ * s ) . By considering such lift, the terms appearing in the Eq. (27) can be reduced to the quantities which can be evaluated with quantum devices. In the following, we explain how to evaluate the terms in the right hand side of Eq. (27).

A. Evaluation of the first term
The first term of Eq. (27) is computed by discretization of the closed loop C R and numerical integration of the integrand. We discretize the value of the systemparameters R on C R as R 0 , . . . , R K−1 appropriately and also define R K = R 0 . The VQE algorithm is performed for all points { R p } K p=0 and the optimal circuitparameters are obtained as { θ * p = θ * ( R p )} K p=0 . We define θ * 0 := θ * s ( R 0 ) and θ * K := θ * t ( R 0 ) and stress again that θ * 0 = θ * K may hold in general due to the redundancy of the ansatz. Here because we choose the phase freedom of the lift C N0 = {|ϕ 0 ( θ * ) } which is independent on θ * , the integrand of the first term can be written as so it is evaluated by measuring the expectation value of P a for the state U * 1:a |ψ 0 , which can be evaluated on quantum devices. Therefore the integral is approximated by

C. Comparison with previous studies
We here compare previous work on calculating Berry's phase on quantum devices with our method. In Ref. [49], Berry's phase is calculated by simulating adiabatic dynamics of the system U C = T e −i T 0 dsH(s) , where T is the time-ordered product and H(s) is a time-dependent Hamiltonian, which varies along the closed loop C in sufficiently long time T . U C is implemented on quantum computers by the Suzuki-Trotter decomposition, and the Hadamard test like Fig. 2 is performed to detect the phase difference between the initial ground state |χ( R 0 ) and the time-evolved state U C |χ( R 0 ) . The phase difference between |χ( R 0 ) and U C |χ( R 0 ) contains the dynamical phase and Berry's phase, but the former phase can be neglected by combining the forward-and backward-time evolutions by assuming the Hamiltonian of the system has the time-reversal symmetry. Compared with this strategy, our proposal for calculating Berry's phase based on the VQE has two features. First, we do not have to assume the time-reversal symmetry in the system to remove the contribution from the dynamical phase like in the previous study because we directly calculate Berry's phase based on the definition Eq. (6). Our method can apply to general quantum systems. Second, the causes of errors are quite different. More concretely, while the errors in the previous methods arise from the Trotterization of the time evolution operator, the errors in our method mainly come from two sources: one is the approximation error of the eigenstates obtained by the VQE, and the other is the numerical error of integration in Eq. (30). These errors can be reduced by deepening the ansatz circuits and taking more discretized points on C, respectively. We comment that further research is needed to conclude the difference in the performance between our method and these previous methods.
Finally, we introduce another method to calculate Berry's phase based on the VQE with a lot of Hadamard tests. Using the formula with taking the principal branch of the complex logarithm, −π ≤ Im(z) < π for z ∈ C, is one of the candidates for avoiding discretization error of the closed loop C and numerical instability [71]. The value of ϕ 0 ( θ * i )|ϕ 0 ( θ * i+1 ) is evaluated by the Hadamard test in Fig. 2 by substituting θ * 0(K) with θ * i(i+1) .

D. Summary
Berry's phase can be calculated based on the VQE as follows: 1. Discretize the closed loop C in the systemparameters space as { R p } K p=0 appropriately and perform the VQE for all points. 3. If necessary, evaluate the phase difference arg( ϕ 0 ( θ * s )|ϕ 0 ( θ * t ) ) by the Hadamard test shown in Fig. 2. 4. Substituting all values obtained in the previous steps into Eq. (27) gives the Berry's phase.

VI. EXPERIMENT ON A REAL QUANTUM DEVICE
In this section, we show an experimental result of our algorithm for the 1-NAC on a real quantum device and the non-adiabatic molecular dynamics (MD) simulation based on the experimental values of the 1-NAC.
We consider a lithium fluoride (LiF) molecule with bond length R and its electronic states under the Born-Oppenheimer approximation. In this system, the potential energy curves, or eigenenergies E as a function of R, of the two lowest 1 Σ + states are known to exhibit the avoided crossing [74][75][76][77], and it plays a crucial role for nonadiabatic dynamics such as photodissociation. Here we focus on this avoided crossing and the resulting nonadiabaitc dynamics by modeling the system with a simple two-state model. Specifically, the electronic Hamiltonian of LiF at bond distance R is constructed by two orbitals obtained by the state-average complete active space selfconsistent field (CASSCF) method [78]. By considering symmetries in the system, one can obtain a two-qubit Hamiltonian H LiF (R) from that electronic Hamiltonian. Further details are described in Appendix A.
We run our algorithm to calculate the 1-NAC (Sec. IV) of H LiF (R) at various bond length R in the IBM Q cloud quantum device (ibmq_valencia) [50]. First, the SSVQE calculation for H LiF (R) is performed on a classical simulator where the exact and noiseless expectation values of observables are obtained. The ansatz for the SSVQE, U ( θ), is depicted in Fig. 4 and we optimize the parameters θ to minimize the cost function (10). The initial states to which the ansatz applied are |ψ 0 = |00 , |ψ 1 = |01 states. To calculate the singlet (S = 0) states 1 Σ + , we penalize the triplet (S = 1) states by modifying the Hamiltonian in the cost function as H = H LiF (R) + βŜ 2 , whereŜ 2 is the spin squared operator and β = 4.0 is a constant [79]. We obtain the (approximate) eigenstates |ϕ 0,1 ( θ * ) = U ( θ * ) |ψ 0,1 and eigenenergiesẼ 0,1 for two 1 Σ + states in this way from the classical simulator ( Fig. 3(a) shows the potential energy curves). Next, the transition amplitude ϕ 0 ( θ * )| dHLiF(R) dR |ϕ 1 ( θ * ) is computed on a quantum device by evaluating each term of the right hand sides of Eq. (13) with 8192 shots. Since the wavefunction is real by definition of the ansatz U ( θ), we consider only the real part of Eq. (13). Finally, the 1-NAC is obtained by plugging the estimate of ϕ 0 ( θ * )| dHLiF(R) dR |ϕ 1 ( θ * ) into the numerator of Eq. (3) and that ofẼ 0 −Ẽ 1 into the denominator. As we mentioned in Sec. II, the 1-NAC is not gauge invariant. Even when considering only real wavefunctions as in this case, the 1-NAC still has indefinite sign, but here the sign is determined by considering the continuity of the 1-NAC with respect to the nuclear coordinate R. We comment that the way of calculating the 1-NAC here is chosen to remove the effect of the noise and errors in the real quantum device during the SSVQE and focus on evaluating the transition amplitude.
The result of the 1-NAC is shown in Fig. 3(a). Because of the noise in the real quantum device, the values of the transition amplitude are smaller than the exact results but still qualitatively consistent with them. This shrinking could be resolved by, for example, using the error mitigation technique [80,81] for near-term quantum devices.
In addition, we perform the trajectory surface hopping (TSH) [82] molecular dynamics calculation using Tully's fewest switches algorithm [38] based on the obtained values of the potential energy curves and the 1-NAC ( Fig. 3(a)). We assume a situation where a LiF molecule is excited by light to the first electronic excited state. In the TSH simulation, the nuclear energy gradient ∂Ẽ0,1 ∂R and the 1-NAC at various bond lengths R are requested by the TSH program code, and we feed it with the values interpolated from the results of the SSVQE and the 1-NAC experiment. The details of the interpolation are described in Appendix B. We prepare a set of 500 molecular geometries and nuclear velocities as a harmonic-oscillator Wigner distribution for the vibrational ground state at the equilibrium geometrical structure in the electronic ground state 1 1 Σ + . We run the trajectories from the first electronic excited state 2 1 Σ + with a time step of 0.1 fs and find that the trajectories hop from 2 1 Σ + to 1 1 Σ + where the dissociation occurs as shown in Fig. 3(b). The dynamics of the populations of the 1 1 Σ + state and 2 1 Σ + state is calculated based on the results of 500 trajectories and shown in Fig. 3(c). Since the experimentally-obtained 1-NAC values are smaller than the noiseless simulation ones (Fig. 3(a)), the decay of the 2 1 Σ + population calculated by using the experimental values of 1-NAC is slightly slower than that is calculated by using the noiseless simulation values of 1-NAC. Nevertheless, the overall dynamics of the populations is similar to each other and this indicates the possibility of performing TSH in a quantum device in the near future.
Although we consider the 1-NAC throughout this section, we perform additional numerical demonstrations of our methods by simulating the quantum circuits to calculate the 1-NAC and 2-NAC of the hydrogen molecules and Berry's phase of a two-site spin model in Appendix D.

VII. DISCUSSION
Our proposed methods presented in the previous sections are based on the analytical derivative of the eigenstates obtained by the SSVQE. This section compares our methods with those using the numerical derivative of the eigenstates by the finite difference method. The comparison will be made from two points of view: (1) the number of distinct Hamiltonians to perform the SSVQE to obtain the optimal circuit-parameters θ * and (2) the total number of measurements required to evaluate the NACs and Berry's phase after performing the SSVQE. We ignore the cost of classical computation throughout analyses in this section.
We recall that N θ and N x are the dimensions of the circuit-parameters and system-parameters, respectively. The number of qubits in the system is denoted as n and the number of Pauli terms in the Hamiltonian as N H . For quantum chemistry problems N H is typically O(n 4 ) [2,3], but several methods for reducing N H are proposed [86,87]. The detailed derivations of the formulas in this section (Eqs. (32), (34), and (37)) are presented in Appendix C.

A. Cost of 1-NAC
Let us consider our proposed method to calculate the 1-NAC d I kl with fixed k, l and all I = 1, . . . N x at some fixed system-parameters R. The evaluation of the 1-NAC is performed by calculating the denominator and numerator in Eq. (3). Both can be obtained by applying a single optimized ansatz circuit U ( θ * ( R)) resulting from the SSVQE to several initial states and measuring appropriate observables, so the number of distinct Hamiltonians to perform the SSVQE is just one. Meanwhile, the number of measurements to calculate the 1-NAC is estimated as follows. When we write the Hamiltonian as H( R) = ∂h i ( R)/∂R I ϕ k |P i |ϕ l are necessary to compute the 1-NAC. The term ϕ k(l) |P i |ϕ k(l) is evaluated as the expectation value of P i , and similarly the term ϕ k |P i |ϕ l (i = 1, . . . , N H ) is evaluated as a sum of four expectation values in the right hand sides of Eq. (13). By taking into account errors in evaluating those expectations values, the number of measurements to estimate the 1-NAC within the error is given by where A I (∆E k,l ) is the numerator (denominator) of Eq. (3), H = i |h i |, ∂H/∂R I = i |∂h i /∂R I |, and I * = argmax I |∆E k,l | ∂H ∂R I + H |A I | . It scales with the number of the Pauli terms in the Hamiltonian but does not depend on the number of circuit-parameters N θ by virtue of Eq. (3). The dependence on I, or the system-parameters, is also absent because it is absorbed into the classical computation of the coefficient ∂h i ( R)/∂R I [29].
To compare with our method, one can consider a method to evaluate the 1-NAC based on numerical differentiation of the eigenstates obtained by the SSVQE. In such approach, the 1-NAC can be evaluated by the following formula, where τ k,l ( R, R ± h e I )= ϕ k ( θ * ( R))|ϕ l ( θ * ( R ± h e I )) , h is a positive number, and e I is the unit vector in I-th direction. For simplicity, we will represent τ k,l ( R, R ± h e I ) as τ ±,I k,l . When evaluating each term of Eq. (33), we assume that τ ±,I k,l is estimated from the overlap , which can be easily evaluated from measurements on near-term quantum devices if |ψ k , |ψ l are computational basis states [25]. We then obtain the value of τ ±,I k,l by taking the square root of the overlaps with a positive sign (real value) [6]. This treatment can be justified when we solve problems in quantum chemistry, where wavefunctions are often real, and we adopt an ansatz which produces a real wavefunction. This approach to evaluate Eq. (33) avoids the costly Hadamard test and is considered to be feasible on nearterm quantum devices. We note that our methods are always applicable without the assumption above.
To evaluate the 1-NAC with Eq. (33), we need optimal parameters θ * ( R), θ * ( R ± h e I ) for all I, so the number of distinct Hamiltonians to perform the SSVQE in the finite difference method is 2N x + 1 = O(N x ). By considering the error in estimating the overlaps, the number of measurements to estimate the 1-NAC with the precision of in the finite difference method is at least where T k,l = min σ=±,I τ σ,I k,l and we assume the condition k,l (s) and τ (I, 3) k,l (s) = d 3 ds 3 τ k,l ( R, R + s e I ). Both our method (32) and the finite difference method (34) scale with −2 , so the prefactors determine the efficiency of them. When N x , or the number of system-parameters (nuclei of the molecule), becomes large, the finite difference method will suffer from a large number of the SSVQE runs and the measurements compared with our method.

B. Cost of 2-NAC
To calculate the 2-NAC D I kl with fixed k, l and all I = 1, . . . N x at R with our method, we require one optimized circuit-parameter θ * ( R), so the number of distinct Hamiltonians to perform the SSVQE is again one. Let us consider the number of measurements. We need the derivatives { is very complicated, we here let the error of the solutions be (see Ref. [29] for similar discussion). The value of ϕ k |∂ a ∂ b ϕ l (a, b = 1, . . . , N θ ) in Eq. (25) is obtained by measuring O(N 2 θ / 2 ) times within error . Similarly, ϕ k |∂ c ϕ l (c = 1, . . . , N θ ) in Eq. (26) is calculated by O(N θ / 2 ) measurements. Therefore, the total number of measurements to evaluate the 2-NAC by our method is roughly given by where the meaning of has to be cared for. It does not depend on the number of system-parameters N x . The finite difference method is based on the following formula: where we used τ k,l ( R, R) = 0. Similarly to the case of the 1-NAC, the number of distinct Hamiltonians to perform the SSVQE is O(N x ). To bound the error of the 2-NAC by , the finite difference method requires at least For calculating Berry's phase, the closed path C is discretized into K points. The integrand (Eq. (27)) is evaluated at all K points and numerically integrated both in our method and the finite difference method. We therefore compare our method with the finite difference method only in terms of the cost to obtain the integrand at all the discretized points.
In our method, the integrand at each discretized point can be evaluated by Eq. (29). The total number of distinct Hamiltonians to perform the VQE is K +1 = O(K). If we bound the error in estimating each integrand by , the total number of measurements is given by In the finite difference, the integrand can be evaluated with the finite difference method by the following formula, where v k ∝ R k+1 − R k is the unit vector along the closed loop C. The minimal number of measurements to obtain all integrands within error can be derived in a similar way for the 1-NAC, and the result is

VIII. CONCLUSION
In this paper, we have proposed methods to calculate the NACs and Berry's phase based on the VQE. We utilize the SSVQE and the MCVQE algorithms, which enable us to evaluate transition amplitudes of observables between approximate eigenstates. We explicitly present quantum circuits and classical post-processings to evaluate the NACs and Berry's phase in the framework of the VQE. For the 1-NAC, the calculations are simplified by taking advantage of the formula (3). The 2-NAC is obtained by combining the projective measurements and the expectation-value measurements of Pauli operators. The evaluation of Berry's phase is also carried out by the measurements of expectation values of Pauli operators with numerical integration of the definition of Berry's phase in addition to performing the Hadamard test once. We note that our method for calculating Berry's phase is applicable for molecular systems which have the conical intersection [88]. To show the potential feasibility of our method for the 1-NAC on a near-term quantum device, we evaluate the value of the 1-NAC of a lithium fluoride molecule on the IBM Q processor. Based on those results, we perform the nonadiabatic molecular dynamics simulation of photodissociation of a lithium fluoride for the first time. The methods given in the present paper contribute to enlarging the usage of the VQE and accelerate further developments to investigate quantum chemistry and quantum many-body problems on nearterm quantum devices.
We lastly comment on the effect of the barren plateau problem [89][90][91][92][93] to our methods. The barren plateau problem states that the gradients of the expectation values of observables with respect to ansatz circuit parameters vanish exponentially with the increase of the number of qubits when the ansatz has enough expressibility. When the barren plateau occurs, it is difficult to obtain the optimal circuit parameters θ * because the gradients become too small to optimize the parameters. Our proposed methods in this paper totally discuss the procedures after the optimal circuit parameters have been obtained (i.e., the VQE has successfully converged). Although our methods do not work when we cannot obtain θ * , several techniques [94][95][96][97][98] to avoid or ameliorate the barren plateau for the VQE and other variational quantum algorithms have been proposed.
/∂R I P i as the same in the main text. The Hoeffding's inequality [104] implies that an expectation value P i can be estimated within the precision P with high probability 1 − δ by measuring P i for O(ln(1/δ)/ 2 P ) times. When we perform O(ln(1/δ)/ 2 P ) measurements for estimating each P i , the total error in estimating the expectation value of H( R) is where · · · is an estimated value of · · · and H = N H i=1 |h i |. In the same way, the error of ∂H( R) where ∂H/∂R I = N H i=1 |∂h i /∂R I |. By using Eqs. (C1)(C2), we can derive the total number of measurements needed to estimate the 1-NAC within error . Equation (3) is evaluated in our method as When we perform O(ln(1/δ)/ 2 P ) measurements to estimate the expectation value of each P i appearing in ϕ ± k,l | ∂H( R) ∂R I |ϕ ± k,l , ϕ k(l) |H( R)|ϕ k(l) , the error propagation follows that the total error is given by Recalling that the 1-NACs for all I = 1, . . . , N x are obtained by the same results of the measurements P i , we conclude that the total number of measurements to estimate all the 1-NACs within the error with probability 1 − δ is given by where I * = argmax I |∆E k,l | ∂H ∂R I + H |A I | .

Cost of the finite difference method for 1-NAC
Let us discuss the number of measurements to calculate the 1-NAC with the finite difference method based on Eq. (33).
When the error of the overlap s ±,I k,l is bounded by s , or s ±,I k,l − s ±,I k,l ≤ s , it follows that The error of the 1-NAC with the finite difference method is then be expressed as .

(C9)
Therefore, the total number of measurements needed to evaluate the 1-NAC for all I = 1, . . . , N x within the precision with the finite difference method is where T 2 k,l = min I τ I k,l and M 3 = max I M (I) 3 . Moreover, to clarify the dependence of , we take h such that s attains the maximum with respect to h and that N 1-NAC total takes the minimum. This is realized for h = 3 T k,l . In such case, it follows that and obtain where M 4 = max I M (I) 4 . Similarly to the analysis for the 1-NAC, when we take h = 3 M4 , we obtain s,max = In this section, we demonstrate our methods for calculating the NACs, the DBOC, and Berry's phase by numerical simulations. Regarding the NACs, we consider the different electronic states of the hydrogen molecules. For the DBOC, we also take the electronic state of the hydrogen molecules. As for Berry's phase, we take a simple two-site spin model with a "twist" parameter where Berry's phase is quantized. In all the cases, numerical simulations of our method exhibit almost perfect agreement with the exact results. In addition, we can reproduce the shift of the equilibrium distance of the hydrogen atom by adding the DBOC to the potential energy curve obtained by the VQE [55]. These results further validates our methods proposed in the main text.

NACs of the hydrogen molecule
In the numerical simulation of the NACs and the DBOC, the electronic Hamiltonians of the hydrogen molecules are prepared in bond lengths from 0.5Å to 2.0Å with the interval of 0.1Å. Furthermore, we arrange the electronic Hamiltonian around the equilibrium point from 0.7320Å to 0.7350Å fine enough to see the shift of the equilibrium distance with the interval of 0.0001Å. We perform the standard Hartree-Fock calculation by employing STO-3G minimal basis set and compute the fermionic second-quantized Hamiltonian [2,3] with opensource libraries PySCF [102] and OpenFermion [103]. The Hamiltonians are mapped to the sum of the Pauli operators (qubit Hamiltonians) by the Jordan-Wigner transformation [105].
The SSVQE algorithm for the qubit Hamiltonians is executed with an ansatz consisting of SO(4) gates [23] FIG. 5. Ansatz quantum circuit for the VQE of the hydrogen molecule [23]. Each θ has six parameters, and RY (θ) = e −i θ 2 Y . The total number of parameters is 36.
shown in Fig. 5. This ansatz gives real-valued wavefunctions for any parameters θ. To obtain charge-neutral and spin-singlet eigenstates, we add penalty terms containing the total particle number operatorN and the total spin squared operatorŜ 2 to the Hamiltonian whose expectation value is to be minimized [79]. The cost function is where N 0 = 2 is the number of electrons and β S = β N = 10 are the penalty coefficients. We choose M = 3 and obtain the singlet ground state for i = 0 and another electronic state for i = 2 which has a non-zero value of NACs between the ground state (i.e., having the same symmetry as the ground state). The reference states and the weights are taken as |ψ 0 = |0011 , |ψ 1 = |0101 , |ψ 2 = |0110 and w 0 = 3, w 1 = 2, w 2 = 1. The circuitparameters θ are optimized by the BFGS algorithm implemented in Scipy library [72]. All simulations are run by the high-speed quantum circuit simulator Qulacs [106]. The results of the numerical calculation are shown in Fig. 6. We calculate the 1-NAC d 02 and 2-NAC D 02 between the ground state i = 0 (S 0 state) and the excited state i = 2 (S 2 state) as well as the DBOC of the ground state (S 0 state). The results are in agreement with the values computed by numerical differentiation of the full configuration interaction (Full-CI) results based on the definition of the NACs (Eq. (1) and Eq. (2) in the main text). In addition, the result shown in Fig. 6(c) ex-hibits the shift of the equilibrium distance from 0.7349Å to 0.7348Å by considering the DBOC based on Eq. (18) when k = l = 0. As mentioned in the main text, the 2-NAC also has indefinite sign, thus here we determine sign by considering the continuity of the 2-NAC with respect to the nuclear coordinate.

Berry's phase of twisted 2-spin model
To demonstrate our method for Berry's phase, we use a two-site spin-1/2 model with a twist. The Hamiltonian is defined as where S ± i = 1 2 (X i ± Y i ), S z i = 1 2 Z i , and ρ is a twist angle. ∆ is the parameter determines type and strength of the interaction between spins. The ground state for −1 < ∆ of this model is while for ∆ < −1 it is degenerate as Since H(ρ = 0) = H(ρ = 2π), we can consider Berry's phase Π C associated to the closed path C from ρ = 0 to ρ = 2π. From the exact expression of the ground state above, the analytical values of Π C can be calculate as Π C = π for −1 < ∆ and Π C = 0 for ∆ < −1. Berry's phase Π C , in this case, is called the local Z 2 Berry's phase and known to detect the topological nature of the ground state of quantum many-body systems [45,46]. We perform the VQE for the model (D2) with the ansatz depicted in Fig. 7. Again, the BFGS algorithm implemented in Scipy library [72] is used and all quantum circuit simulations are run by Qulacs in the noiseless case. We discretize the path from ρ = 0 as ρ = 2π into 100 points uniformly and run the VQE at each point. The first term of Eq. (27) is calculated by the summation (30) and the phase difference arg( ϕ( θ * s )|ϕ( θ * t ) ) in Eq. (27) is evaluated by the Hadamard test in Fig. 2. The result is shown in Fig. 8. The value of Berry's phase Π C exhibits the sharp transition reflecting the change of the ground state. These results illustrate the validity of our method to calculate Berry's phase based on the VQE.  7. (a) Ansatz quantum circuit to find the ground state of Eq. (D2). It contains four parameters (x0, x1, θ, φ) to be optimized. (b) Definition of the particle-number-preserving gate A(θ, φ) [62]. We define R(θ, φ) as R(θ, φ) = RY (θ + π/2)RZ (φ + π), where RY (θ) = e iθY /2 and RZ (φ) = e iφZ/2 .